diff --git a/Catalog/C3Catalog.py b/Catalog/C3Catalog.py
new file mode 100644
index 0000000000000000000000000000000000000000..f12436cf1f2cb50d8ac263322cd7be2089eb21a6
--- /dev/null
+++ b/Catalog/C3Catalog.py
@@ -0,0 +1,198 @@
+import os
+import galsim
+import random
+import numpy as np
+import h5py as h5
+import healpy as hp
+import astropy.constants as cons
+from astropy.coordinates import spherical_to_cartesian
+from astropy.table import Table
+from scipy import interpolate
+
+from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar
+from ObservationSim.MockObject._util import seds, sed_assign, extAv, tag_sed, getObservedSED
+
+NSIDE = 128
+
+class C3Catalog(CatalogBase):
+    def __init__(self, config, **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"]
+        self.normalize_dir = os.path.join(config["data_dir"], config["SLS_path"]["SLS_norm"])
+
+        self.normF_star = Table.read(os.path.join(self.normalize_dir, 'SLOAN_SDSS.g.fits'))
+        self.normF_galaxy = Table.read(os.path.join(self.normalize_dir, 'lsst_throuput_g.fits'))
+        
+        try:
+            self.chip = kwargs["chip"]
+        except:
+            raise ValueError('For Cycle-3 Catalog class, must give a "chip" object to initiate')
+
+        if "star_cat" in config["input_path"] and config["input_path"]["star_cat"]:
+            star_file = config["input_path"]["star_cat"]
+            star_SED_file = config["SED_templates_path"]["star_SED"]
+            self.star_path = os.path.join(self.cat_dir, star_file)
+            self.star_SED_path = os.path.join(config["data_dir"], star_SED_file)
+            self._load_SED_lib_star()
+        if "galaxy_cat" in config["input_path"] and config["input_path"]["galaxy_cat"]:
+            galaxy_file = config["input_path"]["galaxy_cat"]
+            self.galaxy_path = os.path.join(self.cat_dir, galaxy_file)
+            self.galaxy_SED_path = os.path.join(config["data_dir"], config["SED_templates_path"]["galaxy_SED"])
+            self._load_SED_lib_gals()
+        if "rotateEll" in config["shear_setting"]:
+            self.rotation = float(int(config["shear_setting"]["rotateEll"]/45.))
+        else:
+            self.rotation = 0.
+
+        self._get_healpix_list()
+        self._load()
+
+    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]))
+        vertices = spherical_to_cartesian(1., dec, ra)
+        self.pix_list = hp.query_polygon(NSIDE, np.array(vertices).T, inclusive=True)
+        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 self.normF_galaxy
+        else:
+            return None
+
+    def _load_SED_lib_star(self):
+        self.tempSED_star = h5.File(self.star_SED_path,'r')
+
+    def _load_SED_lib_gals(self):
+        self.tempSed_gal, self.tempRed_gal = seds("galaxy.list", seddir=self.galaxy_SED_path)
+
+    def _load_gals(self, gals, pix_id=None):
+        ngals = len(gals['galaxyID'])
+        self.rng_sedGal = random.Random()
+        self.rng_sedGal.seed(pix_id) # Use healpix index as the random seed
+        self.ud = galsim.UniformDeviate(pix_id)
+        for igals in range(ngals):
+            param = self.initialize_param()
+            param['ra'] = gals['ra_true'][igals]
+            param['dec'] = gals['dec_true'][igals]
+            param['z'] = gals['redshift_true'][igals]
+            param['model_tag'] = 'None'
+            param['gamma1'] = 0
+            param['gamma2'] = 0
+            param['kappa'] = 0
+            param['delta_ra'] = 0
+            param['delta_dec'] = 0
+            sersicB = gals['sersic_bulge'][igals]
+            hlrMajB = gals['size_bulge_true'][igals]
+            hlrMinB = gals['size_minor_bulge_true'][igals]
+            sersicD = gals['sersic_disk'][igals]
+            hlrMajD = gals['size_disk_true'][igals]
+            hlrMinD = gals['size_minor_disk_true'][igals]
+            aGal = gals['size_true'][igals]
+            bGal = gals['size_minor_true'][igals]
+            param['bfrac'] = gals['bulge_to_total_ratio_i'][igals]
+            param['theta'] = gals['position_angle_true'][igals]
+            param['hlr_bulge'] = np.sqrt(hlrMajB * hlrMinB)
+            param['hlr_disk'] = np.sqrt(hlrMajD * hlrMinD)
+            param['ell_bulge'] = (hlrMajB - hlrMinB)/(hlrMajB + hlrMinB)
+            param['ell_disk'] = (hlrMajD - hlrMinD)/(hlrMajD + hlrMinD)
+            param['ell_tot'] = (aGal - bGal) / (aGal + bGal)
+
+            # Assign each galaxy a template SED
+            param['sed_type'] = sed_assign(phz=param['z'], btt=param['bfrac'], rng=self.rng_sedGal)
+            param['redden'] = self.tempRed_gal[param['sed_type']]
+            param['av'] = self.avGal[int(self.ud()*self.nav)]
+            if param['sed_type'] <= 5:
+                param['av'] = 0.0
+                param['redden'] = 0
+            param['star'] = 0   # Galaxy
+            if param['sed_type'] >= 29:
+                param['av'] = 0.6 * param['av'] / 3.0 # for quasar, av=[0, 0.2], 3.0=av.max-av.im
+                param['star'] = 2 # Quasar
+
+            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
+                continue
+            param['mag_use_normal'] = gals['mag_true_g_lsst'][igals]
+            if param['mag_use_normal'] >= 26.5:
+                continue
+            self.ids += 1
+            param['id'] = self.ids
+            
+            if param['star'] == 0:
+                obj = Galaxy(param, self.rotation)
+                self.objs.append(obj)
+            if param['star'] == 2:
+                obj = Quasar(param)
+                self.objs.append(obj)
+
+    def _load_stars(self, stars, pix_id=None):
+        nstars = len(stars['sourceID'])
+        for istars in range(nstars):
+            param = self.initialize_param()
+            param['ra'] = stars['RA'][istars]
+            param['dec'] = stars['Dec'][istars]
+            if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
+                continue
+            param['mag_use_normal'] = stars['app_sdss_g'][istars]
+            if param['mag_use_normal'] >= 26.5:
+                continue
+            self.ids += 1
+            param['id'] = self.ids
+            param['sed_type'] = stars['sourceID'][istars]
+            param['model_tag'] = stars['model_tag'][istars]
+            param['teff'] = stars['teff'][istars]
+            param['logg'] = stars['grav'][istars]
+            param['feh'] = stars['feh'][istars]
+            param['z'] = 0.0
+            param['star'] = 1   # Star
+            obj = Star(param)
+            self.objs.append(obj)
+
+    def _load(self, **kwargs):
+        self.nav = 15005
+        self.avGal = extAv(self.nav, seed=self.seed_Av)
+        gals_cat = h5.File(self.galaxy_path, 'r')['galaxies']
+        star_cat = h5.File(self.star_path, 'r')['catalog']
+        self.objs = []
+        self.ids = 0
+        for pix in self.pix_list:
+            gals = gals_cat[str(pix)]
+            stars = star_cat[str(pix)]
+            self._load_gals(gals, pix_id=pix)
+            self._load_stars(stars, pix_id=pix)
+        print("number of objects in catalog: ", len(self.objs))
+        del self.avGal
+
+
+    def load_sed(self, obj, **kwargs):
+        if obj.type == 'star':
+            _, wave, flux = tag_sed(
+                h5file=self.tempSED_star,
+                model_tag=obj.param['model_tag'],
+                teff=obj.param['teff'],
+                logg=obj.param['logg'],
+                feh=obj.param['feh']
+            )
+        elif obj.type == 'galaxy' or obj.type == 'quasar':
+            sed_data = getObservedSED(
+                sedCat=self.tempSed_gal[obj.sed_type],
+                redshift=obj.z,
+                av=obj.param["av"],
+                redden=obj.param["redden"]
+            )
+            wave, flux = sed_data[0], sed_data[1]
+        else:
+            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)
+        y = speci(lamb)
+        # erg/s/cm2/A --> photo/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'))
+        return sed
diff --git a/Catalog/Catalog_example.py b/Catalog/Catalog_example.py
new file mode 100644
index 0000000000000000000000000000000000000000..c8af865d4e43f2344a99d6f34a8771929d96bf28
--- /dev/null
+++ b/Catalog/Catalog_example.py
@@ -0,0 +1,198 @@
+import os
+import numpy as np
+import astropy.constants as cons
+from astropy.table import Table
+from scipy import interpolate
+
+from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar
+
+class Catalog_example(CatalogBase):
+    """An user customizable class for reading in catalog(s) of objects and SEDs.
+    
+    NOTE: must inherit the "CatalogBase" abstract class
+
+    ...
+    
+    Attributes
+    ----------
+    cat_dir : str
+        a directory that contains the catalog file(s)
+    star_path : str
+        path to the star catalog file
+    star_SED_path : str
+        path to the star SED data
+    objs : list
+        a list of ObservationSim.MockObject (Star, Galaxy, or Quasar)
+        NOTE: must have "obj" list when implement your own Catalog
+    
+    Methods
+    ----------
+    load_sed(obj, **kwargs):
+        load the corresponding SED data for one object
+    load_norm_filt(obj):
+        load the filter throughput for the input catalog's photometric system.
+    """
+
+    def __init__(self, config, **kwargs):
+        """Constructor method.
+        
+        Parameters
+        ----------
+        config : dict
+            configuration dictionary which is parsed from the input YAML file
+        **kwargs : dict
+            other needed input parameters (in key-value pairs), please modify corresponding
+            initialization call in "ObservationSim.py" as you need.
+        
+        Returns
+        ----------
+        None
+        """
+        
+        super().__init__()
+        self.cat_dir = os.path.join(config["data_dir"], config["input_path"]["cat_dir"])
+        if "star_cat" in config["input_path"] and config["input_path"]["star_cat"]:
+            star_file = config["input_path"]["star_cat"]
+            star_SED_file = config["SED_templates_path"]["star_SED"]
+            self.star_path = os.path.join(self.cat_dir, star_file)
+            self.star_SED_path = os.path.join(config["data_dir"], star_SED_file)
+        # NOTE: must call _load() method here to read in all objects
+        self.objs = []
+        self._load()
+
+    def _load(self, **kwargs):
+        """Read in all objects in from the catalog file(s).
+
+        This is a must implemented method which is used to read in all objects, and
+        then convert them to ObservationSim.MockObject (Star, Galaxy, or Quasar).
+        
+        Currently, 
+        the model of ObservationSim.MockObject.Star class requires:
+            param["star"] : int
+                specify the object type: 0: galaxy, 1: star, 2: quasar
+            param["id"] : int
+                ID number of the object
+            param["ra"] : float
+                Right ascension (in degrees)
+            param["dec"] : float
+                Declination (in degrees)
+            param["mag_use_normal"]: float
+                the absolute magnitude in a particular filter
+                NOTE: if that filter is not the corresponding CSST filter, the 
+                load_norm_filt(obj) function must be implemented to load the filter
+                throughput of that particular photometric system
+        
+        the model of ObservationSim.MockObject.Galaxy class requires:
+            param["star"] : int
+                specify the object type: 0: galaxy, 1: star, 2: quasar
+            param["id"] : int
+                ID number of the object
+            param["ra"] : float
+                Right ascension (in degrees)
+            param["dec"] : float
+                Declination (in degrees)
+            param["mag_use_normal"]: float
+                the absolute magnitude in a particular filter
+                NOTE: if that filter is not the corresponding CSST filter, the 
+                load_norm_filt(obj) function must be implemented to load the filter
+                throughput of that particular photometric system
+            param["theta"] : float
+                the position angle (in degrees)
+            param["bfrac"] : float
+                the bulge fraction
+            param["hlr_bulge] : float
+                the half-light-radius of the bulge
+            param["hlr_disk] : float
+                the half-light-radius of the disk
+             param["ell_bulge] : float
+                the ellipticity of the bulge
+             param["ell_disk] : float
+                the ellipticity of the disk
+             param["ell_tot] : float
+                the total ellipticity
+        the model of ObservationSim.MockObject.Galaxy class requires:
+            Currently a Quasar is modeled as a point source, just like a Star.
+
+        NOTE: To construct an object, according to its type, just call:
+        Star(param), Galaxy(param), or Quasar(param)
+
+        NOTE: All constructed objects should be appened to "self.objs".
+
+        Parameters
+        ----------
+        **kwargs : dict
+            other needed input parameters (in key-value pairs), please modify corresponding
+            initialization call in "ObservationSim.py" as you need.
+
+        Returns
+        ----------
+        None
+        """
+
+        stars = Table.read(self.star_path)
+        nstars = stars['sourceID'].size
+        for istars in range(nstars):
+            param = self.initialize_param()
+            param['id'] = istars + 1
+            param['ra'] = stars['RA'][istars]
+            param['dec'] = stars['Dec'][istars]
+            param['sed_type'] = stars['sourceID'][istars]
+            param['model_tag'] = stars['model_tag'][istars]
+            param['z'] = 0.0
+            param['star'] = 1   # Star
+            param['mag_use_normal'] = stars['app_sdss_g'][istars]
+            obj = Star(param)
+            self.objs.append(obj)
+
+    def load_sed(self, obj, **kwargs):
+        """Load the corresponding SED data for a particular obj.
+
+        Parameters
+        ----------
+        obj : ObservationSim.MockObject
+            the object to get SED data for
+        **kwargs : dict
+            other needed input parameters (in key-value pairs), please modify corresponding
+            initialization call in "ObservationSim.py" as you need.
+
+        Returns
+        ----------
+        sed : Astropy.Table
+            the SED Table with two columns (namely, "WAVELENGTH", "FLUX"):
+                sed["WAVELENGTH"] : wavelength in Angstroms
+                sed["FLUX"] : fluxes in photons/s/m^2/A
+            NOTE: the range of wavelengthes must at least cover [2450 - 11000] Angstorms
+        """
+
+        if obj.type == 'star':
+            wave = Table.read(self.star_SED_path, path=f"/SED/wave_{obj.model_tag}")
+            flux = Table.read(self.star_SED_path, path=f"/SED/{obj.sed_type}")
+            wave, flux = wave['col0'].data, flux['col0'].data
+        else:
+            raise ValueError("Object type not known")
+        speci = interpolate.interp1d(wave, flux)
+        lamb = np.arange(2400, 11001 + 0.5, 0.5)
+        y = speci(lamb)
+        # erg/s/cm^2/A --> photons/s/m^2/A
+        all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
+        sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
+        return sed
+
+    def load_norm_filt(self, obj):
+        """Load the corresponding thourghput for the input magnitude "param["mag_use_normal"]".
+
+        NOTE: if the input magnitude is already in CSST magnitude, simply return None
+
+        Parameters
+        ----------
+        obj : ObservationSim.MockObject
+            the object to get thourghput data for
+
+        Returns
+        ----------
+        norm_filt : Astropy.Table
+            the throughput Table with two columns (namely, "WAVELENGTH", "SENSITIVITY"):
+                norm_filt["WAVELENGTH"] : wavelengthes in Angstroms
+                norm_filt["SENSITIVITY"] : efficiencies
+        """
+        return None
\ No newline at end of file
diff --git a/ObservationSim/Config/Header/__pycache__/ImageHeader.cpython-37.pyc b/ObservationSim/Config/Header/__pycache__/ImageHeader.cpython-37.pyc
deleted file mode 100644
index 8d4dfb7928373cd54f48c9ac5298e1b4ad691f77..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Config/Header/__pycache__/ImageHeader.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/Config/Header/__pycache__/ImageHeader.cpython-38.pyc b/ObservationSim/Config/Header/__pycache__/ImageHeader.cpython-38.pyc
deleted file mode 100644
index 5fa40c48881f6052fa33c5dcaef8bd21e388b296..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Config/Header/__pycache__/ImageHeader.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/Config/Header/__pycache__/__init__.cpython-37.pyc b/ObservationSim/Config/Header/__pycache__/__init__.cpython-37.pyc
deleted file mode 100644
index 4bfc708af882253a8bf3e3ac92c04b2ee8f3cbee..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Config/Header/__pycache__/__init__.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/Config/Header/__pycache__/__init__.cpython-38.pyc b/ObservationSim/Config/Header/__pycache__/__init__.cpython-38.pyc
deleted file mode 100644
index 1341eda0d8a2615b2a576d9fc08910af4b3b5ddd..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Config/Header/__pycache__/__init__.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/Config/__pycache__/ChipOutput.cpython-38.pyc b/ObservationSim/Config/__pycache__/ChipOutput.cpython-38.pyc
deleted file mode 100644
index f155567189d29c105786f72576465fb0589c6466..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Config/__pycache__/ChipOutput.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/Config/__pycache__/Config.cpython-38.pyc b/ObservationSim/Config/__pycache__/Config.cpython-38.pyc
deleted file mode 100644
index ab25d1fc4ed1628753d8f0d812fd18ce6da7509b..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Config/__pycache__/Config.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/Config/__pycache__/__init__.cpython-38.pyc b/ObservationSim/Config/__pycache__/__init__.cpython-38.pyc
deleted file mode 100644
index c334df4b112cf42398995ec284caf5859e8f31ff..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Config/__pycache__/__init__.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/Instrument/Chip/__pycache__/Chip.cpython-38.pyc b/ObservationSim/Instrument/Chip/__pycache__/Chip.cpython-38.pyc
deleted file mode 100644
index 494be8e7feb822165cd70af27a4429697328180a..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Instrument/Chip/__pycache__/Chip.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/Instrument/Chip/__pycache__/Effects.cpython-38.pyc b/ObservationSim/Instrument/Chip/__pycache__/Effects.cpython-38.pyc
deleted file mode 100644
index 5e2fc0e5926b1ffd19e65dd2b8aebc6b4f963435..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Instrument/Chip/__pycache__/Effects.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/Instrument/Chip/__pycache__/__init__.cpython-38.pyc b/ObservationSim/Instrument/Chip/__pycache__/__init__.cpython-38.pyc
deleted file mode 100644
index 08322b16d722afbc1f14d8449b5c2b3cb2df92ff..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Instrument/Chip/__pycache__/__init__.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/Instrument/Filter.py b/ObservationSim/Instrument/Filter.py
index 082d6a6e75cc37f21837ba741f6ddd8248f6f7ba..5268ddc5119177c8273f73fb73a821af101d5d15 100755
--- a/ObservationSim/Instrument/Filter.py
+++ b/ObservationSim/Instrument/Filter.py
@@ -14,6 +14,13 @@ class Filter(object):
         self._getParam(filter_param, filter_type)
         if self.filter_dir is not None:
             self.bandpass_full, self.bandpass_sub_list = self._get_bandpasses(self.filter_dir)
