Chip.py

import galsim
import os
import numpy as np
import Instrument.Chip.Effects as effects
from Instrument.FocalPlane import FocalPlane
from astropy.table import Table
from numpy.random import Generator, PCG64

class Chip(FocalPlane):
    def __init__(self, chipID, ccdEffCurve_dir, CRdata_dir, normalize_dir=None, sls_dir=None, config=None, treering_func=None):
        # Get focal plane (instance of paraent class) info
        # TODO: use chipID to config individual chip?
        super().__init__()
        # if config is not None:
        #     self.npix_x     = config["npix_x"]
        #     self.npix_y     = config["npix_y"]
        #     self.read_noise = config["read_noise"]
        #     self.dark_noise = config["dark_noise"]
        #     self.pix_scale  = config["pix_scale"]
        #     self.gain       = config["gain"]
        #     self.bias_level = config["bias_level"]
        #     self.overscan   = config["overscan"]
        # else:
        # Default setting
        self.npix_x = 9216
        self.npix_y = 9232
        self.read_noise = 5.0   # e/pix
        self.dark_noise = 0.02  # e/pix/s
        self.pix_scale  = 0.074 # pixel scale
        # self.gain       = 1.0
        # self.bias_level = 1000  # e-/pix
        self.gain = float(config["gain"])
        self.bias_level = float(config["bias_level"])
        self.overscan   = 1000
        self.exptime    = 150   # second

        # A chip ID must be assigned
        self.chipID = int(chipID)
        self._getChipRowCol()

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

        # Get boundary (in pix)
        self.bound = self.getChipLim()
        self.ccdEffCurve_dir = ccdEffCurve_dir
        self.CRdata_dir = CRdata_dir
        self.normalize_dir = normalize_dir
        self.sls_dir=sls_dir
        # self.sls_conf = os.path.join(self.sls_dir, self.getChipSLSConf())
        slsconfs = self.getChipSLSConf()
        if np.size(slsconfs) == 1:
            self.sls_conf = [os.path.join(self.sls_dir, slsconfs)]
        else:
            self.sls_conf = [os.path.join(self.sls_dir, slsconfs[0]), os.path.join(self.sls_dir, slsconfs[1])]
        
        if self.normalize_dir is not None:
            self._getNormF()
        self.effCurve = self._getChipEffCurve(self.filter_type)
        self._getCRdata()

        # Define the sensor
        if config["bright_fatter"].lower() == "y":
            self.sensor = galsim.SiliconSensor(strength=config["df_strength"], treering_func=treering_func)
        else:
            self.sensor = galsim.Sensor()

    # def _getChipRowCol(self):
    #     self.rowID = (self.chipID - 1) // 5 + 1
    #     self.colID = (self.chipID - 1) % 5 + 1
    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 ["GI", "GV", "GU"]:
            return "spectroscopic"
        else:
            return "photometric"

    def _getNormF(self):
        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'))

    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']: filename = 'Astro_MB.txt'
        if filter_type in ['i', 'z', 'y', 'GI']: filename = 'Basic_NIR.txt'
        # Mirror efficiency:
        if filter_type == 'nuv': mirror_eff = 0.54
        if filter_type == 'u': mirror_eff = 0.68
        if filter_type in ['g', 'r', 'i', 'z', 'y']: mirror_eff = 0.8
        if filter_type in ['GU', 'GV', 'GI']: mirror_eff = 1. # Not sure if this is right
        
        path = os.path.join(self.ccdEffCurve_dir, filename)
        table = Table.read(path, format='ascii')
        throughput = galsim.LookupTable(x=table['col1'], f=table['col2']*mirror_eff, interpolant='linear')
        bandpass = galsim.Bandpass(throughput, wave_type='nm')
        return bandpass

    def _getCRdata(self):
        path = os.path.join(self.CRdata_dir, 'wfc-cr-attachpixel.dat')
        self.attachedSizes = np.loadtxt(path)

    def getChipFilter(self, chipID=None, filter_layout=None):
        """Return the filter index and type for a given chip #(chipID)
        """
        filter_type_list = ["nuv","u", "g", "r", "i","z","y","GU", "GV", "GI"]
        # TODO: maybe a more elegent way other than hard coded?
        # e.g. use something like a nested dict:
        if filter_layout is not None:
            return filter_layout[chipID][0], filter_layout[chipID][1]
        if chipID == None:
            chipID = self.chipID

        # updated configurations
        # if chipID>30 or chipID<1: raise ValueError("!!! Chip ID: [1,30]")
        # if chipID in [10, 15, 16, 21]: filter_type = 'y'
        # if chipID in [11, 20]:         filter_type = "z"
        # if chipID in [9, 22]:           filter_type = "i"
        # if chipID in [12, 19]:         filter_type = "u"
        # if chipID in [7, 24]:         filter_type = "r"
        # if chipID in [14, 13, 18, 17]:    filter_type = "nuv"
        # if chipID in [8, 23]:         filter_type = "g"
        # if chipID in [6, 5, 25, 26]:    filter_type = "GI"
        # if chipID in [27, 30, 1, 4]:    filter_type = "GV"
        # if chipID in [28, 29, 2, 3]:    filter_type = "GU"
        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"
        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 observation.
        Parameters:
            chipID:         int
                            the index of the chip
        Returns:
            A galsim BoundsD object
        """
        # if chipID == None:
        #     chipID = self.chipID
        
        # gx = self.npix_gap_x
        # gy1, gy2 = self.npix_gap_y

        # # xlim of a given ccd chip
        # xrem = (chipID-1)%self.nchip_x - self.nchip_x // 2
        # xcen = (self.npix_x + gx) * xrem
        # nx0 = xcen - self.npix_x//2 + 1
        # nx1 = xcen + self.npix_x//2

        # # ylim of a given ccd chip
        # yrem = 2*((chipID-1)//self.nchip_x) - (self.nchip_y-1)
        # ycen = (self.npix_y//2 + gy1//2) * yrem
        # if chipID <= 6: ycen = (self.npix_y//2 + gy1//2) * yrem - (gy2-gy1)
        # if chipID >= 25: ycen = (self.npix_y//2 + gy1//2) * yrem + (gy2-gy1)
        # ny0 = ycen - self.npix_y//2 + 1
        # ny1 = ycen + self.npix_y//2

        if chipID == None:
            chipID = self.chipID
            rowID, colID = self.rowID, self.colID
        else:
            rowID, colID = self.getChipRowCol(chipID)
        gx1, gx2 = self.npix_gap_x
        gy = self.npix_gap_y

        # xlim of a given CCD chip
        xrem = 2*(colID - 1) - (self.nchip_x - 1)
        xcen = (self.npix_x//2 + gx1//2) * xrem
        if chipID >= 26 or chipID == 21:
            xcen = (self.npix_x//2 + gx1//2) * xrem - (gx2-gx1)
        if chipID <= 5 or chipID == 10:
            xcen = (self.npix_x//2 + gx1//2) * xrem + (gx2-gx1)
        nx0 = xcen - self.npix_x//2 + 1
        nx1 = xcen + self.npix_x//2

        # ylim of a given CCD chip
        yrem = (rowID - 1) - self.nchip_y // 2
        ycen = (self.npix_y + gy) * yrem
        ny0 = ycen - self.npix_y//2 + 1
        ny1 = ycen + self.npix_y//2

        return galsim.BoundsD(nx0-1, nx1-1, ny0-1, ny1-1)


    def getSkyCoverage(self, wcs):
        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, dec_obj, wcs=None, margin=1):
        # magin in number of pix
        if wcs is None:
            wcs = self.img.wcs
        pos_obj = wcs.toImage(galsim.CelestialCoord(ra=ra_obj*galsim.degrees, dec=dec_obj*galsim.degrees))
        xmin, xmax = self.bound.xmin - margin, self.bound.xmax + margin
        ymin, ymax = self.bound.ymin - margin, self.bound.ymax + margin
        if (pos_obj.x - xmin) * (pos_obj.x - xmax) > 0.0 or (pos_obj.y - ymin) * (pos_obj.y - 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 addNoise_phot(self, img, exptime=150.0, sky_noise=0., seed=31415):
        rng = galsim.BaseDeviate(seed)
        dark_noise = galsim.DeviateNoise(galsim.PoissonDeviate(rng, self.dark_noise*exptime))
        img.addNoise(dark_noise)
        ccd_noise = galsim.CCDNoise(rng, sky_level=sky_noise, gain=self.gain, read_noise=self.read_noise)
        img.addNoise(ccd_noise)
        return img

    def addNoise_spec(self, config, tel, img, sky_map, exptime=150.0, seed=31415):
        if img.array.shape != sky_map.shape:
            raise ValueError("The shape img and sky_map must be equal.")

        # n_img_arrar = (img.array + sky_map) * tel.pupil_area * exptime
        # Should be the following?
        rng = galsim.BaseDeviate(seed)
        noise_img = galsim.Image(img, copy=True)
        dark_noise = galsim.DeviateNoise(galsim.PoissonDeviate(rng, self.dark_noise*exptime))
        img.addNoise(dark_noise)
        if config["abs_back"].lower() == "y":
            img += (sky_map * tel.pupil_area * exptime)
        ccd_noise = galsim.CCDNoise(rng, gain=self.gain, read_noise=self.read_noise)
        img.addNoise(ccd_noise)
        return img
        # pNoise = self.dark_noise * exptime
        # noise_img = galsim.Image(n_img_arrar, copy=True)
        # dNose = galsim.PoissonNoise(rng=rng, sky_level=pNoise)
        # noise_img.addNoise(dNose)
        # rdNoise = galsim.GaussianNoise(rng=rng, sigma=self.read_noise)
        # noise_img.addNoise(rdNoise)
        # return noise_img

    def addNoise(self, config, tel, filt, img, sky_map=None, exptime=150.0):
        if self.survey_type == "photometric":
            sky_level = filt.getSkyNoise(exptime=exptime, gain=self.gain)
            img = self.addNoise_phot(
                img=img, 
                exptime=exptime, 
                sky_noise=filt.getSkyNoise(exptime=exptime, gain=self.gain), 
                seed=int(config["seed_skynoise"]))
            if config["abs_back"].lower() == "y":
                img += sky_level
        elif self.survey_type == "spectroscopic":
            if sky_map is None:
                raise ValueError("Must provide sky_map for spectroscopic chip")
            img = self.addNoise_spec(
                config=config, 
                tel=tel, 
                img=img, 
                sky_map=sky_map, 
                exptime=exptime, 
                seed=int(config["seed_skynoise"]))
        return img

    def getChipSLSConf(self):
        confFile = ''
        if self.chipID == 1: confFile = ['CSST_GI2.conf', 'CSST_GI1.conf']
        if self.chipID == 2: confFile = ['CSST_GV4.conf', 'CSST_GV3.conf']
        if self.chipID == 3: confFile = ['CSST_GU2.conf', 'CSST_GU1.conf']
        if self.chipID == 4: confFile = ['CSST_GU4.conf', 'CSST_GU3.conf']
        if self.chipID == 5: confFile = ['CSST_GV2.conf', 'CSST_GV1.conf']
        if self.chipID == 10: confFile = ['CSST_GI4.conf', 'CSST_GI3.conf']
        if self.chipID == 21: confFile = ['CSST_GI6.conf', 'CSST_GI5.conf']
        if self.chipID == 26: confFile = ['CSST_GV8.conf', 'CSST_GV7.conf']
        if self.chipID == 27: confFile = ['CSST_GU6.conf', 'CSST_GU5.conf']
        if self.chipID == 28: confFile = ['CSST_GU8.conf', 'CSST_GU7.conf']
        if self.chipID == 29: confFile = ['CSST_GV6.conf', 'CSST_GV5.conf']
        if self.chipID == 30: confFile = ['CSST_GI8.conf', 'CSST_GI7.conf']
        return confFile

    def addEffects(self, config, img, chip_output, filt, exptime=150., pointing_ID=0):
        SeedGainNonuni=int(config["seed_gainNonUniform"])
        SeedBiasNonuni=int(config["seed_biasNonUniform"])
        SeedRnNonuni = int(config["seed_rnNonUniform"])
        SeedBadColumns = int(config["seed_badcolumns"])
        SeedDefective = int(config["seed_defective"])
        SeedCosmicRay = int(config["seed_CR"])
        fullwell = int(config["full_well"])
        if config["add_hotpixels"].lower() == "y":
            BoolHotPix = True
        else:
            BoolHotPix = False
        if config["add_deadpixels"].lower() == "y":
            BoolDeadPix = True
        else:
            BoolDeadPix = False

        # Apply flat-field large scale structure for one chip
        if config["flat_fielding"].lower() == "y":
            print("  Creating and applying Flat-Fielding", flush=True)
            print(img.bounds, flush=True)
            flat_img = effects.MakeFlatSmooth(
                img.bounds, 
                int(config["seed_flat"]))
            flat_normal = flat_img / np.mean(flat_img.array)
            img *= flat_normal
            del flat_normal
            if config["flat_output"].lower() == "n":
                del flat_img

        # Apply Shutter-effect for one chip
        if config["shutter_effect"].lower() == "y":
            print("  Apply shutter effect", flush=True)
            shuttimg = effects.ShutterEffectArr(img, t_shutter=1.3, dist_bearing=735, dt=1E-3)    # shutter effect normalized image for this chip
            img *= shuttimg
            if config["shutter_output"].lower() == "y":    # output 16-bit shutter effect image with pixel value <=65535
                shutt_gsimg = galsim.ImageUS(shuttimg*6E4)
                shutt_gsimg.write("%s/ShutterEffect_%s_1.fits" % (chip_output.subdir, self.chipID))
                del shutt_gsimg
            del shuttimg

        # Apply PRNU effect and output PRNU flat file:
        if config["prnu_effect"].lower() == "y":
            print("  Applying PRNU effect", flush=True)
            prnu_img = effects.PRNU_Img(
                xsize=self.npix_x, 
                ysize=self.npix_y, 
                sigma=0.01, 
                seed=int(config["seed_prnu"]+self.chipID))
            img *= prnu_img
            if config["prnu_output"].lower() == "y":
                prnu_img.write("%s/FlatImg_PRNU_%s.fits" % (chip_output.subdir,self.chipID))
            if config["flat_output"].lower() == "n":
                del prnu_img

        # Add dark current
        img += self.dark_noise*exptime

        # Add cosmic-rays
        if config["cosmic_ray"].lower() == "y":
            print("  Adding Cosmic-Ray", flush=True)
            cr_map = effects.produceCR_Map(
                xLen=self.npix_x, yLen=self.npix_y, 
                exTime=exptime, 
                cr_pixelRatio=0.003, 
                gain=self.gain, 
                attachedSizes=self.attachedSizes,
                seed=SeedCosmicRay+pointing_ID*30+self.chipID)
            img += cr_map
            cr_map[cr_map > 65535] = 65535
            cr_map[cr_map < 0] = 0
            crmap_gsimg = galsim.Image(cr_map, dtype=np.uint16)
            # crmap_gsimg.write("%s/CosmicRay_%s_1.fits" % (chip_output.subdir, self.chipID))
            crmap_gsimg.write("%s/CosmicRay_%s.fits" % (chip_output.subdir, self.chipID))
            del crmap_gsimg

        # Add Bias level
        if config["add_bias"].lower() == "y":
            print("  Adding Bias level and 16-channel non-uniformity")
            # img += float(config["bias_level"])
            img = effects.AddBiasNonUniform16(img, 
                bias_level=float(config["bias_level"]), 
                nsecy = 2, nsecx=8, 
                seed=SeedBiasNonuni+self.chipID)

        # Bias output
        if config["bias_output"].lower() == "y":
            print("  Output N frame Bias files", flush=True)
            NBias = int(config["NBias"])
            for i in range(NBias):
                BiasCombImg, BiasTag = effects.MakeBiasNcomb(
                    self.npix_x, self.npix_y, 
                    bias_level=float(config["bias_level"]), 
                    ncombine=1, read_noise=self.read_noise, gain=1,
                    seed=SeedBiasNonuni+self.chipID)
                
                # Non-Linearity for Bias
                if config["non_linear"].lower() == "y":
                    print("  Applying Non-Linearity on the Bias image", flush=True)
                    BiasCombImg = effects.NonLinearity(GSImage=BiasCombImg, beta1=1.e-7, beta2=1.e-10)

                BiasCombImg = effects.ApplyGainNonUniform16(BiasCombImg, gain=self.gain, 
                    nsecy = 2, nsecx=8, 
                    seed=SeedGainNonuni+self.chipID)
                # BiasCombImg = effects.AddOverscan(
                #     BiasCombImg, 
                #     overscan=float(config["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)
                BiasCombImg.write("%s/BiasImg_%s_%s_%s.fits" % (chip_output.subdir, BiasTag, self.chipID, i+1))
            del BiasCombImg

        # Export combined (ncombine, Vignetting + PRNU) & single vignetting flat-field file
        if config["flat_output"].lower() == "y":
            print("  Output N frame Flat-Field files", flush=True)
            NFlat = int(config["NFlat"])
            if config["add_bias"].lower() == "y":
                biaslevel = self.bias_level
                overscan = biaslevel-2
            elif config["add_bias"].lower() == "n":
                biaslevel = 0
                overscan = 0
            darklevel = self.dark_noise*self.exptime
            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=biaslevel,
                    seed_bias=SeedDefective+self.chipID
                    )
                if config["cosmic_ray"].lower() == "y":
                    FlatCombImg += cr_map

                if config["non_linear"].lower() == "y":
                    print("  Applying Non-Linearity on the Flat image", flush=True)
                    FlatCombImg = effects.NonLinearity(GSImage=FlatCombImg, beta1=1.e-7, beta2=1.e-10)

                if config["cte_trail"].lower() == "y":
                    FlatCombImg = effects.CTE_Effect(GSImage=FlatCombImg, threshold=3)

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

                # 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=self.bias_level)

                FlatCombImg = effects.ApplyGainNonUniform16(FlatCombImg, gain=self.gain, 
                    nsecy = 2, nsecx=8, 
                    seed=SeedGainNonuni+self.chipID)
                # 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)
                FlatCombImg.write("%s/FlatImg_%s_%s_%s.fits" % (chip_output.subdir, FlatTag, self.chipID, i+1))
            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["dark_output"].lower() == "y":
            print("  Output N frame Dark Current files", flush=True)
            NDark = int(config["NDark"])
            if config["add_bias"].lower() == "y":
                biaslevel = self.bias_level
                overscan = biaslevel-2
            elif config["add_bias"].lower() == "n":
                biaslevel = 0
                overscan = 0
            for i in range(NDark):
                DarkCombImg, DarkTag = effects.MakeDarkNcomb(
                    self.npix_x, self.npix_y, 
                    overscan=overscan, bias_level=biaslevel, darkpsec=0.02, exptime=150,
                    ncombine=1, read_noise=self.read_noise, 
                    gain=1, seed_bias=SeedBiasNonuni+self.chipID)
                if config["cosmic_ray"].lower() == "y":
                    DarkCombImg += cr_map

                # Non-Linearity for Dark
                if config["non_linear"].lower() == "y":
                    print("  Applying Non-Linearity on the Dark image", flush=True)
                    DarkCombImg = effects.NonLinearity(GSImage=DarkCombImg, beta1=1.e-7, beta2=1.e-10)

                if config["cte_trail"].lower() == "y":
                    DarkCombImg = effects.CTE_Effect(GSImage=DarkCombImg, threshold=3)

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

                # 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=self.bias_level)

                DarkCombImg = effects.ApplyGainNonUniform16(
                    DarkCombImg, gain=self.gain, 
                    nsecy = 2, nsecx=8, 
                    seed=SeedGainNonuni+self.chipID)
                # DarkCombImg = effects.AddOverscan(
                #     DarkCombImg, 
                #     overscan=overscan, gain=self.gain, 
                #     widthl=27, widthr=27, widtht=8, widthb=8)
                DarkCombImg.replaceNegative(replace_value=0)
                DarkCombImg.quantize()
                DarkCombImg = galsim.ImageUS(DarkCombImg)
                DarkCombImg.write("%s/DarkImg_%s_%s_%s.fits" % (chip_output.subdir, DarkTag, self.chipID, i+1))
            del DarkCombImg

        # garbage collection of cosmic-ray array
        if config["cosmic_ray"].lower() == "y":
            del cr_map

        # Apply Nonlinearity on the chip image
        if config["non_linear"].lower() == "y":
            print("  Applying Non-Linearity on the chip image", flush=True)
            img = effects.NonLinearity(GSImage=img, beta1=1.e-7, beta2=1.e-10)

        # Apply CCD Saturation & Blooming
        if config["saturbloom"].lower() == "y":
            print("  Applying CCD Saturation & Blooming")
            img = effects.SaturBloom(GSImage=img, nsect_x=1, nsect_y=1, fullwell=fullwell)

        # Apply CTE Effect
        if config["cte_trail"].lower() == "y":
            print("  Apply CTE Effect")
            img = effects.CTE_Effect(GSImage=img, threshold=27)

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

        # 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=self.bias_level)
        

        # Apply Gain & Quantization
        print("  Applying Gain (and 16 channel non-uniformity) & Quantization", flush=True)
        img = effects.ApplyGainNonUniform16(
            img, gain=self.gain, 
            nsecy = 2, nsecx=8, 
            seed=SeedGainNonuni+self.chipID)
        img.array[img.array > 65535] = 65535
        img.replaceNegative(replace_value=0)
        img.quantize()
        # img = galsim.ImageUS(img)

        # # 16 output channel, with each a single image file
        # if config["readout16"].lower() == "y":
        #     print("  16 Output Channel simulation")
        #     for coli in [0, 1]:
        #         for rowi in range(8):
        #             sub_img = effects.readout16(
        #                 GSImage=img, 
        #                 rowi=rowi, 
        #                 coli=coli, 
        #                 overscan_value=self.overscan)
        #             rowcoltag = str(rowi) + str(coli)
        #             img_name_root = chip_output.img_name.split(".")[0]
        #             sub_img.write("%s/%s_%s.fits" % (chip_output.subdir, img_name_root, rowcoltag))
        #     del sub_img
        return img