adjust to new C6 catalog

Catalog/C6_Catalog.py

@@ -5,6 +5,7 @@ import numpy as np
 import h5py as h5
 import healpy as hp
 import astropy.constants as cons
+import traceback
 from astropy.coordinates import spherical_to_cartesian
 from astropy.table import Table
 from scipy import interpolate
@@ -14,6 +15,11 @@ from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar
 from ObservationSim.MockObject._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist
 from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position
 
+# (TEST)
+from astropy.cosmology import FlatLambdaCDM
+from astropy import constants
+from astropy import units as U
+
 try:
     import importlib.resources as pkg_resources
 except ImportError:
@@ -22,26 +28,40 @@ except ImportError:
 
 NSIDE = 128
 
+def get_bundleIndex(healpixID_ring, bundleOrder=4, healpixOrder=7):
+    assert NSIDE == 2**healpixOrder
+    shift = healpixOrder - bundleOrder
+    shift = 2*shift
+
+    nside_bundle = 2**bundleOrder
+    nside_healpix= 2**healpixOrder
+
+    healpixID_nest= hp.ring2nest(nside_healpix, healpixID_ring)
+    bundleID_nest = (healpixID_nest >> shift)
+    bundleID_ring = hp.nest2ring(nside_bundle, bundleID_nest)
+
+    return bundleID_ring
+
 class C6_Catalog(CatalogBase):
     def __init__(self, config, chip, pointing, chip_output, filt, **kwargs):
         super().__init__()
         self.cat_dir = os.path.join(config["data_dir"], config["input_path"]["cat_dir"])
         self.seed_Av = config["random_seeds"]["seed_Av"]
 
+        # (TEST)
+        self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111)
+
         self.chip_output = chip_output
         self.filt = filt
         self.logger = chip_output.logger
 
         with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path:
             self.normF_star = Table.read(str(filter_path))
-        # with pkg_resources.path('Catalog.data', 'lsst_throuput_g.fits') as filter_path:
-        #     self.normF_galaxy = Table.read(str(filter_path))
         
         self.config = config
         self.chip = chip
         self.pointing = pointing
 
-        # (DEBUG)
         self.max_size = 0.
 
         if "star_cat" in config["input_path"] and config["input_path"]["star_cat"] and not config["run_option"]["galaxy_only"]:
@@ -80,8 +100,13 @@ class C6_Catalog(CatalogBase):
         ra_min, ra_max, dec_min, dec_max = self.sky_coverage.xmin, self.sky_coverage.xmax, self.sky_coverage.ymin, self.sky_coverage.ymax
         ra = np.deg2rad(np.array([ra_min, ra_max, ra_max, ra_min]))
         dec = np.deg2rad(np.array([dec_max, dec_max, dec_min, dec_min]))
-        vertices = spherical_to_cartesian(1., dec, ra)
-        self.pix_list = hp.query_polygon(NSIDE, np.array(vertices).T, inclusive=True)
+        # vertices = spherical_to_cartesian(1., dec, ra)
+        self.pix_list = hp.query_polygon(
+            NSIDE,
+            hp.ang2vec(np.radians(90.) - dec, ra),
+            inclusive=True
+        )
+        # self.pix_list = hp.query_polygon(NSIDE, np.array(vertices).T, inclusive=True)
         if self.logger is not None:
             msg = str(("HEALPix List: ", self.pix_list))
             self.logger.info(msg)
@@ -144,16 +169,21 @@ class C6_Catalog(CatalogBase):
             )
 
         for igals in range(ngals):
-            # (TEST)
+            # # (TEST)
             # if igals > 100:
             #     break
+            # if igals != 1186:
+            #     continue
             
             param = self.initialize_param()
             param['ra'] = ra_arr[igals]
             param['dec'] = dec_arr[igals]
             param['ra_orig'] = gals['ra'][igals]
             param['dec_orig'] = gals['dec'][igals]
-            # param['mag_use_normal'] = gals['mag_true_g_lsst'][igals]
+            param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals]
+            if self.filt.is_too_dim(mag=param['mag_use_normal'], margin=self.config["obs_setting"]["mag_lim_margin"]):
+                continue
+
             # if param['mag_use_normal'] >= 26.5:
             #     continue
             param['z'] = gals['redshift'][igals]
