Commit 3e9d0f8e authored by xin's avatar xin
Browse files

sls direct image

parent f93152b4
......@@ -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
......
......@@ -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:
......
......@@ -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])
......
......@@ -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(
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment