Chip.py

import galsim
import os
import numpy as np
import pickle
import json
import ObservationSim.Instrument._util as _util
from astropy.table import Table
from numpy.random import Generator, PCG64
from astropy.io import fits
from datetime import datetime

from ObservationSim.Instrument.Chip import Effects as effects
from ObservationSim.Instrument.FocalPlane import FocalPlane
from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader
from ObservationSim.Instrument._util import rotate_conterclockwise
from ObservationSim.Instrument.Chip import ChipUtils as chip_utils

from ObservationSim.Instrument.Chip.libCTI.CTI_modeling import CTI_sim

try:
    import importlib.resources as pkg_resources
except ImportError:
    # Try backported to PY<37 'importlib_resources'
    import importlib_resources as pkg_resources


class Chip(FocalPlane):
    def __init__(self, chipID, ccdEffCurve_dir=None, CRdata_dir=None, sls_dir=None, config=None, treering_func=None, logger=None):
        # Get focal plane (instance of paraent class) info
        super().__init__()
        self.nsecy = 2
        self.nsecx = 8
        self.gain_channel = np.ones(self.nsecy * self.nsecx)
        self._set_attributes_from_config(config)

        self.logger = logger

        # A chip ID must be assigned
        self.chipID = int(chipID)
        self.chip_name = str(chipID).rjust(2, '0')

        # Get corresponding filter info
        self.filter_id, self.filter_type = self.getChipFilter()
        self.survey_type = self._getSurveyType()

        if self.filter_type != "FGS":
            self._getChipRowCol()

        # Set the relavent specs for detectors
        try:
            with pkg_resources.files('ObservationSim.Instrument.data.ccd').joinpath("chip_definition.json") as chip_definition:
                with open(chip_definition, "r") as f:
                    chip_dict = json.load(f)[str(self.chipID)]
        except AttributeError:
            with pkg_resources.path('ObservationSim.Instrument.data.ccd', "chip_definition.json") as chip_definition:
                with open(chip_definition, "r") as f:
                    chip_dict = json.load(f)[str(self.chipID)]
        for key in chip_dict:
            setattr(self, key, chip_dict[key])

        self.fdModel = None
        if self.filter_type == "FGS":
            fgs_name = self.chip_name[0:4]
            try:
                with pkg_resources.files('ObservationSim.Instrument.data.field_distortion').joinpath("FieldDistModelGlobal_pr4_%s.pickle" % (fgs_name.lower())) as field_distortion:
                    with open(field_distortion, "rb") as f:
                        self.fdModel = pickle.load(f)
            except AttributeError:
                with pkg_resources.path('ObservationSim.Instrument.data.field_distortion', "FieldDistModelGlobal_pr4_%s.pickle" % (fgs_name.lower())) as field_distortion:
                    with open(field_distortion, "rb") as f:
                        self.fdModel = pickle.load(f)
        else:
            # Get the corresponding field distortion model
            try:
                with pkg_resources.files('ObservationSim.Instrument.data.field_distortion').joinpath("FieldDistModel_v2.0.pickle") as field_distortion:
                    with open(field_distortion, "rb") as f:
                        self.fdModel = pickle.load(f)
            except AttributeError:
                with pkg_resources.path('ObservationSim.Instrument.data.field_distortion', "FieldDistModelGlobal_mainFP_v1.0.pickle") as field_distortion:
                    with open(field_distortion, "rb") as f:
                        self.fdModel = pickle.load(f)

        # Get boundary (in pix)
        self.bound = self.getChipLim()

        self.ccdEffCurve_dir = ccdEffCurve_dir
        self.CRdata_dir = CRdata_dir

        slsconfs = chip_utils.getChipSLSConf(chipID=self.chipID)
        if np.size(slsconfs) == 1:
            try:
                with pkg_resources.files('ObservationSim.Instrument.data.sls_conf').joinpath(slsconfs) as conf_path:
                    self.sls_conf = str(conf_path)
            except AttributeError:
                with pkg_resources.path('ObservationSim.Instrument.data.sls_conf', slsconfs) as conf_path:
                    self.sls_conf = str(conf_path)
        else:
            # self.sls_conf = [os.path.join(self.sls_dir, slsconfs[0]), os.path.join(self.sls_dir, slsconfs[1])]
            self.sls_conf = []
            try:
                with pkg_resources.files('ObservationSim.Instrument.data.sls_conf').joinpath(slsconfs[0]) as conf_path:
                    self.sls_conf.append(str(conf_path))
            except AttributeError:
                with pkg_resources.path('ObservationSim.Instrument.data.sls_conf', slsconfs[0]) as conf_path:
                    self.sls_conf.append(str(conf_path))
            try:
                with pkg_resources.files('ObservationSim.Instrument.data.sls_conf').joinpath(slsconfs[1]) as conf_path:
                    self.sls_conf.append(str(conf_path))
            except AttributeError:
                with pkg_resources.path('ObservationSim.Instrument.data.sls_conf', slsconfs[1]) as conf_path:
                    self.sls_conf.append(str(conf_path))

        self.effCurve = self._getChipEffCurve(self.filter_type)
        self._getCRdata()

        # # Define the sensor model
        self.sensor = galsim.Sensor()

        self.flat_cube = None  # for spectroscopic flat field cube simulation

    def _set_attributes_from_config(self, config):
        # Default setting
        self.read_noise = 5.0   # e/pix
        self.dark_noise = 0.02  # e/pix/s
        self.rotate_angle = 0.
        self.overscan = 1000
        # Override default values
        # for key in ["gain", "bias_level, dark_exptime", "flat_exptime", "readout_time", "full_well", "read_noise", "dark_noise", "overscan"]:
        #     if key in config["ins_effects"]:
        #         setattr(self, key, config["ins_effects"][key])

    def _getChipRowCol(self):
        self.rowID, self.colID = self.getChipRowCol(self.chipID)

    def getChipRowCol(self, chipID):
        rowID = ((chipID - 1) % 5) + 1
        colID = 6 - ((chipID - 1) // 5)
        return rowID, colID

    def _getSurveyType(self):
        if self.filter_type in _util.SPEC_FILTERS:
            return "spectroscopic"
        elif self.filter_type in _util.PHOT_FILTERS:
            return "photometric"
        # elif self.filter_type in ["FGS"]:
        #     return "FGS"

    def _getChipEffCurve(self, filter_type):
        # CCD efficiency curves
        if filter_type in ['NUV', 'u', 'GU']:
            filename = 'UV0.txt'
        if filter_type in ['g', 'r', 'GV', 'FGS']:
            # TODO, need to switch to the right efficiency curvey for FGS CMOS
            filename = 'Astro_MB.txt'
        if filter_type in ['i', 'z', 'y', 'GI']:
            filename = 'Basic_NIR.txt'
        try:
            with pkg_resources.files('ObservationSim.Instrument.data.ccd').joinpath(filename) as ccd_path:
                table = Table.read(ccd_path, format='ascii')
        except AttributeError:
            with pkg_resources.path('ObservationSim.Instrument.data.ccd', filename) as ccd_path:
                table = Table.read(ccd_path, format='ascii')
        throughput = galsim.LookupTable(
            x=table['col1'], f=table['col2'], interpolant='linear')
        bandpass = galsim.Bandpass(throughput, wave_type='nm')
        return bandpass

    def _getCRdata(self):
        try:
            with pkg_resources.files('ObservationSim.Instrument.data').joinpath("wfc-cr-attachpixel.dat") as cr_path:
                self.attachedSizes = np.loadtxt(cr_path)
        except AttributeError:
            with pkg_resources.path('ObservationSim.Instrument.data', "wfc-cr-attachpixel.dat") as cr_path:
                self.attachedSizes = np.loadtxt(cr_path)

    # def loadSLSFLATCUBE(self, flat_fn='flat_cube.fits'):
    #     try:
    #         with pkg_resources.files('ObservationSim.Instrument.data').joinpath(flat_fn) as data_path:
    #             flat_fits = fits.open(data_path, ignore_missing_simple=True)
    #     except AttributeError:
    #         with pkg_resources.path('ObservationSim.Instrument.data', flat_fn) as data_path:
    #             flat_fits = fits.open(data_path, ignore_missing_simple=True)

    #     fl = len(flat_fits)
    #     fl_sh = flat_fits[0].data.shape
    #     assert fl == 4, 'FLAT Field Cube is Not 4 layess!!!!!!!'
    #     self.flat_cube = np.zeros([fl, fl_sh[0], fl_sh[1]])
    #     for i in np.arange(0, fl, 1):
    #         self.flat_cube[i] = flat_fits[i].data

    def getChipFilter(self, chipID=None):
        """Return the filter index and type for a given chip #(chipID)
        """
        filter_type_list = _util.ALL_FILTERS
        if chipID == None:
            chipID = self.chipID

        # updated configurations
        if chipID > 42 or chipID < 1:
            raise ValueError("!!! Chip ID: [1,42]")
        if chipID in [6, 15, 16, 25]:
            filter_type = "y"
        if chipID in [11, 20]:
            filter_type = "z"
        if chipID in [7, 24]:
            filter_type = "i"
        if chipID in [14, 17]:
            filter_type = "u"
        if chipID in [9, 22]:
            filter_type = "r"
        if chipID in [12, 13, 18, 19]:
            filter_type = "NUV"
        if chipID in [8, 23]:
            filter_type = "g"
        if chipID in [1, 10, 21, 30]:
            filter_type = "GI"
        if chipID in [2, 5, 26, 29]:
            filter_type = "GV"
        if chipID in [3, 4, 27, 28]:
            filter_type = "GU"
        if chipID in range(31, 43):
            filter_type = 'FGS'
        filter_id = filter_type_list.index(filter_type)

        return filter_id, filter_type

    def getChipLim(self, chipID=None):
        """Calculate the edges in pixel for a given CCD chip on the focal plane
        NOTE: There are 5*4 CCD chips in the focus plane for photometric / spectroscopic observation.
        Parameters:
            chipID:         int
                            the index of the chip
        Returns:
            A galsim BoundsD object
        """
        xmin, xmax, ymin, ymax = 1e10, -1e10, 1e10, -1e10
        xcen = self.x_cen / self.pix_size
        ycen = self.y_cen / self.pix_size
        sign_x = [-1., 1., -1., 1.]
        sign_y = [-1., -1., 1., 1.]
        for i in range(4):
            x = xcen + sign_x[i] * self.npix_x / 2.
            y = ycen + sign_y[i] * self.npix_y / 2.
            x, y = _util.rotate_conterclockwise(
                x0=xcen, y0=ycen, x=x, y=y, angle=self.rotate_angle)
            xmin, xmax = min(xmin, x), max(xmax, x)
            ymin, ymax = min(ymin, y), max(ymax, y)
        return galsim.BoundsD(xmin, xmax, ymin, ymax)

    def getSkyCoverage(self, wcs):
        # print("In getSkyCoverage: xmin = %.3f, xmax = %.3f, ymin = %.3f, ymax = %.3f"%(self.bound.xmin, self.bound.xmax, self.bound.ymin, self.bound.ymax))
        return super().getSkyCoverage(wcs, self.bound.xmin, self.bound.xmax, self.bound.ymin, self.bound.ymax)

    def getSkyCoverageEnlarged(self, wcs, margin=0.5):
        """The enlarged sky coverage of the chip
        """
        margin /= 60.0
        bound = self.getSkyCoverage(wcs)
        return galsim.BoundsD(bound.xmin - margin, bound.xmax + margin, bound.ymin - margin, bound.ymax + margin)

    def isContainObj(self, ra_obj=None, dec_obj=None, x_image=None, y_image=None, wcs=None, margin=1):
        # magin in number of pix
        if (ra_obj is not None) and (dec_obj is not None):
            if wcs is None:
                wcs = self.img.wcs
            pos_obj = wcs.toImage(galsim.CelestialCoord(
                ra=ra_obj*galsim.degrees, dec=dec_obj*galsim.degrees))
            x_image, y_image = pos_obj.x, pos_obj.y
        elif (x_image is None) or (y_image is None):
            raise ValueError(
                "Either (ra_obj, dec_obj) or (x_image, y_image) should be given")

        xmin, xmax = self.bound.xmin - margin, self.bound.xmax + margin
        ymin, ymax = self.bound.ymin - margin, self.bound.ymax + margin
        if (x_image - xmin) * (x_image - xmax) > 0.0 or (y_image - ymin) * (y_image - ymax) > 0.0:
            return False
        return True

    def getChipNoise(self, exptime=150.0):
        noise = self.dark_noise * exptime + self.read_noise**2
        return noise

    # def addEffects(self, config, img, chip_output, filt, ra_cen, dec_cen, img_rot, exptime=150., pointing_ID=0, timestamp_obs=1621915200, pointing_type='SCI', sky_map=None, post_flash_map=None, tel=None, logger=None):
    #     # Set random seeds
    #     SeedGainNonuni = int(config["random_seeds"]["seed_gainNonUniform"])
    #     SeedBiasNonuni = int(config["random_seeds"]["seed_biasNonUniform"])
    #     SeedRnNonuni = int(config["random_seeds"]["seed_rnNonUniform"])
    #     SeedBadColumns = int(config["random_seeds"]["seed_badcolumns"])
    #     SeedDefective = int(config["random_seeds"]["seed_defective"])
    #     SeedCosmicRay = int(config["random_seeds"]["seed_CR"])
    #     fullwell = int(self.full_well)
    #     if config["ins_effects"]["add_hotpixels"] == True:
    #         BoolHotPix = True
    #     else:
    #         BoolHotPix = False
    #     if config["ins_effects"]["add_deadpixels"] == True:
    #         BoolDeadPix = True
    #     else:
    #         BoolDeadPix = False
    #     self.logger = logger

    #     # Get Poisson noise generator
    #     rng_poisson, poisson_noise = chip_utils.get_poisson(
    #         seed=int(config["random_seeds"]["seed_poisson"]) + pointing_ID*30 + self.chipID, sky_level=0.)

    #     # Add sky background
    #     if config["ins_effects"]["add_back"] == True:
    #         img, sky_map = chip_utils.add_sky_background(
    #             img=img, filt=filt, exptime=exptime, sky_map=sky_map, tel=tel)
    #         del sky_map

    #     # Apply flat-field large scale structure for one chip
    #     if config["ins_effects"]["flat_fielding"] == True:
    #         chip_utils.log_info(
    #             msg="  Creating and applying Flat-Fielding", logger=self.logger)
    #         chip_utils.log_info(msg=str(img.bounds), logger=self.logger)
    #         flat_img, flat_normal = chip_utils.get_flat(
    #             img=img, seed=int(config["random_seeds"]["seed_flat"]))
    #         if self.survey_type == "photometric":
    #             img *= flat_normal
    #         del flat_normal
    #         if config["output_setting"]["flat_output"] == False:
    #             del flat_img

    #     if post_flash_map is not None:
    #         img = img + post_flash_map

    #     # Apply Shutter-effect for one chip
    #     if config["ins_effects"]["shutter_effect"] == True:
    #         chip_utils.log_info(
    #             msg="  Apply shutter effect", logger=self.logger)
    #         # shutter effect normalized image for this chip
    #         shuttimg = effects.ShutterEffectArr(
    #             img, t_shutter=1.3, dist_bearing=735, dt=1E-3)
    #         if self.survey_type == "photometric":
    #             img *= shuttimg
    #         # output 16-bit shutter effect image with pixel value <=65535
    #         if config["output_setting"]["shutter_output"] == True:
    #             shutt_gsimg = galsim.ImageUS(shuttimg*6E4)
    #             shutt_gsimg.write("%s/ShutterEffect_%s_1.fits" %
    #                               (chip_output.subdir, self.chipID))
    #             del shutt_gsimg
    #         del shuttimg
    #     # # Add Poisson noise to the resulting images
    #     # # (NOTE): this can only applied to the slitless image
    #     # # since it dose not use photon shooting to draw stamps
    #     # if self.survey_type == "spectroscopic":
    #     #     img.addNoise(poisson_noise)

    #     # Add cosmic-rays
    #     if config["ins_effects"]["cosmic_ray"] == True and pointing_type == 'SCI':
    #         chip_utils.log_info(msg="  Adding Cosmic-Ray", logger=self.logger)
    #         img, crmap_gsimg, cr_event_num = chip_utils.add_cosmic_rays(img=img, chip=self, exptime=exptime,
    #                                                                     seed=SeedCosmicRay+pointing_ID*30+self.chipID)
    #         chip_utils.outputCal(
    #             chip=self,
    #             img=crmap_gsimg,
    #             ra_cen=ra_cen,
    #             dec_cen=dec_cen,
    #             img_rot=img_rot,
    #             im_type='CRS',
    #             pointing_ID=pointing_ID,
    #             output_dir=chip_output.subdir,
    #             exptime=exptime,
    #             project_cycle=config["project_cycle"],
    #             run_counter=config["run_counter"],
    #             timestamp=timestamp_obs)
    #         del crmap_gsimg

    #     # Apply PRNU effect and output PRNU flat file:
    #     if config["ins_effects"]["prnu_effect"] == True:
    #         chip_utils.log_info(
    #             msg="  Applying PRNU effect", logger=self.logger)
    #         img, prnu_img = chip_utils.add_PRNU(img=img, chip=self,
    #                                             seed=int(config["random_seeds"]["seed_prnu"]+self.chipID))
    #         if config["output_setting"]["prnu_output"] == True:
    #             prnu_img.write("%s/FlatImg_PRNU_%s.fits" %
    #                            (chip_output.subdir, self.chipID))
    #         if config["output_setting"]["flat_output"] == False:
    #             del prnu_img

    #     # # Add dark current
    #     # if config["ins_effects"]["add_dark"] == True:
    #     #     dark_noise = galsim.DeviateNoise(galsim.PoissonDeviate(rng_poisson, self.dark_noise*(exptime+0.5*self.readout_time)))
    #     #     img.addNoise(dark_noise)

    #     # Add dark current & Poisson noise
    #     InputDark = False
    #     if config["ins_effects"]["add_dark"] == True:
    #         if InputDark:
    #             img = chip_utils.add_inputdark(
    #                 img=img, chip=self, exptime=exptime)
    #         else:
    #             img, _ = chip_utils.add_poisson(
    #                 img=img, chip=self, exptime=exptime, poisson_noise=poisson_noise)
    #     else:
    #         img, _ = chip_utils.add_poisson(
    #             img=img, chip=self, exptime=exptime, poisson_noise=poisson_noise, dark_noise=0.)

    #     # Add diffusion & brighter-fatter effects
    #     if config["ins_effects"]["bright_fatter"] == True:
    #         img = chip_utils.add_brighter_fatter(img=img)

    #     # Add Hot Pixels or/and Dead Pixels
    #     rgbadpix = Generator(PCG64(int(SeedDefective+self.chipID)))
    #     badfraction = 5E-5*(rgbadpix.random()*0.5+0.7)
    #     img = effects.DefectivePixels(img, IfHotPix=BoolHotPix, IfDeadPix=BoolDeadPix,
    #                                   fraction=badfraction, seed=SeedDefective+self.chipID, biaslevel=0)

    #     # Apply Bad lines
    #     if config["ins_effects"]["add_badcolumns"] == True:
    #         img = effects.BadColumns(
    #             img, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger)

    #     # Apply Nonlinearity on the chip image
    #     if config["ins_effects"]["non_linear"] == True:
    #         chip_utils.log_info(
    #             msg="  Applying Non-Linearity on the chip image", logger=self.logger)
    #         img = effects.NonLinearity(GSImage=img, beta1=5.e-7, beta2=0)

    #     # Apply CCD Saturation & Blooming
    #     if config["ins_effects"]["saturbloom"] == True:
    #         chip_utils.log_info(
    #             msg="  Applying CCD Saturation & Blooming", logger=self.logger)
    #         img = effects.SaturBloom(
    #             GSImage=img, nsect_x=1, nsect_y=1, fullwell=fullwell)

    #     # Apply CTE Effect
    #     # if config["ins_effects"]["cte_trail"] == True:
    #     # chip_utils.log_info(msg="  Apply CTE Effect", logger=self.logger)
    #     # img = effects.CTE_Effect(GSImage=img, threshold=27)

    #     pre1 = self.prescan_x  # 27
    #     over1 = self.overscan_x  # 71
    #     pre2 = self.prescan_y  # 0 #4
    #     over2 = self.overscan_y  # 84 #80

    #     if config["ins_effects"]["cte_trail"] == True:
    #         chip_utils.log_info(msg="  Apply CTE Effect", logger=self.logger)
    #         # img = effects.CTE_Effect(GSImage=img, threshold=27)
    #         # CTI_modeling
    #         # 2*8 -> 1*16 img-layout
    #         img = chip_utils.formatOutput(GSImage=img)
    #         self.nsecy = 1
    #         self.nsecx = 16

    #         img_arr = img.array
    #         ny, nx = img_arr.shape
    #         dx = int(nx/self.nsecx)
    #         dy = int(ny/self.nsecy)
    #         newimg = galsim.Image(nx, int(ny+over2), init_value=0)
    #         for ichannel in range(16):
    #             print('\n***add CTI effects: pointing-{:} chip-{:} channel-{:}***'.format(
    #                 pointing_ID, self.chipID, ichannel+1))
    #             # nx,ny,noverscan,nsp,nmax = 4608,4616,84,3,10
    #             noverscan, nsp, nmax = over2, 3, 10
    #             beta, w, c = 0.478, 84700, 0
    #             t = np.array([0.74, 7.7, 37], dtype=np.float32)
    #             rho_trap = np.array([0.6, 1.6, 1.4], dtype=np.float32)
    #             trap_seeds = np.array(
    #                 [0, 1000, 10000], dtype=np.int32) + ichannel + self.chipID*16
    #             release_seed = 50 + ichannel + pointing_ID*30 + self.chipID*16
    #             newimg.array[:, 0+ichannel*dx:dx+ichannel*dx] = CTI_sim(
    #                 img_arr[:, 0+ichannel*dx:dx+ichannel*dx], dx, dy, noverscan, nsp, nmax, beta, w, c, t, rho_trap, trap_seeds, release_seed)
    #         newimg.wcs = img.wcs
    #         del img
    #         img = newimg

    #         # 1*16 -> 2*8 img-layout
    #         img = chip_utils.formatRevert(GSImage=img)
    #         self.nsecy = 2
    #         self.nsecx = 8

    #     # prescan & overscan
    #     if config["ins_effects"]["add_prescan"] == True:
    #         chip_utils.log_info(
    #             msg="  Apply pre/over-scan", logger=self.logger)
    #         if config["ins_effects"]["cte_trail"] == False:
    #             img = chip_utils.AddPreScan(
    #                 GSImage=img, pre1=pre1, pre2=pre2, over1=over1, over2=over2)
    #         if config["ins_effects"]["cte_trail"] == True:
    #             img = chip_utils.AddPreScan(
    #                 GSImage=img, pre1=pre1, pre2=pre2, over1=over1, over2=0)

    #     # 1*16 output
    #     if config["ins_effects"]["format_output"] == True:
    #         chip_utils.log_info(msg="  Apply 1*16 format", logger=self.logger)
    #         img = chip_utils.formatOutput(GSImage=img)
    #         self.nsecy = 1
    #         self.nsecx = 16

    #     # Add Bias level
    #     if config["ins_effects"]["add_bias"] == True:
    #         chip_utils.log_info(
    #             msg="  Adding Bias level and 16-channel non-uniformity", logger=self.logger)
    #         if config["ins_effects"]["bias_16channel"] == True:
    #             img = effects.AddBiasNonUniform16(img,
    #                                               bias_level=float(
    #                                                   self.bias_level),
    #                                               nsecy=self.nsecy, nsecx=self.nsecx,
    #                                               seed=SeedBiasNonuni+self.chipID,
    #                                               logger=self.logger)
    #         elif config["ins_effects"]["bias_16channel"] == False:
    #             img += self.bias_level

    #     # Add Read-out Noise
    #     if config["ins_effects"]["add_readout"] == True:
    #         seed = int(config["random_seeds"]["seed_readout"]
    #                    ) + pointing_ID*30 + self.chipID
    #         rng_readout = galsim.BaseDeviate(seed)
    #         readout_noise = galsim.GaussianNoise(
    #             rng=rng_readout, sigma=self.read_noise)
    #         img.addNoise(readout_noise)

    #     # Apply Gain & Quantization
    #     chip_utils.log_info(
    #         msg="  Applying Gain (and 16 channel non-uniformity) & Quantization", logger=self.logger)
    #     if config["ins_effects"]["gain_16channel"] == True:
    #         img, self.gain_channel = effects.ApplyGainNonUniform16(
    #             img, gain=self.gain,
    #             nsecy=self.nsecy, nsecx=self.nsecx,
    #             seed=SeedGainNonuni+self.chipID,
    #             logger=self.logger)
    #     elif config["ins_effects"]["gain_16channel"] == False:
    #         img /= self.gain

    #     img.array[img.array > 65535] = 65535
    #     img.replaceNegative(replace_value=0)
    #     img.quantize()

    #     ######################################################################################
    #     # Output images for calibration pointing
    #     ######################################################################################
    #     # Bias output
    #     if config["ins_effects"]["add_bias"] == True and config["output_setting"]["bias_output"] == True and pointing_type == 'CAL':
    #         if self.logger is not None:
    #             self.logger.info("  Output N frame Bias files")
    #         else:
    #             print("  Output N frame Bias files", flush=True)
    #         NBias = int(config["output_setting"]["NBias"])
    #         for i in range(NBias):
    #             # BiasCombImg, BiasTag = effects.MakeBiasNcomb(
    #             # self.npix_x, self.npix_y,
    #             # bias_level=float(self.bias_level),
    #             # ncombine=1, read_noise=self.read_noise, gain=1,
    #             # seed=SeedBiasNonuni+self.chipID,
    #             # logger=self.logger)
    #             BiasCombImg = galsim.Image(
    #                 self.npix_x, self.npix_y, init_value=0)
    #             if config["ins_effects"]["add_bias"] == True:
    #                 biaslevel = self.bias_level
    #                 overscan = biaslevel-2
    #             elif config["ins_effects"]["add_bias"] == False:
    #                 biaslevel = 0
    #                 overscan = 0

    #             # Readout noise for Biases is not generated with random seeds. So readout noise for bias images can't be reproduced.
    #             if config["ins_effects"]["cosmic_ray"] == True:
    #                 if config["ins_effects"]["cray_differ"] == True:
    #                     cr_map, cr_event_num = effects.produceCR_Map(
    #                         xLen=self.npix_x, yLen=self.npix_y,
    #                         exTime=0.01,
    #                         cr_pixelRatio=0.003 *
    #                         (0.01+0.5*self.readout_time)/150.,
    #                         gain=self.gain,
    #                         attachedSizes=self.attachedSizes,
    #                         seed=SeedCosmicRay+pointing_ID*30+self.chipID+1)
    #                     # seed: obj-imaging:+0; bias:+1; dark:+2; flat:+3;
    #                 BiasCombImg += cr_map
    #                 del cr_map

    #             # Apply Bad lines
    #             if config["ins_effects"]["add_badcolumns"] == True:
    #                 BiasCombImg = effects.BadColumns(
    #                     BiasCombImg-float(self.bias_level)+5, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger) + float(self.bias_level)-5

    #             # Non-Linearity for Bias
    #             if config["ins_effects"]["non_linear"] == True:
    #                 if self.logger is not None:
    #                     self.logger.info(
    #                         "  Applying Non-Linearity on the Bias image")
    #                 else:
    #                     print(
    #                         "  Applying Non-Linearity on the Bias image", flush=True)
    #                 BiasCombImg = effects.NonLinearity(
    #                     GSImage=BiasCombImg, beta1=5.e-7, beta2=0)

    #             # START
    #             pre1 = self.prescan_x  # 27
    #             over1 = self.overscan_x  # 71
    #             pre2 = self.prescan_y  # 0 #4
    #             over2 = self.overscan_y  # 84 #80

    #             # prescan & overscan
    #             if config["ins_effects"]["add_prescan"] == True:
    #                 chip_utils.log_info(
    #                     msg="  Apply pre/over-scan", logger=self.logger)
    #                 BiasCombImg = chip_utils.AddPreScan(
    #                     GSImage=BiasCombImg, pre1=pre1, pre2=pre2, over1=over1, over2=over2)

    #             # 1*16 output
    #             if config["ins_effects"]["format_output"] == True:
    #                 chip_utils.log_info(
    #                     msg="  Apply 1*16 format", logger=self.logger)
    #                 BiasCombImg = chip_utils.formatOutput(GSImage=BiasCombImg)
    #                 self.nsecy = 1
    #                 self.nsecx = 16
    #             # END

    #             # Add Bias level
    #             if config["ins_effects"]["add_bias"] == True:
    #                 if self.logger is not None:
    #                     self.logger.info(
    #                         "  Adding Bias level and 16-channel non-uniformity")
    #                 else:
    #                     print("  Adding Bias level and 16-channel non-uniformity")
    #                 BiasCombImg = effects.AddBiasNonUniform16(BiasCombImg,
    #                                                           bias_level=biaslevel,
    #                                                           nsecy=self.nsecy, nsecx=self.nsecx,
    #                                                           seed=SeedBiasNonuni+self.chipID,
    #                                                           logger=self.logger)
    #                 rng = galsim.UniformDeviate()
    #                 ncombine = 1
    #                 NoiseBias = galsim.GaussianNoise(
    #                     rng=rng, sigma=self.read_noise*ncombine**0.5)
    #                 BiasCombImg.addNoise(NoiseBias)

    #             BiasCombImg, self.gain_channel = effects.ApplyGainNonUniform16(BiasCombImg, gain=self.gain,
    #                                                                            nsecy=self.nsecy, nsecx=self.nsecx,
    #                                                                            seed=SeedGainNonuni+self.chipID,
    #                                                                            logger=self.logger)
    #             # BiasCombImg = effects.AddOverscan(
    #             #     BiasCombImg,
    #             #     overscan=float(config["ins_effects"]["bias_level"])-2, gain=self.gain,
    #             #     widthl=27, widthr=27, widtht=8, widthb=8)
    #             BiasCombImg.replaceNegative(replace_value=0)
    #             BiasCombImg.quantize()
    #             BiasCombImg = galsim.ImageUS(BiasCombImg)
    #             timestamp_obs += 10 * 60
    #             chip_utils.outputCal(
    #                 chip=self,
    #                 img=BiasCombImg,
    #                 ra_cen=ra_cen,
    #                 dec_cen=dec_cen,
    #                 img_rot=img_rot,
    #                 im_type='BIAS',
    #                 pointing_ID=pointing_ID,
    #                 output_dir=chip_output.subdir,
    #                 exptime=0.0,
    #                 project_cycle=config["project_cycle"],
    #                 run_counter=config["run_counter"],
    #                 timestamp=timestamp_obs)
    #         del BiasCombImg

    #     # Export combined (ncombine, Vignetting + PRNU) & single vignetting flat-field file
    #     if config["ins_effects"]["flat_fielding"] == True and config["output_setting"]["flat_output"] == True and pointing_type == 'CAL':
    #         if self.logger is not None:
    #             self.logger.info("  Output N frame Flat-Field files")
    #         else:
    #             print("  Output N frame Flat-Field files", flush=True)
    #         NFlat = int(config["output_setting"]["NFlat"])
    #         if config["ins_effects"]["add_bias"] == True:
    #             biaslevel = self.bias_level
    #             overscan = biaslevel-2
    #         elif config["ins_effects"]["add_bias"] == False:
    #             biaslevel = 0
    #             overscan = 0
    #         darklevel = self.dark_noise * \
    #             (self.flat_exptime+0.5*self.readout_time)
    #         for i in range(NFlat):
    #             FlatSingle = flat_img * prnu_img + darklevel
    #             FlatCombImg, FlatTag = effects.MakeFlatNcomb(
    #                 flat_single_image=FlatSingle,
    #                 ncombine=1,
    #                 read_noise=self.read_noise,
    #                 gain=1,
    #                 overscan=overscan,
    #                 biaslevel=0,
    #                 seed_bias=SeedDefective+self.chipID,
    #                 logger=self.logger
    #             )
    #             if config["ins_effects"]["cosmic_ray"] == True:
    #                 if config["ins_effects"]["cray_differ"] == True:
    #                     cr_map, cr_event_num = effects.produceCR_Map(
    #                         xLen=self.npix_x, yLen=self.npix_y,
    #                         exTime=self.flat_exptime+0.5*self.readout_time,
    #                         cr_pixelRatio=0.003 *
    #                         (self.flat_exptime+0.5*self.readout_time)/150.,
    #                         gain=self.gain,
    #                         attachedSizes=self.attachedSizes,
    #                         seed=SeedCosmicRay+pointing_ID*30+self.chipID+3)
    #                     # seed: obj-imaging:+0; bias:+1; dark:+2; flat:+3;
    #                 FlatCombImg += cr_map
    #                 del cr_map

    #             # Add Hot Pixels or/and Dead Pixels
    #             rgbadpix = Generator(PCG64(int(SeedDefective+self.chipID)))
    #             badfraction = 5E-5*(rgbadpix.random()*0.5+0.7)
    #             FlatCombImg = effects.DefectivePixels(
    #                 FlatCombImg, IfHotPix=BoolHotPix, IfDeadPix=BoolDeadPix, fraction=badfraction, seed=SeedDefective+self.chipID, biaslevel=0)

    #             # Apply Bad lines
    #             if config["ins_effects"]["add_badcolumns"] == True:
    #                 FlatCombImg = effects.BadColumns(
    #                     FlatCombImg, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger)

    #             if config["ins_effects"]["non_linear"] == True:
    #                 if self.logger is not None:
    #                     self.logger.info(
    #                         "  Applying Non-Linearity on the Flat image")
    #                 else:
    #                     print(
    #                         "  Applying Non-Linearity on the Flat image", flush=True)
    #                 FlatCombImg = effects.NonLinearity(
    #                     GSImage=FlatCombImg, beta1=5.e-7, beta2=0)

    #             # if config["ins_effects"]["cte_trail"] == True:
    #             # FlatCombImg = effects.CTE_Effect(GSImage=FlatCombImg, threshold=3)
    #             # START
    #             pre1 = self.prescan_x  # 27
    #             over1 = self.overscan_x  # 71
    #             pre2 = self.prescan_y  # 0 #4
    #             over2 = self.overscan_y  # 84 #80
    #             if config["ins_effects"]["cte_trail"] == True:
    #                 chip_utils.log_info(
    #                     msg="  Apply CTE Effect", logger=self.logger)
    #                 # img = effects.CTE_Effect(GSImage=img, threshold=27)
    #                 # CTI_modeling
    #                 # 2*8 -> 1*16 img-layout
    #                 FlatCombImg = chip_utils.formatOutput(GSImage=FlatCombImg)
    #                 self.nsecy = 1
    #                 self.nsecx = 16

    #                 img_arr = FlatCombImg.array
    #                 ny, nx = img_arr.shape
    #                 dx = int(nx/self.nsecx)
    #                 dy = int(ny/self.nsecy)
    #                 newimg = galsim.Image(nx, int(ny+over2), init_value=0)
    #                 for ichannel in range(16):
    #                     print('\n***add CTI effects: pointing-{:} chip-{:} channel-{:}***'.format(
    #                         pointing_ID, self.chipID, ichannel+1))
    #                     # nx,ny,noverscan,nsp,nmax = 4608,4616,84,3,10
    #                     noverscan, nsp, nmax = over2, 3, 10
    #                     beta, w, c = 0.478, 84700, 0
    #                     t = np.array([0.74, 7.7, 37], dtype=np.float32)
    #                     rho_trap = np.array([0.6, 1.6, 1.4], dtype=np.float32)
    #                     trap_seeds = np.array(
    #                         [0, 1000, 10000], dtype=np.int32) + ichannel + self.chipID*16
    #                     release_seed = 50 + ichannel + pointing_ID*30 + self.chipID*16
    #                     newimg.array[:, 0+ichannel*dx:dx+ichannel*dx] = CTI_sim(
    #                         img_arr[:, 0+ichannel*dx:dx+ichannel*dx], dx, dy, noverscan, nsp, nmax, beta, w, c, t, rho_trap, trap_seeds, release_seed)
    #                 newimg.wcs = FlatCombImg.wcs
    #                 del FlatCombImg
    #                 FlatCombImg = newimg

    #                 # 1*16 -> 2*8 img-layout
    #                 FlatCombImg = chip_utils.formatRevert(GSImage=FlatCombImg)
    #                 self.nsecy = 2
    #                 self.nsecx = 8

    #             # prescan & overscan
    #             if config["ins_effects"]["add_prescan"] == True:
    #                 chip_utils.log_info(
    #                     msg="  Apply pre/over-scan", logger=self.logger)
    #                 if config["ins_effects"]["cte_trail"] == False:
    #                     FlatCombImg = chip_utils.AddPreScan(
    #                         GSImage=FlatCombImg, pre1=pre1, pre2=pre2, over1=over1, over2=over2)
    #                 if config["ins_effects"]["cte_trail"] == True:
    #                     FlatCombImg = chip_utils.AddPreScan(
    #                         GSImage=FlatCombImg, pre1=pre1, pre2=pre2, over1=over1, over2=0)

    #             # 1*16 output
    #             if config["ins_effects"]["format_output"] == True:
    #                 chip_utils.log_info(
    #                     msg="  Apply 1*16 format", logger=self.logger)
    #                 FlatCombImg = chip_utils.formatOutput(GSImage=FlatCombImg)
    #                 self.nsecy = 1
    #                 self.nsecx = 16
    #             # END

    #             # Add Bias level
    #             if config["ins_effects"]["add_bias"] == True:
    #                 if self.logger is not None:
    #                     self.logger.info(
    #                         "  Adding Bias level and 16-channel non-uniformity")
    #                 else:
    #                     print("  Adding Bias level and 16-channel non-uniformity")
    #                 # img += float(config["ins_effects"]["bias_level"])
    #                 FlatCombImg = effects.AddBiasNonUniform16(FlatCombImg,
    #                                                           bias_level=biaslevel,
    #                                                           nsecy=self.nsecy, nsecx=self.nsecx,
    #                                                           seed=SeedBiasNonuni+self.chipID,
    #                                                           logger=self.logger)

    #             # Add Read-out Noise
    #             if config["ins_effects"]["add_readout"] == True:
    #                 seed = int(config["random_seeds"]["seed_readout"]
    #                            ) + pointing_ID*30 + self.chipID + 3
    #                 rng_readout = galsim.BaseDeviate(seed)
    #                 readout_noise = galsim.GaussianNoise(
    #                     rng=rng_readout, sigma=self.read_noise)
    #                 FlatCombImg.addNoise(readout_noise)

    #             FlatCombImg, self.gain_channel = effects.ApplyGainNonUniform16(FlatCombImg, gain=self.gain,
    #                                                                            nsecy=self.nsecy, nsecx=self.nsecx,
    #                                                                            seed=SeedGainNonuni+self.chipID,
    #                                                                            logger=self.logger)
    #             # FlatCombImg = effects.AddOverscan(FlatCombImg, overscan=overscan, gain=self.gain, widthl=27, widthr=27, widtht=8, widthb=8)
    #             FlatCombImg.replaceNegative(replace_value=0)
    #             FlatCombImg.quantize()
    #             FlatCombImg = galsim.ImageUS(FlatCombImg)
    #             timestamp_obs += 10 * 60
    #             chip_utils.outputCal(
    #                 chip=self,
    #                 img=FlatCombImg,
    #                 ra_cen=ra_cen,
    #                 dec_cen=dec_cen,
    #                 img_rot=img_rot,
    #                 im_type='FLAT',
    #                 pointing_ID=pointing_ID,
    #                 output_dir=chip_output.subdir,
    #                 exptime=self.flat_exptime,
    #                 project_cycle=config["project_cycle"],
    #                 run_counter=config["run_counter"],
    #                 timestamp=timestamp_obs)

    #         del FlatCombImg, FlatSingle, prnu_img
    #         # flat_img.replaceNegative(replace_value=0)
    #         # flat_img.quantize()
    #         # galsim.ImageUS(flat_img).write("%s/FlatImg_Vignette_%s.fits" % (chip_output.subdir, self.chipID))
    #         del flat_img

    #     # Export Dark current images
    #     if config["ins_effects"]["add_dark"] == True and config["output_setting"]["dark_output"] == True and pointing_type == 'CAL':
    #         if self.logger is not None:
    #             self.logger.info("  Output N frame Dark Current files")
    #         else:
    #             print("  Output N frame Dark Current files", flush=True)
    #         NDark = int(config["output_setting"]["NDark"])
    #         if config["ins_effects"]["add_bias"] == True:
    #             biaslevel = self.bias_level
    #             overscan = biaslevel-2
    #         elif config["ins_effects"]["add_bias"] == False:
    #             biaslevel = 0
    #             overscan = 0
    #         for i in range(NDark):
    #             DarkCombImg, DarkTag = effects.MakeDarkNcomb(
    #                 self.npix_x, self.npix_y,
    #                 overscan=overscan, bias_level=0, darkpsec=0.02, exptime=self.dark_exptime+0.5*self.readout_time,
    #                 ncombine=1, read_noise=self.read_noise,
    #                 gain=1, seed_bias=SeedBiasNonuni+self.chipID,
    #                 logger=self.logger)
    #             if config["ins_effects"]["cosmic_ray"] == True:
    #                 if config["ins_effects"]["cray_differ"] == True:
    #                     cr_map, cr_event_num = effects.produceCR_Map(
    #                         xLen=self.npix_x, yLen=self.npix_y,
    #                         exTime=self.dark_exptime+0.5*self.readout_time,
    #                         cr_pixelRatio=0.003 *
    #                         (self.dark_exptime+0.5*self.readout_time)/150.,
    #                         gain=self.gain,
    #                         attachedSizes=self.attachedSizes,
    #                         seed=SeedCosmicRay+pointing_ID*30+self.chipID+2)
    #                     # seed: obj-imaging:+0; bias:+1; dark:+2; flat:+3;
    #                 DarkCombImg += cr_map
    #                 cr_map[cr_map > 65535] = 65535
    #                 cr_map[cr_map < 0] = 0
    #                 crmap_gsimg = galsim.Image(cr_map, dtype=np.uint16)
    #                 del cr_map
    #                 # START
    #                 # prescan & overscan
    #                 if config["ins_effects"]["add_prescan"] == True:
    #                     chip_utils.log_info(
    #                         msg="  Apply pre/over-scan", logger=self.logger)
    #                     crmap_gsimg = chip_utils.AddPreScan(
    #                         GSImage=crmap_gsimg, pre1=pre1, pre2=pre2, over1=over1, over2=over2)

    #                 # 1*16 output
    #                 if config["ins_effects"]["format_output"] == True:
    #                     chip_utils.log_info(
    #                         msg="  Apply 1*16 format", logger=self.logger)
    #                     crmap_gsimg = chip_utils.formatOutput(
    #                         GSImage=crmap_gsimg)
    #                     self.nsecy = 1
    #                     self.nsecx = 16
    #                 # END
    #                 chip_utils.outputCal(
    #                     chip=self,
    #                     img=crmap_gsimg,
    #                     ra_cen=ra_cen,
    #                     dec_cen=dec_cen,
    #                     img_rot=img_rot,
    #                     im_type='CRD',
    #                     pointing_ID=pointing_ID,
    #                     output_dir=chip_output.subdir,
    #                     exptime=self.dark_exptime,
    #                     project_cycle=config["project_cycle"],
    #                     run_counter=config["run_counter"],
    #                     timestamp=timestamp_obs)
    #                 del crmap_gsimg

    #             # Add Hot Pixels or/and Dead Pixels
    #             rgbadpix = Generator(PCG64(int(SeedDefective+self.chipID)))
    #             badfraction = 5E-5*(rgbadpix.random()*0.5+0.7)
    #             DarkCombImg = effects.DefectivePixels(
    #                 DarkCombImg, IfHotPix=BoolHotPix, IfDeadPix=BoolDeadPix, fraction=badfraction, seed=SeedDefective+self.chipID, biaslevel=0)

    #             # Apply Bad lines
    #             if config["ins_effects"]["add_badcolumns"] == True:
    #                 DarkCombImg = effects.BadColumns(
    #                     DarkCombImg, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger)

    #             # Non-Linearity for Dark
    #             if config["ins_effects"]["non_linear"] == True:
    #                 if self.logger is not None:
    #                     self.logger.info(
    #                         "  Applying Non-Linearity on the Dark image")
    #                 else:
    #                     print(
    #                         "  Applying Non-Linearity on the Dark image", flush=True)
    #                 DarkCombImg = effects.NonLinearity(
    #                     GSImage=DarkCombImg, beta1=5.e-7, beta2=0)

    #             # if config["ins_effects"]["cte_trail"] == True:
    #             # DarkCombImg = effects.CTE_Effect(GSImage=DarkCombImg, threshold=3)
    #             # START
    #             pre1 = self.prescan_x  # 27
    #             over1 = self.overscan_x  # 71
    #             pre2 = self.prescan_y  # 0 #4
    #             over2 = self.overscan_y  # 84 #80
    #             if config["ins_effects"]["cte_trail"] == True:
    #                 chip_utils.log_info(
    #                     msg="  Apply CTE Effect", logger=self.logger)
    #                 # img = effects.CTE_Effect(GSImage=img, threshold=27)
    #                 # CTI_modeling
    #                 # 2*8 -> 1*16 img-layout
    #                 DarkCombImg = chip_utils.formatOutput(GSImage=DarkCombImg)
    #                 self.nsecy = 1
    #                 self.nsecx = 16

    #                 img_arr = DarkCombImg.array
    #                 ny, nx = img_arr.shape
    #                 dx = int(nx/self.nsecx)
    #                 dy = int(ny/self.nsecy)
    #                 newimg = galsim.Image(nx, int(ny+over2), init_value=0)
    #                 for ichannel in range(16):
    #                     print('\n***add CTI effects: pointing-{:} chip-{:} channel-{:}***'.format(
    #                         pointing_ID, self.chipID, ichannel+1))
    #                     # nx,ny,noverscan,nsp,nmax = 4608,4616,84,3,10
    #                     noverscan, nsp, nmax = over2, 3, 10
    #                     beta, w, c = 0.478, 84700, 0
    #                     t = np.array([0.74, 7.7, 37], dtype=np.float32)
    #                     rho_trap = np.array([0.6, 1.6, 1.4], dtype=np.float32)
    #                     trap_seeds = np.array(
    #                         [0, 1000, 10000], dtype=np.int32) + ichannel + self.chipID*16
    #                     release_seed = 50 + ichannel + pointing_ID*30 + self.chipID*16
    #                     newimg.array[:, 0+ichannel*dx:dx+ichannel*dx] = CTI_sim(
    #                         img_arr[:, 0+ichannel*dx:dx+ichannel*dx], dx, dy, noverscan, nsp, nmax, beta, w, c, t, rho_trap, trap_seeds, release_seed)
    #                 newimg.wcs = DarkCombImg.wcs
    #                 del DarkCombImg
    #                 DarkCombImg = newimg

    #                 # 1*16 -> 2*8 img-layout
    #                 DarkCombImg = chip_utils.formatRevert(GSImage=DarkCombImg)
    #                 self.nsecy = 2
    #                 self.nsecx = 8

    #             # prescan & overscan
    #             if config["ins_effects"]["add_prescan"] == True:
    #                 chip_utils.log_info(
    #                     msg="  Apply pre/over-scan", logger=self.logger)
    #                 if config["ins_effects"]["cte_trail"] == False:
    #                     DarkCombImg = chip_utils.AddPreScan(
    #                         GSImage=DarkCombImg, pre1=pre1, pre2=pre2, over1=over1, over2=over2)
    #                 if config["ins_effects"]["cte_trail"] == True:
    #                     DarkCombImg = chip_utils.AddPreScan(
    #                         GSImage=DarkCombImg, pre1=pre1, pre2=pre2, over1=over1, over2=0)

    #             # 1*16 output
    #             if config["ins_effects"]["format_output"] == True:
    #                 chip_utils.log_info(
    #                     msg="  Apply 1*16 format", logger=self.logger)
    #                 DarkCombImg = chip_utils.formatOutput(GSImage=DarkCombImg)
    #                 self.nsecy = 1
    #                 self.nsecx = 16
    #             # END

    #             # Add Bias level
    #             if config["ins_effects"]["add_bias"] == True:
    #                 if self.logger is not None:
    #                     self.logger.info(
    #                         "  Adding Bias level and 16-channel non-uniformity")
    #                 else:
    #                     print("  Adding Bias level and 16-channel non-uniformity")
    #                 # img += float(config["ins_effects"]["bias_level"])
    #                 DarkCombImg = effects.AddBiasNonUniform16(DarkCombImg,
    #                                                           bias_level=biaslevel,
    #                                                           nsecy=self.nsecy, nsecx=self.nsecx,
    #                                                           seed=SeedBiasNonuni+self.chipID,
    #                                                           logger=self.logger)