@@ -165,11 +195,16 @@ class C6_Catalog(CatalogBase):
             param['e2'] = gals['ellipticity_true'][igals][1]
             
             # For shape calculation
+            
+            param['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2)
+            if param['ell_total'] > 0.9:
+                continue
             param['e1_disk'] = param['e1']
             param['e2_disk'] = param['e2']
             param['e1_bulge'] = param['e1']
             param['e2_bulge'] = param['e2']
 
+
             param['delta_ra'] = 0
             param['delta_dec'] = 0
 
@@ -211,7 +246,7 @@ class C6_Catalog(CatalogBase):
             # TEMP
             self.ids += 1
             # param['id'] = self.ids
-            param['id'] = str(pix_id) + '%02d'%(cat_id) + '%08d'%(igals)
+            param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
             
             if param['star'] == 0:
                 obj = Galaxy(param, self.rotation, logger=self.logger)
@@ -264,8 +299,8 @@ class C6_Catalog(CatalogBase):
                 input_time_str=time_str
             )
         for istars in range(nstars):
-            # (TEST)
-            # if istars > 10:
+            # # (TEST)
+            # if istars > 100:
             #     break
 
             param = self.initialize_param()
@@ -283,7 +318,6 @@ class C6_Catalog(CatalogBase):
             # if param['mag_use_normal'] >= 26.5:
             #     continue
             self.ids += 1
-            # param['id'] = self.ids
             param['id'] = stars['sourceID'][istars]
             param['sed_type'] = stars['sourceID'][istars]
             param['model_tag'] = stars['model_tag'][istars]
@@ -316,15 +350,28 @@ class C6_Catalog(CatalogBase):
                     print(e)
         
         if "galaxy_cat" in self.config["input_path"] and self.config["input_path"]["galaxy_cat"] and not self.config["run_option"]["star_only"]:
-            for i in range(76):
-                file_path = os.path.join(self.galaxy_path, "galaxies%04d_C6.h5"%(i))
-                gals_cat = h5.File(file_path, 'r')['galaxies']
+            # for i in range(76):
+            #     file_path = os.path.join(self.galaxy_path, "galaxies%04d_C6.h5"%(i))
+            #     gals_cat = h5.File(file_path, 'r')['galaxies']
+            #     for pix in self.pix_list:
+            #         try:
+            #             gals = gals_cat[str(pix)]
+            #             self._load_gals(gals, pix_id=pix, cat_id=i)
+            #             del gals
+            #         except Exception as e:
+            #             traceback.print_exc()
+            #             self.logger.error(str(e))
+            #             print(e)
             for pix in self.pix_list:
                 try:
+                    bundleID  = get_bundleIndex(pix)
+                    file_path = os.path.join(self.galaxy_path, "galaxies_C6_bundle{:06}.h5".format(bundleID))
+                    gals_cat = h5.File(file_path, 'r')['galaxies']
                     gals = gals_cat[str(pix)]
-                        self._load_gals(gals, pix_id=pix, cat_id=i)
+                    self._load_gals(gals, pix_id=pix, cat_id=bundleID)
                     del gals
                 except Exception as e:
+                    traceback.print_exc()
                     self.logger.error(str(e))
                     print(e)
         if self.logger is not None:
@@ -333,6 +380,19 @@ class C6_Catalog(CatalogBase):
         else:
             print("number of objects in catalog: ", len(self.objs))
 
+    # (TEST)
+    # def convert_mags(self, obj):
+    #     spec = np.matmul(obj.coeff, self.pcs.T)
+    #     lamb = self.lamb_gal * U.angstrom
+    #     unit_sed = ((lamb * lamb)/constants.c*U.erg/U.second/U.cm**2/U.angstrom).to(U.jansky)
+    #     unit_sed = unit_sed.value
+    #     lamb = lamb.value
+    #     lamb *= (1 + obj.z)
+    #     spec *= (1 + obj.z)
+    #     mags = -2.5 * np.log10(unit_sed*spec) + 8.9
+    #     mags += self.cosmo.distmod(obj.z).value
+    #     return lamb, mags
+
 
     def load_sed(self, obj, **kwargs):
         if obj.type == 'star':
@@ -344,11 +404,13 @@ class C6_Catalog(CatalogBase):
                 feh=obj.param['feh']
             )
         elif obj.type == 'galaxy' or obj.type == 'quasar':
