update saturation magnitude cut based on selected CSST band

Catalog/NGPCatalog.py

@@ -304,7 +304,8 @@ class NGPCatalog(CatalogBase):
             raise ValueError("Object type not known")
         speci = interpolate.interp1d(wave, flux)
         # lamb = np.arange(2500, 10001 + 0.5, 0.5)
-        lamb = np.arange(2400, 11001 + 0.5, 0.5)
+        # lamb = np.arange(2400, 11001 + 0.5, 0.5)
+        lamb = np.arange(2000, 18001 + 0.5, 0.5)
         y = speci(lamb)
         # erg/s/cm2/A --> photo/s/m2/A
         all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13

ObservationSim/Config/ChipOutput.py

@@ -33,7 +33,9 @@ class ChipOutput(object):
         fh = logging.FileHandler(os.path.join(self.subdir, logger_filename), mode='w+', encoding='utf-8')
         fh.setLevel(logging.DEBUG)
         self.logger.setLevel(logging.DEBUG)
-        formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s')
+	logging.getLogger('numba').setLevel(logging.WARNING)
+        # formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s')
+        formatter = logging.Formatter('%(asctime)s - %(msecs)d - %(levelname)-8s - [%(filename)s:%(lineno)d] - %(message)s')
         fh.setFormatter(formatter)
         self.logger.addHandler(fh)
 

ObservationSim/Instrument/Filter.py

@@ -40,8 +40,8 @@ class Filter(object):
         # self.filter_dir = filter_param.filter_dir
 
     def is_too_bright(self, mag):
-        # return mag <= self.mag_saturation - 1.0
-        return mag <= 14.0
+        return mag <= self.mag_saturation - 1.0
+        # return mag <= 14.0
 
     def is_too_dim(self, mag):
         return mag >= self.mag_limiting + 1.0
@@ -100,6 +100,8 @@ class Filter(object):
         return self.sky_background * exptime / gain
 
     def update_limit_saturation_mags(self, exptime=150., psf_fwhm=0.1969, skyFn='sky_emiss_hubble_50_50_A.dat', chip=None):
+        if self.filter_type in ["GI", "GV", "GU"]:
+            return
         if chip is not None:
             pix_scale = chip.pix_scale
             read_noise = chip.read_noise

ObservationSim/Instrument/_util.py

@@ -99,7 +99,7 @@ def calculateLimitMag(aperture = 2.0, psf_fwhm = 0.1969,pixelSize = 0.074, pmRat
 	yp, xp = np.mgrid[0:m_size, 0:m_size]
 	psfMap = gaussShape(xp, yp)
 	maxRatio = np.amax(psfMap)/np.sum(psfMap)
-	print(maxRatio)
+	# print(maxRatio)
 
 	flux_sat = fw/maxRatio*exNum
 	satMag = -2.5*log10(flux_sat/(54799275581.04437 * trapz(wavey*eff[:,1]/wave,wave, 0.1)*exTime*exNum*pi*(aperture/2)*(aperture/2)));

ObservationSim/MockObject/Galaxy.py

@@ -45,6 +45,7 @@ class Galaxy(MockObject):
             full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
         except Exception as e:
             print(e)
+            self.logger.error(e)
             return -1
         for i in range(len(bandpass_list)):
             bandpass = bandpass_list[i]
@@ -52,6 +53,7 @@ class Galaxy(MockObject):
                 sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
             except Exception as e:
                 print(e)
+                self.logger.error(e)
                 return -1
             
             ratio = sub/full
@@ -86,6 +88,7 @@ class Galaxy(MockObject):
             full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
         except Exception as e:
             print(e)
+            self.logger.error(e)
             return False
 
         nphotons_sum = 0
@@ -111,6 +114,7 @@ class Galaxy(MockObject):
                 sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
             except Exception as e:
                 print(e)
+                self.logger.error(e)
                 # return False
                 continue
             

ObservationSim/MockObject/MockObject.py

@@ -126,6 +126,7 @@ class MockObject(object):
             full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
         except Exception as e:
             print(e)
+            self.logger.error(e)
             return False
 
         nphotons_sum = 0
@@ -145,6 +146,7 @@ class MockObject(object):
                 sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
             except Exception as e:
                 print(e)
+                self.logger.error(e)
                 # return False
                 continue
         

ObservationSim/MockObject/Star.py

@@ -37,6 +37,7 @@ class Star(MockObject):
             full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
         except Exception as e:
             print(e)
+            self.logger.error(e)
             return -1
 
         for i in range(len(bandpass_list)):
@@ -46,6 +47,7 @@ class Star(MockObject):
                 sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
             except Exception as e:
                 print(e)
+                self.logger.error(e)
                 return -1
         
             ratio = sub/full

ObservationSim/ObservationSim.py

@@ -75,9 +75,6 @@ class Observation(object):
         chip_output.logger.info('Chip : %d' % chip.chipID)
         chip_output.logger.info(':::::::::::::::::::::::::::END:::::::::::::::::::::::::::::::::::')
 
-        # Update the limiting magnitude using exposure time in pointing
-        filt.update_limit_saturation_mags(exptime=pointing.exp_time, chip=chip)
-
         if self.config["psf_setting"]["psf_model"] == "Gauss":
             psf_model = PSFGauss(chip=chip)
         elif self.config["psf_setting"]["psf_model"] == "Interp":
@@ -181,7 +178,7 @@ class Observation(object):
                 elif obj.type == 'quasar' and self.config["run_option"]["star_only"]:
                     continue
 
-                # load SED
+                # load and convert SED; also caculate object's magnitude in all CSST bands
                 try:
                     sed_data = self.cat.load_sed(obj)
                     norm_filt = self.cat.load_norm_filt(obj)
@@ -191,26 +188,47 @@ class Observation(object):
                         target_filt=filt, 
                         norm_filt=norm_filt,
                     )
