backup

.gitignore

+*.fits
+*.cat
+*.log
+*.list
+*.png
+*.pyc
+*.so
\ No newline at end of file

config/config_injection_20231203.yaml

+---
+###############################################
+#
+#  Configuration file for CSST object injection
+#      Last modified: 2022/06/19
+#
+###############################################
+# n_objects: 500
+rotate_objs: NO
+use_mpi: YES
+run_name: "test_20231203"
+project_cycle: 6
+run_counter: 1
+
+pos_sampling:
+  type: "HexGrid"
+  # type: "RectGrid"
+  # grid_spacing: 18.5 # arcsec (~250 pixels)
+  grid_spacing: 15 # arcsec (~500 pixels)
+  # type: "uniform"
+  # object_density: 37 # arcmin^-2
+
+output_img_dir: "/share/home/fangyuedong/injection_pipeline/workspace"
+input_img_list: "/share/home/fangyuedong/injection_pipeline/input_L1_IMG_20231203.list"
+
+###############################################
+# PSF setting
+###############################################
+psf_setting:
+  
+  # Which PSF model to use:
+  # "Gauss": simple gaussian profile
+  # "Interp": Interpolated PSF from sampled ray-tracing data
+  psf_model: "Interp"
+
+  # PSF size [arcseconds]
+  # radius of 80% energy encircled
+  # NOTE: only valid for "Gauss" PSF
+  psf_rcont: 0.15
+
+  # path to PSF data
+  # NOTE: only valid for "Interp" PSF
+  psf_dir: "/share/simudata/CSSOSDataProductsSims/data/psfCube1"
+
+###############################################
+# Input path setting
+# (NOTE) Used NGP Catalog for testing
+###############################################
+# Default path settings for NGP footprint simulation
+data_dir: "/share/simudata/CSSOSDataProductsSims/data/"
+
+input_path:
+  cat_dir: "OnOrbitCalibration/CTargets20211231"
+  star_cat: "CT-NGP_r1.8_G28.hdf5"
+  galaxy_cat: "galaxyCats_r_3.0_healpix_shift_192.859500_27.128300.hdf5"
+
+SED_templates_path:
+  star_SED: "Catalog_20210126/SpecLib.hdf5"
+  galaxy_SED: "Templates/Galaxy/"
+
+
+###############################################
+# Instrumental effects setting
+# (NOTE) Here only used to construct
+# ObservationSim.Instrument.Chip object
+# (TODO) Should readout from header
+###############################################
+ins_effects:
+  # switches
+  bright_fatter:  ON # Whether to simulate Brighter-Fatter (also diffusion) effect
+
+  # values
+  # dark_exptime:   300   # Exposure time for dark current frames [seconds]
+  # flat_exptime:   150   # Exposure time for flat-fielding frames [seconds]
+  # readout_time:   40    # The read-out time for each channel [seconds]
+  # df_strength:    2.3   # Sillicon sensor diffusion strength
+  # bias_level:     500   # bias level [e-/pixel]
+  # gain:           1.1   # Gain
+  # full_well:      90000 # Full well depth [e-]
+
+###############################################
+# Random seeds
+###############################################
+random_seeds:
+  seed_Av:              121212    # Seed for generating random intrinsic extinction
+
+###############################################
+# Measurement setting
+###############################################
+measurement_setting:
+  input_img_list: "/share/home/fangyuedong/injection_pipeline/L1_INJECTED_20231203.list"
+  # input_img_list: "/share/home/fangyuedong/injection_pipeline/input_L1_img_MSC_0000000.list"
+  input_wht_list: "/share/home/fangyuedong/injection_pipeline/L1_WHT_20231203.list"
+  input_flg_list: "/share/home/fangyuedong/injection_pipeline/L1_FLG_20231203.list"
+  # input_psf_list: "/share/home/fangyuedong/injection_pipeline/psf_img_MSC_0000000.list"
+  sex_config: "/share/home/fangyuedong/injection_pipeline/config/default.config"
+  sex_param: "/share/home/fangyuedong/injection_pipeline/config/default.param"
+  n_jobs: 18
+  output_dir: "/share/home/fangyuedong/injection_pipeline/workspace"
\ No newline at end of file

evaluation/calculate_completeness_fraction.py

@@ -8,16 +8,30 @@ from ObservationSim.Instrument import Telescope, Filter, FilterParam
 
 VC_A = 2.99792458e+18  # speed of light: A/s
 
+# def define_options():
+#     parser = argparse.ArgumentParser()
+#     parser.add_argument('--TU_catalog', dest='TU_catalog', type=str, required=True,
+#                         help='path to the (injected) truth catalog')
+#     parser.add_argument('--source_catalog', dest='source_catalog', type=str, required=True,
+#                         help='path to the (extracted) injected catalog')
+#     parser.add_argument('--orig_catalog', dest='orig_catalog', type=str, required=True,
+#                         help='path to the (extracted) original catalog')
+#     parser.add_argument('--image', dest='image', type=str, required=True,
+#                         help='path to the image, used to get the header info')
+#     parser.add_argument('--output_dir', dest='output_dir', type=str, required=False,
+#                         default='./workspace', help='output path')
+#     return parser
+
 def define_options():
     parser = argparse.ArgumentParser()