-            dist_L_pc = (1 + obj.z) * comoving_dist(z=obj.z)[0]
-            factor = (10 / dist_L_pc)**2
+            # dist_L_pc = (1 + obj.z) * comoving_dist(z=obj.z)[0]
+            # factor = (10 / dist_L_pc)**2
+            factor = 10**(-.4 * self.cosmo.distmod(obj.z).value)
             flux = np.matmul(self.pcs, obj.coeff) * factor
-            if np.any(flux < 0):
-                raise ValueError("Glaxy %s: negative SED fluxes"%obj.id)
+            # if np.any(flux < 0):
+            #     raise ValueError("Glaxy %s: negative SED fluxes"%obj.id)
+            flux[flux < 0] = 0.
             sedcat = np.vstack((self.lamb_gal, flux)).T
             sed_data = getObservedSED(
                 sedCat=sedcat,
@@ -368,15 +430,16 @@ class C6_Catalog(CatalogBase):
         all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
         sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
         
-        if obj.type == 'galaxy' or obj.type == 'quasar':
-            # integrate to get the magnitudes
-            sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T
-            # print("sed_photon = ", sed_photon)
-            sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest')
-            sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False)
-            interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=self.filt.bandpass_full)
-            obj.param['mag_use_normal'] = getABMAG(interFlux, self.filt.bandpass_full)
-            # print("mag_use_normal = %.3f"%obj.param['mag_use_normal'])
+        # if obj.type == 'galaxy' or obj.type == 'quasar':
+        #     # integrate to get the magnitudes
+        #     sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T
+        #     sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest')
+        #     sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False)
+        #     interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=self.filt.bandpass_full)
+        #     obj.param['mag_use_normal'] = getABMAG(interFlux, self.filt.bandpass_full)
+        #     # print("mag_use_normal = %.3f"%obj.param['mag_use_normal'])
+        #     # mag = getABMAG(interFlux, self.filt.bandpass_full)
+        #     # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal']))
         
         del wave
         del flux

ObservationSim/Instrument/Chip/Chip.py

@@ -229,6 +229,7 @@ class Chip(FocalPlane):
 
 
     def getSkyCoverage(self, wcs):
+        # print("In getSkyCoverage: xmin = %.3f, xmax = %.3f, ymin = %.3f, ymax = %.3f"%(self.bound.xmin, self.bound.xmax, self.bound.ymin, self.bound.ymax))
         return super().getSkyCoverage(wcs, self.bound.xmin, self.bound.xmax, self.bound.ymin, self.bound.ymax)
 
 

ObservationSim/MockObject/Galaxy.py

@@ -137,6 +137,7 @@ class Galaxy(MockObject):
                 # return False
                 continue
             nphotons_sum += nphotons
+
             # print("nphotons_sub-band_%d = %.2f"%(i, nphotons))
 
             psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
@@ -171,8 +172,10 @@ class Galaxy(MockObject):
             # photons_list.append(photons)
 
             stamp = gal.drawImage(wcs=real_wcs_local, method='phot', offset=offset, save_photons=True)
+            
             xmax = max(xmax, stamp.xmax)
             ymax = max(ymax, stamp.ymax)
+
             photons = stamp.photons
             photons.x += x_nominal
             photons.y += y_nominal
@@ -184,6 +187,8 @@ class Galaxy(MockObject):
         stamp.wcs = real_wcs_local
         stamp.setCenter(x_nominal, y_nominal)
         bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1)
+
+        if bounds.area() > 0:
             chip.img.setOrigin(0, 0)
             stamp[bounds] = chip.img[bounds]
 
@@ -237,7 +242,7 @@ class Galaxy(MockObject):
 
     def drawObj_slitless(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0,
                          exptime=150., normFilter=None, grating_split_pos=3685):
-
+        if normFilter is not None:
             norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001
             sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'], spectrum=self.sed,
                                                         norm_thr=normFilter,
@@ -245,6 +250,8 @@ class Galaxy(MockObject):
                                                         eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0]))
             if sedNormFactor == 0:
                 return False
+        else:
+            sedNormFactor = 1.
         normalSED = Table(np.array([self.sed['WAVELENGTH'], self.sed['FLUX'] * sedNormFactor]).T,
                           names=('WAVELENGTH', 'FLUX'))
 

ObservationSim/MockObject/MockObject.py

 import galsim
 import numpy as np
 import astropy.constants as cons