+                    for ifilt in range(len(self.filter_list)):
+                        temp_filter = self.filter_list[ifilt]
+                        _, obj.param["mag_%s"%temp_filter.filter_type] = self.cat.convert_sed(
+                            mag=obj.param["mag_use_normal"],
+                            sed=sed_data,
+                            target_filt=temp_filter, 
+                            norm_filt=norm_filt,
+                        )
+                        # Update the limiting magnitude using exposure time in pointing
+                        temp_filter.update_limit_saturation_mags(exptime=pointing.exp_time, chip=chip)
+
+                        # Select cutting band filter for saturation/limiting magnitude
+                        if temp_filter.filter_type.lower() == self.config["obs_setting"]["cut_in_band"].lower():
+                            cut_filter = temp_filter
+
                 except Exception as e:
-                    # print(e)
+                    print(e)
                     chip_output.logger.error(e)
                     continue
+                
+                # chip_output.logger.info("debug point #1")
 
                 # Exclude very bright/dim objects (for now)
                 # if filt.is_too_bright(mag=obj.getMagFilter(filt)):
-                if filt.is_too_bright(mag=obj.mag_use_normal):
+                # if filt.is_too_bright(mag=obj.mag_use_normal):
+                if cut_filter.is_too_bright(mag=obj.param["mag_%s"%self.config["obs_setting"]["cut_in_band"].lower()]):
                     # print("obj too birght!!", flush=True)
                     if obj.type != 'galaxy':
                         bright_obj += 1
                         obj.unload_SED()
                         continue
                 if filt.is_too_dim(mag=obj.getMagFilter(filt)):
+                # if cut_filter.is_too_dim(mag=obj.param["mag_%s"%self.config["obs_setting"]["cut_in_band"].lower()]):
                     # print("obj too dim!!", flush=True)
                     dim_obj += 1
                     obj.unload_SED()
                     # 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.
@@ -224,11 +242,15 @@ 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")
 
                 pos_img, offset, local_wcs = obj.getPosImg_Offset_WCS(img=chip.img, fdmodel=self.fd_model, chip=chip, verbose=False)
                 if pos_img.x == -1 or pos_img.y == -1:
@@ -238,8 +260,12 @@ class Observation(object):
                     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
                         pos_shear = 0.
@@ -268,6 +294,7 @@ class Observation(object):
                             exptime=pointing.exp_time,
                             normFilter=norm_filt,
                             )
+                    # chip_output.logger.info("debug point #7")
                     if isUpdated:
                         # TODO: add up stats
                         # print("updating output catalog...")
@@ -277,7 +304,7 @@ class Observation(object):
                         # print("object omitted", flush=True)
                         continue
                 except Exception as e:
-                    # print(e)
+                    print(e)
                     chip_output.logger.error(e)
                     pass
                 # Unload SED:

config/config_NGP.yaml

@@ -10,9 +10,9 @@
 # Can add some of the command-line arguments here as well;
 # OK to pass either way or both, as long as they are consistent
 # work_dir: "/public/home/fangyuedong/sim_code_release/CSST/test/"
-work_dir: "/public/home/fangyuedong/temp/CSST/workplace/"
+work_dir: "/public/home/fangyuedong/test/CSST/workplace/"
 data_dir: "/data/simudata/CSSOSDataProductsSims/data/"
-run_name: "NGP_Astrometry-on"
+run_name: "TEST_LIMITING_MAG"
 
 # (Optional) a file of point list 
 # if you just want to run default pointing:
@@ -48,7 +48,7 @@ obs_setting:
   # "Photometric": simulate photometric chips only
   # "Spectroscopic": simulate slitless spectroscopic chips only
   # "All": simulate full focal plane
-  survey_type: "Photometric"
+  survey_type: "All"
   
   # Exposure time [seconds]
   exp_time: 150.
@@ -73,7 +73,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: [ 5,  7, 11, 14, 19, 60, 70, 82, 88]
+  run_pointings: [0, 1]
   
   # Run specific chip(s):
   # - give a list of indexes of chips: [ip_1, ip_2...]
@@ -85,6 +85,9 @@ obs_setting:
   # astrometric_lib: "libshao.so"
   enable_astrometric_model: True
 
+  # Cut by saturation/limiting magnitude in which band?
+  cut_in_band: "g"
+
 ###############################################
 # Input path setting
 ###############################################
@@ -165,7 +168,7 @@ ins_effects:
   non_linear:     ON # Whether to add non-linearity
   cosmic_ray:     ON # Whether to add cosmic-ray
   cray_differ:    ON # Whether to generate different cosmic ray maps CAL and MS output
-  cte_trail:      OFF # Whether to simulate CTE trails
+  cte_trail:      ON # Whether to simulate CTE trails
   saturbloom:     ON # Whether to simulate Saturation & Blooming
   add_badcolumns: ON # Whether to add bad columns
   add_hotpixels:  ON # Whether to add hot pixels

run_NGP.pbs

@@ -12,8 +12,8 @@
 #PBS -q batch
 #PBS -u fangyuedong
 
-NP=30
+NP=80
 date
 echo $NP
 
-mpirun -np $NP python /public/home/fangyuedong/temp/CSST/run_sim.py config_NGP.yaml -c /public/home/fangyuedong/temp/CSST/config
+mpirun -np $NP python /public/home/fangyuedong/test/CSST/run_sim.py config_NGP.yaml -c /public/home/fangyuedong/test/CSST/config
\ No newline at end of file