Skip to content
GitLab
Commit 6928c30c authored by JX's avatar JX 😵
Browse files

Merge remote-tracking branch 'origin/develop'

parents 749253d9 f5c20f83
Pipeline #4249 passed with stage
in 0 seconds
Show whitespace changes
Inline Side-by-side
.gitignore
View file @ 6928c30c
......@@ -6,3 +6,5 @@ dist/*
*disperse.c
*interp.c
!*libshao.so
*.out
pnodes
\ No newline at end of file
Catalog/C6_50sqdeg_ns.py
View file @ 6928c30c
......@@ -39,9 +39,9 @@ bundle_file_list = ['galaxies_C6_bundle000199.h5','galaxies_C6_bundle000200.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.)]
star_center_list = [(170., -23.), (180., 60.), (240., 30.), (300., -60.), (30., -48.),[246.5, 40]]
star_file_list = ['C9_RA170_DECm23_calmag_Nside_128_healpix.hdf5', 'C9_RA180_DECp60_calmag_Nside_128_healpix.hdf5', 'C9_RA240_DECp30_calmag_Nside_128_healpix.hdf5', 'C9_RA300_DECm60_calmag_Nside_128_healpix.hdf5', 'C9_RA30_DECm48_calmag_Nside_128_healpix.hdf5']
star_file_list = ['C9_RA170_DECm23_calmag_Nside_128_healpix.hdf5', 'C9_RA180_DECp60_calmag_Nside_128_healpix.hdf5', 'C9_RA240_DECp30_calmag_Nside_128_healpix.hdf5', 'C9_RA300_DECm60_calmag_Nside_128_healpix.hdf5', 'C9_RA30_DECm48_calmag_Nside_128_healpix.hdf5','trilegal_calMag_mpi_Nside_128_healpix.hdf5']
class StarParm(ctypes.Structure):
_fields_ = [
......@@ -187,7 +187,7 @@ class Catalog(CatalogBase):
rv_c = obj.param['rv']/(atcons.c.value/1000.)
Doppler_factor = np.sqrt((1+rv_c)/(1-rv_c))
wave_RV = wave*Doppler_factor
return wave_RV, spec
return wave_RV, np.power(10,spec[:])
def _load_SED_lib_gals(self):
pcs = h5.File(os.path.join(self.galaxy_SED_path, "pcs.h5"), "r")
......@@ -511,7 +511,6 @@ class Catalog(CatalogBase):
# 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
......
ObservationSim/Config/ChipOutput.py
View file @ 6928c30c
......@@ -3,8 +3,9 @@ import logging
import ObservationSim.Config._util as _util
from ObservationSim.Config.Header import generatePrimaryHeader
class ChipOutput(object):
def __init__(self, config, chip, filt, pointing):
def __init__(self, config, chip, filt, pointing, logger_filename=None):
self.config = config
self.chip = chip
self.filt = filt
......@@ -15,12 +16,12 @@ class ChipOutput(object):
self.h_prim = generatePrimaryHeader(
xlen=chip.npix_x,
ylen=chip.npix_y,
pointing_id = pointing.obs_id,
pointing_type_code = pointing.pointing_type_code,
pointing_id=pointing.obs_id,
pointing_type_code=pointing.pointing_type_code,
ra=pointing.ra,
dec=pointing.dec,
pixel_scale=chip.pix_scale,
time_pt = pointing.timestamp,
time_pt=pointing.timestamp,
exptime=pointing.exp_time,
im_type=pointing.pointing_type,
sat_pos=[pointing.sat_x, pointing.sat_y, pointing.sat_z],
......@@ -29,18 +30,22 @@ class ChipOutput(object):
run_counter=self.config["run_counter"],
chip_name=self.chip_label)
obs_id = _util.get_obs_id(img_type=self.pointing_type, project_cycle=config["project_cycle"], run_counter=config["run_counter"], pointing_id=pointing.obs_id, pointing_type_code = pointing.pointing_type_code)
obs_id = _util.get_obs_id(img_type=self.pointing_type, project_cycle=config["project_cycle"], run_counter=config[
"run_counter"], pointing_id=pointing.obs_id, pointing_type_code=pointing.pointing_type_code)
self.subdir = pointing.output_dir
self.cat_name = self.h_prim['FILENAME'] + '.cat'
if logger_filename is None:
logger_filename = self.h_prim['FILENAME'] + '.log'
self.logger = logging.getLogger()
fh = logging.FileHandler(os.path.join(self.subdir, logger_filename), mode='w+', encoding='utf-8')
fh = logging.FileHandler(os.path.join(
self.subdir, logger_filename), mode='w+', encoding='utf-8')
fh.setLevel(logging.DEBUG)
self.logger.setLevel(logging.DEBUG)
logging.getLogger('numba').setLevel(logging.WARNING)
formatter = logging.Formatter('%(asctime)s - %(msecs)d - %(levelname)-8s - [%(filename)s:%(lineno)d] - %(message)s')
formatter = logging.Formatter(
'%(asctime)s - %(msecs)d - %(levelname)-8s - [%(filename)s:%(lineno)d] - %(message)s')
fh.setFormatter(formatter)
self.logger.addHandler(fh)
......@@ -51,7 +56,7 @@ class ChipOutput(object):
self.hdr = hdr1 + hdr2
self.fmt = fmt1 + fmt2
self.logger.info("pointing_type = %s\n"%(self.pointing_type))
self.logger.info("pointing_type = %s\n" % (self.pointing_type))
def Log_info(self, message):
print(message)
......@@ -67,7 +72,8 @@ class ChipOutput(object):
def create_output_file(self):
if self.pointing_type == 'SCI':
self.cat = open(os.path.join(self.subdir, self.cat_name), "w")
self.logger.info("Creating catalog file %s ...\n"%(os.path.join(self.subdir, self.cat_name)))
self.logger.info("Creating catalog file %s ...\n" %
(os.path.join(self.subdir, self.cat_name)))
if not self.hdr.endswith("\n"):
self.hdr += "\n"
self.cat.write(self.hdr)
......@@ -76,8 +82,9 @@ class ChipOutput(object):
ximg = obj.real_pos.x + 1.0
yimg = obj.real_pos.y + 1.0
line = self.fmt%(
obj.id, int(self.chip_label), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.ra_orig, obj.dec_orig, obj.z, obj.getMagFilter(self.filt), obj.type,
line = self.fmt % (
obj.id, int(self.chip_label), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.ra_orig, obj.dec_orig, obj.z, obj.getMagFilter(
self.filt), obj.type,
obj.pmra, obj.pmdec, obj.rv, obj.parallax)
line += obj.additional_output_str
if not line.endswith("\n"):
......
ObservationSim/MockObject/CatalogBase.py
View file @ 6928c30c
......@@ -27,52 +27,52 @@ class CatalogBase(metaclass=ABCMeta):
@staticmethod
def initialize_param():
param = {
"star":-1,
"id":-1,
"ra":0,
"dec":0.,
"ra_orig":0.,
"dec_orig":0.,
"z":0.,
"sed_type":-1,
"model_tag":"unknown",
"mag_use_normal":100.,
"theta":0.,
"kappa":0.,
"g1":0.,
"g2":0.,
"bfrac":0.,
"av":0.,
"redden":0.,
"hlr_bulge":0.,
"hlr_disk":0.,
"ell_bulge":0.,
"ell_disk":0.,
"ell_tot":0.,
"e1_disk":0.,
"e2_disk":0.,
"e1_bulge":0.,
"e2_bulge":0.,
"teff":0.,
"logg":0.,
"feh":0.,
"DM":0.,
"stellarMass":1.,
"star": -1,
"id": -1,
"ra": 0,
"dec": 0.,
"ra_orig": 0.,
"dec_orig": 0.,
"z": 0.,
"sed_type": -1,
"model_tag": "unknown",
"mag_use_normal": 100.,
"theta": 0.,
"kappa": 0.,
"g1": 0.,
"g2": 0.,
"bfrac": 0.,
"av": 0.,
"redden": 0.,
"hlr_bulge": 0.,
"hlr_disk": 0.,
"ell_bulge": 0.,
"ell_disk": 0.,
"ell_tot": 0.,
"e1_disk": 0.,
"e2_disk": 0.,
"e1_bulge": 0.,
"e2_bulge": 0.,
"teff": 0.,
"logg": 0.,
"feh": 0.,
"DM": 0.,
"stellarMass": 1.,
# C6 galaxies parameters
"e1":0.,
"e2":0.,
"bulgemass":0.,
"diskmass":0.,
"size":0.,
"detA":0.,
"type":0,
"veldisp":0.,
"e1": 0.,
"e2": 0.,
"bulgemass": 0.,
"diskmass": 0.,
"size": 0.,
"detA": 0.,
"type": 0,
"veldisp": 0.,
"coeff": np.zeros(20),
# Astrometry related
"pmra":0.,
"pmdec":0.,
"rv":0.,
"parallax":1e-9
"pmra": 0.,
"pmdec": 0.,
"rv": 0.,
"parallax": 1e-9
}
return param
......@@ -87,27 +87,32 @@ class CatalogBase(metaclass=ABCMeta):
return e1, e2, e_total
@staticmethod
def convert_sed(mag, sed, target_filt, norm_filt=None):
def convert_sed(mag, sed, target_filt, norm_filt=None, mu=1.):
bandpass = target_filt.bandpass_full
if norm_filt is not None:
norm_thr_rang_ids = norm_filt['SENSITIVITY'] > 0.001
else:
norm_filt = Table(
np.array(np.array([bandpass.wave_list*10.0, bandpass.func(bandpass.wave_list)])).T, names=(['WAVELENGTH', 'SENSITIVITY'])
np.array(np.array([bandpass.wave_list*10.0, bandpass.func(
bandpass.wave_list)])).T, names=(['WAVELENGTH', 'SENSITIVITY'])
)
norm_thr_rang_ids = norm_filt['SENSITIVITY'] > 0.001
sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=mag,
spectrum=sed,
norm_thr=norm_filt,
sWave=np.floor(norm_filt[norm_thr_rang_ids][0][0]),
sWave=np.floor(
norm_filt[norm_thr_rang_ids][0][0]),
eWave=np.ceil(norm_filt[norm_thr_rang_ids][-1][0]))
sed_photon = copy.copy(sed)
sed_photon = np.array([sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T
sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest')
sed_photon = np.array(
[sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T
sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(
sed_photon[:, 1] * mu), interpolant='nearest')
# Get magnitude
sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False)
sed_photon = galsim.SED(sed_photon, wave_type='A',
flux_type='1', fast=False)
interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=bandpass)
mag_csst = getABMAG(
interFlux=interFlux,
......
ObservationSim/MockObject/Galaxy.py
View file @ 6928c30c
......@@ -8,6 +8,7 @@ from ObservationSim.MockObject.MockObject import MockObject
# import tracemalloc
class Galaxy(MockObject):
def __init__(self, param, logger=None):
super().__init__(param, logger=logger)
......@@ -16,6 +17,11 @@ class Galaxy(MockObject):
self.disk_sersic_idx = 1.
if not hasattr(self, "bulge_sersic_idx"):
self.bulge_sersic_idx = 4.
if not hasattr(self, "mu"):
if hasattr(self, "detA"):
self.mu = 1./self.detA
else:
self.mu = 1.
def unload_SED(self):
"""(Test) free up SED memory
......@@ -24,14 +30,16 @@ class Galaxy(MockObject):
def getGSObj_multiband(self, tel, psf_list, bandpass_list, filt, nphotons_tot=None, g1=0, g2=0, exptime=150., fd_shear=None):
if len(psf_list) != len(bandpass_list):
raise ValueError("!!!The number of PSF profiles and the number of bandpasses must be equal.")
raise ValueError(
"!!!The number of PSF profiles and the number of bandpasses must be equal.")
objs = []
if nphotons_tot == None:
nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
# print("nphotons_tot = ", nphotons_tot)
try:
full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
full = integrate_sed_bandpass(
sed=self.sed, bandpass=filt.bandpass_full)
except Exception as e:
print(e)
if self.logger:
......@@ -54,10 +62,12 @@ class Galaxy(MockObject):
return -1
psf = psf_list[i]
disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0)
disk = galsim.Sersic(n=self.disk_sersic_idx,
half_light_radius=self.hlr_disk, flux=1.0)
disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
disk = disk.shear(disk_shape)
bulge = galsim.Sersic(n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0)
bulge = galsim.Sersic(n=self.bulge_sersic_idx,
half_light_radius=self.hlr_bulge, flux=1.0)
bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
bulge = bulge.shear(bulge_shape)
......@@ -67,13 +77,15 @@ class Galaxy(MockObject):
gal = bulge
else:
gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk
gal = gal.withFlux(nphotons)
if fd_shear is not None:
g1 += fd_shear.g1
g2 += fd_shear.g2
gal_shear = galsim.Shear(g1=g1, g2=g2)
gal = gal.shear(gal_shear)
# Magnification
gal = gal.magnify(self.mu)
gal = galsim.Convolve(psf, gal)
gal = gal.withFlux(nphotons)
objs.append(gal)
final = galsim.Sum(objs)
......@@ -85,7 +97,8 @@ class Galaxy(MockObject):
# print("nphotons_tot = ", nphotons_tot)
try:
full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
full = integrate_sed_bandpass(
sed=self.sed, bandpass=filt.bandpass_full)
except Exception as e:
print(e)
if self.logger:
......@@ -121,10 +134,12 @@ class Galaxy(MockObject):
is_updated = 0
# Model the galaxy as disk + bulge
disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp)
disk = galsim.Sersic(
n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp)
disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
disk = disk.shear(disk_shape)
bulge = galsim.Sersic(n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp)
bulge = galsim.Sersic(
n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp)
bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
bulge = bulge.shear(bulge_shape)
......@@ -155,7 +170,8 @@ class Galaxy(MockObject):
# print("nphotons_sub-band_%d = %.2f"%(i, nphotons))
# Get PSF model
psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
psf, pos_shear = psf_model.get_PSF(
chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
if self.bfrac == 0:
gal_temp = disk
......@@ -164,10 +180,13 @@ class Galaxy(MockObject):
else:
gal_temp = self.bfrac * bulge + (1.0 - self.bfrac) * disk
gal_temp = gal_temp.shear(gal_shear)
gal_temp = gal_temp.withFlux(nphotons)
# Magnification
gal_temp = gal_temp.magnify(self.mu)
if not big_galaxy: # Not apply PSF for very big galaxy
gal_temp = galsim.Convolve(psf, gal_temp)
gal_temp = gal_temp.withFlux(nphotons)
if i == 0:
gal = gal_temp
else:
......@@ -184,7 +203,8 @@ class Galaxy(MockObject):
# ERROR happens
return 2, pos_shear
stamp.setCenter(x_nominal, y_nominal)
bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1)
bounds = stamp.bounds & galsim.BoundsI(
0, chip.npix_x - 1, 0, chip.npix_y - 1)
if bounds.area() > 0:
chip.img.setOrigin(0, 0)
chip.img[bounds] += stamp[bounds]
......@@ -194,9 +214,9 @@ class Galaxy(MockObject):
if is_updated == 0:
# Return code 0: object photons missed this detector
print("obj %s missed"%(self.id))
print("obj %s missed" % (self.id))
if self.logger:
self.logger.info("obj %s missed"%(self.id))
self.logger.info("obj %s missed" % (self.id))
return 0, pos_shear
# # [C6 TEST]
......@@ -209,7 +229,8 @@ class Galaxy(MockObject):
norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001
sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'], spectrum=self.sed,
norm_thr=normFilter,
sWave=np.floor(normFilter[norm_thr_rang_ids][0][0]),
sWave=np.floor(
normFilter[norm_thr_rang_ids][0][0]),
eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0]))
if sedNormFactor == 0:
return 2, None
......@@ -230,7 +251,6 @@ class Galaxy(MockObject):
chip_wcs_local = self.chip_wcs.local(self.real_pos)
big_galaxy = False
if self.hlr_disk > 3.0 or self.hlr_bulge > 3.0: # Very big galaxy
big_galaxy = True
......@@ -244,7 +264,8 @@ class Galaxy(MockObject):
flat_cube = chip.flat_cube
xOrderSigPlus = {'A':1.3909419820029296,'B':1.4760376591236062,'C':4.035447379743442,'D':5.5684364343742825,'E':16.260021029735388}
xOrderSigPlus = {'A': 1.3909419820029296, 'B': 1.4760376591236062,
'C': 4.035447379743442, 'D': 5.5684364343742825, 'E': 16.260021029735388}
grating_split_pos_chip = 0 + grating_split_pos
branges = np.zeros([len(bandpass_list), 2])
......@@ -267,10 +288,12 @@ class Galaxy(MockObject):
brange = branges[i]
# psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp)
disk = galsim.Sersic(
n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp)
disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
disk = disk.shear(disk_shape)
bulge = galsim.Sersic(n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp)
bulge = galsim.Sersic(
n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp)
bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
bulge = bulge.shear(bulge_shape)
......@@ -286,12 +309,13 @@ class Galaxy(MockObject):
# kfrac = np.random.random()*(1.0 - self.bfrac)
# gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots
gal = gal.withFlux(tel.pupil_area * exptime)
if fd_shear:
g1 += fd_shear.g1
g2 += fd_shear.g2
gal_shear = galsim.Shear(g1=g1, g2=g2)
gal = gal.shear(gal_shear)
gal = gal.magnify(self.mu)
gal = gal.withFlux(tel.pupil_area * exptime)
# gal = galsim.Convolve(psf, gal)
# if not big_galaxy: # Not apply PSF for very big galaxy
......@@ -299,23 +323,25 @@ class Galaxy(MockObject):
# # if fd_shear is not None:
# # gal = gal.shear(fd_shear)
starImg = gal.drawImage(wcs=chip_wcs_local, offset=offset,method = 'real_space')
starImg = gal.drawImage(
wcs=chip_wcs_local, offset=offset, method='real_space')
origin_star = [y_nominal - (starImg.center.y - starImg.ymin),
x_nominal - (starImg.center.x - starImg.xmin)]
starImg.setOrigin(0, 0)
gal_origin = [origin_star[0], origin_star[1]]
gal_end = [origin_star[0] + starImg.array.shape[0] - 1, origin_star[1] + starImg.array.shape[1] - 1]
gal_end = [origin_star[0] + starImg.array.shape[0] -
1, origin_star[1] + starImg.array.shape[1] - 1]
if gal_origin[1] < grating_split_pos_chip < gal_end[1]:
subSlitPos = int(grating_split_pos_chip - gal_origin[1] + 1)
## part img disperse
# part img disperse
subImg_p1 = starImg.array[:, 0:subSlitPos]
star_p1 = galsim.Image(subImg_p1)
star_p1.setOrigin(0, 0)
origin_p1 = origin_star
xcenter_p1 = min(x_nominal,grating_split_pos_chip-1) - 0
xcenter_p1 = min(x_nominal, grating_split_pos_chip-1) - 0
ycenter_p1 = y_nominal-0
sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1,
......@@ -330,9 +356,10 @@ class Galaxy(MockObject):
pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p1, chip=chip, pos_img_local=[xcenter_p1, ycenter_p1],
psf_model=psf_model, bandNo=i + 1,
grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img = pos_img)
local_wcs=chip_wcs_local, pos_img=pos_img)
subImg_p2 = starImg.array[:, subSlitPos+1:starImg.array.shape[1]]
subImg_p2 = starImg.array[:,
subSlitPos+1:starImg.array.shape[1]]
star_p2 = galsim.Image(subImg_p2)
star_p2.setOrigin(0, 0)
origin_p2 = [origin_star[0], grating_split_pos_chip]
......@@ -351,11 +378,11 @@ class Galaxy(MockObject):
pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p2, chip=chip, pos_img_local=[xcenter_p2, ycenter_p2],
psf_model=psf_model, bandNo=i + 1,
grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img = pos_img)
local_wcs=chip_wcs_local, pos_img=pos_img)
del sdp_p1
del sdp_p2
elif grating_split_pos_chip<=gal_origin[1]:
elif grating_split_pos_chip <= gal_origin[1]:
sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED,
......@@ -367,9 +394,9 @@ class Galaxy(MockObject):
pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
psf_model=psf_model, bandNo=i + 1,
grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img = pos_img)
local_wcs=chip_wcs_local, pos_img=pos_img)
del sdp
elif grating_split_pos_chip>=gal_end[1]:
elif grating_split_pos_chip >= gal_end[1]:
sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
ycenter=y_nominal - 0, origin=origin_star,
tar_spec=normalSED,
......@@ -381,7 +408,7 @@ class Galaxy(MockObject):
pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
psf_model=psf_model, bandNo=i + 1,
grating_split_pos=grating_split_pos,
local_wcs=chip_wcs_local, pos_img = pos_img)
local_wcs=chip_wcs_local, pos_img=pos_img)
del sdp
# print(self.y_nominal, starImg.center.y, starImg.ymin)
......@@ -391,11 +418,13 @@ class Galaxy(MockObject):
def getGSObj(self, psf, g1=0, g2=0, flux=None, filt=None, tel=None, exptime=150.):
if flux == None:
flux = self.getElectronFluxFilt(filt, tel, exptime)
disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0)
disk = galsim.Sersic(n=self.disk_sersic_idx,
half_light_radius=self.hlr_disk, flux=1.0)
disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
disk = disk.shear(disk_shape)
bulge = galsim.Sersic(n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0)
bulge = galsim.Sersic(n=self.bulge_sersic_idx,
half_light_radius=self.hlr_bulge, flux=1.0)
bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
bulge = bulge.shear(bulge_shape)
......@@ -407,5 +436,6 @@ class Galaxy(MockObject):
return final
def getObservedEll(self, g1=0, g2=0):
e1_obs, e2_obs, e_obs, theta = eObs(self.e1_total, self.e2_total, g1, g2)
e1_obs, e2_obs, e_obs, theta = eObs(
self.e1_total, self.e2_total, g1, g2)
return self.e1_total, self.e2_total, g1, g2, e1_obs, e2_obs
ObservationSim/MockObject/Quasar.py
View file @ 6928c30c
import galsim
import os, sys
import os
import sys
import numpy as np
import astropy.constants as cons
from astropy.table import Table
......@@ -8,9 +9,15 @@ from scipy import interpolate
from ObservationSim.MockObject.MockObject import MockObject
from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG
class Quasar(MockObject):
def __init__(self, param, logger=None):
super().__init__(param, logger=logger)
if not hasattr(self, "mu"):
if hasattr(self, "detA"):
self.mu = 1./self.detA
else:
self.mu = 1.
def load_SED(self, survey_type, sed_path=None, cosids=None, objtypes=None, sed_templates=None, normFilter=None, target_filt=None):
'''
......@@ -22,8 +29,9 @@ class Quasar(MockObject):
norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001
if sed_templates is None:
# Read SED data directly
itype = objtypes[cosids==self.sed_type][0]
sed_file = os.path.join(sed_path, itype + "_ID%s.sed"%(self.sed_type))
itype = objtypes[cosids == self.sed_type][0]
sed_file = os.path.join(
sed_path, itype + "_ID%s.sed" % (self.sed_type))
if not os.path.exists(sed_file):
raise ValueError("!!! No SED found.")
sed_data = Table.read(sed_file, format="ascii")
......@@ -40,20 +48,26 @@ class Quasar(MockObject):
wave, flux = sed_data[0], sed_data[1]
flux_photon = flux * (wave / (cons.h.value * cons.c.value)) * 1e-13
sed_photon = Table(np.array([wave, flux_photon]).T, names=('WAVELENGTH', 'FLUX'))
sed_photon = Table(
np.array([wave, flux_photon]).T, names=('WAVELENGTH', 'FLUX'))
# Get scaling factor for SED
sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'],
spectrum=sed_photon,
norm_thr=normFilter,
sWave=np.floor(normFilter[norm_thr_rang_ids][0][0]),
sWave=np.floor(
normFilter[norm_thr_rang_ids][0][0]),
eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0]))
sed_photon = np.array([sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T
sed_photon = np.array(
[sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T
# Convert to galsim.SED object
spec = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest')
self.sed = galsim.SED(spec, wave_type='A', flux_type='1', fast=False)
spec = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(
sed_photon[:, 1]), interpolant='nearest')
self.sed = galsim.SED(spec, wave_type='A',
flux_type='1', fast=False)
# Get magnitude
interFlux = integrate_sed_bandpass(sed=self.sed, bandpass=target_filt.bandpass_full)
self.param['mag_%s'%target_filt.filter_type] = getABMAG(
interFlux = integrate_sed_bandpass(
sed=self.sed, bandpass=target_filt.bandpass_full)
self.param['mag_%s' % target_filt.filter_type] = getABMAG(
interFlux=interFlux,
bandpass=target_filt.bandpass_full)
# print('mag_use_normal = ', self.param['mag_use_normal'])
......@@ -61,7 +75,8 @@ class Quasar(MockObject):
elif survey_type == "spectroscopic":
if sed_templates is None:
self.sedPhotons(sed_path=sed_path, cosids=cosids, objtypes=objtypes)
self.sedPhotons(sed_path=sed_path,
cosids=cosids, objtypes=objtypes)
else:
sed_data = sed_templates[self.sed_type]
sed_data = getObservedSED(
......@@ -74,8 +89,8 @@ class Quasar(MockObject):
y = speci(lamb)
# erg/s/cm2/A --> photo/s/m2/A
all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
self.sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
self.sed = Table(
np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
def unload_SED(self):
"""(Test) free up SED memory
......
ObservationSim/MockObject/Star.py
View file @ 6928c30c
import galsim
import os, sys
import os
import sys
import numpy as np
import astropy.constants as cons
from astropy.table import Table
......@@ -8,9 +9,12 @@ from scipy import interpolate
from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG, tag_sed
from ObservationSim.MockObject.MockObject import MockObject
class Star(MockObject):
def __init__(self, param, logger=None):
super().__init__(param, logger=logger)
if not hasattr(self, "mu"):
self.mu = 1.
def unload_SED(self):
"""(Test) free up SED memory
......@@ -28,13 +32,15 @@ class Star(MockObject):
def getGSObj_multiband(self, tel, psf_list, bandpass_list, filt, nphotons_tot=None, g1=0, g2=0, exptime=150.):
if len(psf_list) != len(bandpass_list):
raise ValueError("!!!The number of PSF profiles and the number of bandpasses must be equal.")
raise ValueError(
"!!!The number of PSF profiles and the number of bandpasses must be equal.")
objs = []
if nphotons_tot == None:
nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
try:
full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
full = integrate_sed_bandpass(
sed=self.sed, bandpass=filt.bandpass_full)
except Exception as e:
print(e)
self.logger.error(e)
......
ObservationSim/sim_steps/add_objects.py
View file @ 6928c30c
......@@ -10,6 +10,7 @@ from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp, PSFInterpSL
from astropy.time import Time
from datetime import datetime, timezone
def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
# Get exposure time
......@@ -20,9 +21,12 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
# Load catalogues
if catalog is None:
self.chip_output.Log_error("Catalog interface class must be specified for SCIE-OBS")
raise ValueError("Catalog interface class must be specified for SCIE-OBS")
cat = catalog(config=self.overall_config, chip=chip, pointing=pointing, chip_output=self.chip_output, filt=filt)
self.chip_output.Log_error(
"Catalog interface class must be specified for SCIE-OBS")
raise ValueError(
"Catalog interface class must be specified for SCIE-OBS")
cat = catalog(config=self.overall_config, chip=chip,
pointing=pointing, chip_output=self.chip_output, filt=filt)
# Prepare output file(s) for this chip
# [NOTE] Headers of output .cat file may be updated by Catalog instance
......@@ -31,12 +35,15 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
# Prepare the PSF model
if self.overall_config["psf_setting"]["psf_model"] == "Gauss":
psf_model = PSFGauss(chip=chip, psfRa=self.overall_config["psf_setting"]["psf_rcont"])
psf_model = PSFGauss(
chip=chip, psfRa=self.overall_config["psf_setting"]["psf_rcont"])
elif self.overall_config["psf_setting"]["psf_model"] == "Interp":
if chip.survey_type == "spectroscopic":
psf_model = PSFInterpSLS(chip, filt, PSF_data_prefix=self.overall_config["psf_setting"]["psf_sls_dir"])
psf_model = PSFInterpSLS(
chip, filt, PSF_data_prefix=self.overall_config["psf_setting"]["psf_sls_dir"])
else:
psf_model = PSFInterp(chip=chip, npsf=chip.n_psf_samples, PSF_data_file=self.overall_config["psf_setting"]["psf_pho_dir"])
psf_model = PSFInterp(chip=chip, npsf=chip.n_psf_samples,
PSF_data_file=self.overall_config["psf_setting"]["psf_pho_dir"])
else:
self.chip_output.Log_error("unrecognized PSF model type!!", flush=True)
......@@ -50,20 +57,22 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
# Get the filter which will be used for magnitude cut
for ifilt in range(len(self.all_filters)):
temp_filter = self.all_filters[ifilt]
temp_filter.update_limit_saturation_mags(exptime=pointing.get_full_depth_exptime(temp_filter.filter_type), chip=chip)
temp_filter.update_limit_saturation_mags(
exptime=pointing.get_full_depth_exptime(temp_filter.filter_type), chip=chip)
if temp_filter.filter_type.lower() == self.overall_config["obs_setting"]["cut_in_band"].lower():
cut_filter = temp_filter
# Read in shear values from configuration file if the constant shear type is used
if self.overall_config["shear_setting"]["shear_type"] == "constant":
g1_field, g2_field, _ = get_shear_field(config=self.overall_config)
self.chip_output.Log_info("Use constant shear: g1=%.5f, g2=%.5f"%(g1_field, g2_field))
self.chip_output.Log_info(
"Use constant shear: g1=%.5f, g2=%.5f" % (g1_field, g2_field))
# Get chip WCS
if not hasattr(self, 'h_ext'):
_, _ = self.prepare_headers(chip=chip, pointing=pointing)
chip_wcs = galsim.FitsWCS(header = self.h_ext)
chip_wcs = galsim.FitsWCS(header=self.h_ext)
# Loop over objects
nobj = len(cat.objs)
......@@ -80,17 +89,19 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
try:
sed_data = cat.load_sed(obj)
norm_filt = cat.load_norm_filt(obj)
obj.sed, obj.param["mag_%s"%filt.filter_type.lower()], obj.param["flux_%s"%filt.filter_type.lower()] = cat.convert_sed(
obj.sed, obj.param["mag_%s" % filt.filter_type.lower()], obj.param["flux_%s" % filt.filter_type.lower()] = cat.convert_sed(
mag=obj.param["mag_use_normal"],
sed=sed_data,
target_filt=filt,
norm_filt=norm_filt,
mu=obj.mu
)
_, obj.param["mag_%s"%cut_filter.filter_type.lower()], obj.param["flux_%s"%cut_filter.filter_type.lower()] = cat.convert_sed(
_, obj.param["mag_%s" % cut_filter.filter_type.lower()], obj.param["flux_%s" % cut_filter.filter_type.lower()] = cat.convert_sed(
mag=obj.param["mag_use_normal"],
sed=sed_data,
target_filt=cut_filter,
norm_filt=norm_filt,
mu=obj.mu
)
except Exception as e:
traceback.print_exc()
......@@ -102,16 +113,19 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
# Exclude very bright/dim objects (for now)
if cut_filter.is_too_bright(
mag=obj.param["mag_%s"%self.overall_config["obs_setting"]["cut_in_band"].lower()],
mag=obj.param["mag_%s" %
self.overall_config["obs_setting"]["cut_in_band"].lower()],
margin=self.overall_config["obs_setting"]["mag_sat_margin"]):
self.chip_output.Log_info("obj %s too birght!! mag_%s = %.3f"%(obj.id, cut_filter.filter_type, obj.param["mag_%s"%self.overall_config["obs_setting"]["cut_in_band"].lower()]))
self.chip_output.Log_info("obj %s too birght!! mag_%s = %.3f" % (
obj.id, cut_filter.filter_type, obj.param["mag_%s" % self.overall_config["obs_setting"]["cut_in_band"].lower()]))
bright_obj += 1
obj.unload_SED()
continue
if filt.is_too_dim(
mag=obj.getMagFilter(filt),
margin=self.overall_config["obs_setting"]["mag_lim_margin"]):
self.chip_output.Log_info("obj %s too dim!! mag_%s = %.3f"%(obj.id, filt.filter_type, obj.getMagFilter(filt)))
self.chip_output.Log_info("obj %s too dim!! mag_%s = %.3f" % (
obj.id, filt.filter_type, obj.getMagFilter(filt)))
dim_obj += 1
obj.unload_SED()
continue
......@@ -130,14 +144,16 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
raise ValueError("Unknown shear input")
# Get position of object on the focal plane
pos_img, _, _, _, fd_shear = obj.getPosImg_Offset_WCS(img=chip.img, fdmodel=fd_model, chip=chip, verbose=False, chip_wcs=chip_wcs, img_header=self.h_ext)
pos_img, _, _, _, fd_shear = obj.getPosImg_Offset_WCS(
img=chip.img, fdmodel=fd_model, chip=chip, verbose=False, chip_wcs=chip_wcs, img_header=self.h_ext)
# [TODO] For now, only consider objects which their centers (after field distortion) are projected within the focal plane
# Otherwise they will be considered missed objects
# if pos_img.x == -1 or pos_img.y == -1 or (not chip.isContainObj(x_image=pos_img.x, y_image=pos_img.y, margin=0.)):
if pos_img.x == -1 or pos_img.y == -1:
self.chip_output.Log_info('obj_ra = %.6f, obj_dec = %.6f, obj_ra_orig = %.6f, obj_dec_orig = %.6f'%(obj.ra, obj.dec, obj.ra_orig, obj.dec_orig))
self.chip_output.Log_error("Objected missed: %s"%(obj.id))
self.chip_output.Log_info('obj_ra = %.6f, obj_dec = %.6f, obj_ra_orig = %.6f, obj_dec_orig = %.6f' % (
obj.ra, obj.dec, obj.ra_orig, obj.dec_orig))
self.chip_output.Log_error("Objected missed: %s" % (obj.id))
missed_obj += 1
obj.unload_SED()
continue
......@@ -146,7 +162,8 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
try:
if self.overall_config["run_option"]["out_cat_only"]:
isUpdated = True
obj.real_pos = obj.getRealPos(chip.img, global_x=obj.posImg.x, global_y=obj.posImg.y, img_real_wcs=obj.chip_wcs)
obj.real_pos = obj.getRealPos(
chip.img, global_x=obj.posImg.x, global_y=obj.posImg.y, img_real_wcs=obj.chip_wcs)
pos_shear = 0.
elif chip.survey_type == "photometric" and not self.overall_config["run_option"]["out_cat_only"]:
isUpdated, pos_shear = obj.drawObj_multiband(
......@@ -181,9 +198,10 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
pass
elif isUpdated == 0:
missed_obj += 1
self.chip_output.Log_error("Objected missed: %s"%(obj.id))
self.chip_output.Log_error("Objected missed: %s" % (obj.id))
else:
self.chip_output.Log_error("Draw error, object omitted: %s"%(obj.id))
self.chip_output.Log_error(
"Draw error, object omitted: %s" % (obj.id))
continue
except Exception as e:
traceback.print_exc()
......@@ -196,38 +214,47 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
del psf_model
gc.collect()
self.chip_output.Log_info("Running checkpoint #1 (Object rendering finished): pointing-%d chip-%d pid-%d memory-%6.2fGB"%(pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) ))
self.chip_output.Log_info("Running checkpoint #1 (Object rendering finished): pointing-%d chip-%d pid-%d memory-%6.2fGB" %
(pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024)))
self.chip_output.Log_info("# objects that are too bright %d out of %d"%(bright_obj, nobj))
self.chip_output.Log_info("# objects that are too dim %d out of %d"%(dim_obj, nobj))
self.chip_output.Log_info("# objects that are missed %d out of %d"%(missed_obj, nobj))
self.chip_output.Log_info(
"# objects that are too bright %d out of %d" % (bright_obj, nobj))
self.chip_output.Log_info(
"# objects that are too dim %d out of %d" % (dim_obj, nobj))
self.chip_output.Log_info(
"# objects that are missed %d out of %d" % (missed_obj, nobj))
# Apply flat fielding (with shutter effects)
flat_normal = np.ones_like(chip.img.array)
if obs_param["flat_fielding"] == True:
flat_normal = flat_normal * chip.flat_img.array / np.mean(chip.flat_img.array)
flat_normal = flat_normal * chip.flat_img.array / \
np.mean(chip.flat_img.array)
if obs_param["shutter_effect"] == True:
flat_normal = flat_normal * chip.shutter_img
flat_normal = np.array(flat_normal, dtype='float32')
self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT'], values = [True])
self.updateHeaderInfo(header_flag='ext', keys=[
'SHTSTAT'], values=[True])
else:
self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT','SHTOPEN1','SHTCLOS0'], values = [True,self.h_ext['SHTCLOS1'],self.h_ext['SHTOPEN0']])
self.updateHeaderInfo(header_flag='ext', keys=['SHTSTAT', 'SHTOPEN1', 'SHTCLOS0'], values=[
True, self.h_ext['SHTCLOS1'], self.h_ext['SHTOPEN0']])
chip.img *= flat_normal
del flat_normal
# renew header info
datetime_obs = datetime.utcfromtimestamp(pointing.timestamp)
datetime_obs = datetime_obs.replace(tzinfo=timezone.utc)
t_obs = Time(datetime_obs)
##ccd刷新2s,等待0.s，开始曝光
# ccd刷新2s,等待0.s，开始曝光
t_obs_renew = Time(t_obs.mjd - (2.+0.) / 86400., format="mjd")
t_obs_utc = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_obs_renew.unix).replace(tzinfo=timezone.utc).timestamp(), 1))
self.updateHeaderInfo(header_flag='prim', keys = ['DATE-OBS'], values = [t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]])
t_obs_utc = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(
t_obs_renew.unix).replace(tzinfo=timezone.utc).timestamp(), 1))
self.updateHeaderInfo(header_flag='prim', keys=[
'DATE-OBS'], values=[t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]])
#dark time : 曝光时间+刷新后等带时间0.s+关快门后读出前等待0.s
self.updateHeaderInfo(header_flag='ext', keys = ['DARKTIME'], values = [0.+0.+pointing.exp_time])
# dark time : 曝光时间+刷新后等带时间0.s+关快门后读出前等待0.s
self.updateHeaderInfo(header_flag='ext', keys=[
'DARKTIME'], values=[0.+0.+pointing.exp_time])
return chip, filt, tel, pointing
ObservationSim/sim_steps/add_sky_background.py
View file @ 6928c30c
......@@ -136,7 +136,7 @@ def add_sky_background(self, chip, filt, tel, pointing, catalog, obs_param):
self.updateHeaderInfo(header_flag='prim', keys = ['DATE-OBS'], values = [t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]])
#dark time : 曝光时间+刷新后等带时间0.5s+关闭快门时间1.5s+管快门后读出前等待0.5s
self.updateHeaderInfo(header_flag='ext', keys = ['DARKTIME'], values = [0.5+1.5+0.5+pointing.exp_time])
self.updateHeaderInfo(header_flag='ext', keys = ['DARKTIME'], values = [0.+0.0+0.0+pointing.exp_time])
return chip, filt, tel, pointing
\ No newline at end of file
setup.py
View file @ 6928c30c
......@@ -76,7 +76,7 @@ with open("requirements.txt", "r") as f:
]
setup(name='csst_msc_sim',
version='3.0.0',
version='3.0.0rc',
packages=find_packages(),
# install_requires=[
# # 'numpy>=1.18.5',
......
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