import galsim import os import numpy as np 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 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 # TODO: use chipID to config individual chip? super().__init__() 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 = float(config["ins_effects"]["gain"]) self.bias_level = float(config["ins_effects"]["bias_level"]) self.overscan = 1000 self.exptime = float(config["obs_setting"]["exp_time"]) # second self.dark_exptime = float(config["ins_effects"]['dark_exptime']) self.flat_exptime = float(config["ins_effects"]['flat_exptime']) self.readout_time = float(config["ins_effects"]['readout_time']) self.full_well = int(config["ins_effects"]["full_well"]) self.logger = logger # 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.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)] 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 = [] with pkg_resources.path('ObservationSim.Instrument.data.sls_conf', slsconfs[0]) as conf_path: self.sls_conf.append(str(conf_path)) 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 if config["ins_effects"]["bright_fatter"] == True and self.survey_type == "photometric": self.sensor = galsim.SiliconSensor(strength=config["ins_effects"]["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 _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') 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']*mirror_eff, interpolant='linear') throughput = galsim.LookupTable(x=table['col1'], f=table['col2'], 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) with pkg_resources.path('ObservationSim.Instrument.data', "wfc-cr-attachpixel.dat") as cr_path: self.attachedSizes = np.loadtxt(cr_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 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', sky_map=None, tel=None, logger=None): 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(config["ins_effects"]["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 seed = int(config["random_seeds"]["seed_poisson"]) + pointing_ID*30 + self.chipID rng_poisson = galsim.BaseDeviate(seed) poisson_noise = galsim.PoissonNoise(rng_poisson, sky_level=0.) # Add sky background if sky_map is None: sky_map = filt.getSkyNoise(exptime=self.exptime) sky_map = sky_map * np.ones_like(img.array) sky_map = galsim.Image(array=sky_map) # Apply Poisson noise to the sky map # (NOTE): only for photometric chips # since it utilize the photon shooting # to draw stamps if self.survey_type == "photometric": sky_map.addNoise(poisson_noise) elif img.array.shape != sky_map.shape: raise ValueError("The shape img and sky_map must be equal.") elif tel is not None: # If sky_map is given in flux sky_map = sky_map * tel.pupil_area * self.exptime if config["ins_effects"]["add_back"] == True: img += sky_map del sky_map # Apply flat-field large scale structure for one chip if config["ins_effects"]["flat_fielding"] == True: if self.logger is not None: self.logger.info(" Creating and applying Flat-Fielding") msg = str(img.bounds) self.logger.info(msg) else: print(" Creating and applying Flat-Fielding", flush=True) print(img.bounds, flush=True) flat_img = effects.MakeFlatSmooth( img.bounds, int(config["random_seeds"]["seed_flat"])) flat_normal = flat_img / np.mean(flat_img.array) if self.survey_type == "photometric": img *= flat_normal del flat_normal if config["output_setting"]["flat_output"] == False: del flat_img # Apply Shutter-effect for one chip if config["ins_effects"]["shutter_effect"] == True: if self.logger is not None: self.logger.info(" Apply shutter effect") else: 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 if self.survey_type == "photometric": img *= shuttimg if config["output_setting"]["shutter_output"] == True: # 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 # 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=='MS': if self.logger is not None: self.logger.info((" Adding Cosmic-Ray")) else: print(" Adding Cosmic-Ray", flush=True) cr_map, cr_event_num = effects.produceCR_Map( xLen=self.npix_x, yLen=self.npix_y, exTime=self.exptime+0.5*self.readout_time, cr_pixelRatio=0.003*(self.exptime+0.5*self.readout_time)/600., 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=self.exptime) del crmap_gsimg # Apply PRNU effect and output PRNU flat file: if config["ins_effects"]["prnu_effect"] == True: if self.logger is not None: self.logger.info(" Applying PRNU effect") else: 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["random_seeds"]["seed_prnu"]+self.chipID)) img *= prnu_img 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*(self.exptime+0.5*self.readout_time))) img.addNoise(dark_noise) # 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) # 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") if config["ins_effects"]["bias_16channel"] == True: img = effects.AddBiasNonUniform16(img, bias_level=float(config["ins_effects"]["bias_level"]), nsecy = 2, nsecx=8, seed=SeedBiasNonuni+self.chipID, logger=self.logger) elif config["ins_effects"]["bias_16channel"] == False: img += self.bias_level # Apply Nonlinearity on the chip image if config["ins_effects"]["non_linear"] == True: if self.logger is not None: self.logger.info(" Applying Non-Linearity on the chip image") else: print(" Applying Non-Linearity on the chip image", flush=True) img = effects.NonLinearity(GSImage=img, beta1=5.e-7, beta2=0) # Apply CCD Saturation & Blooming if config["ins_effects"]["saturbloom"] == True: if self.logger is not None: self.logger.info(" Applying CCD Saturation & Blooming") else: print(" Applying CCD Saturation & Blooming") img = effects.SaturBloom(GSImage=img, nsect_x=1, nsect_y=1, fullwell=fullwell) # Apply CTE Effect if config["ins_effects"]["cte_trail"] == True: if self.logger is not None: self.logger.info(" Apply CTE Effect") else: print(" Apply CTE Effect") img = effects.CTE_Effect(GSImage=img, threshold=27) # 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 if self.logger is not None: self.logger.info(" Applying Gain (and 16 channel non-uniformity) & Quantization") else: print(" Applying Gain (and 16 channel non-uniformity) & Quantization", flush=True) if config["ins_effects"]["gain_16channel"] == True: img = effects.ApplyGainNonUniform16( img, gain=self.gain, nsecy = 2, nsecx=8, 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["ins_effects"]["NBias"]) for i in range(NBias): BiasCombImg, BiasTag = effects.MakeBiasNcomb( self.npix_x, self.npix_y, bias_level=float(config["ins_effects"]["bias_level"]), ncombine=1, read_noise=self.read_noise, gain=1, seed=SeedBiasNonuni+self.chipID, logger=self.logger) # 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 # 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) # Apply Bad lines if config["ins_effects"]["add_badcolumns"] == True: BiasCombImg = effects.BadColumns(BiasCombImg-float(config["ins_effects"]["bias_level"])+5, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger) + float(config["ins_effects"]["bias_level"])-5 BiasCombImg = effects.ApplyGainNonUniform16(BiasCombImg, gain=self.gain, nsecy = 2, nsecx=8, 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) # 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 += 10 * 60 self.outputCal( img=BiasCombImg, ra_cen=ra_cen, dec_cen=dec_cen, img_rot=img_rot, im_type='BIAS', 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["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["ins_effects"]["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 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) # 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) # 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 = 2, nsecx=8, 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 = effects.ApplyGainNonUniform16(FlatCombImg, gain=self.gain, nsecy = 2, nsecx=8, 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) # 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 += 10 * 60 self.outputCal( img=FlatCombImg, ra_cen=ra_cen, dec_cen=dec_cen, img_rot=img_rot, im_type='FLAT', 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["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["ins_effects"]["NDark"]) if config["ins_effects"]["add_bias"] == True: biaslevel = self.bias_level overscan = biaslevel-2 elif config["ins_effects"]["add_bias"] == True: 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 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["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) # 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) # 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 = 2, nsecx=8, 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 + 2 rng_readout = galsim.BaseDeviate(seed) readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise) DarkCombImg.addNoise(readout_noise) DarkCombImg = effects.ApplyGainNonUniform16( DarkCombImg, gain=self.gain, nsecy = 2, nsecx=8, seed=SeedGainNonuni+self.chipID, logger=self.logger) # 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 += 10 * 60 self.outputCal( img=DarkCombImg, ra_cen=ra_cen, dec_cen=dec_cen, img_rot=img_rot, im_type='DARK', pointing_ID=pointing_ID, date_obs=date_obs, time_obs=time_obs, output_dir=chip_output.subdir, exptime=self.dark_exptime) del DarkCombImg # img = galsim.ImageUS(img) # # 16 output channel, with each a single image file # if config["ins_effects"]["readout16"] == True: # 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