-    parser.add_argument('--TU_catalog', dest='TU_catalog', type=str, required=True,
-                        help='path to the (injected) truth catalog')
-    parser.add_argument('--source_catalog', dest='source_catalog', type=str, required=True,
+    parser.add_argument('--TU_catalog_list', dest='TU_catalog_list', type=str, required=True,
+                        help='path to the list of (injected) truth catalog')
+    parser.add_argument('--source_catalog_list', dest='source_catalog_list', type=str, required=True,
                         help='path to the (extracted) injected catalog')
-    parser.add_argument('--orig_catalog', dest='orig_catalog', type=str, required=True,
-                        help='path to the (extracted) original catalog')
-    parser.add_argument('--image', dest='image', type=str, required=True,
-                        help='path to the image, used to get the header info')
+    # parser.add_argument('--orig_catalog', dest='orig_catalog', type=str, required=True,
+    #                     help='path to the list of (extracted) original catalog')
+    # parser.add_argument('--image', dest='image', type=str, required=True,
+    #                     help='path to the image, used to get the header info')
     parser.add_argument('--output_dir', dest='output_dir', type=str, required=False,
                         default='./workspace', help='output path')
     return parser
@@ -72,13 +86,24 @@ def convert_catalog(catname):
     fits_filename = os.path.join(data_dir, base_name + '.fits')
     text_file.write(fits_filename, overwrite=True)
 
-def validation_hist(val, idx, name="val", nbins=10, fig_name='detected_counts.png', output_dir='./'):
+def validation_hist(val, idx, name="val", nbins=10, bins=None, fig_name='detected_counts.png', output_dir='./', create_figure=True):
+    if bins is None:
         counts, bins = np.histogram(val, bins=nbins)
+    else:
+        counts, bins = np.histogram(val, bins=bins)
     is_empty = np.full(len(val), False)
     for i in range(len(idx)):
         if idx[i].size == 0:
             is_empty[i] = True
+    if bins is None:
         counts_detected, _ = np.histogram(val[~is_empty], bins=nbins)
+    else:
+        counts_detected, _ = np.histogram(val[~is_empty], bins=bins)
+    if create_figure:
+        create_hist_figure(counts, counts_detected, bins, name, output_dir, fig_name)
+    return counts, counts_detected, bins
+
+def create_hist_figure(counts, counts_detected, bins, name="val", output_dir='./', fig_name='detected_counts.png'):
     plt.figure()
     plt.stairs(counts, bins, color='r', label='TU objects')
     plt.stairs(counts_detected, bins, color='g', label='Detected')
@@ -87,22 +112,38 @@ def validation_hist(val, idx, name="val", nbins=10, fig_name='detected_counts.pn
     plt.legend(loc='upper right', fancybox=True)
     fig_name = os.path.join(output_dir, fig_name)
     plt.savefig(fig_name)
-    return counts, bins
 
-def hist_fraction(val, idx, name='val', nbins=10, normed=False, output_dir='./'):
+def hist_fraction(val, idx, name='val', nbins=10, bins=None, output_dir='./', fig_name="completeness_fraction.png"):
+    if bins is None:
         counts, bins = np.histogram(val, bins=nbins)
+    else:
+        counts, bins = np.histogram(val, bins=bins)
     is_empty = np.full(len(val), False)
     for i in range(len(idx)):
         if idx[i].size == 0:
             is_empty[i] = True
-    counts_detected, _ = np.histogram(val[~is_empty], bins=nbins, density=normed)
+    if bins is None:
+        counts_detected, _ = np.histogram(val[~is_empty], bins=nbins)
+    else:
+        counts_detected, _ = np.histogram(val[~is_empty], bins=bins)
     fraction = counts_detected / counts
     fraction[np.where(np.isnan(fraction))[0]] = 0.
     plt.figure()
     plt.stairs(fraction, bins, color='r', label='completeness fraction')
     plt.xlabel(name, size='x-large')
     plt.title("Completeness Fraction")
-    fig_name = os.path.join(output_dir, "completeness_fraction_%s.png"%(name))
+    fig_name = os.path.join(output_dir, fig_name)
+    plt.savefig(fig_name)
+    return fraction
+
+def create_fraction_figure(counts, counts_detected, bins, name='val', output_dir='./', fig_name="completeness_fraction.png"):
+    fraction = counts_detected / counts
+    fraction[np.where(np.isnan(fraction))[0]] = 0.
+    plt.figure()
+    plt.stairs(fraction, bins, color='r', label='completeness fraction')
+    plt.xlabel(name, size='x-large')
+    plt.title("Completeness Fraction")
+    fig_name = os.path.join(output_dir, fig_name)
     plt.savefig(fig_name)
     return fraction
 
@@ -116,6 +157,25 @@ def calculate_fraction(TU_catalog, source_catalog, output_dir, nbins=10):
     fraction = hist_fraction(val=mag_TU, idx=idx1, name="mag_injected", nbins=10, output_dir=output_dir)
     return counts, bins, fraction
 
+def calculate_fraction_multi_cats(TU_catalog_list, source_catalog_list, output_dir, nbins=10):
+    counts = np.zeros(nbins)
+    counts_detected = np.zeros(nbins)
+    bins = np.linspace(18, 26, num=(nbins+1))
+    for i in range(len(TU_catalog_list)):
+        TU_catalog = TU_catalog_list[i]
+        source_catalog = source_catalog_list[i]
+        convert_catalog(TU_catalog)
+        x_TU_temp, y_TU_temp, col_list = read_catalog(TU_catalog + '.fits', ext_num=1, ra_name="xImage", dec_name="yImage", col_list=["mag"])
+        mag_TU_temp = col_list[0]
+        x_source_temp, y_source_temp, _ = read_catalog(source_catalog, ext_num=1, ra_name="X_IMAGE", dec_name="Y_IMAGE")
+        idx1, _, = match_catalogs_img(x1=x_TU_temp, y1=y_TU_temp, x2=x_source_temp, y2=y_source_temp)
+        counts_temp, counts_detected_temp, _ = validation_hist(val=mag_TU_temp, idx=idx1, name="mag_injected", bins=bins, output_dir=output_dir, create_figure=False)
+        counts += counts_temp
+        counts_detected += counts_detected_temp
+    create_hist_figure(counts, counts_detected, bins, "mag_injected", output_dir)
+    fraction = create_fraction_figure(counts, counts_detected, bins, 'mag_injected', output_dir)
+    return counts, counts_detected, bins, fraction
+
 def calculate_undetected_flux(orig_cat, mag_bins, fraction, mag_low=20.0, mag_high=26.0, image=None,  output_dir='./'):
     # Get info from original image
     hdu = fits.open(image)
@@ -166,19 +226,32 @@ def calculate_undetected_flux(orig_cat, mag_bins, fraction, mag_low=20.0, mag_hi
     undetected_flux /= (float(nx_pix) * float(ny_pix))
     return undetected_flux
 
+# if __name__ == "__main__":
+#     args = define_options().parse_args()
+#     counts, bins, fraction = calculate_fraction(
+#         TU_catalog=args.TU_catalog,
+#         source_catalog=args.source_catalog,
+#         output_dir=args.output_dir,
+#         nbins=20
+#     )
+#     undetected_flux = calculate_undetected_flux(
+#         orig_cat=args.orig_catalog,
+#         mag_bins=bins,
+#         fraction=fraction,
+#         image=args.image,
+#         output_dir=args.output_dir,
+#     )
+#     print(undetected_flux)
+
 if __name__ == "__main__":
     args = define_options().parse_args()
-    counts, bins, fraction = calculate_fraction(
-        TU_catalog=args.TU_catalog,
-        source_catalog=args.source_catalog,
+    with open(args.TU_catalog_list) as file:
+        TU_catalog_list = [line.rstrip() for line in file]
+    with open(args.source_catalog_list) as file:
+        source_catalog_list = [line.rstrip() for line in file]
+    counts, counts_detected, bins, fraction = calculate_fraction_multi_cats(
+        TU_catalog_list=TU_catalog_list,
+        source_catalog_list=source_catalog_list,
         output_dir=args.output_dir,
         nbins=20
     )
\ No newline at end of file
-    undetected_flux = calculate_undetected_flux(
-        orig_cat=args.orig_catalog,
-        mag_bins=bins,
-        fraction=fraction,
-        image=args.image,
-        output_dir=args.output_dir,
-    )
-    print(undetected_flux)
\ No newline at end of file

evaluation/cross_match_catalogs.py

 import numpy as np
 import astropy.units as u
-import matplotlib.pyplot as plt
 from astropy.coordinates import SkyCoord
 from astropy.io import fits
 from astropy.io import ascii
 from sklearn.neighbors import BallTree
 
-# TU_catalog = "test_RectGrid_20220628.cat"
-# source_catalog = "extracted_test_RectGrid_20220628.fits"
 TU_catalog = "injected_bkgsub_img.cat"
 source_catalog = "extracted_injected_bkgsub_img.fits"
 
@@ -26,7 +23,6 @@ def read_catalog(catname, ext_num=1, ra_name='ra', dec_name='dec', col_list=[]):
     if len(col_list) > 0:
         for col in col_list:
             col_other.append(data[col])
-    # print(ra, dec)
     return ra, dec, col_other
 
 def match_catalogs_sky(ra1, dec1, ra2, dec2, max_dist=0.6, others1=[], others2=[], thresh=[]):
@@ -37,33 +33,24 @@ def match_catalogs_sky(ra1, dec1, ra2, dec2, max_dist=0.6, others1=[], others2=[
     # print(idx2)
     # print(np.shape(idx1))
     # print(np.shape(idx2))
+    # TODO
 
-def match_catalogs_img(x1, y1, x2, y2, max_dist=0.5, others1=[], others2=[], thresh=[]):
+def match_catalogs_img(x1, y1, x2, y2, max_dist=2, others1=[], others2=[], thresh=[]):
     cat1 = np.array([(x, y) for x,y in zip(x1, y1)])
     cat2 = np.array([(x, y) for x,y in zip(x2, y2)])
-    # print(np.shape(cat1))
-    # print(np.shape(cat2))
     tree = BallTree(cat2)
     idx1 = tree.query_radius(cat1, r = max_dist)
     tree = BallTree(cat1)
     idx2 = tree.query_radius(cat2, r = max_dist)
-    # print(np.shape(idx1))
     tot = 0
-    print(idx1)
     for idx in idx1:
         if len(idx) == 0:
             continue
+        if len(idx) > 1: print(len(idx))
         tot += 1
-    print(tot)
+    print("number of matched sources = ", tot)
     return idx1, idx2
 
-def validation_hist(val1, idx1, name="val1", nbins=10):
-    counts, bins = np.histogram(val1)
-    plt.stairs(counts, bins)
-    plt.xlabel(name, size='x-large')
-    plt.savefig("detection_completeness.png")
-    # plt.show()
-
 if __name__=="__main__":
     convert_catalog(TU_catalog)
     ra_TU, dec_TU, _ = read_catalog('test_ascii_to_fits.fits', ext_num=1, ra_name="ra", dec_name="dec")
@@ -74,4 +61,3 @@ if __name__=="__main__":
     # match_catalogs_sky(ra1=ra_TU, dec1=dec_TU, ra2=ra_source, dec2=dec_source)
     idx1, idx2, = match_catalogs_img(x1=x_TU, y1=y_TU, x2=x_source, y2=y_source)
     # print(ra_TU, dec_TU)
\ No newline at end of file
-    validation_hist(mag_TU, idx1, name="mag_injected")
\ No newline at end of file

injection_pipeline/Catalog/C6_SimCat.py

+import os
+import galsim
+import random
+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
+from datetime import datetime
+
+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:
+    # Try backported to PY<37 'importlib_resources'
+    import importlib_resources as pkg_resources
+
+# CONSTANTS
+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 SimCat(CatalogBase):
+    def __init__(self, config, chip, nobjects=None):
+        super().__init__()
+        self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
+        self.seed_Av = config["catalog_options"]["seed_Av"]
+        self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111)
+
+        with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path:
+            self.normF_star = Table.read(str(filter_path))
+
+        self.config = config
+        self.chip = chip
+
+        galaxy_dir = config["catalog_options"]["input_path"]["galaxy_cat"]
+        self.galaxy_path = os.path.join(self.cat_dir, galaxy_dir)
+        self.galaxy_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["galaxy_SED"])
+        self._load_SED_lib_gals()
+
+        if "rotateEll" in config["catalog_options"]:
+            self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.))
+        else:
+            self.rotation = 0.
+        
+        self._get_healpix_list()
+        self._load(nobjects=nobjects)
+
+    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
+        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]))
+        self.pix_list = hp.query_polygon(
+            NSIDE,
+            hp.ang2vec(np.radians(90.) - dec, ra),
+            inclusive=True
+        )
+        if self.logger is not None:
+            msg = str(("HEALPix List: ", self.pix_list))
+            self.logger.info(msg)
+        else:
+            print("HEALPix List: ", self.pix_list)
+
+    def load_norm_filt(self, obj):
+        if obj.type == "star":
+            return self.normF_star
+        elif obj.type == "galaxy" or obj.type == "quasar":
+            return None
+        else:
+            return None
+
+    def _load_SED_lib_gals(self):
+        pcs = h5.File(os.path.join(self.galaxy_SED_path, "pcs.h5"), "r")
+        lamb = h5.File(os.path.join(self.galaxy_SED_path, "lamb.h5"), "r")
+        self.lamb_gal = lamb['lamb'][()]
+        self.pcs = pcs['pcs'][()]
+
+    def _load_gals(self, gals, pix_id=None, cat_id=0, nobjects=1):
+        # Load how mnay objects?
+        max_ngals = len(gals['ra'])
+        remain = nobjects 
+        for igals in range(max_ngals):
+            if remain == 0:
+                break
+            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]
+            # [TODO]
+            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
+
+            param['z'] = gals['redshift'][igals]
+            param['model_tag'] = 'None'
+            param['g1'] = gals['shear'][igals][0]
+            param['g2'] = gals['shear'][igals][1]
+            param['kappa'] = gals['kappa'][igals]
+            param['e1'] = gals['ellipticity_true'][igals][0]
+            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
+            
+            remain -= 1
+            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
+
+            # Masses
+            param['bulgemass'] = gals['bulgemass'][igals]
+            param['diskmass'] = gals['diskmass'][igals]
+
+            param['size'] = gals['size'][igals]
+            if param['size'] > self.max_size:
+                self.max_size = param['size']
+
+            # Sersic index
+            param['disk_sersic_idx'] = 1.
+            param['bulge_sersic_idx'] = 4.
+
+            # Sizes
+            param['bfrac'] = param['bulgemass']/(param['bulgemass'] + param['diskmass'])
+            if param['bfrac'] >= 0.6:
+                param['hlr_bulge'] = param['size']
+                param['hlr_disk'] = param['size'] * (1. - param['bfrac'])
+            else:
+                param['hlr_disk'] = param['size']
+                param['hlr_bulge'] = param['size'] * param['bfrac']
+
+            # SED coefficients
+            param['coeff'] = gals['coeff'][igals]
+            param['detA'] = gals['detA'][igals]
+
+            # Others
+            param['galType'] = gals['type'][igals]
+            param['veldisp'] = gals['veldisp'][igals]
+            
+            # TEST no redening and no extinction
+            param['av'] = 0.0
+            param['redden'] = 0
+
+            param['star'] = 0   # Galaxy
+
+            # TEMP
+            self.ids += 1
+            # param['id'] = self.ids
+            param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
+            
+            if param['star'] == 0:
+                obj = Galaxy(param, self.rotation, logger=self.logger)
+
+            self.objs.append(obj)
+        return remain
+    
+    def _load(self, nobjects=1):
+        from itertools import cycle
+        self.objs = []
+        self.ids = 0
+        to_be_read_in = nobjects
+        pool = cycle(self.pix_list)
+        for pix in pool:
+            try:
+                if to_be_read_in == 0:
+                    break
+                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)]
+                to_be_read_in = self._load_gals(gals, pix_id=pix, cat_id=bundleID, n_objects=to_be_read_in)
+                del gals
+            except Exception as e:
+                traceback.print_exc()
+                print(e)
+
+    def load_sed(self, obj, **kwargs):
+        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)
+        flux[flux < 0] = 0.
+        sedcat = np.vstack((self.lamb_gal, flux)).T
+        sed_data = getObservedSED(
+            sedCat=sedcat,
+            redshift=obj.z,
+            av=obj.param["av"],
+            redden=obj.param["redden"]
+        )
+        wave, flux = sed_data[0], sed_data[1]
+        speci = interpolate.interp1d(wave, flux)
+        lamb = np.arange(2000, 11001+0.5, 0.5)
+        y = speci(lamb)
+        # erg/s/cm2/A --> photon/s/m2/A
+        all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
+        sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
+        
+        del wave
+        del flux
+        return sed
\ No newline at end of file

