diff --git a/Catalog/NGPCatalog.py b/Catalog/NGPCatalog.py index c854440228adbc5112f7d62ed6a0901fef79718b..6611501c2b2165b9a49bcea19c11c074097b06e6 100644 --- a/Catalog/NGPCatalog.py +++ b/Catalog/NGPCatalog.py @@ -136,7 +136,7 @@ class NGPCatalog(CatalogBase): if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): continue param['mag_use_normal'] = gals['mag_true_g_lsst'][igals] - if param['mag_use_normal'] >= 26.5: + if param['mag_use_normal'] >= 24: continue param['z'] = gals['redshift_true'][igals] param['model_tag'] = 'None' @@ -234,7 +234,7 @@ class NGPCatalog(CatalogBase): if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): continue param['mag_use_normal'] = stars['app_sdss_g'][istars] - if param['mag_use_normal'] >= 26.5: + if param['mag_use_normal'] >= 26: continue self.ids += 1 # param['id'] = self.ids diff --git a/ObservationSim/Instrument/Chip/Chip.py b/ObservationSim/Instrument/Chip/Chip.py index f9fa4f6afa076e52a2f45a42af8775425b09bb34..c06904bdd88521aecc5e9286da6376513dd00758 100755 --- a/ObservationSim/Instrument/Chip/Chip.py +++ b/ObservationSim/Instrument/Chip/Chip.py @@ -17,7 +17,7 @@ except ImportError: 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): + def __init__(self, chipID, ccdEffCurve_dir=None, CRdata_dir=None, sls_dir=None, config=None, treering_func=None, logger=None, direct_Img_sls=False): # Get focal plane (instance of paraent class) info # TODO: use chipID to config individual chip? super().__init__() @@ -35,6 +35,7 @@ class Chip(FocalPlane): self.readout_time = float(config["ins_effects"]['readout_time']) self.logger = logger + self.direct_Img_sls = direct_Img_sls # A chip ID must be assigned self.chipID = int(chipID) @@ -144,9 +145,14 @@ class Chip(FocalPlane): 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 self.direct_Img_sls: + if chipID in [1, 10, 21, 30]: filter_type = "i" + if chipID in [2, 5, 26, 29]: filter_type = "g" + if chipID in [3, 4, 27, 28]: filter_type = "u" + else: + 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 @@ -297,13 +303,46 @@ class Chip(FocalPlane): hdu1.writeto(fname, output_verify='ignore', overwrite=True) def addEffects(self, config, img, chip_output, filt, ra_cen, dec_cen, img_rot, exptime=150., pointing_ID=0, timestamp_obs=1621915200, pointing_type='MS', sky_map=None, tel=None, logger=None): - SeedGainNonuni=int(config["random_seeds"]["seed_gainNonUniform"]) - SeedBiasNonuni=int(config["random_seeds"]["seed_biasNonUniform"]) + + SeedGainNonuni = int(config["random_seeds"]["seed_gainNonUniform"]) + SeedBiasNonuni = int(config["random_seeds"]["seed_biasNonUniform"]) SeedRnNonuni = int(config["random_seeds"]["seed_rnNonUniform"]) SeedBadColumns = int(config["random_seeds"]["seed_badcolumns"]) SeedDefective = int(config["random_seeds"]["seed_defective"]) SeedCosmicRay = int(config["random_seeds"]["seed_CR"]) fullwell = int(config["ins_effects"]["full_well"]) + + if config["output_setting"]["output_SCI"]: + 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["ins_effects"]["add_back"] == True: + img += sky_map + # del sky_map + + # Add Poisson noise + seed = int(config["random_seeds"]["seed_poisson"]) + pointing_ID * 30 + self.chipID + rng_poisson = galsim.BaseDeviate(seed) + poisson_noise = galsim.PoissonNoise(rng_poisson, sky_level=0.) + img.addNoise(poisson_noise) + + dark_noise = galsim.DeviateNoise( + galsim.PoissonDeviate(rng_poisson, self.dark_noise * (self.exptime + 0.5 * self.readout_time))) + img.addNoise(dark_noise) + + 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) + + img.array[img.array > fullwell] = fullwell + img = img - sky_map - self.dark_noise * (self.exptime + 0.5 * self.readout_time) + return img + + if config["ins_effects"]["add_hotpixels"] == True: BoolHotPix = True else: diff --git a/ObservationSim/ObservationSim.py b/ObservationSim/ObservationSim.py index dd7739851b2ecce01293f515d267efc87aa76255..7825799618e2bf18ee1bfb94588201abe6a63384 100755 --- a/ObservationSim/ObservationSim.py +++ b/ObservationSim/ObservationSim.py @@ -17,7 +17,7 @@ from ObservationSim._util import get_shear_field, makeSubDir_PointingList from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position class Observation(object): - def __init__(self, config, Catalog, work_dir=None, data_dir=None): + def __init__(self, config, Catalog, work_dir=None, data_dir=None, direct_Img_sls=False): self.path_dict = config_dir(config=config, work_dir=work_dir, data_dir=data_dir) self.config = config self.tel = Telescope() @@ -43,7 +43,7 @@ class Observation(object): # Make Chip & Filter lists chip = Chip( chipID=chipID, - config=self.config) + config=self.config, direct_Img_sls=direct_Img_sls) filter_id, filter_type = chip.getChipFilter() filt = Filter(filter_id=filter_id, filter_type=filter_type, @@ -339,7 +339,11 @@ class Observation(object): row_num=chip.rowID, col_num=chip.colID, extName='raw') - chip.img = galsim.Image(chip.img.array, dtype=np.uint16) + out_type = np.uint16 + if self.config["output_setting"]["output_SCI"]: + out_type = np.float32 + + chip.img = galsim.Image(chip.img.array, dtype=out_type) hdu1 = fits.PrimaryHDU(header=h_prim) hdu2 = fits.ImageHDU(chip.img.array, header=h_ext) hdu1 = fits.HDUList([hdu1, hdu2]) diff --git a/run_sim.py b/run_sim.py index ee4ee57d887a4b4f12e3fa71089a544bf7f78849..6d2a0754ee8e4a8ce65aae7c874803801a61ca42 100755 --- a/run_sim.py +++ b/run_sim.py @@ -52,7 +52,7 @@ def run_sim(Catalog): shutil.copy(args.config_file, run_dir) # Initialize the simulation - obs = Observation(config=config, Catalog=Catalog, work_dir=config['work_dir'], data_dir=config['data_dir']) + obs = Observation(config=config, Catalog=Catalog, work_dir=config['work_dir'], data_dir=config['data_dir'], direct_Img_sls=config["obs_setting"]['direct_img_sls']) # Run simulation obs.runExposure_MPI_PointingList(