+        self.survey_type = self._getSurveyType()
+
+    def _getSurveyType(self):
+        if self.filter_type in ["GI", "GV", "GU"]:
+            return "spectroscopic"
+        else:
+            return "photometric"
 
     def _getParam(self, filter_param, filter_type, filter_id=None):
         self.effective_wavelength = filter_param.param[filter_type][0]
diff --git a/ObservationSim/Instrument/__pycache__/Chip.cpython-37.pyc b/ObservationSim/Instrument/__pycache__/Chip.cpython-37.pyc
deleted file mode 100755
index 183388db7049dbd6fa8f3abd15c9cc2ffa283b6e..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Instrument/__pycache__/Chip.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/Instrument/__pycache__/Chip.cpython-38.pyc b/ObservationSim/Instrument/__pycache__/Chip.cpython-38.pyc
deleted file mode 100644
index cae0f4b8b7eccbcc562c8e138024b1c4d41d460f..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Instrument/__pycache__/Chip.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/Instrument/__pycache__/Filter.cpython-37.pyc b/ObservationSim/Instrument/__pycache__/Filter.cpython-37.pyc
deleted file mode 100755
index beccb3c2b88e83b44bac98da550f1f0b719d3a84..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Instrument/__pycache__/Filter.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/Instrument/__pycache__/Filter.cpython-38.pyc b/ObservationSim/Instrument/__pycache__/Filter.cpython-38.pyc
deleted file mode 100644
index 3a3a9a6609e559f8145ea9a1a0883548e44037d0..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Instrument/__pycache__/Filter.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/Instrument/__pycache__/FilterParam.cpython-37.pyc b/ObservationSim/Instrument/__pycache__/FilterParam.cpython-37.pyc
deleted file mode 100755
index db20be7c339424b182f606e30679905cd09fced7..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Instrument/__pycache__/FilterParam.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/Instrument/__pycache__/FilterParam.cpython-38.pyc b/ObservationSim/Instrument/__pycache__/FilterParam.cpython-38.pyc
deleted file mode 100644
index 881c1cf95e5e6c1a90ac9b8e38cc20e4e57f83e8..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Instrument/__pycache__/FilterParam.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/Instrument/__pycache__/FocalPlane.cpython-37.pyc b/ObservationSim/Instrument/__pycache__/FocalPlane.cpython-37.pyc
deleted file mode 100755
index 39d9b6fd4de46101b087f4e50222b5021f88cf09..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Instrument/__pycache__/FocalPlane.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/Instrument/__pycache__/FocalPlane.cpython-38.pyc b/ObservationSim/Instrument/__pycache__/FocalPlane.cpython-38.pyc
deleted file mode 100644
index 7997be3e03e48f50d9ffaf6c74af880e7f256906..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Instrument/__pycache__/FocalPlane.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/Instrument/__pycache__/Telescope.cpython-37.pyc b/ObservationSim/Instrument/__pycache__/Telescope.cpython-37.pyc
deleted file mode 100755
index 369f398d22bbecf3ad259e6a3ef0de54c9378ccc..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Instrument/__pycache__/Telescope.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/Instrument/__pycache__/Telescope.cpython-38.pyc b/ObservationSim/Instrument/__pycache__/Telescope.cpython-38.pyc
deleted file mode 100644
index dae2816701485d0eebfc28213839a1e2a6ac0086..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Instrument/__pycache__/Telescope.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/Instrument/__pycache__/__init__.cpython-37.pyc b/ObservationSim/Instrument/__pycache__/__init__.cpython-37.pyc
deleted file mode 100755
index 274d5de542491de8323c9f2b94a2b8c2677f187b..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Instrument/__pycache__/__init__.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/Instrument/__pycache__/__init__.cpython-38.pyc b/ObservationSim/Instrument/__pycache__/__init__.cpython-38.pyc
deleted file mode 100644
index 92ea9cd16903a0ace81a3b6a93432ec42993d90e..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Instrument/__pycache__/__init__.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/Instrument/__pycache__/_util.cpython-37.pyc b/ObservationSim/Instrument/__pycache__/_util.cpython-37.pyc
deleted file mode 100755
index e99a51d4a89c8a907d883e8d45cfe2714f1e5c62..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Instrument/__pycache__/_util.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/Instrument/__pycache__/_util.cpython-38.pyc b/ObservationSim/Instrument/__pycache__/_util.cpython-38.pyc
deleted file mode 100644
index 59c3215c5df8d78daeb96409c9ab3b54c20e1fbd..0000000000000000000000000000000000000000
Binary files a/ObservationSim/Instrument/__pycache__/_util.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/MockObject/Catalog.py b/ObservationSim/MockObject/Catalog.py
index d99f9cec29e76dcd1f429c12a0668e857a2dc404..a86d36fe2f0c115b1e38e2b37e03b07479772b7b 100644
--- a/ObservationSim/MockObject/Catalog.py
+++ b/ObservationSim/MockObject/Catalog.py
@@ -33,14 +33,14 @@ class Catalog(object):
         # star_file = 'stars_ccd' + chip.getChipLabel(chip.chipID) + '_p_RA'+str(self.pRa) + '_DE' + str(self.pDec) + '.hdf5'
         # galaxy_file = 'galaxies_ccd' + chip.getChipLabel(chip.chipID) + '_p_RA'+str(self.pRa) + '_DE' + str(self.pDec) + '.hdf5'
 