injection_pipeline/SingleEpochImage.py

@@ -57,23 +57,6 @@ class SingleEpochImage(object):
         
         self.setup_image_for_injection()
 
-        # [TODO] Header should come from input image
-        # self.header_ext = generateExtensionHeader(
-        #     chip=self.chip,
-        #     xlen=self.chip.npix_x,
-        #     ylen=self.chip.npix_y,
-        #     ra=self.ra_cen,
-        #     dec=self.dec_cen,
-        #     pa=self.pos_ang,
-        #     gain=self.chip.gain, # {TODO}
-        #     readout=self.chip.read_noise,
-        #     dark=self.chip.dark_noise,
-        #     saturation=90000,
-        #     pixel_scale=self.chip.pix_scale, 
-        #     pixel_size=self.chip.pix_size,
-        #     xcen=self.chip.x_cen,
-        #     ycen=self.chip.y_cen,
-        #     extName = 'raw')
         self.header_ext = self.header_img
 
         # (TODO) specify sub-directory
@@ -95,10 +78,12 @@ class SingleEpochImage(object):
     def setup_image_for_injection(self):
         self.ra_cen = self.wcs.wcs.crval[0]
         self.dec_cen = self.wcs.wcs.crval[1]
-        self.wcs_fp = self.focal_plane.getTanWCS(self.ra_cen, self.dec_cen, self.pos_ang*galsim.degrees, self.pixel_scale)
+        # self.wcs_fp = self.focal_plane.getTanWCS(self.ra_cen, self.dec_cen, self.pos_ang*galsim.degrees, self.pixel_scale)
         self.chip.img = galsim.Image(self.image, copy=True)
         self.chip.img.setOrigin(self.chip.bound.xmin, self.chip.bound.ymin)
