fix the crashe when use a filter that is not in the list of run_chips for mag...

fix the crashe when use a filter that is not in the list of run_chips for mag cut. add new feature: add additional output columns to the .cat file, see NGPCatalog.py as an example

Catalog/NGPCatalog.py

@@ -23,15 +23,13 @@ except ImportError:
 NSIDE = 128
 
 class NGPCatalog(CatalogBase):
-    def __init__(self, config, chip, pointing, **kwargs):
+    def __init__(self, config, chip, pointing, chip_output, **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"]
 
-        if "logger" in kwargs:
-            self.logger = kwargs["logger"]
-        else:
-            self.logger = None
+        self.chip_output = chip_output
+        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))
@@ -58,9 +56,17 @@ class NGPCatalog(CatalogBase):
         else:
             self.rotation = 0.
 
+        # Update output .cat header with catalog specific output columns
+        self._add_output_columns_header()
+
         self._get_healpix_list()
         self._load()
     
+    def _add_output_columns_header(self):
+        self.add_hdr = " model_tag teff logg feh"
+        self.add_fmt = " %10s %8.4f %8.4f %8.4f"
+        self.chip_output.update_ouptut_header(additional_column_names=self.add_hdr)
+
     def _get_healpix_list(self):
         self.sky_coverage = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.2)
         ra_min, ra_max, dec_min, dec_max = self.sky_coverage.xmin, self.sky_coverage.xmax, self.sky_coverage.ymin, self.sky_coverage.ymax
@@ -181,9 +187,12 @@ class NGPCatalog(CatalogBase):
             
             if param['star'] == 0:
                 obj = Galaxy(param, self.rotation, logger=self.logger)
-                self.objs.append(obj)
             if param['star'] == 2:
                 obj = Quasar(param, logger=self.logger)
+            
+            # Need to deal with additional output columns
+            obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.)
+            
             self.objs.append(obj)
 
     def _load_stars(self, stars, pix_id=None):
@@ -249,6 +258,10 @@ class NGPCatalog(CatalogBase):
             param['z'] = 0.0
             param['star'] = 1   # Star
             obj = Star(param, logger=self.logger)
+
+            # Append additional output columns to the .cat file
+            obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh'])
+
             self.objs.append(obj)
 
     def _load(self, **kwargs):

ObservationSim/Config/ChipOutput.py

@@ -11,6 +11,7 @@ class ChipOutput(object):
         self.imgKey2 = imgKey2
         self.exptime = exptime
         self.mjdTime = mjdTime
+        self.pointing_type = pointing_type
         if (ra_cen is not None) and (dec_cen is not None):
             self.ra_cen = ra_cen
             self.dec_cen = dec_cen
@@ -42,78 +43,40 @@ class ChipOutput(object):
         hdr1  = "# obj_ID ID_chip filter xImage yImage ra dec ra_orig dec_orig z mag obj_type "
         hdr2  = "thetaR bfrac hlr_disk hlr_bulge e1_disk e2_disk e1_bulge e2_bulge g1 g2 "
         hdr3  = "sed_type av redden "
-        hdr4 = "pm_ra pm_dec RV parallax\n"
+        hdr4 = "pm_ra pm_dec RV parallax"
 
         fmt1  = "%10d %4d %5s %10.3f %10.3f %15.8f %15.8f %15.8f %15.8f %7.4f %8.4f %15s "
         fmt2  = "%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f "
         fmt3 = "%2d %8.4f %8.4f "
-        fmt4 = "%15.8f %15.8f %15.8f %15.8f\n"
+        fmt4 = "%15.8f %15.8f %15.8f %15.8f"
 
         self.hdr = hdr1 + hdr2 + hdr3 + hdr4
         self.fmt = fmt1 + fmt2 + fmt3 + fmt4
 
         self.logger.info("pointing_type = %s\n"%(pointing_type))
-        if pointing_type == 'MS':
+
+    
+    def update_ouptut_header(self, additional_column_names=""):
+        self.hdr += additional_column_names
+    
+    def create_output_file(self):
+        if self.pointing_type == 'MS':
             self.cat = open(os.path.join(self.subdir, self.cat_name), "w")
             self.logger.info("Creating catalog file %s ...\n"%(os.path.join(self.subdir, self.cat_name)))
+            if not self.hdr.endswith("\n"):
+                self.hdr += "\n"
             self.cat.write(self.hdr)
 
