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 = float(config["gain"]) self.bias_level = float(config["bias_level"]) self.overscan = 1000 self.exptime = float(config["exp_time"]) # second self.dark_exptime = float(config['dark_exptime']) self.flat_exptime = float(config['flat_exptime']) self.readout_time = float(config['readout_time']) # 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 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): 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 # Add sky background if sky_map is None: sky_map = filt.getSkyNoise(exptime=self.exptime) 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["add_back"].lower() == "y": img += sky_map del sky_map # 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 # Add Poisson noise seed = int(config["seed_poisson"]) + pointing_ID*30 + self.chipID rng_poisson = galsim.BaseDeviate(seed) poisson_noise = galsim.PoissonNoise(rng_poisson, sky_level=0.) img.addNoise(poisson_noise) # Add cosmic-rays if config["cosmic_ray"].lower() == "y" and pointing_type=='MS': 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)/150., 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 # 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 if config["add_dark"].lower() == "y": 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["add_badcolumns"].lower() == "y": img = effects.BadColumns(img, seed=SeedBadColumns, chipid=self.chipID) # Add Bias level if config["add_bias"].lower() == "y": print(" Adding Bias level and 16-channel non-uniformity") img = effects.AddBiasNonUniform16(img, bias_level=float(config["bias_level"]), nsecy = 2, nsecx=8, seed=SeedBiasNonuni+self.chipID) # 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=5.e-7, beta2=0) # 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) # Add Read-out Noise if config["add_readout"].lower() == "y": seed = int(config["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 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() ###################################################################################### # Output images for calibration pointing ###################################################################################### # 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, 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)/(self.exptime+0.5*self.readout_time), 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=5.e-7, beta2=0) # Apply Bad lines if config["add_badcolumns"].lower() == "y": BiasCombImg = effects.BadColumns(BiasCombImg-float(config["bias_level"])+5, seed=SeedBadColumns, chipid=self.chipID) + float(config["bias_level"])-5 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 += 10 * 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+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 ) if config["cosmic_ray"].lower() == "y": if config["cray_differ"].lower() == "y": 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)/(self.exptime+0.5*self.readout_time), 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=5.e-7, beta2=0) if config["cte_trail"].lower() == "y": 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["add_badcolumns"].lower() == "y": FlatCombImg = effects.BadColumns(FlatCombImg, seed=SeedBadColumns, chipid=self.chipID) # Add Bias level if config["add_bias"].lower() == "y": print(" Adding Bias level and 16-channel non-uniformity") # img += float(config["bias_level"]) FlatCombImg = effects.AddBiasNonUniform16(FlatCombImg, bias_level=biaslevel, nsecy = 2, nsecx=8, seed=SeedBiasNonuni+self.chipID) # Add Read-out Noise if config["add_readout"].lower() == "y": seed = int(config["seed_readout"]) + pointing_ID*30 + self.chipID 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) # 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='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=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) if config["cosmic_ray"].lower() == "y": if config["cray_differ"].lower() == "y": 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)/(self.exptime+0.5*self.readout_time), 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=5.e-7, beta2=0) if config["cte_trail"].lower() == "y": 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["add_badcolumns"].lower() == "y": DarkCombImg = effects.BadColumns(DarkCombImg, seed=SeedBadColumns, chipid=self.chipID) # Add Bias level if config["add_bias"].lower() == "y": print(" Adding Bias level and 16-channel non-uniformity") # img += float(config["bias_level"]) DarkCombImg = effects.AddBiasNonUniform16(DarkCombImg, bias_level=biaslevel, nsecy = 2, nsecx=8, seed=SeedBiasNonuni+self.chipID) # Add Read-out Noise if config["add_readout"].lower() == "y": seed = int(config["seed_readout"]) + pointing_ID*30 + self.chipID 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) # 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='CLD', 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["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