-        self.chip.img.wcs = self.wcs_fp
+        # self.chip.img.wcs = self.wcs_fp
+        # self.chip.img.setOrigin(0, 0)
+        self.chip.img.wcs = galsim.AstropyWCS(wcs=self.wcs)
 
     def read_initial_image(self, filepath):
         data_dir = os.path.dirname(filepath)
@@ -110,15 +95,20 @@ class SingleEpochImage(object):
         hdu.close()
         
         # Determine which CCD
-        self.chip_ID = int(self.header0['DETECTOR'][-2:])
+        # self.chip_ID = int(self.header0['DETECTOR'][-2:])
+        self.chip_ID = int(self.header_img["CHIPID"])
         # Construnct Chip object
         self.chip = Chip(chipID=self.chip_ID, config=self.config)
         self.exp_time = float(self.header0 ['EXPTIME'])
-        self.chip.gain = float(self.header_img["GAIN1"])
+        # self.chip.gain = float(self.header_img["GAIN1"])
+        self.chip.gain = float(self.header_img["GAIN01"])
 
         # Process L1 image
-        # self.image = self.image * self.exp_time / self.chip.gain
-        self.image = self.image * self.chip.gain
+        self.image = self.image * self.exp_time / self.chip.gain
+        print(self.chip.gain, self.exp_time)
+        print(self.image)
+        print(np.sum(self.image < 0), np.sum(self.image > 0))
+        # self.image = self.image * self.chip.gain
 
         # Process L1 SKY image
         # [TODO]