-    # def updateHDR(self, hdr):
-    #     hdrNew = [{"name":"RDNOISE", "value":self.chip.read_noise,        "comment":"read noise in e-/pixel"},
-    #             {"name":"DARK",    "value":self.chip.dark_noise,        "comment":"Dark noise (e-/pixel/s)"},
-    #             {"name":"EXPTIME", "value":self.exptime,          "comment":"exposure time in second"},
-    #             {"name":"GAIN",    "value":self.chip.gain,             "comment":"CCD gain in e-/ADU"},
-    #             {"name":"SATURATE","value":65535.0,          "comment":"saturation level"},
-    #             {"name":"CCDCHIP",  "value":int(self.chipLabel), "comment":"chip ID in the CCD mosaic"},
-    #             {"name":"FILTER",  "value":self.filt.filter_type,          "comment":"filter name"},
-    #             {"name":"MJD-OBS", "value":self.mjdTime,          "comment":"Modified Julian Date (MJD) of observation"},
-    #             {"name":"DATE-OBS","value":self.imgKey1,          "comment":"Date of observation"},
-    #             {"name":"EQUINOX", "value":2000.0},
-    #             {"name":"RADECSYS","value":"ICRS"},
-    #             {"name":"RA",      "value":self.ra_cen,           "comment":"telescope pointing center"},
-    #             {"name":"DEC",     "value":self.dec_cen,          "comment":"telescope pointing center"},
-    #             {"name":"OBJECT",  "value":"CSS-OS"},
-    #             {"name":"WCSDIM",  "value":2.0,              "comment":"WCS Dimensionality"},
-    #             {"name":"EXTNAME", "value":"IM1",            "comment":"Extension name"},
-    #             {"name":"BSCALE",  "value":1.0},
-    #             {"name":"BZERO",   "value":0.0},
-    #             {"name":"OBSID",   "value":self.imgKey0,          "comment":"Observation ID"},
-    #             {"name":"CCDNAME", "value":"ccd"+self.chipLabel,"comment":"CCD name"},
-    #             {"name":"RSPEED",  "value":10.0,             "comment":"Read speed"},
-    #             {"name":"CHIPTEMP","value":-100.0,           "comment":"Chip temperature"},
-    #             {"name":"DATASEC", "value":"1:%d,1:%d"%(self.chip.npix_x,self.chip.npix_y), "comment":"Data section"},
-    #             {"name":"CCDSUM",  "value":self.chip.npix_x*self.chip.npix_y,      "comment":"CCD pixel summing"},
-    #             {"name":"NSUM",    "value":self.chip.npix_x*self.chip.npix_y,      "comment":"CCD pixel summing"},
-    #             {"name":"AUTHOR",  "value":"CSST-Sim Group"},
-    #             {"name":"GROUP",   "value":"Weak Lensing Working Group for CSST"}]
-    #     for item in hdrNew:
-    #         hdr.add_record(item)
-    #     return hdr
-
-    # def cat_add_obj(self, obj, pos_img, snr, pos_shear, g1, g2):
     def cat_add_obj(self, obj, pos_img, pos_shear):
         ximg = pos_img.x - self.chip.bound.xmin + 1.0
         yimg = pos_img.y - self.chip.bound.ymin + 1.0
-        # if obj.type == 'galaxy':
-        #     line = self.fmt%(obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z, obj.getMagFilter(self.filt), obj.param["star"], obj.thetaR, obj.bfrac, obj.hlr_disk, obj.hlr_bulge,
-        #         obj.e1_disk, obj.e2_disk, obj.e1_bulge, obj.e2_bulge,
-        #         pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS, obj.sed_type, obj.param['av'], obj.param['redden'], 'n', 0, 0, 0)
-        # elif obj.type == "quasar":
-        #     line = self.fmt % (obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z,
-        #                        obj.getMagFilter(self.filt), obj.param["star"], 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
-        #                        pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS, obj.sed_type, obj.param['av'], obj.param['redden'], 'n', 0.0, 0.0, 0.0)
-        # else:
-        #     line = self.fmt%(obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z, obj.getMagFilter(self.filt), obj.param["star"], 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
-        #         pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS, 0, 0.0, 0.0, obj.param['model_tag'], obj.param['teff'], obj.param['logg'],obj.param['feh'])
-        # print(
-        #     obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.ra_orig, obj.dec_orig, obj.z, obj.getMagFilter(self.filt), obj.type, 
-        #     obj.thetaR, obj.bfrac, obj.hlr_disk, obj.hlr_bulge, obj.e1_disk, obj.e2_disk, obj.e1_bulge, obj.e2_bulge, obj.g1, obj.g2,
-        #     obj.sed_type, obj.av, obj.redden,
-        #     obj.pmra, obj.pmdec, obj.rv, obj.parallax)
         
         line = self.fmt%(
             obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.ra_orig, obj.dec_orig, obj.z, obj.getMagFilter(self.filt), obj.type, 
             obj.thetaR, obj.bfrac, obj.hlr_disk, obj.hlr_bulge, obj.e1_disk, obj.e2_disk, obj.e1_bulge, obj.e2_bulge, obj.g1, obj.g2,
             obj.sed_type, obj.av, obj.redden,
             obj.pmra, obj.pmdec, obj.rv, obj.parallax)
