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