'''Sources simulation module
Created on 14 févr. 2019
author: catalano camille , APC/IN2P3/CNRS
'''
import numpy as np
import scipy.integrate as integrate
from common import add_path_data_ref_eclairs
from common.io.fits_tools import read_fits
import common.sky.catalog as cat
from common.sky.cat_builder import catalog_with_channel
from simu.lib.InstruX import SimuECLAIRsMaskProjection
from gammapy.astro.source.tests.test_snr import snr
from common.instru.model_effect import ECLAIRsDetectorEffectDefault
[docs]class ModelSrcInterface(object):
"""
Base object for the creation of the shadowgram model of different sources
"""
idx_chan = 0
e_max = 0
e_min = 0
"""energy upper limit in keV"""
attitude = [0, 0, 0]
"""attitude [RA, DEC, ORI] of ECLAIRs in degrees"""
position = [0, 0, 0]
"""position [X, Y, Z] of SVOM in km in J2000"""
velocity = [0, 0, 0]
"""velocity [Vx, Vy, Vz] of SVOM in km/s"""
time = None
"""start time in UTC?? (TBD)"""
duration = 10
"""time duration in s"""
sim_geom = SimuECLAIRsMaskProjection()
"""instrument simulator for geometry and have coherent convention between source"""
mdl_effect = ECLAIRsDetectorEffectDefault()
def __init__(self):
self._snr = -1
self._name = "Unknown"
self._name_model = "Not defined"
self.always_noisy = False
[docs] def use_always_noisy_model(self):
self.always_noisy = True
[docs] def get_model(self, obs_time=None):
"""should return the shadowgram model
"""
raise
[docs] def get_noisy_model(self, obs_time=None):
"""Return the shadowgram model with poisson noise
"""
return np.random.poisson(self.get_model(obs_time))
[docs] @classmethod
def set_idx_channel(cls, idx_chan):
"""
"""
if idx_chan >= 0:
cls.idx_chan = idx_chan
cls.set_energy_band(cls.mdl_effect.chan_boundary[idx_chan],
cls.mdl_effect.chan_boundary[idx_chan+1])
[docs] @classmethod
def set_energy_band(cls, e_min, e_max):
"""
set the energy lower and upper limit for sources models classes
:param e_min: energy lower limit in keV
:type e_min: float
:param e_max: energy upper limit in keV
:type e_max: float
"""
cls.e_min = e_min
cls.e_max = e_max
[docs] @classmethod
def set_sim_dpix(cls, dpix):
"""
"""
cls.mdl_effect = dpix
[docs] @classmethod
def set_satellite_context(cls, attitude, position, velocity, time):
"""
set the satellite information: attitude, position, velocity and time
:param attitude: [ra, dec, ori] in degrees
:type attitude: [float, float, float]
:param position: [X, Y, Z] in km in J2000
:type position: [float, float, float]
:param velocity: [Vx, Vy, Vz] in km/s
:type velocity: [float, float, float]
:param time: PPS time in s from mjdref
:type time: float
"""
cls.attitude = attitude
cls.position = position
cls.velocity = velocity
cls.time = time
@property
def snr(self, idx_chan):
return self._snr[idx_chan]
@snr.setter
def snr(self, snr, idx_chan):
self._snr[idx_chan] = snr
@property
def name(self):
return self._name
[docs]class ModelFlat(ModelSrcInterface):
def __init__(self, level):
"""
**constructor**
:param level: intensity detected
:type level: float [ph/cm2/s]
"""
super().__init__()
self.pix_area = ModelSrcInterface.sim_geom.get_det_pixel_surface()
shape = ModelSrcInterface.sim_geom.inst.detect.get_shape()
self.detect_ones = np.ones(shape, dtype=np.float64)
self.noise_cm2_s = level
self._name_model = "Flat"
[docs] def get_model(self, obs_time=None):
"""return the shadowgram model
:param obs_time: observation time in s
:type obs_time: float
:return: noise detector image in ph/pix
:rtype: array(float)
"""
if not obs_time :
obs_time = ModelSrcInterface.duration
noise_pixel = self.noise_cm2_s * self.pix_area * obs_time
return self.detect_ones*noise_pixel
[docs]class ModelInternalNoise(ModelFlat):
"""
Model for a flat internal noise, impact on detector of particles in the terrestrial environment
The internal noise is independent from wavelength
"""
def __init__(self, noise_level=0.003):
"""
**constructor**
:param noise_level: internal noise level in ph/s/cm2/keV (default is 0.003)
:type noise_level: float
"""
super().__init__(noise_level)
self._name = "Internal noise"
self.interne_noise = noise_level
[docs] def get_model(self, obs_time=None):
"""return the shadowgram model
:param obs_time: observation time in s
:type obs_time: floatl
:return: noise detector image in ph/pix
:rtype: array(float)
"""
if not obs_time :
obs_time = ModelSrcInterface.duration
self.noise_cm2_s = (self.e_max - self.e_min) * self.interne_noise
return super().get_model(obs_time)
[docs]class ModelCxbFlat(ModelFlat):
"""model for a flat CXB without spectrum
"""
def __init__(self, cxb_level=1):
"""
**constructor**
:param cxb_level: intensity of CXB
:type cxb_level: float [ph/cm2/s]
"""
super().__init__(cxb_level)
self._name = "CXB"
self._name_model = "Flat"
[docs]class ModelCxbShapeBased(ModelSrcInterface):
"""model for a CXB without spectrum based on the pixel solid angle shape
"""
def __init__(self, cxb_level=3035, solidangle_filename=add_path_data_ref_eclairs('instru/pixels_solid_angle.fits')):
"""
**constructor**
:param cxb_level: intensity of CXB
:type cxb_level: float [ph/cm2/s]
:param solidangle_filename: PATH/filename of the solid angle image
:type solidangle_filename: string
"""
super().__init__()
self.cxb_level = cxb_level
self.solidangle_filename = solidangle_filename
solid_angle = read_fits(solidangle_filename, ugts_standard=1)
self.rate = cxb_level * solid_angle / solid_angle.sum()
self._name = "CXB"
self._name_model = "Shape"
[docs] def get_model(self, obs_time=None):
"""return the shadowgram model
:param obs_time: observation time in s
:type obs_time: float
:return: cxb detector image in ph/pix
:rtype: array(float)
"""
if not obs_time :
obs_time = ModelSrcInterface.duration
ret = self.rate*obs_time
if self.always_noisy:
ret = np.random.poisson(ret)
return ret
[docs]class ModelCxbShapeBasedEnergy(ModelSrcInterface):
"""model for a energy dependent CXB based on the pixel solid angle shape
"""
def __init__(self, solidangle_filename=add_path_data_ref_eclairs('instru/pixels_solid_angle.fits')):
"""
**constructor**
:param solidangle_filename: PATH/filename of the solid angle image
:type solidangle_filename: string
"""
super().__init__()
self.norm = 0.109
self.gamma_1 = 1.4
self.gamma_2 = 2.88
self.e_b = 29
self.solidangle_filename = solidangle_filename
self.solid_angle = read_fits(solidangle_filename, ugts_standard=1)
self._name = "energy CXB shape"
[docs] def get_model(self, obs_time=None):
"""return the shadowgram model
:param obs_time: observation time in s
:type obs_time: float
:return: cxb detector image in ph/pix
:rtype: array(float)
"""
if not obs_time :
obs_time = ModelSrcInterface.duration
# flux = integral between e_min and e_max of eq. 4 of https://arxiv.org/pdf/0811.1444.pdf
flux, err = integrate.quad(lambda e: self.norm / ((e/self.e_b)**self.gamma_1 + (e/self.e_b)**self.gamma_2), self.e_min, self.e_max)
ret = obs_time * self.solid_angle * self.sim_geom.get_det_pixel_surface() * flux
if self.always_noisy:
ret = np.random.poisson(ret)
return ret
[docs]class ModelPointSrcIRF(ModelSrcInterface):
"""model for point source
"""
S_cpt = 0
def __init__(self, info_src, verbose=False):
"""
**constructor**
info_src must have at least the same number of intensities as idx_chan
:param info_src: source information {'elev','dir','intensity','name'}
:type info_src: astropy.table.Row or dict
"""
self.info_src = info_src
if ModelPointSrcIRF.S_cpt == 0 and verbose:
print(self.info_src)
print(type(self.info_src))
print(self.info_src.dtype)
ModelPointSrcIRF.S_cpt += 1
self.simu = ModelSrcInterface.sim_geom
self._name = self.info_src["name"]
self._name_model = "with IRF"
self._shadow_cm2 = None
self._last_shadow = None
self.coef_hors_axis = self.set_coef_hors_axis(verbose)
[docs] def set_coef_hors_axis(self, verbose=False):
pix_y, pix_z = ModelSrcInterface.sim_geom.inst.elevdir_to_skypix(np.array(self.info_src['elev']),
np.array(self.info_src['dir']))
return ModelSrcInterface.mdl_effect.get_irf_val_ecllos(pix_y, pix_z, verbose)
[docs] def get_model(self, obs_time=None):
"""return the shadowgram model
:param obs_time: observation time in s
:type obs_time: float
:return: cxb detector image in ph/pix
:rtype: array(float)
"""
if not obs_time :
obs_time = ModelSrcInterface.duration
if self._shadow_cm2 is None:
self._shadow_cm2 = self.simu.get_shadow_surface(self.info_src['elev'], self.info_src['dir']).copy()
try:
intensity = self.info_src['intensity'][ModelSrcInterface.idx_chan]
except:
intensity = self.info_src['intensity']
self._last_shadow = intensity*obs_time*self._shadow_cm2 * self.coef_hors_axis
return self._last_shadow.copy()
[docs]class ModelPointSrcSinElev(ModelPointSrcIRF):
"""model for point source
"""
def __init__(self, info_src):
self.sin_elev = np.sin(np.deg2rad(info_src['elev']))
super().__init__(info_src)
self._name_model = "with sin(elev)"
[docs] def set_coef_hors_axis(self, verbose=False):
return self.sin_elev
[docs]class ModelGRBfromPhotonFile(ModelSrcInterface):
def __init__(self, file_grb):
# read list
# uniform random for each photon
# define ray tracing simulator
pass
[docs] def get_model(self, obs_time=None):
# read list
# uniform random for each photon
# define ray tracing simulator
#
pass
[docs]class ListModelPointSrcFromCatalog(object):
def __init__(self, o_dpix, use_irf=True):
"""
o dpix used to define dpix channel
"""
assert isinstance(o_dpix, ECLAIRsDetectorEffectDefault)
self._dpix = o_dpix
self.use_irf = use_irf
def _get_model_point_src(self, info_src, verbose):
if self.use_irf:
mdl = ModelPointSrcIRF(info_src, verbose)
else:
mdl = ModelPointSrcSinElev(info_src)
return mdl
[docs] def init_with_cat_swift_2012(self, attitude, verbose=False):
""" return a list of ShadowModelPointSrcSinElev from a catalog
attitude used to define fov
"""
self._attitude = attitude
path_cat = add_path_data_ref_eclairs("cat/BAT_70m_catalog_20nov2012.fits")
cat_bat = catalog_with_channel(path_cat, self._dpix, verbose)
assert isinstance(cat_bat, cat.CatalogAstroWithEnergySpecSampling)
# select only sources in fov with CatalogFovEclairs
self.cat_fov = cat.CatalogFovEclairs()
self.cat_fov.set_ptg_instru(self._attitude)
self.cat_fov.from_astro_catalog_with_spectrum(cat_bat)
d_src = {}
for isrc in range(self.cat_fov.get_nb_element()):
info_src = self.cat_fov._catalog[isrc]
if verbose:
print('add ',info_src["name"])
try:
mdl = self._get_model_point_src(info_src, verbose)
d_src[mdl.name] = mdl
except:
pass
if verbose:
print(f"Select {len(d_src)} sources in FOV with SWIFT/BAT catalog")
return d_src
[docs] def init_with_cat_fov(self, cat_fov, verbose=False):
"""
"""
assert isinstance(cat_fov, cat.CatalogFovEclairs)
#assert self._dpix.chan_nb() == cat_fov[0]["intensity"].size
d_src = {}
for isrc in range(cat_fov.get_nb_element()):
info_src = cat_fov._catalog[isrc]
if verbose:
print('add ',info_src["name"])
mdl = self._get_model_point_src(info_src, verbose)
d_src[mdl.name] = mdl
if verbose:
print(f"add {len(d_src)} sources in FOV with user catalog")
return d_src