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 from Config.Header import generatePrimaryHeader, generateExtensionHeader from astropy.io import fits from datetime import datetime 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 # self.dark_exptime = 300 # second # self.flat_exptime = 150 # second self.dark_exptime = float(config['dark_exptime']) self.flat_exptime = float(config['flat_exptime']) # 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 generateHeader(self, ra_cen, dec_cen, img_rot, im_type, pointing_ID, date_obs, time_obs, exptime=150.): h_prim = generatePrimaryHeader( xlen=self.npix_x, ylen=self.npix_y, pointNum = str(pointing_ID), ra=ra_cen, dec=dec_cen, psize=self.pix_scale, row_num=self.rowID, col_num=self.colID, date=date_obs, time_obs=time_obs, im_type = im_type, exptime=exptime ) h_ext = generateExtensionHeader( xlen=self.npix_x, ylen=self.npix_y, ra=ra_cen, dec=dec_cen, pa=img_rot.deg, gain=self.gain, readout=self.read_noise, dark=self.dark_noise, saturation=90000, psize=self.pix_scale, row_num=self.rowID, col_num=self.colID, extName='raw') return h_prim, h_ext def outputCal(self, img, ra_cen, dec_cen, img_rot, im_type, pointing_ID, date_obs, time_obs, output_dir, exptime=150.): h_prim, h_ext = self.generateHeader( ra_cen=ra_cen, dec_cen=dec_cen, img_rot=img_rot, im_type=im_type, pointing_ID=pointing_ID, date_obs=date_obs, time_obs=time_obs, exptime=exptime) hdu1 = fits.PrimaryHDU(header=h_prim) hdu2 = fits.ImageHDU(img.array, header=h_ext) hdu1 = fits.HDUList([hdu1, hdu2]) fname = os.path.join(output_dir, h_prim['FILENAME']+'.fits') hdu1.writeto(fname, output_verify='ignore', overwrite=True) def addEffects(self, config, img, chip_output, filt, ra_cen, dec_cen, img_rot, exptime=150., pointing_ID=0, timestamp_obs=1621915200, pointing_type='MS'): 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" and pointing_type=='MS': 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) # seed: obj-imaging:+0; bias:+1; dark:+2; flat:+3; img += 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 # 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)) datetime_obs = datetime.fromtimestamp(timestamp_obs) date_obs = datetime_obs.strftime("%y%m%d") time_obs = datetime_obs.strftime("%H%M%S") self.outputCal( img=crmap_gsimg, ra_cen=ra_cen, dec_cen=dec_cen, img_rot=img_rot, im_type='CRS', pointing_ID=pointing_ID, date_obs=date_obs, time_obs=time_obs, output_dir=chip_output.subdir, exptime=150.) 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" and pointing_type=='CAL': 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) if config["cosmic_ray"].lower() == "y": if config["cray_differ"].lower() == "y": cr_map = effects.produceCR_Map( xLen=self.npix_x, yLen=self.npix_y, exTime=0.01, cr_pixelRatio=0.003*0.01/self.exptime, 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 # 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)) datetime_obs = datetime.fromtimestamp(timestamp_obs) date_obs = datetime_obs.strftime("%y%m%d") time_obs = datetime_obs.strftime("%H%M%S") timestamp_obs += 5 * 60 self.outputCal( img=BiasCombImg, ra_cen=ra_cen, dec_cen=dec_cen, img_rot=img_rot, im_type='CLB', pointing_ID=pointing_ID, date_obs=date_obs, time_obs=time_obs, output_dir=chip_output.subdir, exptime=0.0) del BiasCombImg # Export combined (ncombine, Vignetting + PRNU) & single vignetting flat-field file if config["flat_output"].lower() == "y" and pointing_type=='CAL': 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.flat_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": if config["cray_differ"].lower() == "y": cr_map = effects.produceCR_Map( xLen=self.npix_x, yLen=self.npix_y, exTime=self.flat_exptime, cr_pixelRatio=0.003*self.flat_exptime/self.exptime, 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 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)) datetime_obs = datetime.fromtimestamp(timestamp_obs) date_obs = datetime_obs.strftime("%y%m%d") time_obs = datetime_obs.strftime("%H%M%S") timestamp_obs += 5 * 60 self.outputCal( img=FlatCombImg, ra_cen=ra_cen, dec_cen=dec_cen, img_rot=img_rot, im_type='CLF', pointing_ID=pointing_ID, date_obs=date_obs, time_obs=time_obs, output_dir=chip_output.subdir, exptime=self.flat_exptime) 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" and pointing_type=='CAL': 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=self.dark_exptime, ncombine=1, read_noise=self.read_noise, gain=1, seed_bias=SeedBiasNonuni+self.chipID) if config["cosmic_ray"].lower() == "y": if config["cray_differ"].lower() == "y": cr_map = effects.produceCR_Map( xLen=self.npix_x, yLen=self.npix_y, exTime=self.dark_exptime, cr_pixelRatio=0.003*self.dark_exptime/self.exptime, 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 datetime_obs = datetime.fromtimestamp(timestamp_obs) date_obs = datetime_obs.strftime("%y%m%d") time_obs = datetime_obs.strftime("%H%M%S") self.outputCal( img=crmap_gsimg, ra_cen=ra_cen, dec_cen=dec_cen, img_rot=img_rot, im_type='CRD', pointing_ID=pointing_ID, date_obs=date_obs, time_obs=time_obs, output_dir=chip_output.subdir, exptime=self.dark_exptime) del crmap_gsimg # 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)) datetime_obs = datetime.fromtimestamp(timestamp_obs) date_obs = datetime_obs.strftime("%y%m%d") time_obs = datetime_obs.strftime("%H%M%S") timestamp_obs += 5 * 60 self.outputCal( img=DarkCombImg, ra_cen=ra_cen, dec_cen=dec_cen, img_rot=img_rot, im_type='CLD', pointing_ID=pointing_ID, date_obs=date_obs, time_obs=time_obs, output_dir=chip_output.subdir, exptime=self.dark_exptime) del DarkCombImg # 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