+from astropy import wcs
 from astropy.table import Table
 
 from ObservationSim.MockObject._util import magToFlux, VC_A, convolveGaussXorders
@@ -99,8 +100,6 @@ class MockObject(object):
         dy = y - self.y_nominal
         self.offset = galsim.PositionD(dx, dy)
 
-        from astropy import wcs
-
         if img_header is not None:
             self.real_wcs = galsim.FitsWCS(header=img_header)
         else:
@@ -224,6 +223,11 @@ class MockObject(object):
         stamp.wcs = real_wcs_local
         stamp.setCenter(x_nominal, y_nominal)
         bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1)
+        
+        # # (DEBUG)
+        # print("stamp bounds: ", stamp.bounds)
+        # print(bounds)
+        if bounds.area() > 0:
             chip.img.setOrigin(0, 0)
             stamp[bounds] = chip.img[bounds]
             for i in range(len(photons_list)):

ObservationSim/MockObject/_util.py

@@ -455,7 +455,10 @@ def getObservedSED(sedCat, redshift=0.0, av=0.0, redden=0):
     sw, sf = sedCat[:,0], sedCat[:,1]
     # reddening
     sf = reddening(sw, sf, av=av, model=redden)
-    sw, sf = sw*z, sf*(z**3)
+    # sw, sf = sw*z, sf*(z**3)
+    sw, sf = sw*z, sf/z
+    # sw, sf = sw*z, sf
+
     # lyman forest correction
     sf = lyman_forest(sw, sf, redshift)
     isedObs = (sw.copy(), sf.copy())

ObservationSim/ObservationSim.py

@@ -109,7 +109,7 @@ class Observation(object):
                 input_vx=pointing.sat_vx,
                 input_vy=pointing.sat_vy,
                 input_vz=pointing.sat_vz,
-                input_epoch="J2015.5",
+                input_epoch="J2000",
                 input_date_str=date_str,
                 input_time_str=time_str
             )
@@ -187,6 +187,7 @@ class Observation(object):
                 #     break
 
                 obj = self.cat.objs[j]
+
                 if obj.type == 'star' and self.config["run_option"]["galaxy_only"]:
                     continue
                 elif obj.type == 'galaxy' and self.config["run_option"]["star_only"]:
@@ -210,14 +211,13 @@ class Observation(object):
                         target_filt=cut_filter, 
                         norm_filt=norm_filt,
                     )
-
                 except Exception as e:
                     # print(e)
                     traceback.print_exc()
                     chip_output.logger.error(e)
                     continue
-                
-                # chip_output.logger.info("debug point #1")
+                # print("mag_%s = %.3f"%(filt.filter_type, obj.param["mag_%s"%filt.filter_type]))
+                # chip_output.logger.info("mag_%s = %.3f"%(filt.filter_type, obj.param["mag_%s"%filt.filter_type]))
 
                 # Exclude very bright/dim objects (for now)
                 # if filt.is_too_bright(mag=obj.getMagFilter(filt)):
@@ -243,8 +243,6 @@ class Observation(object):
                     # print(obj.getMagFilter(filt))
                     continue
 
-                # chip_output.logger.info("debug point #2")
-
                 if self.config["shear_setting"]["shear_type"] == "constant":
                     if obj.type == 'star':
                         obj.g1, obj.g2 = 0., 0.
@@ -258,15 +256,12 @@ class Observation(object):
                         # print("failed to load external shear.")
                         chip_output.logger.error("failed to load external shear.")
                         pass
-
-                    # chip_output.logger.info("debug point #3")
                 elif self.config["shear_setting"]["shear_type"] == "catalog":
                     pass
                 else:
                     chip_output.logger.error("Unknown shear input")
                     raise ValueError("Unknown shear input")
 
-                # chip_output.logger.info("debug point #4")
                 header_wcs = generateExtensionHeader(
                     xlen=chip.npix_x,
                     ylen=chip.npix_y,
@@ -283,18 +278,18 @@ class Observation(object):
                     extName='raw')
 
                 pos_img, offset, local_wcs, real_wcs = obj.getPosImg_Offset_WCS(img=chip.img, fdmodel=self.fd_model, chip=chip, verbose=False, img_header=header_wcs)
+                # pos_img, offset, local_wcs, real_wcs = obj.getPosImg_Offset_WCS(img=chip.img, fdmodel=self.fd_model, chip=chip, verbose=True, img_header=header_wcs)
                 if pos_img.x == -1 or pos_img.y == -1:
                     # Exclude object which is outside the chip area (after field distortion)