-        if "star_cat" in config["input_path"]:
+        if "star_cat" in config["input_path"] and config["input_path"]["star_cat"]:
             star_file = config["input_path"]["star_cat"]
             star_SED_file = config["SED_templates_path"]["star_SED"]
             self.star_path = os.path.join(self.cat_dir, star_file)
             self.star_SED_path = os.path.join(config["data_dir"], star_SED_file)
             self._load_SED_lib_star()
-        if "galaxy_cat" in config["input_path"]:
 
+        if "galaxy_cat" in config["input_path"] and config["input_path"]["galaxy_cat"]:
             galaxy_file = config["input_path"]["galaxy_cat"]
             self.galaxy_path = os.path.join(self.cat_dir, galaxy_file)
             self.galaxy_SED_path = os.path.join(config["data_dir"], config["SED_templates_path"]["galaxy_SED"])
@@ -86,7 +86,6 @@ class Catalog(object):
         self.rng_sedGal = random.Random()
         self.rng_sedGal.seed(pix_id) # Use healpix index as the random seed
         self.ud = galsim.UniformDeviate(pix_id)
-        # print(ngals)
         for igals in range(ngals):
             param = {}
             param['ra'] = gals['ra_true'][igals]
@@ -143,7 +142,6 @@ class Catalog(object):
 
     def _load_stars(self, stars, pix_id=None):
         nstars = len(stars['sourceID'])
-        # print(nstars)
         for istars in range(nstars):
             param = {}
             param['ra'] = stars['RA'][istars]
@@ -168,10 +166,6 @@ class Catalog(object):
     def _load(self):
         self.nav = 15005
         self.avGal = extAv(self.nav, seed=self.seed_Av)
-        # gals = Table.read(self.galaxy_path)
-        # stars = Table.read(self.star_path)
-        # self.ngals = gals['galaxyID'].size
-        # self.nstars = stars['sourceID'].size
         gals_cat = h5.File(self.galaxy_path, 'r')['galaxies']
         star_cat = h5.File(self.star_path, 'r')['catalog']
         self.objs = []
@@ -182,66 +176,4 @@ class Catalog(object):
             self._load_gals(gals, pix_id=pix)
             self._load_stars(stars, pix_id=pix)
         print("number of objects in catalog: ", len(self.objs))
-        del self.avGal
-        # for igals in range(self.ngals):
-        #     param = {}
-        #     param['id'] = igals + 1
-        #     param['ra'] = gals['ra_true'][igals]
-        #     param['dec'] = gals['dec_true'][igals]
-        #     param['z'] = gals['redshift_true'][igals]
-        #     param['model_tag'] = 'None'
-        #     param['gamma1'] = 0
-        #     param['gamma2'] = 0
-        #     param['kappa'] = 0
-        #     param['delta_ra'] = 0
-        #     param['delta_dec'] = 0
-        #     sersicB = gals['sersic_bulge'][igals]
-        #     hlrMajB = gals['size_bulge_true'][igals]
-        #     hlrMinB = gals['size_minor_bulge_true'][igals]
-        #     sersicD = gals['sersic_disk'][igals]
-        #     hlrMajD = gals['size_disk_true'][igals]
-        #     hlrMinD = gals['size_minor_disk_true'][igals]
-        #     aGal = gals['size_true'][igals]
-        #     bGal = gals['size_minor_true'][igals]
-        #     param['bfrac'] = gals['bulge_to_total_ratio_i'][igals]
-        #     param['theta'] = gals['position_angle_true'][igals]
-        #     param['hlr_bulge'] = np.sqrt(hlrMajB * hlrMinB)
-        #     param['hlr_disk'] = np.sqrt(hlrMajD * hlrMinD)
-        #     param['ell_bulge'] = (hlrMajB - hlrMinB)/(hlrMajB + hlrMinB)
-        #     param['ell_disk'] = (hlrMajD - hlrMinD)/(hlrMajD + hlrMinD)
-        #     param['ell_tot'] = (aGal - bGal) / (aGal + bGal)
-
-        #     # Assign each galaxy a template SED
-        #     param['sed_type'] = sed_assign(phz=param['z'], btt=param['bfrac'], rng=self.rng_sedGal)
-        #     param['redden'] = self.tempRed_gal[param['sed_type']]
-        #     param['av'] = avGal[int(self.ud()*nav)]
-        #     if param['sed_type'] <= 5:
-        #         param['av'] = 0.0
-        #         param['redden'] = 0
-        #     param['star'] = 0   # Galaxy
-        #     if param['sed_type'] >= 29:
-        #         param['av'] = 0.6 * param['av'] / 3.0 # for quasar, av=[0, 0.2], 3.0=av.max-av.im
-        #         param['star'] = 2 # Quasar
-        #     param['mag_use_normal'] = gals['mag_true_g_lsst'][igals]
-        #     if param['star'] == 0:
-        #         obj = Galaxy(param, self.rotation)
-        #         self.objs.append(obj)
-        #     if param['star'] == 2:
-        #         obj = Quasar(param)
-        #         self.objs.append(obj)
-
-        # for istars in range(self.nstars):
-        #     param = {}
-        #     param['id'] = self.ngals + istars + 1
-        #     param['ra'] = stars['RA'][istars]
-        #     param['dec'] = stars['Dec'][istars]
-        #     param['sed_type'] = stars['sourceID'][istars]
-        #     param['model_tag'] = stars['model_tag'][istars]
-        #     param['z'] = 0.0
-        #     param['star'] = 1   # Star
-        #     param['mag_use_normal'] = stars['app_sdss_g'][istars]
-        #     obj = Star(param)
-        #     self.objs.append(obj)
-        # del gals
-        # del stars
-        # del avGal
\ No newline at end of file
+        del self.avGal
\ No newline at end of file
diff --git a/ObservationSim/MockObject/CatalogBase.py b/ObservationSim/MockObject/CatalogBase.py
new file mode 100644
index 0000000000000000000000000000000000000000..bcb287e4a19232b98689090438c4d4dcd980e205
--- /dev/null
+++ b/ObservationSim/MockObject/CatalogBase.py
@@ -0,0 +1,85 @@
+import numpy as np
+import galsim
+import copy
+
+from astropy.table import Table
+from abc import abstractmethod, ABCMeta
+from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getABMAG
+
+
+class CatalogBase(metaclass=ABCMeta):
+    def __init__(self):
+        pass
+
+    @abstractmethod
+    def load_sed(self, obj, **kward):
+        pass
+
+    @abstractmethod
+    def _load(self, **kward):
+        pass
+
+    @abstractmethod
+    def load_norm_filt(self, obj):
+        return None
+    
+    @staticmethod
+    def initialize_param():
+        param = {
+            "star":-1,
+            "id":-1,
+            "ra":0,
+            "dec":0,
+            "z":0,
+            "sed_type":-1,
+            "model_tag":"unknown",
+            "mag_use_normal":100,
+            "theta":0,
+            "kappa":0,
+            "gamma1":0,
+            "gamma2":0,
+            "bfrac":0,
+            "hlr_bulge":0,
+            "hlr_disk":0,
+            "ell_bulge":0,
+            "ell_disk":0,
+            "ell_tot":0,
+            "teff":0,
+            "logg":0,
+            "feh":0
+        }
+        return param
+    
+    @staticmethod
+    def convert_sed(mag, sed, target_filt, norm_filt=None):
+        bandpass = target_filt.bandpass_full
+        
+        if norm_filt is not None:
+            norm_thr_rang_ids = norm_filt['SENSITIVITY'] > 0.001
+        else:
+            norm_filt = Table(
+                np.array(np.array([bandpass.wave_list*10.0, bandpass.func(bandpass.wave_list)])).T, names=(['WAVELENGTH', 'SENSITIVITY'])
+            )
+            norm_thr_rang_ids = norm_filt['SENSITIVITY'] > 0.001
+        
+        sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=mag,
+                spectrum=sed,
+                norm_thr=norm_filt,
+                sWave=np.floor(norm_filt[norm_thr_rang_ids][0][0]),
+                eWave=np.ceil(norm_filt[norm_thr_rang_ids][-1][0]))
+        sed_photon = copy.copy(sed)
+        sed_photon = np.array([sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T
+        sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest')
+        # Get magnitude
+        sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False)
+        interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=bandpass)
+        mag_csst = getABMAG(
+            interFlux=interFlux,
+            bandpass=bandpass
+        )
+        if target_filt.survey_type == "photometric":
+            return sed_photon, mag_csst
+        elif target_filt.survey_type == "spectroscopic":
+            del sed_photon
+            return sed, mag_csst
+        
\ No newline at end of file
diff --git a/ObservationSim/MockObject/Galaxy.py b/ObservationSim/MockObject/Galaxy.py
index fdaa80f91a33d1e1b97f43fd6859a001dc273adf..8462114c424aeca4853617995b4657071f3115c1 100755
--- a/ObservationSim/MockObject/Galaxy.py
+++ b/ObservationSim/MockObject/Galaxy.py
@@ -28,65 +28,65 @@ class Galaxy(MockObject):
         if rotation is not None:
             self.rotateEllipticity(rotation)
 
-    def load_SED(self, survey_type, sed_path=None, cosids=None, objtypes=None, sed_templates=None, normFilter=None, target_filt=None):
-        if survey_type == "photometric":
-            norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001
-            if sed_templates is None:
-                # Read SED data directly
-                itype = objtypes[cosids==self.sed_type][0]
-                sed_file = os.path.join(sed_path, itype + "_ID%s.sed"%(self.sed_type))
-                if not os.path.exists(sed_file):
-                    raise ValueError("!!! No SED found.")
-                sed_data = Table.read(sed_file, format="ascii")
-                wave, flux = sed_data["observedLambda"].data, sed_data["observedFlux"].data
-            else:
-                # Load SED from templates
-                sed_data = sed_templates[self.sed_type]
-                # redshift, intrinsic extinction
-                sed_data = getObservedSED(
-                    sedCat=sed_data, 
-                    redshift=self.z, 
-                    av=self.param['av'], 
-                    redden=self.param['redden'])
-                wave, flux = sed_data[0], sed_data[1] 
-            flux_photon = flux * (wave / (cons.h.value * cons.c.value)) * 1e-13
-            sed_photon = Table(np.array([wave, flux_photon]).T, names=('WAVELENGTH', 'FLUX'))
-            # Get scaling factor for SED
-            sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'],
-                spectrum=sed_photon,
-                norm_thr=normFilter,
-                sWave=np.floor(normFilter[norm_thr_rang_ids][0][0]),
-                eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0]))
-            sed_photon = np.array([sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T
-            # Convert to galsim.SED object
-            spec = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest')
-            self.sed = galsim.SED(spec, wave_type='A', flux_type='1', fast=False)
-            # Get magnitude
-            interFlux = integrate_sed_bandpass(sed=self.sed, bandpass=target_filt.bandpass_full)
-            self.param['mag_%s'%target_filt.filter_type] = getABMAG(
-                interFlux=interFlux, 
-                bandpass=target_filt.bandpass_full)
-            # print('mag_use_normal = ', self.param['mag_use_normal'])
-            # print('mag_%s = '%target_filt.filter_type, self.param['mag_%s'%target_filt.filter_type])
-            # print('redshift = %.3f'%(self.z))
-            # print('sed_type = %d, av = %.2f, redden = %d'%(self.sed_type, self.param['av'], self.param['redden']))
-
-        elif survey_type == "spectroscopic":
-            if sed_templates is None:
-                self.sedPhotons(sed_path=sed_path, cosids=cosids, objtypes=objtypes)
-            else:
-                sed_data = sed_templates[self.sed_type]
-                sed_data = getObservedSED(
-                    sedCat=sed_data, 
-                    redshift=self.z, 
-                    av=self.param['av'], 
-                    redden=self.param['redden'])
-                speci = interpolate.interp1d(sed_data[0], sed_data[1])
-                lamb = np.arange(2500, 10001 + 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
-                self.sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
+    # def load_SED(self, survey_type, sed_path=None, cosids=None, objtypes=None, sed_templates=None, normFilter=None, target_filt=None):
+    #     if survey_type == "photometric":
+    #         norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001
+    #         if sed_templates is None:
+    #             # Read SED data directly
+    #             itype = objtypes[cosids==self.sed_type][0]
+    #             sed_file = os.path.join(sed_path, itype + "_ID%s.sed"%(self.sed_type))
+    #             if not os.path.exists(sed_file):
+    #                 raise ValueError("!!! No SED found.")
+    #             sed_data = Table.read(sed_file, format="ascii")
+    #             wave, flux = sed_data["observedLambda"].data, sed_data["observedFlux"].data
+    #         else:
+    #             # Load SED from templates
+    #             sed_data = sed_templates[self.sed_type]
+    #             # redshift, intrinsic extinction
+    #             sed_data = getObservedSED(
+    #                 sedCat=sed_data, 
+    #                 redshift=self.z, 
+    #                 av=self.param['av'], 
+    #                 redden=self.param['redden'])
+    #             wave, flux = sed_data[0], sed_data[1] 
+    #         flux_photon = flux * (wave / (cons.h.value * cons.c.value)) * 1e-13
+    #         sed_photon = Table(np.array([wave, flux_photon]).T, names=('WAVELENGTH', 'FLUX'))
+    #         # Get scaling factor for SED
+    #         sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'],
+    #             spectrum=sed_photon,
+    #             norm_thr=normFilter,
+    #             sWave=np.floor(normFilter[norm_thr_rang_ids][0][0]),
+    #             eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0]))
+    #         sed_photon = np.array([sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T
+    #         # Convert to galsim.SED object
+    #         spec = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest')
+    #         self.sed = galsim.SED(spec, wave_type='A', flux_type='1', fast=False)
+    #         # Get magnitude
+    #         interFlux = integrate_sed_bandpass(sed=self.sed, bandpass=target_filt.bandpass_full)
+    #         self.param['mag_%s'%target_filt.filter_type] = getABMAG(
+    #             interFlux=interFlux, 
+    #             bandpass=target_filt.bandpass_full)
+    #         # print('mag_use_normal = ', self.param['mag_use_normal'])
+    #         # print('mag_%s = '%target_filt.filter_type, self.param['mag_%s'%target_filt.filter_type])
+    #         # print('redshift = %.3f'%(self.z))
+    #         # print('sed_type = %d, av = %.2f, redden = %d'%(self.sed_type, self.param['av'], self.param['redden']))
+
+    #     elif survey_type == "spectroscopic":
+    #         if sed_templates is None:
+    #             self.sedPhotons(sed_path=sed_path, cosids=cosids, objtypes=objtypes)
+    #         else:
+    #             sed_data = sed_templates[self.sed_type]
+    #             sed_data = getObservedSED(
+    #                 sedCat=sed_data, 
+    #                 redshift=self.z, 
+    #                 av=self.param['av'], 
+    #                 redden=self.param['redden'])
+    #             speci = interpolate.interp1d(sed_data[0], sed_data[1])
+    #             lamb = np.arange(2500, 10001 + 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
+    #             self.sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
 
 
     def unload_SED(self):
@@ -94,23 +94,23 @@ class Galaxy(MockObject):
         """
         del self.sed
 
-    def sedPhotons(self, sed_path, cosids, objtypes):
-
-        itype = objtypes[cosids == self.sed_type][0]
-        sed_file = os.path.join(sed_path, itype + "_ID%s.sed" % (self.sed_type))
-        if not os.path.exists(sed_file):
-            raise ValueError("!!! No SED found.")
-        sed = Table.read(sed_file, format="ascii")
-        spec_data = {}
-        f_orig = sed["observedFlux"].data
-        w_orig = sed["observedLambda"].data
-
-        speci = interpolate.interp1d(w_orig, f_orig)
-        lamb = np.arange(2500, 10001 + 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
-        self.sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
+    # def sedPhotons(self, sed_path, cosids, objtypes):
+
+    #     itype = objtypes[cosids == self.sed_type][0]
+    #     sed_file = os.path.join(sed_path, itype + "_ID%s.sed" % (self.sed_type))
+    #     if not os.path.exists(sed_file):
+    #         raise ValueError("!!! No SED found.")
+    #     sed = Table.read(sed_file, format="ascii")
+    #     spec_data = {}
+    #     f_orig = sed["observedFlux"].data
+    #     w_orig = sed["observedLambda"].data
+
+    #     speci = interpolate.interp1d(w_orig, f_orig)
+    #     lamb = np.arange(2500, 10001 + 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
+    #     self.sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
 
     def getGSObj_multiband(self, tel, psf_list, bandpass_list, filt, nphotons_tot=None, g1=0, g2=0, exptime=150.):
         if len(psf_list) != len(bandpass_list):
diff --git a/ObservationSim/MockObject/MockObject.py b/ObservationSim/MockObject/MockObject.py
index d9fa2b43cf2e00c29ed76c1266c7d1ca4a932be4..a0ebdc4b58d973bef8df8285971da88d55dd16fb 100755
--- a/ObservationSim/MockObject/MockObject.py
+++ b/ObservationSim/MockObject/MockObject.py
@@ -3,7 +3,7 @@ import numpy as np
 import astropy.constants as cons
 from astropy.table import Table
 
-from ObservationSim.MockObject._util import magToFlux, vc_A, convolveGaussXorders
+from ObservationSim.MockObject._util import magToFlux, VC_A, convolveGaussXorders
 from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG
 from ObservationSim.MockObject.SpecDisperser import SpecDisperser
 
@@ -24,6 +24,7 @@ class MockObject(object):
         self.sed_type = self.param["sed_type"]
         self.model_tag = self.param["model_tag"]
         self.mag_use_normal = self.param["mag_use_normal"]
+        self.sed = None
 
     def getMagFilter(self, filt):
         if filt.filter_type in ["GI", "GV", "GU"]:
@@ -39,7 +40,7 @@ class MockObject(object):
     def getElectronFluxFilt(self, filt, tel, exptime=150.):
         photonEnergy = filt.getPhotonE()
         flux = magToFlux(self.getMagFilter(filt))
-        factor = 1.0e4 * flux/photonEnergy * vc_A * (1.0/filt.blue_limit - 1.0/filt.red_limit)
+        factor = 1.0e4 * flux/photonEnergy * VC_A * (1.0/filt.blue_limit - 1.0/filt.red_limit)
         return factor * filt.efficiency * tel.pupil_area * exptime
 
     def getPosWorld(self):
diff --git a/ObservationSim/MockObject/SpecDisperser/__pycache__/SpecDisperser.cpython-38.pyc b/ObservationSim/MockObject/SpecDisperser/__pycache__/SpecDisperser.cpython-38.pyc
deleted file mode 100644
index 4b3396d630ade0fecffa135e8888c359181ca56f..0000000000000000000000000000000000000000
Binary files a/ObservationSim/MockObject/SpecDisperser/__pycache__/SpecDisperser.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/MockObject/SpecDisperser/__pycache__/__init__.cpython-38.pyc b/ObservationSim/MockObject/SpecDisperser/__pycache__/__init__.cpython-38.pyc
deleted file mode 100644
index 13c4d80fd308e0ffa85ad98ed98efc58c7b55605..0000000000000000000000000000000000000000
Binary files a/ObservationSim/MockObject/SpecDisperser/__pycache__/__init__.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/MockObject/SpecDisperser/disperse_c/__pycache__/__init__.cpython-37.pyc b/ObservationSim/MockObject/SpecDisperser/disperse_c/__pycache__/__init__.cpython-37.pyc
deleted file mode 100644
index ec4022c1b4dda433c8525033f33c0941d6be771d..0000000000000000000000000000000000000000
Binary files a/ObservationSim/MockObject/SpecDisperser/disperse_c/__pycache__/__init__.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/MockObject/SpecDisperser/disperse_c/__pycache__/__init__.cpython-38.pyc b/ObservationSim/MockObject/SpecDisperser/disperse_c/__pycache__/__init__.cpython-38.pyc
deleted file mode 100644
index 60a53d5729681648bf79e00b9e0a2b4cad6ace70..0000000000000000000000000000000000000000
Binary files a/ObservationSim/MockObject/SpecDisperser/disperse_c/__pycache__/__init__.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/MockObject/SpecDisperser/disperse_c/disperse.cpython-38-x86_64-linux-gnu.so b/ObservationSim/MockObject/SpecDisperser/disperse_c/disperse.cpython-38-x86_64-linux-gnu.so
index 2ed89c7e7d08aa2449d0d8aad1f2add6cb954021..8a04d59950f5f94aa3a95828315b6ab73590d233 100755
Binary files a/ObservationSim/MockObject/SpecDisperser/disperse_c/disperse.cpython-38-x86_64-linux-gnu.so and b/ObservationSim/MockObject/SpecDisperser/disperse_c/disperse.cpython-38-x86_64-linux-gnu.so differ
diff --git a/ObservationSim/MockObject/SpecDisperser/disperse_c/interp.cpython-38-x86_64-linux-gnu.so b/ObservationSim/MockObject/SpecDisperser/disperse_c/interp.cpython-38-x86_64-linux-gnu.so
index 355c81e7e32f056a1241cd702cc401c41ae88800..815f112f1053f0b64e610550fb677e2cc434b517 100755
Binary files a/ObservationSim/MockObject/SpecDisperser/disperse_c/interp.cpython-38-x86_64-linux-gnu.so and b/ObservationSim/MockObject/SpecDisperser/disperse_c/interp.cpython-38-x86_64-linux-gnu.so differ
diff --git a/ObservationSim/MockObject/Star.py b/ObservationSim/MockObject/Star.py
index b1509130d0145582644a6569cfa43d35982cd235..31cba083ddcf255eff3d323248861f1a7d26be2c 100755
--- a/ObservationSim/MockObject/Star.py
+++ b/ObservationSim/MockObject/Star.py
@@ -17,60 +17,60 @@ class Star(MockObject):
         """
         del self.sed
 
-    def load_SED(self, survey_type, normFilter=None, target_filt=None, sed_lib=None, sed_path=None):
-        if survey_type == "photometric":
-            norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001
-            # spec_lambda = Table.read(sed_path, path=f"/SED/wave_{self.model_tag}")
-            # spec_flux = Table.read(sed_path, path=f"/SED/{self.sed_type}")
-            # wave, flux = spec_lambda['col0'].data, spec_flux['col0'].data
-            _, wave, flux = tag_sed(
-                h5file=sed_lib, 
-                model_tag=self.param['model_tag'], 
-                teff=self.param['teff'],
-                logg=self.param['logg'],
-                feh=self.param['feh'])
-            flux_photon = flux * (wave / (cons.h.value * cons.c.value)) * 1e-13
-            sed_photon = Table(np.array([wave, flux_photon]).T, names=('WAVELENGTH', 'FLUX'))
-            # Get scaling factor for SED
-            sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'],
-                spectrum=sed_photon,
-                norm_thr=normFilter,
-                sWave=np.floor(normFilter[norm_thr_rang_ids][0][0]),
-                eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0]))
-            sed_photon = np.array([sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T
-            # Convert to galsim.SED object
-            spec = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest')
-            self.sed = galsim.SED(spec, wave_type='A', flux_type='1', fast=False)
-            # Get magnitude
-            interFlux = integrate_sed_bandpass(sed=self.sed, bandpass=target_filt.bandpass_full)
-            self.param['mag_%s'%target_filt.filter_type] = getABMAG(
-                interFlux=interFlux, 
-                bandpass=target_filt.bandpass_full)
-            # print('mag_use_normal = ', self.param['mag_use_normal'])
-            # print('mag_%s = '%target_filt.filter_type, self.param['mag_%s'%target_filt.filter_type])
+    # def load_SED(self, survey_type, normFilter=None, target_filt=None, sed_lib=None, sed_path=None):
+    #     if survey_type == "photometric":
+    #         norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001
+    #         # spec_lambda = Table.read(sed_path, path=f"/SED/wave_{self.model_tag}")
+    #         # spec_flux = Table.read(sed_path, path=f"/SED/{self.sed_type}")
+    #         # wave, flux = spec_lambda['col0'].data, spec_flux['col0'].data
+    #         _, wave, flux = tag_sed(
+    #             h5file=sed_lib, 
+    #             model_tag=self.param['model_tag'], 
+    #             teff=self.param['teff'],
+    #             logg=self.param['logg'],
+    #             feh=self.param['feh'])
+    #         flux_photon = flux * (wave / (cons.h.value * cons.c.value)) * 1e-13
+    #         sed_photon = Table(np.array([wave, flux_photon]).T, names=('WAVELENGTH', 'FLUX'))
+    #         # Get scaling factor for SED
+    #         sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'],
+    #             spectrum=sed_photon,
+    #             norm_thr=normFilter,
+    #             sWave=np.floor(normFilter[norm_thr_rang_ids][0][0]),
+    #             eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0]))
+    #         sed_photon = np.array([sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T
+    #         # Convert to galsim.SED object
+    #         spec = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest')
+    #         self.sed = galsim.SED(spec, wave_type='A', flux_type='1', fast=False)
+    #         # Get magnitude
+    #         interFlux = integrate_sed_bandpass(sed=self.sed, bandpass=target_filt.bandpass_full)
+    #         self.param['mag_%s'%target_filt.filter_type] = getABMAG(
+    #             interFlux=interFlux, 
+    #             bandpass=target_filt.bandpass_full)
+    #         # print('mag_use_normal = ', self.param['mag_use_normal'])
+    #         # print('mag_%s = '%target_filt.filter_type, self.param['mag_%s'%target_filt.filter_type])
 
-        elif survey_type == "spectroscopic":
-            # self.sedPhotons(sed_path=sed_path)
-            self.sedPhotons(sed_lib=sed_lib)
+    #     elif survey_type == "spectroscopic":
+    #         # self.sedPhotons(sed_path=sed_path)
+    #         self.sedPhotons(sed_lib=sed_lib)
 
-    def sedPhotons(self, sed_path=None, sed_lib=None):
-        # spec_lambda = Table.read(sed_path, path=f"/SED/wave_{self.model_tag}")
-        # spec_flux = Table.read(sed_path, path=f"/SED/{self.sed_type}")
-        _, w_orig, f_orig = tag_sed(
-                h5file=sed_lib, 
-                model_tag=self.param['model_tag'], 
-                teff=self.param['teff'],
-                logg=self.param['logg'],
-                feh=self.param['feh'])
-        # spec_data = {}
-        # f_orig = spec_flux["col0"].data
-        # w_orig = spec_lambda["col0"].data
-        speci = interpolate.interp1d(w_orig, f_orig)
-        lamb = np.arange(2500, 10001 + 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
-        self.sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
+    # def sedPhotons(self, sed_path=None, sed_lib=None):
+    #     # spec_lambda = Table.read(sed_path, path=f"/SED/wave_{self.model_tag}")
+    #     # spec_flux = Table.read(sed_path, path=f"/SED/{self.sed_type}")
+    #     _, w_orig, f_orig = tag_sed(
+    #             h5file=sed_lib, 
+    #             model_tag=self.param['model_tag'], 
+    #             teff=self.param['teff'],
+    #             logg=self.param['logg'],
+    #             feh=self.param['feh'])
+    #     # spec_data = {}
+    #     # f_orig = spec_flux["col0"].data
+    #     # w_orig = spec_lambda["col0"].data
+    #     speci = interpolate.interp1d(w_orig, f_orig)
+    #     lamb = np.arange(2500, 10001 + 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
+    #     self.sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
 
     def getGSObj(self, psf, g1=0, g2=0, flux=None, filt=None, tel=None, exptime=150.):
         if flux == None:
diff --git a/ObservationSim/MockObject/__init__.py b/ObservationSim/MockObject/__init__.py
index 7ef59a144c3ab44c21f3eb31889c40eb73a47546..4868239200c55c30ee405cad6d0b341a647a1859 100755
--- a/ObservationSim/MockObject/__init__.py
+++ b/ObservationSim/MockObject/__init__.py
@@ -1,5 +1,6 @@
 from .MockObject import MockObject
 from .Galaxy import Galaxy
+from .CatalogBase import CatalogBase
 from .Quasar import Quasar
 from .Star import Star
 from .Catalog import Catalog
diff --git a/ObservationSim/MockObject/__pycache__/Catalog.cpython-38.pyc b/ObservationSim/MockObject/__pycache__/Catalog.cpython-38.pyc
deleted file mode 100644
index 2aec03940114bfd44f9bbf11f84bca836deb02e1..0000000000000000000000000000000000000000
Binary files a/ObservationSim/MockObject/__pycache__/Catalog.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/MockObject/__pycache__/Galaxy.cpython-37.pyc b/ObservationSim/MockObject/__pycache__/Galaxy.cpython-37.pyc
deleted file mode 100755
index a7a55ed33245023046f50479936affc0e79ded97..0000000000000000000000000000000000000000
Binary files a/ObservationSim/MockObject/__pycache__/Galaxy.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/MockObject/__pycache__/Galaxy.cpython-38.pyc b/ObservationSim/MockObject/__pycache__/Galaxy.cpython-38.pyc
deleted file mode 100644
index a663442ad59b1faa9182ac376dddb0907d44c9b6..0000000000000000000000000000000000000000
Binary files a/ObservationSim/MockObject/__pycache__/Galaxy.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/MockObject/__pycache__/MockObject.cpython-37.pyc b/ObservationSim/MockObject/__pycache__/MockObject.cpython-37.pyc
deleted file mode 100755
index 6ef3686c2770f13abf0b85b10fd1b605790d656f..0000000000000000000000000000000000000000
Binary files a/ObservationSim/MockObject/__pycache__/MockObject.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/MockObject/__pycache__/MockObject.cpython-38.pyc b/ObservationSim/MockObject/__pycache__/MockObject.cpython-38.pyc
deleted file mode 100644
index 4fe9639e24557f0c661e2900a958163b47b1359b..0000000000000000000000000000000000000000
Binary files a/ObservationSim/MockObject/__pycache__/MockObject.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/MockObject/__pycache__/Quasar.cpython-37.pyc b/ObservationSim/MockObject/__pycache__/Quasar.cpython-37.pyc
deleted file mode 100755
index 91b82a6f9172601338d05093b0621898b6aaf7a9..0000000000000000000000000000000000000000
Binary files a/ObservationSim/MockObject/__pycache__/Quasar.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/MockObject/__pycache__/Quasar.cpython-38.pyc b/ObservationSim/MockObject/__pycache__/Quasar.cpython-38.pyc
deleted file mode 100644
index e3f189462efdcd7f72c2e6525c890b62a279f505..0000000000000000000000000000000000000000
Binary files a/ObservationSim/MockObject/__pycache__/Quasar.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/MockObject/__pycache__/SkybackgroundMap.cpython-38.pyc b/ObservationSim/MockObject/__pycache__/SkybackgroundMap.cpython-38.pyc
deleted file mode 100644
index f1c7e805f9c62cd7bfabda603c47a82abf5b13dd..0000000000000000000000000000000000000000
Binary files a/ObservationSim/MockObject/__pycache__/SkybackgroundMap.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/MockObject/__pycache__/Star.cpython-37.pyc b/ObservationSim/MockObject/__pycache__/Star.cpython-37.pyc
deleted file mode 100755
index 09bcdb35790e13dcf608c19011e8a8f799c3b58b..0000000000000000000000000000000000000000
Binary files a/ObservationSim/MockObject/__pycache__/Star.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/MockObject/__pycache__/Star.cpython-38.pyc b/ObservationSim/MockObject/__pycache__/Star.cpython-38.pyc
deleted file mode 100644
index 512a0c3e954127c7ac3736af3e1e13775937b48a..0000000000000000000000000000000000000000
Binary files a/ObservationSim/MockObject/__pycache__/Star.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/MockObject/__pycache__/__init__.cpython-37.pyc b/ObservationSim/MockObject/__pycache__/__init__.cpython-37.pyc
deleted file mode 100755
index 6444bf499819d696f87598e3ab7e1953726e7e70..0000000000000000000000000000000000000000
Binary files a/ObservationSim/MockObject/__pycache__/__init__.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/MockObject/__pycache__/__init__.cpython-38.pyc b/ObservationSim/MockObject/__pycache__/__init__.cpython-38.pyc
deleted file mode 100644
index 556670e94d9879c41d266edc520fa4d6e51685b4..0000000000000000000000000000000000000000
Binary files a/ObservationSim/MockObject/__pycache__/__init__.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/MockObject/__pycache__/_util.cpython-37.pyc b/ObservationSim/MockObject/__pycache__/_util.cpython-37.pyc
deleted file mode 100755
index d14ec0d5fc96276488f7f6e25c92500cbe51155c..0000000000000000000000000000000000000000
Binary files a/ObservationSim/MockObject/__pycache__/_util.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/MockObject/__pycache__/_util.cpython-38.pyc b/ObservationSim/MockObject/__pycache__/_util.cpython-38.pyc
deleted file mode 100644
index d1e72047f5abebeb3fa7f57bcda249428a676637..0000000000000000000000000000000000000000
Binary files a/ObservationSim/MockObject/__pycache__/_util.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/MockObject/_util.py b/ObservationSim/MockObject/_util.py
index f957e7b0b9367cc1e6178cbe43d259d9818c3cda..ffb978b910cc63d1bf281148554af29c54292d51 100755
--- a/ObservationSim/MockObject/_util.py
+++ b/ObservationSim/MockObject/_util.py
@@ -1,14 +1,13 @@
 import numpy as np
 import os
-from datetime import datetime
-from scipy.interpolate import InterpolatedUnivariateSpline, UnivariateSpline, interp1d
+from scipy.interpolate import InterpolatedUnivariateSpline, interp1d
 from scipy import interpolate
 from astropy.table import Table
 import galsim
 
-vc_A = 2.99792458e+18  # speed of light: A/s
-vc_m = 2.99792458e+8   # speed of light: m/s
-h_Plank = 6.626196e-27 # Plank constant: erg s
+VC_A = 2.99792458e+18  # speed of light: A/s
+VC_M = 2.99792458e+8   # speed of light: m/s
+H_PLANK = 6.626196e-27 # Plank constant: erg s
 
 def magToFlux(mag):
     """
@@ -132,13 +131,13 @@ def tflux(filt, sed, redshift=0.0, av=0.0, redden=0):
     
     sedxx = (sw.copy(), sf.copy())
     
-    sw = vc_A/sw
-    sf = sf*(vc_A/sw**2) # convert flux unit to erg/s/cm^s/Hz
+    sw = VC_A/sw
+    sf = sf*(VC_A/sw**2) # convert flux unit to erg/s/cm^s/Hz
     sw, sf = sw[::-1], sf[::-1]
     sfun = interp1d(sw, sf, kind='linear')
 
     fwave, fresp = filt[:,0], filt[:,1]
-    fwave = vc_A/fwave
+    fwave = VC_A/fwave
     fwave, fresp = fwave[::-1], fresp[::-1]
     tflux = sfun(fwave)
 
@@ -337,8 +336,7 @@ def reddening(sw, sf, av=0.0, model=0):
         flux  = sf*10**(-0.4*A_lambda)
 
     else:
-        print("!!! Please select a proper reddening model")
-        sys.exit(0)
+        raise ValueError("!!! Please select a proper reddening model")
 
     return flux
 
diff --git a/ObservationSim/ObservationSim.py b/ObservationSim/ObservationSim.py
index 563feaeee3c348efded5e3322d748b7be741060d..ece5dd7bee91451bc478b502d92c9aca9313be42 100755
--- a/ObservationSim/ObservationSim.py
+++ b/ObservationSim/ObservationSim.py
@@ -4,7 +4,7 @@ import os
 from ObservationSim.Config import ConfigDir, ChipOutput
 from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader
 from ObservationSim.Instrument import Telescope, Filter, FilterParam, FocalPlane, Chip
-from ObservationSim.MockObject import Catalog, MockObject, Star, Galaxy, Quasar, calculateSkyMap_split_g
+from ObservationSim.MockObject import calculateSkyMap_split_g
 from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp
 from ObservationSim._util import getShearFiled, makeSubDir_PointingList
 
@@ -17,7 +17,7 @@ import logging
 import psutil
 
 class Observation(object):
-    def __init__(self, config, work_dir=None, data_dir=None):
+    def __init__(self, config, Catalog, work_dir=None, data_dir=None):
         self.path_dict = ConfigDir(config=config, work_dir=work_dir, data_dir=data_dir)
         # self.config = ReadConfig(self.path_dict["config_file"])
         self.config = config
@@ -26,6 +26,7 @@ class Observation(object):
         self.filter_param = FilterParam(filter_dir=self.path_dict["filter_dir"]) # Currently the default values are hard coded in
         self.chip_list = []
         self.filter_list = []
+        self.Catalog = Catalog
 
         # if we want to apply field distortion?
         if self.config["ins_effects"]["field_dist"] == True:
@@ -85,7 +86,7 @@ class Observation(object):
 
         if pointing_type == 'MS':
             # Load catalogues and templates
-            self.cat = Catalog(config=self.config, chip=chip, cat_dir=cat_dir, sed_dir=sed_dir)
+            self.cat = self.Catalog(config=self.config, chip=chip, cat_dir=cat_dir, sed_dir=sed_dir)
             self.nobj = len(self.cat.objs)
 
             # Loop over objects
@@ -93,42 +94,29 @@ class Observation(object):
             bright_obj = 0
             dim_obj = 0
             for j in range(self.nobj):
-                # if j >= 20:
+                # if j >= 100:
                 #     break
                 obj = self.cat.objs[j]
+                if obj.type == 'star' and self.config["galaxy_only"]:
+                    continue
+                elif obj.type == 'galaxy' and self.config["star_only"]:
+                    continue
+                elif obj.type == 'quasar' and self.config["star_only"]:
+                    continue
 
-                # Load SED
-                if obj.type == 'star':
-                    normF = chip.normF_star
-                    if self.config["galaxy_only"]:
-                        continue
-                    try:
-                        obj.load_SED(
-                            survey_type=chip.survey_type, 
-                            normFilter=normF, 
-                            target_filt=filt,
-                            sed_lib=self.cat.tempSED_star)
-                    except Exception as e:
-                        print(e)
-                        continue
-                elif obj.type == 'galaxy': # or obj.type == quasar
-                    normF = chip.normF_galaxy
-                    if self.config["star_only"]:
-                        continue
-                    obj.load_SED(
-                        survey_type=chip.survey_type, 
-                        sed_templates=self.cat.tempSed_gal, 
-                        normFilter=normF, 
-                        target_filt=filt)
-                elif obj.type == 'quasar':
-                    normF = chip.normF_galaxy
-                    if self.config["star_only"]:
-                        continue
-                    obj.load_SED(
-                        survey_type=chip.survey_type, 
-                        sed_templates=self.cat.tempSed_gal, 
-                        normFilter=normF, 
-                        target_filt=filt)
+                # load SED
+                try:
+                    sed_data = self.cat.load_sed(obj)
+                    norm_filt = self.cat.load_norm_filt(obj)
+                    obj.sed, obj.param["mag_%s"%filt.filter_type] = self.cat.convert_sed(
+                        mag=obj.param["mag_use_normal"],
+                        sed=sed_data,
+                        target_filt=filt, 
+                        norm_filt=norm_filt,
+                    )
+                except Exception as e:
+                    print(e)
+                    continue
 
                 # Exclude very bright/dim objects (for now)
                 if filt.is_too_bright(mag=obj.getMagFilter(filt)):
@@ -156,6 +144,10 @@ class Observation(object):
                     except:
                         print("failed to load external shear.")
                         pass
+                elif self.config["shear_setting"]["shear_type"] == "catalog":
+                    pass
+                else:
+                    raise ValueError("Unknown shear input")
 
                 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:
@@ -191,7 +183,9 @@ class Observation(object):
                             g1=g1, 
                             g2=g2, 
                             exptime=exptime, 
-                            normFilter=normF)
+                            # normFilter=normF,
+                            normFilter=norm_filt,
+                            )
                     if isUpdated:
                         # TODO: add up stats
                         chip_output.cat_add_obj(obj, pos_img, pos_shear, g1, g2)
@@ -274,53 +268,6 @@ class Observation(object):
 
         print("check running:2: pointing-{:} chip-{:} pid-{:} memory-{:6.2}GB".format(pointing_ID, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) ), flush=True)
 
-
-    def runExposure(self, ra_cen=None, dec_cen=None, timestamp_obs=1621915200, pointing_ID=0, pointing_type='MS', img_rot=None, exptime=150., shear_cat_file=None, chips=None):
-
-        if (ra_cen == None) or (dec_cen == None):
-            ra_cen = self.config["obs_setting"]["ra_center"]
-            dec_cen = self.config["obs_setting"]["dec_center"]
-        if img_rot == None:
-            img_rot = self.config["obs_setting"]["image_rot"]
-
-        sub_img_dir, prefix = makeSubDir_PointingList(path_dict=self.path_dict, config=self.config, pointing_ID=pointing_ID)
-
-        # Loop over chips
-        for i in range(len(self.chip_list)):
-            chip = self.chip_list[i]
-            filt = self.filter_list[i]
-
-            # Only run a particular set of chips
-            if chips is not None:
-                if chip.chipID not in chips:
-                    continue
-
-            # Prepare output files
-            chip_output = ChipOutput(
-                config=self.config, 
-                focal_plane=self.focal_plane, 
-                chip=chip, 
-                filt=filt,  
-                exptime=exptime,
-                pointing_type=pointing_type,
-                pointing_ID=pointing_ID,
-                subdir=sub_img_dir,
-                prefix=prefix)
-            
-            self.runOneChip(
-                chip=chip, 
-                filt=filt, 
-                chip_output=chip_output, 
-                pointing_ID = pointing_ID,
-                ra_cen=ra_cen, 
-                dec_cen=dec_cen, 
-                img_rot=img_rot, 
-                exptime=exptime, 
-                timestamp_obs=timestamp_obs,
-                pointing_type=pointing_type,
-                cat_dir=self.path_dict["cat_dir"])
-            print("finished running chip#%d..."%(chip.chipID), flush=True)
-
     def runExposure_MPI_PointingList(self, ra_cen=None, dec_cen=None, pRange=None, timestamp_obs=np.array([1621915200]), pointing_type=np.array(['MS']), img_rot=None, exptime=150., shear_cat_file=None, chips=None, use_mpi=False):
         if use_mpi:
             comm = MPI.COMM_WORLD
diff --git a/ObservationSim/PSF/PSFInterp/__pycache__/PSFConfig.cpython-37.pyc b/ObservationSim/PSF/PSFInterp/__pycache__/PSFConfig.cpython-37.pyc
deleted file mode 100644
index cc0e9bfe56d95a888f852ac6dde70ffc950929b7..0000000000000000000000000000000000000000
Binary files a/ObservationSim/PSF/PSFInterp/__pycache__/PSFConfig.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/PSF/PSFInterp/__pycache__/PSFConfig.cpython-38.pyc b/ObservationSim/PSF/PSFInterp/__pycache__/PSFConfig.cpython-38.pyc
deleted file mode 100644
index 468f0294e213a98c00d3ca280a3712c028f70af0..0000000000000000000000000000000000000000
Binary files a/ObservationSim/PSF/PSFInterp/__pycache__/PSFConfig.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/PSF/PSFInterp/__pycache__/PSFInterp.cpython-38.pyc b/ObservationSim/PSF/PSFInterp/__pycache__/PSFInterp.cpython-38.pyc
deleted file mode 100644
index 324cb82e4d50e76b23b6a794ffb6f267f8534f8a..0000000000000000000000000000000000000000
Binary files a/ObservationSim/PSF/PSFInterp/__pycache__/PSFInterp.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/PSF/PSFInterp/__pycache__/PSFUtil.cpython-37.pyc b/ObservationSim/PSF/PSFInterp/__pycache__/PSFUtil.cpython-37.pyc
deleted file mode 100644
index b4bcdf22aeeba74c04062d529edd97f0259f77da..0000000000000000000000000000000000000000
Binary files a/ObservationSim/PSF/PSFInterp/__pycache__/PSFUtil.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/PSF/PSFInterp/__pycache__/PSFUtil.cpython-38.pyc b/ObservationSim/PSF/PSFInterp/__pycache__/PSFUtil.cpython-38.pyc
deleted file mode 100644
index 14727bc555a3d3e9ba3958ceafb9cf9c05d5b2aa..0000000000000000000000000000000000000000
Binary files a/ObservationSim/PSF/PSFInterp/__pycache__/PSFUtil.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/PSF/PSFInterp/__pycache__/__init__.cpython-38.pyc b/ObservationSim/PSF/PSFInterp/__pycache__/__init__.cpython-38.pyc
deleted file mode 100644
index 8fac220f537665f2e6018313955a178434ff72fa..0000000000000000000000000000000000000000
Binary files a/ObservationSim/PSF/PSFInterp/__pycache__/__init__.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/PSF/__pycache__/FieldDistortion.cpython-38.pyc b/ObservationSim/PSF/__pycache__/FieldDistortion.cpython-38.pyc
deleted file mode 100644
index 574f7adb3acb80c5dfcddcb139e383443c9a990a..0000000000000000000000000000000000000000
Binary files a/ObservationSim/PSF/__pycache__/FieldDistortion.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/PSF/__pycache__/PSFGauss.cpython-38.pyc b/ObservationSim/PSF/__pycache__/PSFGauss.cpython-38.pyc
deleted file mode 100644
index 526253816d6c36de939653768ad8c5a8825cc110..0000000000000000000000000000000000000000
Binary files a/ObservationSim/PSF/__pycache__/PSFGauss.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/PSF/__pycache__/PSFModel.cpython-37.pyc b/ObservationSim/PSF/__pycache__/PSFModel.cpython-37.pyc
deleted file mode 100755
index a1422b6d92c11f25aa9368e1322fe36c8bd0a4c2..0000000000000000000000000000000000000000
Binary files a/ObservationSim/PSF/__pycache__/PSFModel.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/PSF/__pycache__/PSFModel.cpython-38.pyc b/ObservationSim/PSF/__pycache__/PSFModel.cpython-38.pyc
deleted file mode 100644
index c7adf737f573d401e0b954effab184402bfb5686..0000000000000000000000000000000000000000
Binary files a/ObservationSim/PSF/__pycache__/PSFModel.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/PSF/__pycache__/__init__.cpython-37.pyc b/ObservationSim/PSF/__pycache__/__init__.cpython-37.pyc
deleted file mode 100755
index 0f542f8ae312d144aaf2a355417f41ccbde96425..0000000000000000000000000000000000000000
Binary files a/ObservationSim/PSF/__pycache__/__init__.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/PSF/__pycache__/__init__.cpython-38.pyc b/ObservationSim/PSF/__pycache__/__init__.cpython-38.pyc
deleted file mode 100644
index 953fc965bdbc4b947834ddd9befdfa08787b2090..0000000000000000000000000000000000000000
Binary files a/ObservationSim/PSF/__pycache__/__init__.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/__pycache__/ChipOutput.cpython-37.pyc b/ObservationSim/__pycache__/ChipOutput.cpython-37.pyc
deleted file mode 100755
index 7a0ebcf2761553dc61ad0ccc8bace2364705e729..0000000000000000000000000000000000000000
Binary files a/ObservationSim/__pycache__/ChipOutput.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/__pycache__/ChipOutput.cpython-38.pyc b/ObservationSim/__pycache__/ChipOutput.cpython-38.pyc
deleted file mode 100644
index 283aa63fcd34841cb43c2406210e8f7cd9f58007..0000000000000000000000000000000000000000
Binary files a/ObservationSim/__pycache__/ChipOutput.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/__pycache__/Config.cpython-37.pyc b/ObservationSim/__pycache__/Config.cpython-37.pyc
deleted file mode 100644
index 926b5309c212df628860deb69dec52b024174e62..0000000000000000000000000000000000000000
Binary files a/ObservationSim/__pycache__/Config.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/__pycache__/Config.cpython-38.pyc b/ObservationSim/__pycache__/Config.cpython-38.pyc
deleted file mode 100644
index 6088fea448ee536eb8272290efe5ab5d9cf18d60..0000000000000000000000000000000000000000
Binary files a/ObservationSim/__pycache__/Config.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/__pycache__/ObservationSim.cpython-37.pyc b/ObservationSim/__pycache__/ObservationSim.cpython-37.pyc
deleted file mode 100644
index ef6e9f154c0aecd31c992adbb277a8a90f007611..0000000000000000000000000000000000000000
Binary files a/ObservationSim/__pycache__/ObservationSim.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/__pycache__/ObservationSim.cpython-38.pyc b/ObservationSim/__pycache__/ObservationSim.cpython-38.pyc
deleted file mode 100644
index 4c27078d55d9396bd26e6cd82fa92d62f2c05050..0000000000000000000000000000000000000000
Binary files a/ObservationSim/__pycache__/ObservationSim.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/__pycache__/Pointing.cpython-38.pyc b/ObservationSim/__pycache__/Pointing.cpython-38.pyc
deleted file mode 100644
index 49945ad93f669fe1b8b86ad00e905644ced0c89d..0000000000000000000000000000000000000000
Binary files a/ObservationSim/__pycache__/Pointing.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/__pycache__/__init__.cpython-38.pyc b/ObservationSim/__pycache__/__init__.cpython-38.pyc
deleted file mode 100644
index 611fe47a01b30a6df8913fdc8beaf5b0bf349d99..0000000000000000000000000000000000000000
Binary files a/ObservationSim/__pycache__/__init__.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/__pycache__/_util.cpython-37.pyc b/ObservationSim/__pycache__/_util.cpython-37.pyc
deleted file mode 100755
index 85a261115d1398f4d2166edce8a1bd9d73bf7493..0000000000000000000000000000000000000000
Binary files a/ObservationSim/__pycache__/_util.cpython-37.pyc and /dev/null differ
diff --git a/ObservationSim/__pycache__/_util.cpython-38.pyc b/ObservationSim/__pycache__/_util.cpython-38.pyc
deleted file mode 100644
index 3f3b03b1491efd33880b8588ae9f3f65d17b535c..0000000000000000000000000000000000000000
Binary files a/ObservationSim/__pycache__/_util.cpython-38.pyc and /dev/null differ
diff --git a/ObservationSim/_util.py b/ObservationSim/_util.py
index 79bda9344db54331b41e01d0ee9e59ef74e40248..2e69c3c053c5a1e0b85a2d56c3991fc9cb803d2e 100755
--- a/ObservationSim/_util.py
+++ b/ObservationSim/_util.py
@@ -15,6 +15,83 @@ def parse_args():
     parser.add_argument('-w', '--work_dir', help='The path for output.')
     return parser.parse_args()
 
+def generate_pointings(config, pointing_filename=None, data_dir=None):
+    pRA = []
+    pDEC = []
+    if pointing_filename and data_dir:
+        pointing_file = os.path.join(data_dir, pointing_filename)
+        f = open(pointing_file, 'r')
+        for _ in range(1):
+            header = f.readline()
+        iline = 0
+        for line in f:
+            line = line.strip()
+            columns = line.split()
+            pRA.append(float(columns[0]))
+            pDEC.append(float(columns[1]))
+        f.close()
+    else:
+        pRA.append(config["obs_setting"]["ra_center"])
+        pDEC.append(config["obs_setting"]["dec_center"])
+    pRA = np.array(pRA)
+    pDEC = np.array(pDEC)
+
+    # Create calibration pointings
+    # NOTE: temporary implementation
+    ncal = config['obs_setting']['np_cal']
+    pointing_type = ['MS']*len(pRA)
+    pRA = np.append([pRA[0]]*ncal, pRA)
+    pDEC = np.append([pDEC[0]]*ncal, pDEC)
+    pointing_type = ['CAL']*ncal + pointing_type
+
+    # Calculate starting time(s)
+    # NOTE: temporary implementation
+    t0 = datetime(2021, 5, 25, 12, 0, 0)
+    t = datetime.timestamp(t0)
+    timestamp_obs = []
+    delta_t = 10 # Time elapsed between exposures (minutes)
+    for i in range(len(pointing_type)):
+        timestamp_obs.append(t)
+        if pointing_type[i] == 'CAL':
+            t += 3 * delta_t * 60 # 3 calibration exposures for each pointing
+        elif pointing_type[i] == 'MS':
+            t += delta_t * 60
+    timestamp_obs = np.array(timestamp_obs)
+    pointing_type = np.array(pointing_type)
+
+    if config['obs_setting']['run_pointings'] is None:
+        pRange = list(range(len(pRA)))
+    else:
+        ncal = config['obs_setting']['np_cal']
+        plist = config['obs_setting']['run_pointings']
+        pRange = list(range(ncal)) + [x + ncal for x in plist]
+    return pRA, pDEC, timestamp_obs, pointing_type, pRange
+
+def make_run_dirs(work_dir, run_name, nPointings, pRange=None):
+    if not os.path.exists(work_dir):
+        try:
+            os.makedirs(work_dir, exist_ok=True)
+        except OSError:
+            pass
+    imgDir = os.path.join(work_dir, run_name)
+    if not os.path.exists(imgDir):
+        try:
+            os.makedirs(imgDir, exist_ok=True)
+        except OSError:
+            pass
+    prefix = "MSC_"
+    for pointing_ID in range(nPointings):
+        if pRange is not None:
+            if pointing_ID not in pRange:
+                continue
+        fname=prefix + str(pointing_ID).rjust(7, '0')
+        subImgDir = os.path.join(imgDir, fname)
+        if not os.path.exists(subImgDir):
+            try:
+                os.makedirs(subImgDir, exist_ok=True)
+            except OSError:
+                pass
+
 def imgName(tt=0):
     ut = datetime.utcnow()
     eye, emo, eda, eho, emi, ese = str(ut.year), str(ut.month), str(ut.day), str(ut.hour), str(ut.minute), str(ut.second)
diff --git a/config/config_sim.yaml b/config/config_C3.yaml
similarity index 95%
rename from config/config_sim.yaml
rename to config/config_C3.yaml
index 595c11133bb6f12a45639e2b850a53a06b6e6354..f4bd7fe300e0267b324d27f2b89e456607e6771c 100644
--- a/config/config_sim.yaml
+++ b/config/config_C3.yaml
@@ -12,7 +12,7 @@
 # work_dir: "/public/home/fangyuedong/sim_code_release/CSST/test/"
 work_dir: "/public/home/fangyuedong/20211203/CSST/workplace/"
 data_dir: "/data/simudata/CSSOSDataProductsSims/data/"
-run_name: "TEST"
+run_name: "C3_20211216"
 
 # (Optional) a file of point list 
 # if you just want to run default pointing:
@@ -24,7 +24,10 @@ pointing_file: "pointing10_20210202.dat"
 # Whether to use MPI
 run_option:
   use_mpi: YES
-  n_cores: 40
+  # NOTE: "n_threads" paramters is currently not used in the backend
+  # simulation codes. It should be implemented later in the web frontend
+  # in order to config the number of threads to request from NAOC cluster
+  n_threads: 40
 
 # Output catalog only?
 # If yes, no imaging simulation will run
@@ -65,19 +68,19 @@ obs_setting:
 
   # Number of calibration pointings
   # Note: only valid when a pointing list is specified
-  np_cal: 0
+  np_cal: 3
 
   # Run specific pointing(s):
   # - 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, 6]
+  run_pointings: [1, 2, 3, 4, 5]
   
   # 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: [1, 25]
+  run_chips: [1, 26, 29, 6, 7, 22, 23, 19, 20]
 
 ###############################################
 # Input path setting
diff --git a/config/config_example.yaml b/config/config_example.yaml
new file mode 100644
index 0000000000000000000000000000000000000000..723a9d115d3e81b3d2a8ec3f69ad8b0e3bb7008d
--- /dev/null
+++ b/config/config_example.yaml
@@ -0,0 +1,223 @@
+---
+###############################################
+#
+#  Configuration file for CSST simulation
+#  CSST-Sim Group, 2021/10/07, version 0.3
+#
+###############################################
+
+# Base diretories and naming setup
+# 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/20211203/CSST/workplace/"
+data_dir: "/data/simudata/CSSOSDataProductsSims/data/"
+run_name: "example_20211217"
+
+# (Optional) a file of point list 
+# if you just want to run default pointing:
+# - pointing_dir: null
+# - pointing_file: null
+pointing_dir: null
+pointing_file: null
+
+# Whether to use MPI
+run_option:
+  use_mpi: YES
+  # NOTE: "n_threads" paramters is currently not used in the backend
+  # simulation codes. It should be implemented later in the web frontend
+  # in order to config the number of threads to request from NAOC cluster
+  n_threads: 40
+
+# Output catalog only?
+# If yes, no imaging simulation will run
+out_cat_only: NO
+
+# Only simulate stars?
+star_only: NO
+
+# Only simulate galaxies?
+galaxy_only: NO
+
+###############################################
+# Observation setting
+###############################################
+obs_setting:
+
+  # Options for survey types:
+  # "Photometric": simulate photometric chips only
+  # "Spectroscopic": simulate slitless spectroscopic chips only
+  # "All": simulate full focal plane
+  survey_type: "All"
+  
+  # Exposure time [seconds]
+  exp_time: 150.
+
+  # Observation starting date & time
+  # (Subject to change)
+  date_obs: "210525" # [yymmdd]
+  time_obs: "120000" # [hhmmss]
+
+  # Default Pointing [degrees]
+  # Note: NOT valid when a pointing list file is specified
+  ra_center: 60.9624
+  dec_center:  -41.5032
+
+  # Image rotation [degree]
+  image_rot: -113.4333
+
+  # Number of calibration pointings
+  np_cal: 2
+
+  # Run specific pointing(s):
+  # - 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: null
+  
+  # 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: [18]
+
+###############################################
+# Input path setting
+###############################################
+
+# Default path settings for WIDE survey simulation
+input_path:
+  cat_dir: "catalog_points_7degree2/point_RA60.9624_DE-41.5032/"
+  star_cat: "stars_ccd13_p_RA60.9624_DE-41.5032.hdf5"
+  galaxy_cat: null
+
+SED_templates_path:
+  star_SED: "SED_MMW_Gaia_Cluster_D20_SS.hdf5"
+  galaxy_SED: null
+
+# TODO: should the following path settings be made hidden from user?
+Efficiency_curve_path:
+  filter_eff: "Filters/"
+  ccd_eff: "Filter_CCD_Mirror/ccd/"
+  mirror_eff: "Filter_CCD_Mirror/mirror_ccdnote.txt"
+
+CR_data_path: "CRdata/"
+sky_data_path: "skybackground/sky_emiss_hubble_50_50_A.dat"
+
+SLS_path:
+  SLS_conf: "CONF"
+  SLS_norm: "normalize_filter"
+
+
+###############################################
+# 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: "csstPSFdata/CSSOS_psf_20210108/CSST_psf_ciomp_2p5um_cycle3_ccr90"
+
+  # path to field-distortion model
+  # Note: only valid when ins_effects: field_dist is "ON"
+  fd_path: "FieldDistModelv2.0.pickle"
+  
+  # sigma_spin: 0.0 # psf spin?
+
+###############################################
+# Shear setting
+###############################################
+
+shear_setting:
+  # Options to generate mock shear field:
+  # "constant": all galaxies are assigned a constant reduced shear
+  # "catalog": from catalog (not available yet)
+  # "extra": from seprate file
+  shear_type: "constant"
+
+  # For constant shear filed
+  reduced_g1: 0.026
+  reduced_g2: 0.015
+
+  # Representation of the shear vector?
+  reShear: "E"
+
+  # rotate galaxy ellipticity
+  rotateEll: 0. # [degree]
+
+  # Extra shear catalog
+  # (currently not used)
+  # shear_cat: "mockShear.cat"
+
+###############################################
+# Instrumental effects setting
+###############################################
+ins_effects:
+  # switches
+  field_dist:     ON  # Whether to add field distortions
+  add_back:       ON  # Whether to add sky background
+  add_dark:       ON  # Whether to add dark noise
+  add_readout:    ON  # Whether to add read-out (Gaussian) noise
+  add_bias:       ON  # Whether to add bias-level to images
+  shutter_effect: ON  # Whether to add shutter effect
+  flat_fielding:  ON  # Whether to add flat-fielding effect
+  prnu_effect:    ON  # Whether to add PRNU effect
+  non_linear:     OFF # 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:      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
+  add_deadpixels: ON  # Whether to add dead(dark) pixels
+  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-]
+  NBias:          1     # Number of bias frames to be exported for each exposure
+  NDark:          1     # Number of dark frames to be exported for each exposure
+  NFlat:          1     # Number of flat frames to be exported for each exposure
+
+###############################################
+# Output options
+###############################################
+output_setting:
+  readout16:      OFF # Whether to export as 16 channels (subimages) with pre- and over-scan
+  shutter_output: OFF # Whether to export shutter effect 16-bit image
+  bias_output:    ON  # Whether to export bias frames
+  dark_output:    ON  # Whether to export the combined dark current files
+  flat_output:    ON  # Whether to export the combined flat-fielding files
+  prnu_output:    OFF # Whether to export the PRNU (pixel-to-pixel flat-fielding) files
+
+###############################################
+# Random seeds
+###############################################
+random_seeds:
+  seed_Av:              121212    # Seed for generating random intrinsic extinction
+  seed_poisson:         20210601  # Seed for Poisson noise
+  seed_CR:              20210317  # Seed for generating random cosmic ray maps
+  seed_flat:            20210101  # Seed for generating random flat fields
+  seed_prnu:            20210102  # Seed for photo-response non-uniformity
+  seed_gainNonUniform:  20210202  # Seed for gain nonuniformity
+  seed_biasNonUniform:  20210203  # Seed for bias nonuniformity
+  seed_rnNonUniform:    20210204  # Seed for readout-noise nonuniformity
+  seed_badcolumns:      20240309  # Seed for bad columns
+  seed_defective:       20210304  # Seed for defective (bad) pixels
+  seed_readout:         20210601  # Seed for read-out gaussian noise
+...
\ No newline at end of file
diff --git a/example.sh b/example.sh
new file mode 100755
index 0000000000000000000000000000000000000000..c37eecc113aac3e94f30ab7beecd80f0a59e7713
--- /dev/null
+++ b/example.sh
@@ -0,0 +1,5 @@
+#!/bin/bash
+
+date
+
+python /public/home/fangyuedong/20211203/CSST/run_example.py config_example.yaml -c /public/home/fangyuedong/20211203/CSST/config
diff --git a/run.pbs b/mpi_example.pbs
similarity index 82%
rename from run.pbs
rename to mpi_example.pbs
index 497e0c6a07c9c321d418161c8b2374b53a0bae8c..23730fb3147273e3c0e6e174103f720aa912b880 100755
--- a/run.pbs
+++ b/mpi_example.pbs
@@ -16,4 +16,4 @@ NP=40
 date
 echo $NP
 
-mpirun -np $NP python /public/home/fangyuedong/20211203/CSST/runExposure.py config_sim.yaml -c /public/home/fangyuedong/20211203/CSST/config
+mpirun -np $NP python /public/home/fangyuedong/20211203/CSST/run_sim.py config_example.yaml -c /public/home/fangyuedong/20211203/CSST/config
diff --git a/runExposure.py b/runExposure.py
deleted file mode 100755
index c77e1cd59efbdc8681154fe4ab2b5a8d4a970482..0000000000000000000000000000000000000000
--- a/runExposure.py
+++ /dev/null
@@ -1,125 +0,0 @@
-from ObservationSim.ObservationSim import Observation
-from ObservationSim._util import parse_args
-from datetime import datetime
-import os
-import numpy as np
-import galsim
-import yaml
-
-import gc
-gc.enable()
-
-def Pointing(config, pointing_filename=None, data_dir=None):
-    pRA = []
-    pDEC = []
-    if pointing_filename and data_dir:
-        pointing_file = os.path.join(data_dir, pointing_filename)
-        f = open(pointing_file, 'r')
-        for _ in range(1):
-            header = f.readline()
-        iline = 0
-        for line in f:
-            line = line.strip()
-            columns = line.split()
-            pRA.append(float(columns[0]))
-            pDEC.append(float(columns[1]))
-        f.close()
-    else:
-        pRa.append(config["obs_setting"]["ra_center"])
-        pDec.append(config["obs_setting"]["dec_center"])
-    pRA = np.array(pRA)
-    pDEC = np.array(pDEC)
-
-    # Create calibration pointings
-    # NOTE: temporary implementation
-    ncal = config['obs_setting']['np_cal']
-    pointing_type = ['MS']*len(pRA)
-    pRA = np.append([pRA[0]]*ncal, pRA)
-    pDEC = np.append([pDEC[0]]*ncal, pDEC)
-    pointing_type = ['CAL']*ncal + pointing_type
-
-    # Calculate starting time(s)
-    # NOTE: temporary implementation
-    t0 = datetime(2021, 5, 25, 12, 0, 0)
-    t = datetime.timestamp(t0)
-    timestamp_obs = []
-    delta_t = 10 # Time elapsed between exposures (minutes)
-    for i in range(len(pointing_type)):
-        timestamp_obs.append(t)
-        if pointing_type[i] == 'CAL':
-            t += 3 * delta_t * 60 # 3 calibration exposures for each pointing
-        elif pointing_type[i] == 'MS':
-            t += delta_t * 60
-    timestamp_obs = np.array(timestamp_obs)
-    pointing_type = np.array(pointing_type)
-    return pRA, pDEC, timestamp_obs, pointing_type
-
-def MakeDirectories(work_dir, run_name, nPointings, pRange=None):
-    if not os.path.exists(work_dir):
-        try:
-            os.makedirs(work_dir, exist_ok=True)
-        except OSError:
-            pass
-    imgDir = os.path.join(work_dir, run_name)
-    if not os.path.exists(imgDir):
-        try:
-            os.makedirs(imgDir, exist_ok=True)
-        except OSError:
-            pass
-    prefix = "MSC_"
-    for pointing_ID in range(nPointings):
-        if pRange is not None:
-            if pointing_ID not in pRange:
-                continue
-        fname=prefix + str(pointing_ID).rjust(7, '0')
-        subImgDir = os.path.join(imgDir, fname)
-        if not os.path.exists(subImgDir):
-            try:
-                os.makedirs(subImgDir, exist_ok=True)
-            except OSError:
-                pass
-
-def runSim():
-    """
-    Main simulation call.
-    """
-    args = parse_args()
-    if args.config_dir is None:
-        args.config_dir = ''
-    args.config_dir = os.path.abspath(args.config_dir)
-    args.config_file = os.path.join(args.config_dir, args.config_file)
-    with open(args.config_file, "r") as stream:
-        try:
-            config = yaml.safe_load(stream)
-            for key, value in config.items():
-                print (key + " : " + str(value))
-        except yaml.YAMLError as exc:
-            print(exc)
-    config["obs_setting"]["image_rot"] = float(config["obs_setting"]["image_rot"])*galsim.degrees
-
-    if args.data_dir is not None:
-        config['data_dir'] = args.data_dir
-    if args.work_dir is not None:
-        config['work_dir'] = args.work_dir
-
-    pRA, pDEC, timestamp_obs, pointing_type = Pointing(config=config, pointing_filename=config['pointing_file'], data_dir=config['pointing_dir'])
-
-    MakeDirectories(work_dir=config['work_dir'], run_name=config['run_name'], nPointings=len(pRA), pRange=config['obs_setting']['run_pointings'])
-
-    obs = Observation(config=config, work_dir=config['work_dir'], data_dir=config['data_dir'])
-    if config["pointing_file"] is None:
-        obs.runExposure(chips=config["obs_setting"]["run_chips"])
-    else:
-        obs.runExposure_MPI_PointingList(
-            ra_cen=pRA, 
-            dec_cen=pDEC, 
-            pRange=config["obs_setting"]["run_pointings"], 
-            timestamp_obs=timestamp_obs,
-            pointing_type=pointing_type,
-            exptime=config["obs_setting"]["exp_time"],
-            use_mpi=config["run_option"]["use_mpi"],
-            chips=config["obs_setting"]["run_chips"]
-            )
-
-if __name__=='__main__':
-    runSim()
diff --git a/run_C3.pbs b/run_C3.pbs
new file mode 100755
index 0000000000000000000000000000000000000000..d421f6ab63863cab061294fcecee27654d597626
--- /dev/null
+++ b/run_C3.pbs
@@ -0,0 +1,19 @@
+#!/bin/bash
+
+#PBS -N SIMS
+##PBS -l walltime=70:00:00
+
+##mpdallexit
+##mpdboot -n 10 -f  ./mpd.hosts
+
+##PBS -j oe
+#PBS -l nodes=comput108:ppn=80
+#####PBS -q longq
+#PBS -q batch
+#PBS -u fangyuedong
+
+NP=80
+date
+echo $NP
+
+mpirun -np $NP python /public/home/fangyuedong/20211203/CSST/run_sim.py config_C3.yaml -c /public/home/fangyuedong/20211203/CSST/config
diff --git a/run_sim.py b/run_sim.py
new file mode 100755
index 0000000000000000000000000000000000000000..796e7350bc3da8b8e72a226b1644ba2f79413cc4
--- /dev/null
+++ b/run_sim.py
@@ -0,0 +1,65 @@
+from ObservationSim.ObservationSim import Observation
+from ObservationSim._util import parse_args, generate_pointings, make_run_dirs
+from Catalog.Catalog_example import Catalog_example
+import os
+import galsim
+import yaml
+
+import gc
+gc.enable()
+
+def runSim(Catalog):
+    """
+    Main method for simulation call.
+    """
+    args = parse_args()
+    if args.config_dir is None:
+        args.config_dir = ''
+    args.config_dir = os.path.abspath(args.config_dir)
+    args.config_file = os.path.join(args.config_dir, args.config_file)
+    with open(args.config_file, "r") as stream:
+        try:
+            config = yaml.safe_load(stream)
+            for key, value in config.items():
+                print (key + " : " + str(value))
+        except yaml.YAMLError as exc:
+            print(exc)
+    config["obs_setting"]["image_rot"] = float(config["obs_setting"]["image_rot"])*galsim.degrees
+
+    # Overwrite the data and working directories
+    # if they are specified by the command line inputs
+    if args.data_dir is not None:
+        config['data_dir'] = args.data_dir
+    if args.work_dir is not None:
+        config['work_dir'] = args.work_dir
+
+    # Generate lists of RA, DEC, time_stamps, pointing_type (MS or CAL), and pRange 
+    # (selections of pointings to run) based on the input pointing list (or default 
+    # pointing RA, DEC) and "config["obs_setting"]["run_pointings"]".
+    # "config['obs_setting']['np_cal']"" is the number of CAL pointings which will be 
+    # appended to the front of all above lists.
+    # NOTE: the implementation of gerenating time_stamps is temporary.
+    pRA, pDEC, timestamp_obs, pointing_type, pRange = generate_pointings(config=config, pointing_filename=config['pointing_file'], data_dir=config['pointing_dir'])
+
+    # Make the main output directories
+    make_run_dirs(work_dir=config['work_dir'], run_name=config['run_name'], nPointings=len(pRA), pRange=pRange)
+
+    # Initialize the simulation
+    obs = Observation(config=config, Catalog=Catalog, work_dir=config['work_dir'], data_dir=config['data_dir'])
+    
+    # Run simulation
+    obs.runExposure_MPI_PointingList(
+        ra_cen=pRA, 
+        dec_cen=pDEC, 
+        pRange=pRange, 
+        timestamp_obs=timestamp_obs,
+        pointing_type=pointing_type,
+        exptime=config["obs_setting"]["exp_time"],
+        use_mpi=config["run_option"]["use_mpi"],
+        chips=config["obs_setting"]["run_chips"])
+
+if __name__=='__main__':
+    runSim(Catalog=Catalog_example)
+    # To run cycle-3 simulation
+    # from Catalog.C3Catalog import C3Catalog
+    # runSim(Catalog=C3Catalog)