@@ -128,7 +118,8 @@ class SingleEpochImage(object):
         # image -= sky_img
 
         # Process L1 FLAG image
-        flag_img_path = os.path.join(data_dir, img_name.replace("img_L1", "flg_L1"))
+        # flag_img_path = os.path.join(data_dir, img_name.replace("img_L1", "flg_L1"))
+        flag_img_path = os.path.join(data_dir, img_name.replace("IMG", "FLG"))
         flag_img = fits.getdata(flag_img_path)
         self.image[flag_img > 0] = 0.
 
@@ -143,7 +134,7 @@ class SingleEpochImage(object):
         # hdu1.writeto(fname, output_verify='ignore', overwrite=True)
 
     def _get_wcs(self, header):
-        self.pos_ang = float(header['POS_ANG'])
+        # self.pos_ang = float(header['POS_ANG'])
         self.wcs = wcs.WCS(header)
 
         self.pixel_scale = 0.074
@@ -186,7 +177,7 @@ class SingleEpochImage(object):
             try:
                 sed_data = cat.load_sed(obj)
                 norm_filt = cat.load_norm_filt(obj)
-                obj.sed, obj.param["mag_%s"%self.filt.filter_type], obj.param["flux_%s"%self.filt.filter_type] = cat.convert_sed(
+                obj.sed, obj.param["mag_%s"%self.filt.filter_type.lower()], obj.param["flux_%s"%self.filt.filter_type.lower()] = cat.convert_sed(
                         mag=obj.param["mag_use_normal"],
                         sed=sed_data,
                         target_filt=self.filt, 
@@ -233,9 +224,9 @@ class SingleEpochImage(object):
         self.chip.img = galsim.Image(self.chip.img.array, dtype=np.float32)
         
         # [TODO] TEST
-        # self.chip.img *= self.chip.gain
-        # self.chip.img /= self.exp_time
-        self.chip.img /= self.chip.gain
+        self.chip.img *= self.chip.gain
+        self.chip.img /= self.exp_time
+        # self.chip.img /= self.chip.gain
 
         hdu1 = fits.PrimaryHDU(header=self.header0)
         hdu2 = fits.ImageHDU(self.chip.img.array, header=self.header_img)

measurement_pipeline/L1_pipeline/csst_msc_mbi_instrument.py

+from datetime import datetime
+from typing import Optional
+
+import numpy as np
+from itertools import product
+from astropy.io import fits
+from astropy.stats import sigma_clipped_stats
+
+SATURATE = 50000
+
+def core_msc_l1_mbi_instrument(
+    image_path: str = "/path/to/image",
+    bias_path: str = "/path/to/bias",
+    dark_path: str = "/path/to/dark",
+    flat_path: str = "/path/to/flat",
+    shutter_path: str = "/path/to/shutter",
+    image_output_path: str = "/path/to/image_output",
+    weight_output_path: str = "/path/to/weight_output",
+    flag_output_path: str = "/path/to/flag_output",
+    deepcr_model_path: Optional[str] = None,
+    config_ccd_info: Optional[str] = None,
+    config_bad_pixel: Optional[str] = None,
+    device: Optional[str] = "CPU",
+):
+    """
+    Make the instrument correction for one chip of CSST data.
+
+    This function corrects the instrument effect, such as bias, dark, flat, shutter, cosmicrays etc.
+
+    Parameters
+    ----------
+    image_path : str
+        The sci file path.
+    bias_path : str
+        The bias file path.
+    dark_path : str
+        The dark file path.
+    flat_path : str
+        The flat file path.
+    shutter_path : str
+        The shutter file path.
+    image_output_path : str
+        The image output file path.
+    weight_output_path : str
+        The weight output file path.
+    flag_output_path : str
+        The flag output file path.
+    deepcr_model_path : str, optional
+        The deepcr model file path.
+    config_ccd_info : str, optional
+        The config for ccd file path.
+    config_bad_pixel : str, optional
+        The config for bad pixel file path.
+    device : str
+        The deepcr device input 'CPU' or 'GPU', by default 'CPU'.
+
+    Returns
+    -------
+    CsstResult
+        Result containing `status` and `files`.
+    """
+    status_header = _set_status_header()
+
+    # CCD image processing
+    with fits.open(image_path) as raw:
+        hdu = fits.HDUList([
+            fits.PrimaryHDU(
+                header=raw[0].header.copy()),
+            fits.ImageHDU(
+                header=raw[1].header.copy(),
+                data=raw[1].data.copy())
+        ])
+        # 创建wht[1].data hdulist
+        wht = fits.HDUList([
+            fits.PrimaryHDU(
+                header=raw[0].header.copy()),
+            fits.ImageHDU(
+                header=raw[1].header.copy(),
+                data=np.zeros_like(hdu[1].data, dtype=np.float32))
+        ])
+        # 创建flg[1].data hdulist
+        flg = fits.HDUList([
+            fits.PrimaryHDU(
+                header=raw[0].header.copy()),
+            fits.ImageHDU(
+                header=raw[1].header.copy(),
+                data=np.zeros_like(hdu[1].data, dtype=np.uint32))])
+
+        img_err = np.zeros_like(hdu[1].data)
+        # shutter
+        shutter = fits.getdata(shutter_path, 1)
+        bias, bias_error = _read_ccds(bias_path)
+        dark, dark_error = _read_ccds(dark_path)
+        flat, flat_error = _read_ccds(flat_path)
+
+        hdu[1].data, img_err = subtract_bias(
+            hdu[1].data,
+            img_err,
+            bias,
+            bias_error,
+        )
+        status_header.set('STA_BIAS', 0)
+
+        hdu[1].data, img_err = subtract_dark(
+            hdu[1].data,
+            img_err,
+            dark,
+            dark_error,
+            hdu[0].header["EXPTIME"]
+        )
+        status_header.set('STA_DARK', 0)
+
+        hdu[1].data, img_err = _divide_flat(
+            hdu[1].data,
+            img_err,
+            flat,
+            flat_error,
+        )
+        status_header.set('STA_FLAT', 0)
+
+        # 非线性改正
+        hdu[1].data = _fix_nonlinear(hdu[1].data)
+        status_header.set('STA_NLIN', 0)
+        # 坏像素检测
+        flg[1].data = _check_badpix(
+            flg[1].data,
+            flat
+        )
+
+        # 热像素和暖像素检测
+        flg[1].data = _check_hot_and_warm_pix(
+            flg[1].data,
+            dark,
+            hdu[0].header["EXPTIME"],
+            hdu[1].header["RON01"]
+        )
+        # 过饱和检测
+        saturate = SATURATE
+        flg[1].data = _check_over_saturation(
+            flg[1].data,
+            hdu[1].data,
+            saturate
+        )
+        saturate = (saturate
+                    * hdu[1].header["GAIN01"]
+                    / hdu[0].header["EXPTIME"])
+        status_header.set('SATURATE', saturate)
+
+        # 时间归一
+        hdu[1].data = hdu[1].data/hdu[0].header["EXPTIME"]
+        # 单位转换
+        hdu[1].data = _multiply_gain_map(hdu[1].data, hdu[1].header)
+        status_header.set('BUNIT', 'e/s')
+        
+        # # 检测宇宙线
+        # hdu[1].data, flg[1].data, cr_count = remove_cr_deepcr(
+        #     hdu[1].data,
+        #     flg[1].data,
+        #     device,
+        #     hdu[1].header['CHIPID'],
+        # )
+        # status_header.set('STA_CRS', 0)
+        # status_header.set('CRCOUNT', cr_count)
+
+        # 记录BKG与RMS
+        _, bkg, rms = sigma_clipped_stats(
+            data=hdu[1].data,
+            mask=(flg[1].data == 16)
+        )
+        status_header.set('SKY_BKG', bkg)
+        status_header.set('SKY_RMS', rms)
+        # 检测卫星拖尾
+        flg[1].data = _check_satellitetrail(
+            flg[1].data,
+            hdu[1].data
+        )
+        status_header.set('STA_SAT', 0)
+
+        # 生成wht[1].data
+        wht[1].data = _create_weight(
+            hdu[1].data,
+            flg[1].data,
+            hdu[1].header["RON01"],
+            hdu[0].header["EXPTIME"]
+        )
+        # 数据类型
+        hdu[1].data = hdu[1].data.astype(np.float32)
+        
+        # 快门改正
+        hdu[1].data = _fix_shutter(
+            hdu[1].data,
+            shutter,
+            hdu[0].header["EXPTIME"]
+        )
+        status_header.set('STA_SHUT', 0)
+
+        if all([status_header["STA_BIAS"] == 0,
+                status_header["STA_DARK"] == 0,
+                status_header["STA_FLAT"] == 0,
+                status_header["STA_CRS"] == 0,
+                status_header["STA_NLIN"] == 0,
+                # status_header["STA_CTE"] == 0,  # 这个功能还未添加
+                status_header["STA_SAT"] == 0,
+                status_header["STA_SHUT"] == 0,
+                status_header["SKY_BKG"] != -9999,
+                status_header["SKY_RMS"] != -9999,
+                status_header["SATURATE"] != -9999,
+                status_header["CRCOUNT"] != -9999]):
+            status_header.set('STA_INST', 0)
+        
+        # hdu[1].header.extend(status_header, bottom=True)
+        # 修正wht和flg数据类型
+        wht[1].header = hdu[1].header.copy()
+        wht[1].data = wht[1].data.astype(np.float32)
+        wht[1].header.remove('BUNIT')
+
+        flg[1].header = hdu[1].header.copy()
+        flg[1].data = flg[1].data.astype(np.int32)
+        flg[1].header.remove('BUNIT')
+
+        hdu.writeto(image_output_path, overwrite=True)
+        wht.writeto(weight_output_path, overwrite=True)
+        flg.writeto(flag_output_path, overwrite=True)
+        
+
+def _set_status_header():
+    status_header = fits.Header()  # template Header
+    status_header.append(
+        ('COMMENT', '=' * 66), bottom=True)
+    status_header.append(
+        ('COMMENT', 'instrumental correction information'), bottom=True)
+    status_header.append(
+        ('COMMENT', '=' * 66), bottom=True)
+    # status_header.append(
+    #     ('VER_INST', csst_common.__version__, 'version of instrument'), bottom=True)
+    # status_header.append(
+    #     ('STM_INST', csst_common.now(), 'time stamp of instrument processing'), bottom=True)
+    status_header.append(
+        ('STA_INST', 1, '0=done 1=wrong'), bottom=True)
+    status_header.append(
+        ('STA_BIAS', 1, 'status flag for bias frame correction'), bottom=True)
+    status_header.append(
+        ('STA_DARK', 1, 'status flag for dark frame correction'), bottom=True)
+    status_header.append(
+        ('STA_FLAT', 1, 'status flag for flat frame correction'), bottom=True)
+    status_header.append(
+        ('SKY_BKG', -9999, 'estimated sky background (e-/s per pixel)'), bottom=True)
+    status_header.append(
+        ('SKY_RMS', -9999, 'standard dev of frame background (ADU)-> e-/s'), bottom=True)
+    status_header.append(
+        ('SATURATE', -9999, 'the flux limit of saturated pixel (e-/s)'), bottom=True)
+    status_header.append(
+        ('STA_CTE', 1, 'status flag for CTE correction'), bottom=True)
+    status_header.append(
+        ('STA_SAT', 1, 'status flag for satellite correction'), bottom=True)
+    status_header.append(
+        ('STA_CRS', 1, 'status flag for cosmic rays mask'), bottom=True)
+    status_header.append(
+        ('CRCOUNT', -9999, 'cosmic rays pixel counts'), bottom=True)
+    status_header.append(
+        ('STA_NLIN', 1, 'status flag for non-linear correction'), bottom=True)
+    status_header.append(
+        ('STA_SHUT', 1, 'status flag for shutter effect correction'), bottom=True)
+    return status_header
+
+def _read_ccds(
+    ref_path: str,
+) -> tuple[np.array, np.array]:
+    ref = fits.getdata(ref_path)
+    return ref, np.zeros_like(ref)
+
+def _divide_flat(
+    image: np.ndarray,
+    image_error: np.ndarray,
+    flat: np.ndarray,
+    flat_error: np.ndarray,
+) -> tuple[np.ndarray, np.ndarray]:
+    def divide(a, b):
+        return np.divide(a, b, out=np.zeros_like(a, float), where=(b != 0))
+    result = divide(image, flat)
+    result_error = np.abs(result) * ((divide(image_error, image)) ** 2 +
+                                     (divide(flat_error, flat)) ** 2)**0.5
+    return result, result_error
+
+def _fix_nonlinear(
+    image: np.ndarray,
+    beta1: float = 5e-7
+) -> np.ndarray:
+    # 计算非线性系数
+    beta = 5e-7
+    f1 = [1]
+    fnlin0 = 1+beta
+    fnlin = [fnlin0]
+    for i in range(1000, 65000, 1000):
+        f1.append(i)
+        fnlin.append(1+beta*i)
+    # 插值函数
+    from scipy.interpolate import PchipInterpolator as PCHIP
+    # PCHIP stands for Piecewise Cubic Hermite Interpolating Polynomial
+    interp_func = PCHIP(f1, fnlin)
+    # 非线性系数矩阵
+    imgnlin = interp_func(image)
+    # 图像改正
+    image = imgnlin*image
+    return image
+
+def _check_badpix(
+    flag: np.ndarray,
+    flat: np.ndarray,
+) -> np.ndarray:
+    med = np.median(flat)
+    _flag = (flat < 0.5 * med) | (1.5 * med < flat)
+    flag = flag | (_flag * 1)
+    return flag
+
+def _check_hot_and_warm_pix(
+    flag: np.ndarray,
+    dark: np.ndarray,
+    exptime: float,
+    rdnoise: float,
+) -> np.ndarray:
+    _dark = dark * exptime
+    _dark[_dark < 0] = 0
+    _flag = 1 * rdnoise ** 2 <= _dark  # 不确定是否包含 暂定包含
+    flag = flag | (_flag * 2)
+    _flag = (0.5 * rdnoise ** 2 < _dark) & (_dark < 1 * rdnoise ** 2)
+    flag = flag | (_flag * 4)
+    return flag
+
+def _check_over_saturation(
+    flag: np.ndarray,
+    image: np.ndarray,
+    saturate: int = 65535,
+    iterations: int = 0,
+    flag_value: int = 8
+) -> np.ndarray:
+    _flag = image >= saturate
+    flag_dilated = _dilation(_flag, iterations=iterations)
+    flag = flag | (flag_dilated * flag_value)
+    return flag
+
+def _dilation(array_orig: np.ndarray, iterations: int = 1):
+    """
+    Make a dilation for array
+
+    This function makes a dilation for the 2D mask array (saturated pixels)
+
+    Parameters
+    ----------
+    array_orig : np.ndarray
+        The mask array without dilation.
+    structure : int
+        The number of dilations performed on the structure with itself.
+
+    Returns
+    -------
+    array_out : np.ndarray
+        The mask array after dilation.
+    """
+    from scipy import ndimage
+    struct1 = ndimage.generate_binary_structure(
+        2, 1)  # give a saturate structure
+    struct_ext = ndimage.iterate_structure(
+        struct1, iterations)  # make a dilation
+    if iterations == 0:
+        array_out = array_orig
+    else:
+        array_out = ndimage.binary_dilation(
+            array_orig, structure=struct_ext, iterations=1).astype(array_orig.dtype)
+
+    return array_out
+
+def _check_satellitetrail(
+    flag: np.ndarray,
+    image: np.ndarray,
+) -> np.ndarray:
+    from skimage import exposure, transform
+    from skimage.draw import line
+    from skimage.feature import canny
+    # 调整图像
+    p1, p2 = np.percentile(image, (0.1, 99.9))
+    image = exposure.rescale_intensity(image, in_range=(p1, p2))
+    # 边界识别，转化为二值图像
+    immax = np.max(image)
+    # canny边界
+    edge = canny(image=image,
+                 sigma=2.0,
+                 low_threshold=immax * 0.1,
+                 high_threshold=immax * 0.5)
+
+    # 概率霍夫变换
+    angle = np.radians(np.arange(2, 178, 0.5, dtype=float))
+    # 直线端点
+    result = transform.probabilistic_hough_line(image=edge,
+                                                threshold=210,
+                                                line_length=400,
+                                                line_gap=75,
+                                                theta=angle)
+    result = np.asarray(result)
+    # 绘制mask
+    mask = np.zeros(image.shape)
+    # 遍历每一条线段
+    for (x1, y1), (x2, y2) in result:
+        xx, yy = line(x1, y1, x2, y2)
+        mask[yy, xx] = 1
+    mask = mask.astype(np.int32)
+    flag = flag | (mask * 32)
+    return flag
+
+def _create_weight(
+    img: np.ndarray,
+    flg: np.ndarray,
+    rdnoise: float,
+    exptime: float,
+) -> np.ndarray:
+    data = img.copy()
+    data[img < 0] = 0
+    var_raw = data * exptime + rdnoise ** 2
+    var_bias = 0.0
+    weight = 1. / (var_raw + var_bias) * exptime ** 2
+    weight[flg > 0] = 0
+    wht = weight
+    return wht
+
+def _fix_shutter(
+    img: np.ndarray,
+    shutter: np.ndarray,
+    exptime: float = 150.
+) -> np.ndarray:
+    img_cor = np.float32(img * exptime*1000. / (exptime*1000. + shutter))
+    return img_cor
+
+def _multiply_gain_map(
+        image: np.ndarray,
+        header: fits.header,
+) -> np.ndarray:
+    gain_map = np.zeros_like(image)
+    h = header['NAXIS1'] // 8
+    w = header['NAXIS2'] // 2
+    for y, x in product(range(8), (0, 1)):
+        y = y if x == 0 else 7-y
+        gain = header['GAIN' + str(x*8+y+1).zfill(2)]
+        gain_map[w*x:w*(x+1), h*y:h*(y+1)] = gain
+    return image * gain_map
+
+def subtract_bias(image: np.ndarray,
+                  image_err: np.ndarray,
+                  bias: np.ndarray,
+                  bias_err: np.ndarray,
+                  ) -> tuple[np.ndarray, np.ndarray]:
+    """
+    Subtract bias.
+
+    Subtract bias.
+
+    Parameters
+    ----------
+    image : numpy.ndarray
+        The input image to be corrected.
+    image_err : numpy.ndarray
+        The uncertainty of input image.
+    bias : numpy.ndarray
+        The input bias to be subtracted.
+    bias_err : numpy.ndarray
+        The uncertainty of input bias.
+
+    Returns
+    -------
+    output : np.ndarray
+        Output corrected image.
+    output_err : np.ndarray
+        Output uncertainty map.
+    """
+
+    output = image - bias
+    output_err = np.sqrt(image_err ** 2 + bias_err ** 2)
+
+    return output, output_err
+
+
+def subtract_dark(image: np.ndarray,
+                  image_err: np.ndarray,
+                  dark: np.ndarray,
+                  dark_err: np.ndarray,
+                  tdark_image: float = 1.0,
+                  ) -> tuple[np.ndarray, np.ndarray]:
+    """
+    Subtract dark current.
+
+    Subtract dark current.
+
+    Parameters
+    ----------
+    image : numpy.ndarray
+        The input image to be corrected.
+    image_err : numpy.ndarray
+        The uncertainty of input image.
+    dark : numpy.ndarray
+        The input dark current image to be subtracted.
+    dark_err : numpy.ndarray
+        The uncertainty of input dark current.
+    tdark_image : float, optional
+        The effective dark current cumulation time of input image. Default value is 1.0.
+
+    Returns
+    -------
+    output : np.ndarray
+        Output corrected image.
+    output_err : np.ndarray
+        Output uncertainty map.
+    """
+
+    output = image - dark * tdark_image
+    output_err = np.sqrt(image_err ** 2 + (dark_err * tdark_image) ** 2)
+
+    return output, output_err
+
+

measurement_pipeline/run_measurement.py

@@ -36,7 +36,7 @@ class MeasurementPipeline(object):
             self.wht_list = [line.rstrip() for line in input_list]
         with open(self.config["measurement_setting"]["input_flg_list"]) as input_list:
             self.flg_list = [line.rstrip() for line in input_list]
-        if self.config["measurement_setting"]["input_psf_list"]:
+        if "input_psf_list" in self.config["measurement_setting"] and self.config["measurement_setting"]["input_psf_list"]:
             with open(self.config["measurement_setting"]["input_psf_list"]) as input_list:
                 self.psf_list = [line.rstrip() for line in input_list]
         else:

run_calculate_completeness_fraction.sh

 #!/bin/bash
 
 python /share/home/fangyuedong/injection_pipeline/evaluation/calculate_completeness_fraction.py \
-    --TU_catalog /share/home/fangyuedong/injection_pipeline/workspace/test_20230517/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_08_img_L1_injected.cat \
-    --source_catalog /share/home/fangyuedong/injection_pipeline/workspace/test_20230517/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_08_img_L1_injected_extracted.cat \
-    --orig_catalog /share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_08_cat.fits \
-    --image /share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_08_img_L1.fits \
-    --output_dir /share/home/fangyuedong/injection_pipeline/workspace/test_20230517
\ No newline at end of file
+    --TU_catalog_list /share/home/fangyuedong/injection_pipeline/cat_INJECTED_20231203.list \
+    --source_catalog_list /share/home/fangyuedong/injection_pipeline/cat_EXTRACTED_20231203.list \
+    --output_dir /share/home/fangyuedong/injection_pipeline/workspace/test_20231203
+    # --TU_catalog /share/home/fangyuedong/injection_pipeline/workspace/test_20230517/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_08_img_L1_injected.cat \
+    # --source_catalog /share/home/fangyuedong/injection_pipeline/workspace/test_20230517/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_08_img_L1_injected_extracted.cat \
+    # --orig_catalog /share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_08_cat.fits \
+    # --image /share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_08_img_L1.fits \
+    # --output_dir /share/home/fangyuedong/injection_pipeline/workspace/test_20230517
\ No newline at end of file

run_injection.sh

 #!/bin/bash
 
-python /share/home/fangyuedong/injection_pipeline/injection_pipeline/injection_pipeline.py config_injection.yaml -c /share/home/fangyuedong/injection_pipeline/config
\ No newline at end of file
+# python /share/home/fangyuedong/injection_pipeline/injection_pipeline/injection_pipeline.py config_injection.yaml -c /share/home/fangyuedong/injection_pipeline/config
+python /share/home/fangyuedong/injection_pipeline/injection_pipeline/injection_pipeline.py config_injection_20231203.yaml -c /share/home/fangyuedong/injection_pipeline/config
\ No newline at end of file

run_measurement.sh

 #!/bin/bash
 
-python /share/home/fangyuedong/injection_pipeline/measurement_pipeline/run_measurement.py config_injection.yaml -c /share/home/fangyuedong/injection_pipeline/config
\ No newline at end of file
+# python /share/home/fangyuedong/injection_pipeline/measurement_pipeline/run_measurement.py config_injection.yaml -c /share/home/fangyuedong/injection_pipeline/config
+python /share/home/fangyuedong/injection_pipeline/measurement_pipeline/run_measurement.py config_injection_20231203.yaml -c /share/home/fangyuedong/injection_pipeline/config
\ No newline at end of file

run_injection.pbs → submit_injection.pbs

@@ -10,4 +10,5 @@ cd $PBS_O_WORKDIR
 NP=20
 date
 
-mpirun -np $NP python /share/home/fangyuedong/injection_pipeline/injection_pipeline/injection_pipeline.py config_injection.yaml -c /share/home/fangyuedong/injection_pipeline/config
\ No newline at end of file
+# mpirun -np $NP python /share/home/fangyuedong/injection_pipeline/injection_pipeline/injection_pipeline.py config_injection.yaml -c /share/home/fangyuedong/injection_pipeline/config
+mpirun -np $NP python /share/home/fangyuedong/injection_pipeline/injection_pipeline/injection_pipeline.py config_injection_20231203.yaml -c /share/home/fangyuedong/injection_pipeline/config
\ No newline at end of file