-                    # print("obj missed!!")
+                    # print('obj_ra = ', obj.ra, 'obj_dec = ', obj.dec, 'obj_ra_orig = ', obj.ra_orig, 'obj_dec_orig = ',obj.dec_orig)
+                    # print("obj missed: %s"%(obj.id))
+                    chip_output.logger.error("obj missed: %s"%(obj.id))
                     missed_obj += 1
                     obj.unload_SED()
                     continue
 
-                # chip_output.logger.info("debug point #5")
-
                 # Draw object & update output catalog
                 try:
-                    # chip_output.logger.info("debug point #6")
                     # chip_output.logger.info("current filter type: %s"%filt.filter_type)
                     if self.config["run_option"]["out_cat_only"]:
                         isUpdated = True
@@ -313,6 +308,7 @@ class Observation(object):
                             g2=obj.g2, 
                             exptime=pointing.exp_time
                             )
+
                     elif chip.survey_type == "spectroscopic" and not self.config["run_option"]["out_cat_only"]:
                         isUpdated, pos_shear = obj.drawObj_slitless(
                             tel=self.tel, 
@@ -325,7 +321,6 @@ class Observation(object):
                             g2=obj.g2, 
                             exptime=pointing.exp_time,
                             normFilter=norm_filt)
-                    # chip_output.logger.info("debug point #7")
                     if isUpdated:
                         # TODO: add up stats
                         # print("updating output catalog...")

config/config_C6.yaml

@@ -12,7 +12,7 @@
 work_dir: "/share/home/fangyuedong/test_psf_rot/csst-simulation/workplace/"
 # work_dir: "/share/"
 data_dir: "/share/simudata/CSSOSDataProductsSims/data/"
-run_name: "C6_test"
+run_name: "C6_photometry"
 
 # (Optional) a file of point list 
 # if you just want to run default pointing:
@@ -49,7 +49,7 @@ obs_setting:
   # "Photometric": simulate photometric chips only
   # "Spectroscopic": simulate slitless spectroscopic chips only
   # "All": simulate full focal plane
-  survey_type: "All"
+  survey_type: "Photometric"
   
   # Exposure time [seconds]
   exp_time: 150.
@@ -74,7 +74,7 @@ obs_setting:
   # - give a list of indexes of pointings: [ip_1, ip_2...]
   # - run all pointings: null
   # Note: only valid when a pointing list is specified
-  run_pointings: [0]
+  run_pointings: null
   
   # Run specific chip(s):
   # - give a list of indexes of chips: [ip_1, ip_2...]
@@ -86,7 +86,7 @@ obs_setting:
   # astrometric_lib: "libshao.so"
   enable_astrometric_model: True
 
-  # Cut by saturation/limiting magnitude in which band?
+  # Cut by saturation magnitude in which band?
   cut_in_band: "z"
 
   # saturation magnitude margin
@@ -103,7 +103,7 @@ obs_setting:
 input_path:
   cat_dir: "Catalog_C6_20221212"
   star_cat: "C6_MMW_GGC_Astrometry_healpix.hdf5"
-  galaxy_cat: "cat2CSSTSim/"
+  galaxy_cat: "cat2CSSTSim_bundle/"
 
 SED_templates_path:
   star_SED: "Catalog_20210126/SpecLib.hdf5"

run_C6.pbs

 #!/bin/bash
 
 #PBS -N SIMS
-#PBS -l nodes=wcl-1:ppn=16+wcl-2:ppn=16+wcl-3:ppn=16+wcl-4:ppn=16+wcl-5:ppn=16+wcl-6:ppn=16
+#PBS -l nodes=wcl-1:ppn=60
+###PBS -l nodes=wcl-1:ppn=24+wcl-2:ppn=24+wcl-3:ppn=24+wcl-4:ppn=24+wcl-5:ppn=24+wcl-6:ppn=24
 #PBS -u fangyuedong
 ###PBS -j oe
 
 cd $PBS_O_WORKDIR
 
-NP=96
+NP=60
 date
 
 mpirun -np $NP python3 /share/home/fangyuedong/test_psf_rot/csst-simulation/run_sim.py config_C6.yaml -c /share/home/fangyuedong/test_psf_rot/csst-simulation/config