diff --git a/ObservationSim/Config/ChipOutput.py b/ObservationSim/Config/ChipOutput.py index d3d8ff55e434aeee51c173959026c270fcb4e052..685e12d6c4f684553d10e1bafa6f5191207af808 100755 --- a/ObservationSim/Config/ChipOutput.py +++ b/ObservationSim/Config/ChipOutput.py @@ -37,8 +37,10 @@ class ChipOutput(object): fmt5 = "%10s %8.4f %8.4f %8.4f\n" self.hdr = hdr1 + hdr2 + hdr3 + hdr4 + hdr5 self.fmt = fmt1 + fmt2 + fmt3 + fmt4 + fmt5 + print("pointing_type = %s\n"%(pointing_type)) if pointing_type == 'MS': self.cat = open(os.path.join(self.subdir, self.cat_name), "w") + print("Creating catalog file %s ...\n"%(os.path.join(self.subdir, self.cat_name))) self.cat.write(self.hdr) def updateHDR(self, hdr): diff --git a/ObservationSim/Config/__pycache__/ChipOutput.cpython-38.pyc b/ObservationSim/Config/__pycache__/ChipOutput.cpython-38.pyc index 8629f5acd8af472c3d00ab6b178614b1411845fb..8d63bc483452a94c3210f9a5899706ba9d74c33d 100644 Binary files a/ObservationSim/Config/__pycache__/ChipOutput.cpython-38.pyc and b/ObservationSim/Config/__pycache__/ChipOutput.cpython-38.pyc differ diff --git a/ObservationSim/Instrument/Chip/Chip.py b/ObservationSim/Instrument/Chip/Chip.py index f3d518f310956b92bc2ca239c60ba1f2dbadc591..d77e6c79730b74b6a411096ef11e5db95048a37d 100755 --- a/ObservationSim/Instrument/Chip/Chip.py +++ b/ObservationSim/Instrument/Chip/Chip.py @@ -30,16 +30,13 @@ class Chip(FocalPlane): self.read_noise = 5.0 # e/pix self.dark_noise = 0.02 # e/pix/s self.pix_scale = 0.074 # pixel scale - # self.gain = 1.0 - # self.bias_level = 1000 # e-/pix self.gain = float(config["gain"]) self.bias_level = float(config["bias_level"]) self.overscan = 1000 - self.exptime = 150 # second - # self.dark_exptime = 300 # second - # self.flat_exptime = 150 # second + self.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) @@ -234,60 +231,6 @@ class Chip(FocalPlane): 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'] @@ -351,7 +294,7 @@ class Chip(FocalPlane): 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'): + 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"]) @@ -368,6 +311,17 @@ class Chip(FocalPlane): else: BoolDeadPix = False + # Add sky background + if sky_map == 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) @@ -392,30 +346,19 @@ class Chip(FocalPlane): 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 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 = effects.produceCR_Map( xLen=self.npix_x, yLen=self.npix_y, - exTime=exptime, - cr_pixelRatio=0.003, + exTime=self.exptime+0.5*self.readout_time, + cr_pixelRatio=0.003*(1+0.5*self.readout_time/self.exptime), gain=self.gain, attachedSizes=self.attachedSizes, seed=SeedCosmicRay+pointing_ID*30+self.chipID) # seed: obj-imaging:+0; bias:+1; dark:+2; flat:+3; @@ -442,15 +385,78 @@ class Chip(FocalPlane): 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 += float(config["bias_level"]) 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) @@ -466,7 +472,7 @@ class Chip(FocalPlane): cr_map = effects.produceCR_Map( xLen=self.npix_x, yLen=self.npix_y, exTime=0.01, - cr_pixelRatio=0.003*0.01/self.exptime, + 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) @@ -477,7 +483,11 @@ class Chip(FocalPlane): # 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.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, @@ -493,7 +503,7 @@ class Chip(FocalPlane): datetime_obs = datetime.fromtimestamp(timestamp_obs) date_obs = datetime_obs.strftime("%y%m%d") time_obs = datetime_obs.strftime("%H%M%S") - timestamp_obs += 5 * 60 + timestamp_obs += 10 * 60 self.outputCal( img=BiasCombImg, ra_cen=ra_cen, @@ -517,7 +527,7 @@ class Chip(FocalPlane): elif config["add_bias"].lower() == "n": biaslevel = 0 overscan = 0 - darklevel = self.dark_noise*self.flat_exptime + 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( @@ -526,15 +536,15 @@ class Chip(FocalPlane): read_noise=self.read_noise, gain=1, overscan=overscan, - biaslevel=biaslevel, + biaslevel=0, seed_bias=SeedDefective+self.chipID ) if config["cosmic_ray"].lower() == "y": if config["cray_differ"].lower() == "y": cr_map = effects.produceCR_Map( xLen=self.npix_x, yLen=self.npix_y, - exTime=self.flat_exptime, - cr_pixelRatio=0.003*self.flat_exptime/self.exptime, + 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) @@ -544,19 +554,35 @@ class Chip(FocalPlane): 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) + 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 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) + # 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, @@ -569,7 +595,7 @@ class Chip(FocalPlane): datetime_obs = datetime.fromtimestamp(timestamp_obs) date_obs = datetime_obs.strftime("%y%m%d") time_obs = datetime_obs.strftime("%H%M%S") - timestamp_obs += 5 * 60 + timestamp_obs += 10 * 60 self.outputCal( img=FlatCombImg, ra_cen=ra_cen, @@ -601,15 +627,15 @@ class Chip(FocalPlane): for i in range(NDark): DarkCombImg, DarkTag = effects.MakeDarkNcomb( self.npix_x, self.npix_y, - overscan=overscan, bias_level=biaslevel, darkpsec=0.02, exptime=self.dark_exptime, + 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 = effects.produceCR_Map( xLen=self.npix_x, yLen=self.npix_y, - exTime=self.dark_exptime, - cr_pixelRatio=0.003*self.dark_exptime/self.exptime, + 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) @@ -638,19 +664,35 @@ class Chip(FocalPlane): # 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) + 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 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) + # 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, @@ -667,7 +709,7 @@ class Chip(FocalPlane): datetime_obs = datetime.fromtimestamp(timestamp_obs) date_obs = datetime_obs.strftime("%y%m%d") time_obs = datetime_obs.strftime("%H%M%S") - timestamp_obs += 5 * 60 + timestamp_obs += 10 * 60 self.outputCal( img=DarkCombImg, ra_cen=ra_cen, @@ -680,41 +722,6 @@ class Chip(FocalPlane): output_dir=chip_output.subdir, exptime=self.dark_exptime) del DarkCombImg - - # Apply Nonlinearity on the chip image - if config["non_linear"].lower() == "y": - print(" Applying Non-Linearity on the chip image", flush=True) - img = effects.NonLinearity(GSImage=img, beta1=1.e-7, beta2=1.e-10) - - # Apply CCD Saturation & Blooming - if config["saturbloom"].lower() == "y": - print(" Applying CCD Saturation & Blooming") - img = effects.SaturBloom(GSImage=img, nsect_x=1, nsect_y=1, fullwell=fullwell) - - # Apply CTE Effect - if config["cte_trail"].lower() == "y": - print(" Apply CTE Effect") - img = effects.CTE_Effect(GSImage=img, threshold=27) - - # Apply Bad lines - if config["add_badcolumns"].lower() == "y": - img = effects.BadColumns(img, seed=SeedBadColumns, chipid=self.chipID) - - # Add Hot Pixels or/and Dead Pixels - rgbadpix = Generator(PCG64(int(SeedDefective+self.chipID))) - badfraction = 5E-5*(rgbadpix.random()*0.5+0.7) - img = effects.DefectivePixels(img, IfHotPix=BoolHotPix, IfDeadPix=BoolDeadPix, fraction=badfraction, seed=SeedDefective+self.chipID, biaslevel=self.bias_level) - - - # Apply Gain & Quantization - print(" Applying Gain (and 16 channel non-uniformity) & Quantization", flush=True) - img = effects.ApplyGainNonUniform16( - img, gain=self.gain, - nsecy = 2, nsecx=8, - seed=SeedGainNonuni+self.chipID) - img.array[img.array > 65535] = 65535 - img.replaceNegative(replace_value=0) - img.quantize() # img = galsim.ImageUS(img) # # 16 output channel, with each a single image file diff --git a/ObservationSim/Instrument/Chip/Effects.py b/ObservationSim/Instrument/Chip/Effects.py index 63844062857161749651c81bf9e89e917a0d7ac5..441359ef9fa82e90bcb0f6995fb3a02c70d7e21f 100644 --- a/ObservationSim/Instrument/Chip/Effects.py +++ b/ObservationSim/Instrument/Chip/Effects.py @@ -5,6 +5,7 @@ from numpy.core.fromnumeric import mean, size from numpy.random import Generator, PCG64 import math from numba import jit +from astropy import stats def AddOverscan(GSImage, overscan=1000, gain=1, widthl=27, widthr=27, widtht=8, widthb=8, read_noise=5): @@ -31,7 +32,7 @@ def AddOverscan(GSImage, overscan=1000, gain=1, widthl=27, widthr=27, widtht=8, return newimg -def DefectivePixels(GSImage, IfHotPix=True, IfDeadPix=True, fraction=1E-4, seed=20210304, biaslevel=500): +def DefectivePixels(GSImage, IfHotPix=True, IfDeadPix=True, fraction=1E-4, seed=20210304, biaslevel=0): # Also called bad pixels, including hot pixels and dead pixels # Hot Pixel > 20e-/s # Dead Pixel < 70%*Mean @@ -68,10 +69,11 @@ def DefectivePixels(GSImage, IfHotPix=True, IfDeadPix=True, fraction=1E-4, seed= return GSImage -def BadColumns(GSImage, seed=20210309, chipid=1): +def BadColumns(GSImage, seed=20240309, chipid=1): # Set bad column values ysize,xsize = GSImage.array.shape subarr = GSImage.array[int(ysize*0.1):int(ysize*0.12), int(xsize*0.1):int(xsize*0.12)] + subarr = stats.sigma_clip(subarr, sigma=4, cenfunc='median', maxiters=3, masked=False) meanimg = np.median(subarr) stdimg = np.std(subarr) seed += chipid @@ -80,15 +82,19 @@ def BadColumns(GSImage, seed=20210309, chipid=1): rgxpos = Generator(PCG64(int(seed*1.2))) rgdn = Generator(PCG64(int(seed*1.3))) - nbadsecA,nbadsecD = rgn.integers(low=1, high=5, size=2) + nbadsecA,nbadsecD = rgn.integers(low=1, high=5, size=2) collen = rgcollen.integers(low=int(ysize*0.1), high=int(ysize*0.7), size=(nbadsecA+nbadsecD)) xposit = rgxpos.integers(low=int(xsize*0.05), high=int(xsize*0.95), size=(nbadsecA+nbadsecD)) - signs = 2*rgdn.integers(0,2,size=(nbadsecA+nbadsecD))-1 - dn = rgdn.integers(low=meanimg*0.4, high=meanimg*0.8, size=(nbadsecA+nbadsecD))*signs + print(xposit+1) + # signs = 2*rgdn.integers(0,2,size=(nbadsecA+nbadsecD))-1 + if meanimg>0: + dn = rgdn.integers(low=meanimg*1.3+50, high=meanimg*2+150, size=(nbadsecA+nbadsecD)) #*signs + elif meanimg<0: + dn = rgdn.integers(low=meanimg*2-150, high=meanimg*1.3-50, size=(nbadsecA+nbadsecD)) #*signs for badcoli in range(nbadsecA): - GSImage.array[(ysize-collen[badcoli]):ysize,xposit[badcoli]:(xposit[badcoli]+1)] += np.random.random((collen[badcoli],1))*stdimg*3+dn[badcoli] + GSImage.array[(ysize-collen[badcoli]):ysize,xposit[badcoli]:(xposit[badcoli]+1)] = (np.abs(np.random.normal(0, stdimg*2, (collen[badcoli],1)))+dn[badcoli]) for badcoli in range(nbadsecD): - GSImage.array[0:collen[badcoli+nbadsecA],xposit[badcoli+nbadsecA]:(xposit[badcoli+nbadsecA]+1)] += np.random.random((collen[badcoli+nbadsecA],1))*stdimg*3+dn[badcoli+nbadsecA] + GSImage.array[0:collen[badcoli+nbadsecA],xposit[badcoli+nbadsecA]:(xposit[badcoli+nbadsecA]+1)] = (np.abs(np.random.normal(0, stdimg*2, (collen[badcoli+nbadsecA],1)))+dn[badcoli+nbadsecA]) return GSImage @@ -223,7 +229,7 @@ def PRNU_Img(xsize, ysize, sigma=0.01, seed=202101): return prnuimg -def NonLinearity(GSImage, beta1=1E-7, beta2=1E-10): +def NonLinearity(GSImage, beta1=5E-7, beta2=0): NonLinear_f = lambda x, beta_1, beta_2: x - beta_1*x*x + beta_2*x*x*x GSImage.applyNonlinearity(NonLinear_f, beta1, beta2) return GSImage diff --git a/ObservationSim/Instrument/Chip/__pycache__/Chip.cpython-38.pyc b/ObservationSim/Instrument/Chip/__pycache__/Chip.cpython-38.pyc index 714df698948c5c51688fad876725832c1f66b762..04dcf2fffaf8dbb23429b6504ae597b6bc5e9b96 100644 Binary files a/ObservationSim/Instrument/Chip/__pycache__/Chip.cpython-38.pyc and b/ObservationSim/Instrument/Chip/__pycache__/Chip.cpython-38.pyc differ diff --git a/ObservationSim/Instrument/Chip/__pycache__/Effects.cpython-38.pyc b/ObservationSim/Instrument/Chip/__pycache__/Effects.cpython-38.pyc index 44f89fdced2fa9a27886b382325bbac16bb0f113..5b8b9b9c4e96791918411bdc73a6fe5ccf036c3a 100644 Binary files a/ObservationSim/Instrument/Chip/__pycache__/Effects.cpython-38.pyc and b/ObservationSim/Instrument/Chip/__pycache__/Effects.cpython-38.pyc differ diff --git a/ObservationSim/Instrument/Filter.py b/ObservationSim/Instrument/Filter.py index 0ce2cb74e823b48c7be594e55e3ef759a8af5ca3..abad4059c3a902a68f0fb700f20f1d60ce5d2a7f 100755 --- a/ObservationSim/Instrument/Filter.py +++ b/ObservationSim/Instrument/Filter.py @@ -73,5 +73,5 @@ class Filter(object): def getPhotonE(self): return photonEnergy(self.effective_wavelength) - def getSkyNoise(self, exptime, gain): + def getSkyNoise(self, exptime, gain=1.): return self.sky_background * exptime / gain diff --git a/ObservationSim/Instrument/__pycache__/Filter.cpython-38.pyc b/ObservationSim/Instrument/__pycache__/Filter.cpython-38.pyc index 25130bb497d27f287065bb9c8502433ea6743b78..fcebf5167a91b96e35555e3025f1b38a8f86e165 100644 Binary files a/ObservationSim/Instrument/__pycache__/Filter.cpython-38.pyc and b/ObservationSim/Instrument/__pycache__/Filter.cpython-38.pyc differ diff --git a/ObservationSim/ObservationSim.py b/ObservationSim/ObservationSim.py index 117e95ae82b22e6549d72fa9d09f15f16f0fa3b0..ecbb76ef8ac75e3c8b1c25f732079ad1da571a9b 100755 --- a/ObservationSim/ObservationSim.py +++ b/ObservationSim/ObservationSim.py @@ -97,6 +97,7 @@ class Observation(object): # Load SED if obj.type == 'star': normF = chip.normF_star + #continue try: obj.load_SED( survey_type=chip.survey_type, @@ -108,6 +109,7 @@ class Observation(object): continue elif obj.type == 'galaxy': # or obj.type == quasar normF = chip.normF_galaxy + # continue obj.load_SED( sed_path=sed_dir, survey_type=chip.survey_type, @@ -116,6 +118,7 @@ class Observation(object): target_filt=filt) elif obj.type == 'quasar': normF = chip.normF_galaxy + # continue obj.load_SED( sed_path=sed_dir, survey_type=chip.survey_type, @@ -200,7 +203,7 @@ class Observation(object): # Detector Effects # =========================================================== - chip.img = chip.addNoise(config=self.config, tel=self.tel, filt=filt, img=chip.img, sky_map=sky_map) + # chip.img = chip.addNoise(config=self.config, tel=self.tel, filt=filt, img=chip.img, sky_map=sky_map) # chip.img = chip.addEffects(config=self.config, img=chip.img, chip_output=chip_output, filt=filt, pointing_ID=pointing_ID) # whether to output zero, dark, flat calibration images. @@ -214,7 +217,8 @@ class Observation(object): img_rot=img_rot, pointing_ID=pointing_ID, timestamp_obs=timestamp_obs, - pointing_type=pointing_type) + pointing_type=pointing_type, + sky_map=sky_map) if pointing_type == 'MS': datetime_obs = datetime.fromtimestamp(timestamp_obs) @@ -309,13 +313,20 @@ class Observation(object): sed_dir=self.path_dict["SED_dir"]) print("finished running chip#%d..."%(chip.chipID), flush=True) - def runExposure_MPI_PointingList(self, ra_cen=None, dec_cen=None, pRange=None, timestamp_obs=[1621915200], pointing_type=['MS'], img_rot=None, exptime=150., input_cat_name=None, shear_cat_file=None): + def runExposure_MPI_PointingList(self, ra_cen=None, dec_cen=None, pRange=None, timestamp_obs=np.array([1621915200]), pointing_type=np.array(['MS']), img_rot=None, exptime=150., input_cat_name=None, shear_cat_file=None): comm = MPI.COMM_WORLD ind_thread = comm.Get_rank() num_thread = comm.Get_size() nchips_per_fp = len(self.chip_list) + # TEMP + if len(timestamp_obs) == 1: + timestamp_obs = np.tile(timestamp_obs, len(ra_cen)) + pointing_type = np.tile(pointing_type, len(ra_cen)) + + timestamp_obs = timestamp_obs[pRange] + pointing_type = pointing_type[pRange] ra_cen = ra_cen[pRange] dec_cen = dec_cen[pRange] @@ -325,11 +336,6 @@ class Observation(object): else: pStart = 0 - # TEMP - if len(timestamp_obs) == 1: - timestamp_obs *= len(pRange) - pointing_type *= len(pRange) - for ipoint in range(len(ra_cen)): for ichip in range(nchips_per_fp): i = ipoint*nchips_per_fp + ichip diff --git a/ObservationSim/PSF/PSFInterp/PSFInterp.py b/ObservationSim/PSF/PSFInterp/PSFInterp.py index 0c4380317ca0a41bb2a1071a45cd6de9bc5b6cbd..92d79ca7efa517f25d641fdfe0833a85d3c53be8 100644 --- a/ObservationSim/PSF/PSFInterp/PSFInterp.py +++ b/ObservationSim/PSF/PSFInterp/PSFInterp.py @@ -171,8 +171,8 @@ class PSFInterp(PSFModel): img = galsim.ImageF(imPSF, scale=pixSize) gsp = galsim.GSParams(folding_threshold=folding_threshold) self.psf = galsim.InterpolatedImage(img, gsparams=gsp) - gaussPSF = galsim.Gaussian(sigma=0.04, gsparams=gsp) - self.psf = galsim.Convolve(gaussPSF, self.psf) + # gaussPSF = galsim.Gaussian(sigma=0.04, gsparams=gsp) + # self.psf = galsim.Convolve(gaussPSF, self.psf) return self.PSFspin(x=px/0.01, y=py/0.01) return imPSF diff --git a/ObservationSim/PSF/PSFInterp/__pycache__/PSFInterp.cpython-38.pyc b/ObservationSim/PSF/PSFInterp/__pycache__/PSFInterp.cpython-38.pyc index 85d2ea33d95a0877e0ad35d1afd592187a96f0c2..61560ce79bc8f050debdcea91cf431befc1037c6 100644 Binary files a/ObservationSim/PSF/PSFInterp/__pycache__/PSFInterp.cpython-38.pyc and b/ObservationSim/PSF/PSFInterp/__pycache__/PSFInterp.cpython-38.pyc differ diff --git a/ObservationSim/Pointing.py b/ObservationSim/Pointing.py index f38ebfb94f241f92166205faf059c50ca6cebcf4..138af6681c2977cee2bab84957b1b448e020d016 100644 --- a/ObservationSim/Pointing.py +++ b/ObservationSim/Pointing.py @@ -9,6 +9,8 @@ data_dir = "/data/simudata/CSSOSDataProductsSims/data/" config_file = os.path.join(work_dir, "ObservationSim.cfg") config = ReadConfig(config_file) +exp_time = float(config['exp_time']) + # Read Pointing list pointing_file = os.path.join(data_dir, "pointing10_20210202.dat") f = open(pointing_file, 'r') @@ -37,12 +39,16 @@ pointing_type = ['CAL']*ncal + pointing_type t0 = datetime(2021, 5, 25, 12, 0, 0) t = datetime.timestamp(t0) timestamp_obs = [] +delta_t = 10 # Time elapsed between exposures (minutes) for i in range(len(pointing_type)): timestamp_obs.append(t) if pointing_type[i] == 'CAL': - t += 3 * 5 * 60 # 3 calibration exposure + t += 3 * delta_t * 60 # 3 calibration exposure elif pointing_type[i] == 'MS': - t += 5 * 60 + t += delta_t * 60 + +timestamp_obs = np.array(timestamp_obs) +pointing_type = np.array(pointing_type) # Define the range of pointing list pRange = range(0, 2) \ No newline at end of file diff --git a/ObservationSim/__pycache__/ObservationSim.cpython-38.pyc b/ObservationSim/__pycache__/ObservationSim.cpython-38.pyc index 8927854296f4b745d0d52cb32028896a9b0f8857..77c1d72f8838d139cc4a15cbbd231ce09587a820 100644 Binary files a/ObservationSim/__pycache__/ObservationSim.cpython-38.pyc and b/ObservationSim/__pycache__/ObservationSim.cpython-38.pyc differ diff --git a/ObservationSim/__pycache__/Pointing.cpython-38.pyc b/ObservationSim/__pycache__/Pointing.cpython-38.pyc index 57ecf0ddccba0c3b01b99904dbffed5c19a22a82..22442f6ccad9d75cc79c4e5e4a811f12344183ad 100644 Binary files a/ObservationSim/__pycache__/Pointing.cpython-38.pyc and b/ObservationSim/__pycache__/Pointing.cpython-38.pyc differ diff --git a/ObservationSim/run.pbs b/ObservationSim/run.pbs index 4502bc834f784b43984657ab8ad111e5bb5e4511..395472a421badab7104fe5d545e10f7fff386703 100755 --- a/ObservationSim/run.pbs +++ b/ObservationSim/run.pbs @@ -7,7 +7,7 @@ ##mpdboot -n 10 -f ./mpd.hosts ##PBS -j oe -#PBS -l nodes=comput110:ppn=80 +#PBS -l nodes=comput108:ppn=80 #####PBS -q longq #PBS -q batch #PBS -u fangyuedong diff --git a/ObservationSim/runExposure.py b/ObservationSim/runExposure.py index 4f2eeff9a86c0b65600bed294ae8a0655f76c5d1..a90d30e0cd79e736c4be0c01b7c8929dde640ee5 100755 --- a/ObservationSim/runExposure.py +++ b/ObservationSim/runExposure.py @@ -20,4 +20,4 @@ gc.enable() # Testing run pointing list (using MPI) obs = Observation(work_dir=work_dir, data_dir=data_dir) # obs.runExposure_MPI_PointingList(ra_cen=pRA, dec_cen=pDEC, pRange=pRange, timestamp_obs=timestamp_obs, pointing_type=pointing_type) -obs.runExposure_MPI_PointingList(ra_cen=pRA, dec_cen=pDEC, pRange=pRange, timestamp_obs=timestamp_obs, pointing_type=pointing_type, exptime=150.) +obs.runExposure_MPI_PointingList(ra_cen=pRA, dec_cen=pDEC, pRange=pRange, timestamp_obs=timestamp_obs, pointing_type=pointing_type, exptime=exp_time) diff --git a/test/ObservationSim.cfg b/test/ObservationSim.cfg index ea2c2cef154274a0f95065d3362a12bd697cb43b..eaf3ca674a54f8e1d77df319298b5aa55926f380 100755 --- a/test/ObservationSim.cfg +++ b/test/ObservationSim.cfg @@ -7,6 +7,7 @@ # Observation setting date_obs 210525 # Observation date [yymmdd] time_obs 120000 # Observation time [hhmmss] +exp_time 150. # Exposure time [seconds] ra_center 60.0 # Telesscope pointing center [degree] (Default) dec_center -40.0 psf_rcont 0.15,0.8 # Radius of the 80% flux concentration (for a Gaussian PSF) @@ -18,7 +19,7 @@ survey_type Photometric # Survey simulation option: # 'Photometric': only run sims for photometric chips # 'Spectroscopic': only run sims for spectroscopic chps # 'All': run sims for all chips -psf_model Interp # which PSF model to use: +psf_model Interp # Which PSF model to use: # 'Gauss': simple gaussin profile # 'Interp': Interpolated PSF from measured data @@ -26,8 +27,9 @@ psf_model Interp # which PSF model to use: # Instrument effects #=========================================================== field_dist y # Y/N: Whether to add field distortion -abs_back y # Y: add sky background + sky noise & dark + dark noise - # N: only add sky noise and dark noise +add_back y # Y/N: Whether to add sky background +add_dark y # Y/N: Whether to add dark noise +add_readout y # Y/N: Whether to add read-out (Gaussian) noise add_bias y # Y/N: Whether to add bias-level to image shutter_effect y # Y/N: Whether to apply shutter effect shutter_output n # Y/N: Whether to export shutter effect 16-bit image (<=65535) @@ -56,6 +58,7 @@ NDark 1 # Number of frames to be exported NFlat 1 # Number of frames to be exported dark_exptime 300 # Exposure time for dark current frame flat_exptime 150 # Exposure time for flat field frame +readout_time 40 # Time for read-out by each channel #=========================================================== # Shear method @@ -85,6 +88,7 @@ seed_mock 12345678 # Seed for generating random values of ra, seed_star 121212121 # Seed for generating random redshift, sed_type for stars seed_gal 212121212 # Seed for generating random redshift, sed_type for galaxies seed_Av 121212 # Seed for generating random intrinsic extinction +seed_poisson 20210601 # Seed for Poisson noise seed_CR 20210317 # Seed for generating random cosmic ray map seed_flat 20210101 # Seed for generating random flat field seed_prnu 20210102 # Seed for photo-response non-uniformity @@ -93,4 +97,5 @@ seed_gainNonUniform 20210202 # Seed for Gain nonuniformity seed_biasNonUniform 20210203 # Seed for Bias nonuniformity seed_rnNonUniform 20210204 # Seed for ReadNoise nonuniformity seed_badcolumns 20240309 # Initial Seed for Bad Columns -seed_defective 20210304 # Seed for Defective(bad) pixels \ No newline at end of file +seed_defective 20210304 # Seed for Defective(bad) pixels +seed_readout 20210601 # Seed for Read-out gaussian noise \ No newline at end of file