diff --git a/Catalog/C6_50sqdeg.py b/Catalog/C6_50sqdeg.py new file mode 100644 index 0000000000000000000000000000000000000000..f8ded9cdbb9ce7709d0ea3e6e9e4d69d195fcf5c --- /dev/null +++ b/Catalog/C6_50sqdeg.py @@ -0,0 +1,462 @@ +import os +import galsim +import random +import numpy as np +import h5py as h5 +import healpy as hp +import astropy.constants as cons +import traceback +from astropy.coordinates import spherical_to_cartesian +from astropy.table import Table +from scipy import interpolate +from datetime import datetime + +from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar +from ObservationSim.MockObject._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist +from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position + +# (TEST) +from astropy.cosmology import FlatLambdaCDM +from astropy import constants +from astropy import units as U +from astropy.coordinates import SkyCoord +from astropy.io import fits + +try: + import importlib.resources as pkg_resources +except ImportError: + # Try backported to PY<37 'importlib_resources' + import importlib_resources as pkg_resources + +NSIDE = 128 + +bundle_file_list = ['galaxies_C6_bundle000199.h5','galaxies_C6_bundle000200.h5','galaxies_C6_bundle000241.h5','galaxies_C6_bundle000242.h5','galaxies_C6_bundle000287.h5','galaxies_C6_bundle000288.h5','galaxies_C6_bundle000714.h5','galaxies_C6_bundle000715.h5','galaxies_C6_bundle000778.h5','galaxies_C6_bundle000779.h5','galaxies_C6_bundle000842.h5','galaxies_C6_bundle000843.h5','galaxies_C6_bundle002046.h5','galaxies_C6_bundle002110.h5','galaxies_C6_bundle002111.h5','galaxies_C6_bundle002173.h5','galaxies_C6_bundle002174.h5','galaxies_C6_bundle002238.h5','galaxies_C6_bundle002596.h5','galaxies_C6_bundle002597.h5','galaxies_C6_bundle002656.h5','galaxies_C6_bundle002657.h5','galaxies_C6_bundle002711.h5','galaxies_C6_bundle002712.h5','galaxies_C6_bundle002844.h5','galaxies_C6_bundle002845.h5','galaxies_C6_bundle002884.h5','galaxies_C6_bundle002885.h5','galaxies_C6_bundle002921.h5','galaxies_C6_bundle002922.h5'] + +qsosed_file_list = ['quickspeclib_interp1d_run1.fits','quickspeclib_interp1d_run2.fits','quickspeclib_interp1d_run3.fits','quickspeclib_interp1d_run4.fits','quickspeclib_interp1d_run5.fits','quickspeclib_interp1d_run6.fits','quickspeclib_interp1d_run7.fits','quickspeclib_interp1d_run8.fits','quickspeclib_interp1d_run9.fits','quickspeclib_interp1d_run10.fits','quickspeclib_interp1d_run11.fits','quickspeclib_interp1d_run12.fits','quickspeclib_interp1d_run13.fits','quickspeclib_interp1d_run14.fits','quickspeclib_interp1d_run15.fits','quickspeclib_interp1d_run16.fits','quickspeclib_interp1d_run17.fits','quickspeclib_interp1d_run18.fits','quickspeclib_interp1d_run19.fits','quickspeclib_interp1d_run20.fits','quickspeclib_interp1d_run21.fits','quickspeclib_interp1d_run22.fits','quickspeclib_interp1d_run23.fits','quickspeclib_interp1d_run24.fits','quickspeclib_interp1d_run25.fits','quickspeclib_interp1d_run26.fits','quickspeclib_interp1d_run27.fits','quickspeclib_interp1d_run28.fits','quickspeclib_interp1d_run29.fits','quickspeclib_interp1d_run30.fits'] + +star_file_list = ['C7_Gaia_Galaxia_RA170DECm23_healpix.hdf5', 'C7_Gaia_Galaxia_RA180DECp60_healpix.hdf5', 'C7_Gaia_Galaxia_RA240DECp30_healpix.hdf5', 'C7_Gaia_Galaxia_RA300DECm60_healpix.hdf5', 'C7_Gaia_Galaxia_RA30DECm48_healpix.hdf5'] +star_center_list = [(170., -23.), (180., 60.), (240., 30.), (300., -60.), (30., -48.)] + +def get_bundleIndex(healpixID_ring, bundleOrder=4, healpixOrder=7): + assert NSIDE == 2**healpixOrder + shift = healpixOrder - bundleOrder + shift = 2*shift + + nside_bundle = 2**bundleOrder + nside_healpix= 2**healpixOrder + + healpixID_nest= hp.ring2nest(nside_healpix, healpixID_ring) + bundleID_nest = (healpixID_nest >> shift) + bundleID_ring = hp.nest2ring(nside_bundle, bundleID_nest) + + return bundleID_ring + +def get_agnsed_file(bundle_file_name): + return qsosed_file_list[bundle_file_list.index(bundle_file_name)] + +def get_star_cat(ra_pointing, dec_pointing): + pointing_c = SkyCoord(ra=ra_pointing*U.deg, dec=dec_pointing*U.deg) + max_dist = 10 + return_star_path = None + for star_file, center in zip(star_file_list, star_center_list): + center_c = SkyCoord(ra=center[0]*U.deg, dec=center[1]*U.deg) + dist = pointing_c.separation(center_c).to(U.deg).value + if dist < max_dist: + return_star_path = star_file + max_dist = dist + return return_star_path + +class Catalog(CatalogBase): + def __init__(self, config, chip, pointing, chip_output, filt, **kwargs): + super().__init__() + self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"]) + self.seed_Av = config["catalog_options"]["seed_Av"] + + self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111) + + self.chip_output = chip_output + self.filt = filt + self.logger = chip_output.logger + + with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path: + self.normF_star = Table.read(str(filter_path)) + + self.config = config + self.chip = chip + self.pointing = pointing + + self.max_size = 0. + + if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]: + # Get the cloest star catalog file + star_file_name = get_star_cat(ra_pointing=self.pointing.ra, dec_pointing=self.pointing.dec) + star_path = os.path.join(config["catalog_options"]["input_path"]["star_cat"], star_file_name) + star_SED_file = config["catalog_options"]["SED_templates_path"]["star_SED"] + self.star_path = os.path.join(self.cat_dir, star_path) + self.star_SED_path = os.path.join(config["data_dir"], star_SED_file) + self._load_SED_lib_star() + + if "galaxy_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["galaxy_cat"] and not config["catalog_options"]["star_only"]: + galaxy_dir = config["catalog_options"]["input_path"]["galaxy_cat"] + self.galaxy_path = os.path.join(self.cat_dir, galaxy_dir) + self.galaxy_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["galaxy_SED"]) + self._load_SED_lib_gals() + self.agn_seds = {} + + if "AGN_SED" in config["catalog_options"]["SED_templates_path"] and not config["catalog_options"]["star_only"]: + self.AGN_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["AGN_SED"]) + + if "rotateEll" in config["catalog_options"]: + self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.)) + else: + self.rotation = 0. + + # Update output .cat header with catalog specific output columns + self._add_output_columns_header() + + self._get_healpix_list() + self._load() + + def _add_output_columns_header(self): + self.add_hdr = " model_tag teff logg feh" + self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp " + + self.add_fmt = " %10s %8.4f %8.4f %8.4f" + self.add_fmt += " %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8.4f " + self.chip_output.update_ouptut_header(additional_column_names=self.add_hdr) + + def _get_healpix_list(self): + self.sky_coverage = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.2) + ra_min, ra_max, dec_min, dec_max = self.sky_coverage.xmin, self.sky_coverage.xmax, self.sky_coverage.ymin, self.sky_coverage.ymax + ra = np.deg2rad(np.array([ra_min, ra_max, ra_max, ra_min])) + dec = np.deg2rad(np.array([dec_max, dec_max, dec_min, dec_min])) + self.pix_list = hp.query_polygon( + NSIDE, + hp.ang2vec(np.radians(90.) - dec, ra), + inclusive=True + ) + if self.logger is not None: + msg = str(("HEALPix List: ", self.pix_list)) + self.logger.info(msg) + else: + print("HEALPix List: ", self.pix_list) + + def load_norm_filt(self, obj): + if obj.type == "star": + return self.normF_star + elif obj.type == "galaxy" or obj.type == "quasar": + # return self.normF_galaxy + return None + else: + return None + + def _load_SED_lib_star(self): + self.tempSED_star = h5.File(self.star_SED_path,'r') + + def _load_SED_lib_gals(self): + pcs = h5.File(os.path.join(self.galaxy_SED_path, "pcs.h5"), "r") + lamb = h5.File(os.path.join(self.galaxy_SED_path, "lamb.h5"), "r") + self.lamb_gal = lamb['lamb'][()] + self.pcs = pcs['pcs'][()] + + def _load_gals(self, gals, pix_id=None, cat_id=0, agnsed_file=""): + ngals = len(gals['ra']) + + # Apply astrometric modeling + ra_arr = gals['ra'][:] + dec_arr = gals['dec'][:] + if self.config["obs_setting"]["enable_astrometric_model"]: + ra_list = ra_arr.tolist() + dec_list = dec_arr.tolist() + pmra_list = np.zeros(ngals).tolist() + pmdec_list = np.zeros(ngals).tolist() + rv_list = np.zeros(ngals).tolist() + parallax_list = [1e-9] * ngals + dt = datetime.utcfromtimestamp(self.pointing.timestamp) + date_str = dt.date().isoformat() + time_str = dt.time().isoformat() + ra_arr, dec_arr = on_orbit_obs_position( + input_ra_list=ra_list, + input_dec_list=dec_list, + input_pmra_list=pmra_list, + input_pmdec_list=pmdec_list, + input_rv_list=rv_list, + input_parallax_list=parallax_list, + input_nstars=ngals, + input_x=self.pointing.sat_x, + input_y=self.pointing.sat_y, + input_z=self.pointing.sat_z, + input_vx=self.pointing.sat_vx, + input_vy=self.pointing.sat_vy, + input_vz=self.pointing.sat_vz, + input_epoch="J2000", + input_date_str=date_str, + input_time_str=time_str + ) + + for igals in range(ngals): + # # (TEST) + # if igals > 100: + # break + + param = self.initialize_param() + param['ra'] = ra_arr[igals] + param['dec'] = dec_arr[igals] + param['ra_orig'] = gals['ra'][igals] + param['dec_orig'] = gals['dec'][igals] + # param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals] + if self.filt.filter_type == 'NUV': + param['mag_use_normal'] = gals['mag_csst_nuv'][igals] + else: + param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals] + if self.filt.is_too_dim(mag=param['mag_use_normal'], margin=self.config["obs_setting"]["mag_lim_margin"]): + continue + + param['z'] = gals['redshift'][igals] + param['model_tag'] = 'None' + param['g1'] = gals['shear'][igals][0] + param['g2'] = gals['shear'][igals][1] + param['kappa'] = gals['kappa'][igals] + param['e1'] = gals['ellipticity_true'][igals][0] + param['e2'] = gals['ellipticity_true'][igals][1] + + # For shape calculation + + param['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2) + if param['ell_total'] > 0.9: + continue + param['e1_disk'] = param['e1'] + param['e2_disk'] = param['e2'] + param['e1_bulge'] = param['e1'] + param['e2_bulge'] = param['e2'] + + + param['delta_ra'] = 0 + param['delta_dec'] = 0 + + # Masses + param['bulgemass'] = gals['bulgemass'][igals] + param['diskmass'] = gals['diskmass'][igals] + + param['size'] = gals['size'][igals] + if param['size'] > self.max_size: + self.max_size = param['size'] + + # Sersic index + param['disk_sersic_idx'] = 1. + param['bulge_sersic_idx'] = 4. + + # Sizes + param['bfrac'] = param['bulgemass']/(param['bulgemass'] + param['diskmass']) + if param['bfrac'] >= 0.6: + param['hlr_bulge'] = param['size'] + param['hlr_disk'] = param['size'] * (1. - param['bfrac']) + else: + param['hlr_disk'] = param['size'] + param['hlr_bulge'] = param['size'] * param['bfrac'] + + # SED coefficients + param['coeff'] = gals['coeff'][igals] + param['detA'] = gals['detA'][igals] + + # Others + param['galType'] = gals['type'][igals] + param['veldisp'] = gals['veldisp'][igals] + + # TEST no redening and no extinction + param['av'] = 0.0 + param['redden'] = 0 + + # Is this an Quasar? + param['qsoindex'] = gals['qsoindex'][igals] + if param['qsoindex'] == -1: + param['star'] = 0 # Galaxy + param['agnsed_file'] = "" + else: + param['star'] = 2 # Quasar + param['agnsed_file'] = agnsed_file + + # NOTE: this cut cannot be put before the SED type has been assigned + if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): + continue + + # TEMP + self.ids += 1 + param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals) + + if param['star'] == 0: + obj = Galaxy(param, self.rotation, logger=self.logger) + elif param['star'] == 2: + obj = Quasar(param, logger=self.logger) + + # Need to deal with additional output columns + obj.additional_output_str = self.add_fmt%("n", 0., 0., 0., + param['bulgemass'], param['diskmass'], param['detA'], + param['e1'], param['e2'], param['kappa'], param['g1'], param['g2'], param['size'], + param['galType'], param['veldisp']) + + self.objs.append(obj) + + def _load_stars(self, stars, pix_id=None): + nstars = len(stars['sourceID']) + # Apply astrometric modeling + ra_arr = stars["RA"][:] + dec_arr = stars["Dec"][:] + pmra_arr = stars['pmra'][:] + pmdec_arr = stars['pmdec'][:] + rv_arr = stars['RV'][:] + parallax_arr = stars['parallax'][:] + if self.config["obs_setting"]["enable_astrometric_model"]: + ra_list = ra_arr.tolist() + dec_list = dec_arr.tolist() + pmra_list = pmra_arr.tolist() + pmdec_list = pmdec_arr.tolist() + rv_list = rv_arr.tolist() + parallax_list = parallax_arr.tolist() + dt = datetime.utcfromtimestamp(self.pointing.timestamp) + date_str = dt.date().isoformat() + time_str = dt.time().isoformat() + ra_arr, dec_arr = on_orbit_obs_position( + input_ra_list=ra_list, + input_dec_list=dec_list, + input_pmra_list=pmra_list, + input_pmdec_list=pmdec_list, + input_rv_list=rv_list, + input_parallax_list=parallax_list, + input_nstars=nstars, + input_x=self.pointing.sat_x, + input_y=self.pointing.sat_y, + input_z=self.pointing.sat_z, + input_vx=self.pointing.sat_vx, + input_vy=self.pointing.sat_vy, + input_vz=self.pointing.sat_vz, + input_epoch="J2000", + input_date_str=date_str, + input_time_str=time_str + ) + for istars in range(nstars): + # # (TEST) + # if istars > 100: + # break + + param = self.initialize_param() + param['ra'] = ra_arr[istars] + param['dec'] = dec_arr[istars] + param['ra_orig'] = stars["RA"][istars] + param['dec_orig'] = stars["Dec"][istars] + param['pmra'] = pmra_arr[istars] + param['pmdec'] = pmdec_arr[istars] + param['rv'] = rv_arr[istars] + param['parallax'] = parallax_arr[istars] + 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] + self.ids += 1 + param['id'] = stars['sourceID'][istars] + param['sed_type'] = stars['sourceID'][istars] + param['model_tag'] = stars['model_tag'][istars] + param['teff'] = stars['teff'][istars] + param['logg'] = stars['grav'][istars] + param['feh'] = stars['feh'][istars] + param['z'] = 0.0 + param['star'] = 1 # Star + obj = Star(param, logger=self.logger) + + # Append additional output columns to the .cat file + obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh'], + 0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.) + + self.objs.append(obj) + + def _load(self, **kwargs): + self.objs = [] + self.ids = 0 + + if "star_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["star_cat"] and not self.config["catalog_options"]["galaxy_only"]: + star_cat = h5.File(self.star_path, 'r')['catalog'] + for pix in self.pix_list: + try: + stars = star_cat[str(pix)] + self._load_stars(stars, pix_id=pix) + del stars + except Exception as e: + self.logger.error(str(e)) + print(e) + + if "galaxy_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["galaxy_cat"] and not self.config["catalog_options"]["star_only"]: + for pix in self.pix_list: + try: + bundleID = get_bundleIndex(pix) + bundle_file = "galaxies_C6_bundle{:06}.h5".format(bundleID) + file_path = os.path.join(self.galaxy_path, bundle_file) + gals_cat = h5.File(file_path, 'r')['galaxies'] + gals = gals_cat[str(pix)] + + # Get corresponding AGN SED file + agnsed_file = get_agnsed_file(bundle_file) + agnsed_path = os.path.join(self.AGN_SED_path, agnsed_file) + self.agn_seds[agnsed_file] = fits.open(agnsed_path)[0].data + + self._load_gals(gals, pix_id=pix, cat_id=bundleID, agnsed_file=agnsed_file) + + del gals + except Exception as e: + traceback.print_exc() + self.logger.error(str(e)) + print(e) + + if self.logger is not None: + self.logger.info("maximum galaxy size: %.4f"%(self.max_size)) + self.logger.info("number of objects in catalog: %d"%(len(self.objs))) + else: + print("number of objects in catalog: ", len(self.objs)) + + def load_sed(self, obj, **kwargs): + if obj.type == 'star': + _, wave, flux = tag_sed( + h5file=self.tempSED_star, + model_tag=obj.param['model_tag'], + teff=obj.param['teff'], + logg=obj.param['logg'], + feh=obj.param['feh'] + ) + elif obj.type == 'galaxy' or obj.type == 'quasar': + factor = 10**(-.4 * self.cosmo.distmod(obj.z).value) + if obj.type == 'galaxy': + flux = np.matmul(self.pcs, obj.coeff) * factor + # if np.any(flux < 0): + # raise ValueError("Glaxy %s: negative SED fluxes"%obj.id) + flux[flux < 0] = 0. + sedcat = np.vstack((self.lamb_gal, flux)).T + sed_data = getObservedSED( + sedCat=sedcat, + redshift=obj.z, + av=obj.param["av"], + redden=obj.param["redden"] + ) + wave, flux = sed_data[0], sed_data[1] + elif obj.type == 'quasar': + flux = self.agn_seds[obj.agnsed_file][int(obj.qsoindex)] * 1e-17 + flux[flux < 0] = 0. + wave = self.lamb_gal + else: + raise ValueError("Object type not known") + speci = interpolate.interp1d(wave, flux) + lamb = np.arange(2000, 11001+0.5, 0.5) + y = speci(lamb) + # erg/s/cm2/A --> photon/s/m2/A + all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 + sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX')) + + if obj.type == 'quasar': + # integrate to get the magnitudes + sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T + sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest') + sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False) + interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=self.filt.bandpass_full) + obj.param['mag_use_normal'] = getABMAG(interFlux, self.filt.bandpass_full) + # mag = getABMAG(interFlux, self.filt.bandpass_full) + # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal'])) + del wave + del flux + return sed diff --git a/ObservationSim/Instrument/Chip/Chip.py b/ObservationSim/Instrument/Chip/Chip.py index 9810f11c521c77a088e852891e481983f670beb9..cc0b4ce81623d1b551952b951eda084c43a6aca7 100755 --- a/ObservationSim/Instrument/Chip/Chip.py +++ b/ObservationSim/Instrument/Chip/Chip.py @@ -147,18 +147,18 @@ class Chip(FocalPlane): def _getSurveyType(self): if self.filter_type in ["GI", "GV", "GU"]: return "spectroscopic" - elif self.filter_type in ["nuv", "u", "g", 'r', 'i', 'z', 'y', 'FGS']: + elif self.filter_type in ["NUV", "u", "g", 'r', 'i', 'z', 'y', 'FGS']: 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 ['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' # Mirror efficiency: - # if filter_type == 'nuv': mirror_eff = 0.54 + # 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 @@ -185,7 +185,7 @@ class Chip(FocalPlane): 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", "FGS"] + filter_type_list = ["NUV","u", "g", "r", "i","z","y","GU", "GV", "GI", "FGS"] if filter_layout is not None: return filter_layout[chipID][0], filter_layout[chipID][1] if chipID == None: @@ -198,7 +198,7 @@ class Chip(FocalPlane): 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 [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" diff --git a/ObservationSim/Instrument/FilterParam.py b/ObservationSim/Instrument/FilterParam.py index 80057a52a157963efa260632cd77c48bd4a530e2..8f87b8fffd1f661b7c1eb2ef3eec352a5e24f94f 100755 --- a/ObservationSim/Instrument/FilterParam.py +++ b/ObservationSim/Instrument/FilterParam.py @@ -26,7 +26,7 @@ class FilterParam(object): # 8) dim end magnitude if filter_param == None: filtP = { - "nuv": [2867.7, 705.4, 2470.0, 3270.0, 0.1404, 0.004, 15.7, 25.4], + "NUV": [2867.7, 705.4, 2470.0, 3270.0, 0.1404, 0.004, 15.7, 25.4], "u": [3601.1, 852.1, 3120.0, 4090.0, 0.2176, 0.021, 16.1, 25.4], "g": [4754.5, 1569.8, 3900.0, 5620.0, 0.4640, 0.164, 17.2, 26.3], "r": [6199.8, 1481.2, 5370.0, 7030.0, 0.5040, 0.207, 17.0, 26.0], diff --git a/config/config_50sqdeg.yaml b/config/config_50sqdeg.yaml new file mode 100644 index 0000000000000000000000000000000000000000..1a2b8554a45f1331223fd315da2757d7b0392cf2 --- /dev/null +++ b/config/config_50sqdeg.yaml @@ -0,0 +1,221 @@ +--- +############################################### +# +# Configuration file for CSST simulation +# CSST-Sim Group, 2023/04/25 +# +############################################### + +# Base diretories and naming setup +# Can add some of the command-line arguments here as well; +# OK to pass either way or both, as long as they are consistent +work_dir: "/share/home/fangyuedong/csst-simulation/workplace/" +data_dir: "/share/simudata/CSSOSDataProductsSims/data/" +run_name: "50sqdeg_test" + +# Whether to use MPI +run_option: + use_mpi: NO + # NOTE: "n_threads" paramters is currently not used in the backend + # simulation codes. It should be implemented later in the web frontend + # in order to config the number of threads to request from NAOC cluster + n_threads: 80 + + # Output catalog only? + # If yes, no imaging simulation will run + out_cat_only: NO + +############################################### +# Catalog setting +############################################### +# Configure your catalog: options to be implemented +# in the corresponding (user defined) 'Catalog' class +catalog_options: + input_path: + cat_dir: "" + # star_cat: "C6_MMW_GGC_Astrometry_healpix.hdf5" + star_cat: "starcat/" + galaxy_cat: "qsocat/cat2CSSTSim_bundle-50sqDeg/" + # AGN_cat: "AGN_C6_ross13_rand_pos_rmax-1.3.fits" + + SED_templates_path: + star_SED: "Catalog_20210126/SpecLib.hdf5" + galaxy_SED: "Catalog_C6_20221212/sedlibs/" + AGN_SED: "qsocat/qsosed/" + # AGN_SED_WAVE: "wave_ross13.npy" + + # Only simulate stars? + star_only: NO + + # Only simulate galaxies? + galaxy_only: NO + + # rotate galaxy ellipticity + rotateEll: 0. # [degree] + + seed_Av: 121212 # Seed for generating random intrinsic extinction + +############################################### +# Observation setting +############################################### +obs_setting: + + # Options for survey types: + # "Photometric": simulate photometric chips only + # "Spectroscopic": simulate slitless spectroscopic chips only + # "FGS": simulate FGS chips only (31-42) + # "All": simulate full focal plane + survey_type: "Photometric" + + # Exposure time [seconds] + exp_time: 150. + + # Observation starting date & time + date_obs: "210525" # [yymmdd] + time_obs: "120000" # [hhmmss] + + # Default Pointing [degrees] + # Note: NOT valid when a pointing list file is specified + ra_center: 192.8595 + dec_center: 27.1283 + # Image rotation [degree] + image_rot: -113.4333 + + # (Optional) a file of point list + # if you just want to run default pointing: + # - pointing_dir: null + # - pointing_file: null + pointing_dir: "/share/home/fangyuedong/50sqdeg_pointings/" + pointing_file: "pointing_50_combined.dat" + + # Number of calibration pointings + np_cal: 0 + + # Run specific pointing(s): + # - give a list of indexes of pointings: [ip_1, ip_2...] + # - run all pointings: null + # Note: only valid when a pointing list is specified + run_pointings: [0] + + # Run specific chip(s): + # - give a list of indexes of chips: [ip_1, ip_2...] + # - run all chips: null + # Note: for all pointings + run_chips: [8] + + # Whether to enable astrometric modeling + enable_astrometric_model: True + + # Whether to enable straylight model + enable_straylight_model: True + + # Cut by saturation magnitude in which band? + cut_in_band: "z" + + # saturation magnitude margin + mag_sat_margin: -2.5 + + # limiting magnitude margin + mag_lim_margin: +1.0 + +############################################### +# PSF setting +############################################### +psf_setting: + + # Which PSF model to use: + # "Gauss": simple gaussian profile + # "Interp": Interpolated PSF from sampled ray-tracing data + psf_model: "Interp" + + # PSF size [arcseconds] + # radius of 80% energy encircled + # NOTE: only valid for "Gauss" PSF + psf_rcont: 0.15 + + # path to PSF data + # NOTE: only valid for "Interp" PSF + psf_dir: "/share/simudata/CSSOSDataProductsSims/data/psfCube1" + +############################################### +# Shear setting +############################################### + +shear_setting: + # Options to generate mock shear field: + # "constant": all galaxies are assigned a constant reduced shear + # "catalog": from catalog + shear_type: "catalog" + + # For constant shear filed + reduced_g1: 0. + reduced_g2: 0. + +############################################### +# Instrumental effects setting +############################################### +ins_effects: + # switches + # Note: bias_16channel, gain_16channel, and shutter_effect + # is currently not applicable to "FGS" observations + field_dist: ON # Whether to add field distortions + add_back: ON # Whether to add sky background + add_dark: ON # Whether to add dark noise + add_readout: ON # Whether to add read-out (Gaussian) noise + add_bias: ON # Whether to add bias-level to images + bias_16channel: OFF # Whether to add different biases for 16 channels + gain_16channel: OFF # Whether to make different gains for 16 channels + shutter_effect: OFF # Whether to add shutter effect + flat_fielding: OFF # Whether to add flat-fielding effect + prnu_effect: OFF # Whether to add PRNU effect + non_linear: OFF # Whether to add non-linearity + cosmic_ray: OFF # Whether to add cosmic-ray + cray_differ: OFF # Whether to generate different cosmic ray maps CAL and MS output + cte_trail: OFF # Whether to simulate CTE trails + saturbloom: OFF # Whether to simulate Saturation & Blooming + add_badcolumns: OFF # Whether to add bad columns + add_hotpixels: OFF # Whether to add hot pixels + add_deadpixels: OFF # Whether to add dead(dark) pixels + bright_fatter: OFF # Whether to simulate Brighter-Fatter (also diffusion) effect + + # Values: + # default values have been defined individually for each chip in: + # ObservationSim/Instrument/data/ccd/chip_definition.json + # Set them here will override the default values + # dark_exptime: 300 # Exposure time for dark current frames [seconds] + # flat_exptime: 150 # Exposure time for flat-fielding frames [seconds] + # readout_time: 40 # The read-out time for each channel [seconds] + # df_strength: 2.3 # Sillicon sensor diffusion strength + # bias_level: 500 # bias level [e-/pixel] + # gain: 1.1 # Gain + # full_well: 90000 # Full well depth [e-] + +############################################### +# Output options (for calibration pointings only) +############################################### +output_setting: + readout16: OFF # Whether to export as 16 channels (subimages) with pre- and over-scan + shutter_output: OFF # Whether to export shutter effect 16-bit image + bias_output: ON # Whether to export bias frames + dark_output: ON # Whether to export the combined dark current files + flat_output: ON # Whether to export the combined flat-fielding files + prnu_output: OFF # Whether to export the PRNU (pixel-to-pixel flat-fielding) files + NBias: 1 # Number of bias frames to be exported for each exposure + NDark: 1 # Number of dark frames to be exported for each exposure + NFlat: 1 # Number of flat frames to be exported for each exposure + +############################################### +# Random seeds +############################################### +random_seeds: + seed_poisson: 20210601 # Seed for Poisson noise + seed_CR: 20210317 # Seed for generating random cosmic ray maps + seed_flat: 20210101 # Seed for generating random flat fields + seed_prnu: 20210102 # Seed for photo-response non-uniformity + seed_gainNonUniform: 20210202 # Seed for gain nonuniformity + seed_biasNonUniform: 20210203 # Seed for bias nonuniformity + seed_rnNonUniform: 20210204 # Seed for readout-noise nonuniformity + seed_badcolumns: 20240309 # Seed for bad columns + seed_defective: 20210304 # Seed for defective (bad) pixels + seed_readout: 20210601 # Seed for read-out gaussian noise +... \ No newline at end of file diff --git a/profile_C6.sh b/profile_C6.sh index a64e6a3d1ea434bc68a3d5f78cf4a2d31baefb2d..49c79694a11d0082b595dfecf68e03b5ccca03a2 100755 --- a/profile_C6.sh +++ b/profile_C6.sh @@ -3,9 +3,14 @@ date python -m cProfile -o C6_profiler_test.pstats /share/home/fangyuedong/csst-simulation/run_sim.py \ - --config_file config_C6.yaml \ - --catalog C6_Catalog \ + --config_file config_50sqdeg.yaml \ + --catalog C6_50sqdeg \ -c /share/home/fangyuedong/csst-simulation/config + + # --config_file config_C6.yaml \ + # --catalog C6_Catalog \ + # -c /share/home/fangyuedong/csst-simulation/config + # --config_file test_fd_C6.yaml \ # --catalog fd_test_C6 \ # --config_file config_C6_test_wcs.yaml \