+        line += obj.additional_output_str
+        if not line.endswith("\n"):
+            line += "\n"
         self.cat.write(line)

ObservationSim/MockObject/CatalogBase.py

@@ -92,4 +92,3 @@ class CatalogBase(metaclass=ABCMeta):
         elif target_filt.survey_type == "spectroscopic":
             del sed_photon
             return sed, mag_csst
-        
\ No newline at end of file

ObservationSim/MockObject/MockObject.py

@@ -44,6 +44,7 @@ class MockObject(object):
         self.hlr_bulge = self.param["hlr_bulge"]
         self.e1_disk, self.e2_disk = 0., 0.
         self.e1_bulge, self.e2_bulge = 0., 0.
+        self.additional_output_str = ""
 
         self.logger = logger
 

ObservationSim/ObservationSim.py

@@ -25,6 +25,7 @@ class Observation(object):
         self.filter_param = FilterParam() 
         self.chip_list = []
         self.filter_list = []
+        self.all_filter = []
         self.Catalog = Catalog
 
         # if we want to apply field distortion?
@@ -37,8 +38,6 @@ class Observation(object):
         nchips = self.focal_plane.nchip_x*self.focal_plane.nchip_y
         for i in range(nchips):
             chipID = i + 1
-            if self.focal_plane.isIgnored(chipID=chipID):
-                continue
 
             # Make Chip & Filter lists
             chip = Chip(
@@ -48,8 +47,10 @@ class Observation(object):
             filt = Filter(filter_id=filter_id, 
                 filter_type=filter_type, 
                 filter_param=self.filter_param)
+            if not self.focal_plane.isIgnored(chipID=chipID):
                 self.chip_list.append(chip)
                 self.filter_list.append(filt)
+            self.all_filter.append(filt)
 
         # Read catalog and shear(s)
         self.g1_field, self.g2_field, self.nshear = get_shear_field(config=self.config)
@@ -157,7 +158,8 @@ class Observation(object):
 
         if pointing.pointing_type == 'MS':
             # Load catalogues and templates
-            self.cat = self.Catalog(config=self.config, chip=chip, pointing=pointing, cat_dir=cat_dir, sed_dir=sed_dir, logger=chip_output.logger)
+            self.cat = self.Catalog(config=self.config, chip=chip, pointing=pointing, cat_dir=cat_dir, sed_dir=sed_dir, chip_output=chip_output)
+            chip_output.create_output_file()
             self.nobj = len(self.cat.objs)
 
             # Loop over objects
@@ -188,8 +190,8 @@ class Observation(object):
                         target_filt=filt, 
                         norm_filt=norm_filt,
                     )
-                    for ifilt in range(len(self.filter_list)):
-                        temp_filter = self.filter_list[ifilt]
+                    for ifilt in range(len(self.all_filter)):
+                        temp_filter = self.all_filter[ifilt]
                         _, obj.param["mag_%s"%temp_filter.filter_type] = self.cat.convert_sed(
                             mag=obj.param["mag_use_normal"],
                             sed=sed_data,

config/config_C3.yaml

@@ -84,6 +84,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
 ###############################################

config/config_NGP.yaml

@@ -73,13 +73,13 @@ 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, 1]
+  run_pointings: [0]
   
   # Run specific chip(s):
   # - give a list of indexes of chips: [ip_1, ip_2...]
   # - run all chips: null
   # Note: for all pointings
-  run_chips: null
+  run_chips: [24, 19, 14, 9]
 
   # Whether to enable astrometric modeling
   # astrometric_lib: "libshao.so"

config/config_example.yaml

@@ -80,6 +80,9 @@ obs_setting:
   # Note: for all pointings
   run_chips: [18]
 
+  # Cut by saturation/limiting magnitude in which band?
+  cut_in_band: "g"
+
 ###############################################
 # Input path setting
 ###############################################