import galsim import os import numpy as np import pickle import json import ObservationSim.Instrument._util as _util 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 from ObservationSim.Instrument._util import rotate_conterclockwise from ObservationSim.Instrument.Chip import ChipUtils as chip_utils from ObservationSim.Instrument.Chip.libCTI.CTI_modeling import CTI_sim 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 super().__init__() self.nsecy = 2 self.nsecx = 8 self.gain_channel = np.ones(self.nsecy* self.nsecx) self._set_attributes_from_config(config) self.logger = logger # A chip ID must be assigned self.chipID = int(chipID) self.chip_name = str(chipID).rjust(2, '0') # Get corresponding filter info self.filter_id, self.filter_type = self.getChipFilter() self.survey_type = self._getSurveyType() if self.filter_type != "FGS": self._getChipRowCol() # Set the relavent specs for detectors try: with pkg_resources.files('ObservationSim.Instrument.data.ccd').joinpath("chip_definition.json") as chip_definition: with open(chip_definition, "r") as f: chip_dict = json.load(f)[str(self.chipID)] except AttributeError: with pkg_resources.path('ObservationSim.Instrument.data.ccd', "chip_definition.json") as chip_definition: with open(chip_definition, "r") as f: chip_dict = json.load(f)[str(self.chipID)] for key in chip_dict: setattr(self, key, chip_dict[key]) if self.filter_type == "FGS": if ("field_dist" in config) and (config["ins_effects"]["field_dist"]) == False: self.fdModel = None else: fgs_name = self.chip_name[0:4] try: with pkg_resources.files('ObservationSim.Instrument.data.field_distortion').joinpath("FieldDistModelGlobal_pr4_%s.pickle"%(fgs_name.lower())) as field_distortion: with open(field_distortion, "rb") as f: self.fdModel = pickle.load(f) except AttributeError: with pkg_resources.path('ObservationSim.Instrument.data.field_distortion', "FieldDistModelGlobal_pr4_%s.pickle"%(fgs_name.lower())) as field_distortion: with open(field_distortion, "rb") as f: self.fdModel = pickle.load(f) else: # Get the corresponding field distortion model if ("field_dist" in config) and (config["ins_effects"]["field_dist"] == False): self.fdModel = None else: try: with pkg_resources.files('ObservationSim.Instrument.data.field_distortion').joinpath("FieldDistModel_v2.0.pickle") as field_distortion: with open(field_distortion, "rb") as f: self.fdModel = pickle.load(f) except AttributeError: with pkg_resources.path('ObservationSim.Instrument.data.field_distortion', "FieldDistModelGlobal_mainFP_v1.0.pickle") as field_distortion: with open(field_distortion, "rb") as f: self.fdModel = pickle.load(f) # Get boundary (in pix) self.bound = self.getChipLim() self.ccdEffCurve_dir = ccdEffCurve_dir self.CRdata_dir = CRdata_dir slsconfs = chip_utils.getChipSLSConf(chipID=self.chipID) if np.size(slsconfs) == 1: try: with pkg_resources.files('ObservationSim.Instrument.data.sls_conf').joinpath(slsconfs) as conf_path: self.sls_conf = str(conf_path) except AttributeError: 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 = [] try: with pkg_resources.files('ObservationSim.Instrument.data.sls_conf').joinpath(slsconfs[0]) as conf_path: self.sls_conf.append(str(conf_path)) except AttributeError: with pkg_resources.path('ObservationSim.Instrument.data.sls_conf', slsconfs[0]) as conf_path: self.sls_conf.append(str(conf_path)) try: with pkg_resources.files('ObservationSim.Instrument.data.sls_conf').joinpath(slsconfs[1]) as conf_path: self.sls_conf.append(str(conf_path)) except AttributeError: 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 model if "bright_fatter" in config["ins_effects"] and config["ins_effects"]["bright_fatter"] == True and self.survey_type == "photometric": self.sensor = galsim.SiliconSensor(strength=self.df_strength, treering_func=treering_func) else: self.sensor = galsim.Sensor() self.flat_cube = None # for spectroscopic flat field cube simulation def _set_attributes_from_config(self, config): # Default setting self.read_noise = 5.0 # e/pix self.dark_noise = 0.02 # e/pix/s self.rotate_angle = 0. self.overscan = 1000 # Override default values for key in ["gain", "bias_level, dark_exptime", "flat_exptime", "readout_time", "full_well", "read_noise", "dark_noise", "overscan"]: if key in config["ins_effects"]: setattr(self, key, config["ins_effects"][key]) 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 _util.SPEC_FILTERS: return "spectroscopic" elif self.filter_type in _util.PHOT_FILTERS: return "photometric" # elif self.filter_type in ["FGS"]: # return "FGS" 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', 'FGS']: filename = 'Astro_MB.txt' # TODO, need to switch to the right efficiency curvey for FGS CMOS if filter_type in ['i', 'z', 'y', 'GI']: filename = 'Basic_NIR.txt' try: with pkg_resources.files('ObservationSim.Instrument.data.ccd').joinpath(filename) as ccd_path: table = Table.read(ccd_path, format='ascii') except AttributeError: 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'], interpolant='linear') bandpass = galsim.Bandpass(throughput, wave_type='nm') return bandpass def _getCRdata(self): try: with pkg_resources.files('ObservationSim.Instrument.data').joinpath("wfc-cr-attachpixel.dat") as cr_path: self.attachedSizes = np.loadtxt(cr_path) except AttributeError: with pkg_resources.path('ObservationSim.Instrument.data', "wfc-cr-attachpixel.dat") as cr_path: self.attachedSizes = np.loadtxt(cr_path) def loadSLSFLATCUBE(self, flat_fn='flat_cube.fits'): try: with pkg_resources.files('ObservationSim.Instrument.data').joinpath(flat_fn) as data_path: flat_fits = fits.open(data_path, ignore_missing_simple=True) except AttributeError: with pkg_resources.path('ObservationSim.Instrument.data', flat_fn) as data_path: flat_fits = fits.open(data_path, ignore_missing_simple=True) fl = len(flat_fits) fl_sh = flat_fits[0].data.shape assert fl == 4, 'FLAT Field Cube is Not 4 layess!!!!!!!' self.flat_cube = np.zeros([fl, fl_sh[0], fl_sh[1]]) for i in np.arange(0, fl, 1): self.flat_cube[i] = flat_fits[i].data def getChipFilter(self, chipID=None): """Return the filter index and type for a given chip #(chipID) """ filter_type_list = _util.ALL_FILTERS if chipID == None: chipID = self.chipID # updated configurations if chipID>42 or chipID<1: raise ValueError("!!! Chip ID: [1,42]") 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" if chipID in range(31, 43): filter_type = 'FGS' 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 / spectroscopic observation. Parameters: chipID: int the index of the chip Returns: A galsim BoundsD object """ xmin, xmax, ymin, ymax = 1e10, -1e10, 1e10, -1e10 xcen = self.x_cen / self.pix_size ycen = self.y_cen / self.pix_size sign_x = [-1., 1., -1., 1.] sign_y = [-1., -1., 1., 1.] for i in range(4): x = xcen + sign_x[i] * self.npix_x / 2. y = ycen + sign_y[i] * self.npix_y / 2. x, y = _util.rotate_conterclockwise(x0=xcen, y0=ycen, x=x, y=y, angle=self.rotate_angle) xmin, xmax = min(xmin, x), max(xmax, x) ymin, ymax = min(ymin, y), max(ymax, y) return galsim.BoundsD(xmin, xmax, ymin, ymax) def getSkyCoverage(self, wcs): # print("In getSkyCoverage: xmin = %.3f, xmax = %.3f, ymin = %.3f, ymax = %.3f"%(self.bound.xmin, self.bound.xmax, self.bound.ymin, self.bound.ymax)) 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=None, dec_obj=None, x_image=None, y_image=None, wcs=None, margin=1): # magin in number of pix if (ra_obj is not None) and (dec_obj is not None): if wcs is None: wcs = self.img.wcs pos_obj = wcs.toImage(galsim.CelestialCoord(ra=ra_obj*galsim.degrees, dec=dec_obj*galsim.degrees)) x_image, y_image = pos_obj.x, pos_obj.y elif (x_image is None) or (y_image is None): raise ValueError("Either (ra_obj, dec_obj) or (x_image, y_image) should be given") xmin, xmax = self.bound.xmin - margin, self.bound.xmax + margin ymin, ymax = self.bound.ymin - margin, self.bound.ymax + margin if (x_image - xmin) * (x_image - xmax) > 0.0 or (y_image - ymin) * (y_image - 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 addEffects(self, config, img, chip_output, filt, ra_cen, dec_cen, img_rot, exptime=150., pointing_ID=0, timestamp_obs=1621915200, pointing_type='SCI', sky_map=None, post_flash_map=None, tel=None, logger=None): # Set random seeds 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(self.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 rng_poisson, poisson_noise = chip_utils.get_poisson( seed=int(config["random_seeds"]["seed_poisson"]) + pointing_ID*30 + self.chipID, sky_level=0.) # Add sky background if config["ins_effects"]["add_back"] == True: img, sky_map = chip_utils.add_sky_background(img=img, filt=filt, exptime=exptime, sky_map=sky_map, tel=tel) del sky_map # Apply flat-field large scale structure for one chip if config["ins_effects"]["flat_fielding"] == True: chip_utils.log_info(msg=" Creating and applying Flat-Fielding", logger=self.logger) chip_utils.log_info(msg=str(img.bounds), logger=self.logger) flat_img, flat_normal = chip_utils.get_flat(img=img, seed=int(config["random_seeds"]["seed_flat"])) if self.survey_type == "photometric": img *= flat_normal del flat_normal if config["output_setting"]["flat_output"] == False: del flat_img if post_flash_map is not None: img = img + post_flash_map # Apply Shutter-effect for one chip if config["ins_effects"]["shutter_effect"] == True: chip_utils.log_info(msg=" Apply shutter effect", logger=self.logger) 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=='SCI': chip_utils.log_info(msg=" Adding Cosmic-Ray", logger=self.logger) img, crmap_gsimg, cr_event_num = chip_utils.add_cosmic_rays(img=img, chip=self, exptime=exptime, seed=SeedCosmicRay+pointing_ID*30+self.chipID) chip_utils.outputCal( chip=self, img=crmap_gsimg, ra_cen=ra_cen, dec_cen=dec_cen, img_rot=img_rot, im_type='CRS', pointing_ID=pointing_ID, output_dir=chip_output.subdir, exptime=exptime, project_cycle=config["project_cycle"], run_counter=config["run_counter"], timestamp=timestamp_obs) del crmap_gsimg # Apply PRNU effect and output PRNU flat file: if config["ins_effects"]["prnu_effect"] == True: chip_utils.log_info(msg=" Applying PRNU effect", logger=self.logger) img, prnu_img = chip_utils.add_PRNU(img=img, chip=self, seed=int(config["random_seeds"]["seed_prnu"]+self.chipID)) 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*(exptime+0.5*self.readout_time))) # img.addNoise(dark_noise) # Add dark current & Poisson noise InputDark = False if config["ins_effects"]["add_dark"] == True: if InputDark: img = chip_utils.add_inputdark(img=img, chip=self, exptime=exptime) else: img, _ = chip_utils.add_poisson(img=img, chip=self, exptime=exptime, poisson_noise=poisson_noise) else: img, _ = chip_utils.add_poisson(img=img, chip=self, exptime=exptime, poisson_noise=poisson_noise, dark_noise=0.) # Add diffusion & brighter-fatter effects if config["ins_effects"]["bright_fatter"] == True: img = chip_utils.add_brighter_fatter(img=img) # 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) # Apply Nonlinearity on the chip image if config["ins_effects"]["non_linear"] == True: chip_utils.log_info(msg=" Applying Non-Linearity on the chip image", logger=self.logger) img = effects.NonLinearity(GSImage=img, beta1=5.e-7, beta2=0) # Apply CCD Saturation & Blooming if config["ins_effects"]["saturbloom"] == True: chip_utils.log_info(msg=" Applying CCD Saturation & Blooming", logger=self.logger) img = effects.SaturBloom(GSImage=img, nsect_x=1, nsect_y=1, fullwell=fullwell) # Apply CTE Effect ###if config["ins_effects"]["cte_trail"] == True: ### chip_utils.log_info(msg=" Apply CTE Effect", logger=self.logger) ### img = effects.CTE_Effect(GSImage=img, threshold=27) pre1 = self.prescan_x #27 over1= self.overscan_x #71 pre2 = self.prescan_y #0 #4 over2= self.overscan_y #84 #80 if config["ins_effects"]["cte_trail"] == True: chip_utils.log_info(msg=" Apply CTE Effect", logger=self.logger) ###img = effects.CTE_Effect(GSImage=img, threshold=27) ###CTI_modeling ### 2*8 -> 1*16 img-layout img = chip_utils.formatOutput(GSImage=img) self.nsecy = 1 self.nsecx = 16 img_arr = img.array ny, nx = img_arr.shape dx = int(nx/self.nsecx) dy = int(ny/self.nsecy) newimg = galsim.Image(nx, int(ny+over2), init_value=0) for ichannel in range(16): print('\n***add CTI effects: pointing-{:} chip-{:} channel-{:}***'.format(pointing_ID, self.chipID, ichannel+1)) #nx,ny,noverscan,nsp,nmax = 4608,4616,84,3,10 noverscan,nsp,nmax = over2,3,10 beta,w,c = 0.478,84700,0 t = np.array([0.74,7.7,37],dtype=np.float32) rho_trap = np.array([0.6,1.6,1.4],dtype=np.float32) trap_seeds = np.array([0,1000,10000],dtype=np.int32) + ichannel + self.chipID*16 release_seed = 50 + ichannel + pointing_ID*30 + self.chipID*16 newimg.array[:, 0+ichannel*dx:dx+ichannel*dx] = CTI_sim(img_arr[:, 0+ichannel*dx:dx+ichannel*dx],dx,dy,noverscan,nsp,nmax,beta,w,c,t,rho_trap,trap_seeds,release_seed) newimg.wcs = img.wcs del img img = newimg ### 1*16 -> 2*8 img-layout img = chip_utils.formatRevert(GSImage=img) self.nsecy = 2 self.nsecx = 8 ### prescan & overscan if config["ins_effects"]["add_prescan"] == True: chip_utils.log_info(msg=" Apply pre/over-scan", logger=self.logger) if config["ins_effects"]["cte_trail"] == False: img = chip_utils.AddPreScan(GSImage=img, pre1=pre1, pre2=pre2, over1=over1, over2=over2) if config["ins_effects"]["cte_trail"] == True: img = chip_utils.AddPreScan(GSImage=img, pre1=pre1, pre2=pre2, over1=over1, over2=0) ### 1*16 output if config["ins_effects"]["format_output"] == True: chip_utils.log_info(msg=" Apply 1*16 format", logger=self.logger) img = chip_utils.formatOutput(GSImage=img) self.nsecy = 1 self.nsecx = 16 # Add Bias level if config["ins_effects"]["add_bias"] == True: chip_utils.log_info(msg=" Adding Bias level and 16-channel non-uniformity", logger=self.logger) if config["ins_effects"]["bias_16channel"] == True: img = effects.AddBiasNonUniform16(img, bias_level=float(self.bias_level), nsecy = self.nsecy, nsecx=self.nsecx, seed=SeedBiasNonuni+self.chipID, logger=self.logger) elif config["ins_effects"]["bias_16channel"] == False: img += self.bias_level # 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 chip_utils.log_info(msg=" Applying Gain (and 16 channel non-uniformity) & Quantization", logger=self.logger) if config["ins_effects"]["gain_16channel"] == True: img, self.gain_channel = effects.ApplyGainNonUniform16( img, gain=self.gain, nsecy = self.nsecy, nsecx=self.nsecx, 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["output_setting"]["NBias"]) for i in range(NBias): BiasCombImg, BiasTag = effects.MakeBiasNcomb( self.npix_x, self.npix_y, bias_level=float(self.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(self.bias_level)+5, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger) + float(self.bias_level)-5 BiasCombImg, self.gain_channel = 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) timestamp_obs += 10 * 60 chip_utils.outputCal( chip=self, img=BiasCombImg, ra_cen=ra_cen, dec_cen=dec_cen, img_rot=img_rot, im_type='BIAS', pointing_ID=pointing_ID, output_dir=chip_output.subdir, exptime=0.0, project_cycle=config["project_cycle"], run_counter=config["run_counter"], timestamp=timestamp_obs) 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["output_setting"]["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, self.gain_channel = 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) timestamp_obs += 10 * 60 chip_utils.outputCal( chip=self, img=FlatCombImg, ra_cen=ra_cen, dec_cen=dec_cen, img_rot=img_rot, im_type='FLAT', pointing_ID=pointing_ID, output_dir=chip_output.subdir, exptime=self.flat_exptime, project_cycle=config["project_cycle"], run_counter=config["run_counter"], timestamp=timestamp_obs) 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["output_setting"]["NDark"]) 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 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 chip_utils.outputCal( chip=self, img=crmap_gsimg, ra_cen=ra_cen, dec_cen=dec_cen, img_rot=img_rot, im_type='CRD', pointing_ID=pointing_ID, output_dir=chip_output.subdir, exptime=self.dark_exptime, project_cycle=config["project_cycle"], run_counter=config["run_counter"], timestamp=timestamp_obs) 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, self.gain_channel = 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) timestamp_obs += 10 * 60 chip_utils.outputCal( chip=chip, img=DarkCombImg, ra_cen=ra_cen, dec_cen=dec_cen, img_rot=img_rot, im_type='DARK', pointing_ID=pointing_ID, output_dir=chip_output.subdir, exptime=self.dark_exptime, project_cycle=config["project_cycle"], run_counter=config["run_counter"], timestamp = timestamp_obs) 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