'''Sources simulation module
Created on 14 févr. 2019
author: catalano camille , APC/IN2P3/CNRS
'''
import logging
import numpy as np
from scipy.integrate import quad
from astropy.table import Row
import ecpi.common.sky.catalog as cat
from ecpi.common import add_path_data_ref_eclairs
from ecpi.common.io.fits_tools import read_fits
from ecpi.common.sky.cat_builder import catalog_with_channel
from ecpi.simu.lib.instru_x import SimuECLAIRsMaskProjection
from ecpi.common.instru.model_effect import ECLAIRsDetectorEffect, \
ECLAIRsDetectorEffectDefault
[docs]class ModelSrcInterface(object):
"""
Base object for the creation of the shadowgram model of different sources
"""
static_idx_chan = 0
static_e_max = 0
static_e_min = 0
"""energy upper limit in keV"""
static_attitude = [0, 0, 0]
"""attitude [RA, DEC, ORI] of ECLAIRs in degrees"""
static_position = [0, 0, 0]
"""position [X, Y, Z] of SVOM in km in J2000"""
static_velocity = [0, 0, 0]
"""velocity [Vx, Vy, Vz] of SVOM in km/s"""
static_time = None
"""start time in UTC?? (TBD)"""
static_duration = 10
"""time duration in s"""
static_sim_geom = None
"""instrument simulator for geometry and have coherent convention between source"""
static_mdl_effect = None
def __init__(self):
"""**Constructor**
"""
self._snr = -1
self._name = "Unknown"
self._name_model = "Not defined"
self._init_static_object()
self.logger = logging.getLogger(__name__)
[docs] def get_noisy_model(self, obs_time=None):
"""Return the shadowgram model with poisson noise.
:param obs_time: observation time in seconds
:type obs_time: float
"""
return np.random.poisson(self.get_model(obs_time))
def _init_static_object(self):
"""instrument simulator for geometry and have coherent convention between source
"""
if ModelSrcInterface.static_sim_geom is None:
ModelSrcInterface.static_sim_geom = SimuECLAIRsMaskProjection()
if ModelSrcInterface.static_mdl_effect is None:
ModelSrcInterface.static_mdl_effect = ECLAIRsDetectorEffectDefault()
[docs] @staticmethod
def set_idx_channel(idx_chan):
"""Set channel index and associated energy band.
:param idx_chan: channel index
:type idx_chan: int
"""
if idx_chan >= 0:
ModelSrcInterface.static_idx_chan = idx_chan
ModelSrcInterface.set_energy_band(
ModelSrcInterface.static_mdl_effect.chan_boundary[idx_chan],
ModelSrcInterface.static_mdl_effect.chan_boundary[idx_chan + 1])
[docs] @staticmethod
def set_energy_band(e_min, e_max):
"""Set the energy lower and upper limit for sources models classes
.. warning:: must be renamed as _set_energy_band
:param e_min: energy lower limit in keV
:type e_min: float
:param e_max: energy upper limit in keV
:type e_max: float
"""
ModelSrcInterface.static_e_min = e_min
ModelSrcInterface.static_e_max = e_max
[docs] @staticmethod
def set_sim_dpix(dpix):
ModelSrcInterface.static_mdl_effect = dpix
[docs] @staticmethod
def set_satellite_context(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
"""
ModelSrcInterface.static_attitude = attitude
ModelSrcInterface.static_position = position
ModelSrcInterface.static_velocity = velocity
ModelSrcInterface.static_time = time
@property
def name(self):
"""Get name.
"""
return self._name
[docs]class ModelFlat(ModelSrcInterface):
"""Generic class for flat model
"""
def __init__(self, level):
"""**constructor**
:param level: intensity detected
:type level: float [ph/cm2/s]
"""
super().__init__()
self.pix_area = ModelSrcInterface.static_sim_geom.get_det_pixel_surface()
shape = ModelSrcInterface.static_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.static_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: float
:return: noise detector image in ph/pix
:rtype: array(float)
"""
if not obs_time:
obs_time = ModelSrcInterface.static_duration
size_band = ModelSrcInterface.static_e_max - ModelSrcInterface.static_e_min
self.noise_cm2_s = size_band * 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 ModelCxbFlatBasedEnergy(ModelSrcInterface):
"""Model for a geometrically flat CXB with energy dependence
as specified with the Moretti spectrum.
"""
def __init__(self,
solidangle_filename=\
add_path_data_ref_eclairs('instru/pixels_solid_angle_withcostheta.fits')):
"""**constructor**
:param solidangle_filename: PATH/filename of the solid angle image
:type solidangle_filename: string
"""
super().__init__()
norm = 0.109
gamma_1 = 1.4
gamma_2 = 2.88
e_b = 29
self.solidangle_filename = solidangle_filename
shape_based_solid_angle = read_fits(solidangle_filename, ugts_standard=1)
mean_solid_angle = np.mean(shape_based_solid_angle)
solid_angle = np.ones(shape=shape_based_solid_angle.shape) * mean_solid_angle
self._name = "flat CXB Moretti spectrum"
self.s_angle_pixel = solid_angle * ModelSrcInterface.static_sim_geom.get_det_pixel_surface()
self.func = lambda e: norm / ((e / e_b) ** gamma_1 + (e / e_b) ** gamma_2)
self.solid_angle = solid_angle
[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)
.. warning :: a factor 0.8 is added to this model to
make the count rate equal to the CXB shape case. TBI.
"""
# ## initial function get_model as implemented by CC/PB
val = self.solid_angle[0, 0]
assert (self.solid_angle == val).all()
if not obs_time:
obs_time = ModelSrcInterface.static_duration
# flux = integral between e_min and e_max of eq. 4 of https://arxiv.org/pdf/0811.1444.pdf
flux, _ = quad(self.func,
ModelSrcInterface.static_e_min,
ModelSrcInterface.static_e_max)
ret = (obs_time * flux) * self.s_angle_pixel
self.logger.info(f'ADDED 0.8 FACTOR in get_model cxb flat !!! \
[eband=[{ModelSrcInterface.static_e_min, ModelSrcInterface.static_e_max}] keV]')
ret *= 0.8
return ret
[docs]class ModelCxbShapeBased(ModelSrcInterface):
"""Model for a CXB without spectrum based on the pixel solid angle shape
..warning :: dead code ?
"""
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.static_duration
ret = self.rate * obs_time
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_withcostheta.fits')):
"""**constructor**
:param solidangle_filename: PATH/filename of the solid angle image
:type solidangle_filename: string
"""
super().__init__()
norm = 0.109
gamma_1 = 1.4
gamma_2 = 2.88
e_b = 29 # keV
self.mask_aperture = ModelSrcInterface.static_sim_geom.inst.mask_aperture
self.npix_det = ModelSrcInterface.static_sim_geom.inst.p_nb_pix
self.solidangle_filename = solidangle_filename
self.solid_angle = read_fits(solidangle_filename, ugts_standard=1)
self._name = "energy CXB shape"
self.func = lambda e: norm / ((e / e_b) ** gamma_1 + (e / e_b) ** gamma_2)
det_pix_surf = ModelSrcInterface.static_sim_geom.inst.det_pix_size ** 2
self.s_angle_pix = 0.4 * det_pix_surf * self.solid_angle * (1. / self.mask_aperture)
[docs] def get_model(self, obs_time=None):
"""return the shadowgram model
..note: RMF must be handled here at some point !
..note: doc to Moretti's model:
flux = integral between e_min and e_max of eq. 4
of https://arxiv.org/pdf/0811.1444.pdf
: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.static_duration
# flux = integral between e_min and e_max of eq. 4 of https://arxiv.org/pdf/0811.1444.pdf
flux, _ = quad(self.func,
ModelSrcInterface.static_e_min,
ModelSrcInterface.static_e_max)
ret = (flux * obs_time) * self.s_angle_pix
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: single source information {'elev','dir','intensity','name'}.
Elev and dir in deg.
:type info_src: astropy.table.Row or dict
"""
super().__init__()
self.info_src = info_src
if (ModelPointSrcIRF.S_cpt == 0) and \
isinstance(info_src, Row) and \
verbose:
self.logger.info(f'source info: {self.info_src}')
self.logger.info(f'source info dtype: {self.info_src.dtype}')
ModelPointSrcIRF.S_cpt += 1
self.simu = ModelSrcInterface.static_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()
# super().__init__()
[docs] def set_coef_hors_axis(self):
"""Return coefficient associated to cos theta or sin elev off axis effect.
:return: coefficient (units ?)
:rtype: array(float)
"""
if self.simu is None:
self.simu = SimuECLAIRsMaskProjection()
pix_y, pix_z = self.simu.inst.elevdir_to_skypix(np.array(self.info_src['elev']),
np.array(self.info_src['dir']))
return ModelSrcInterface.static_mdl_effect.get_irf_val_ecllos(pix_y, pix_z)[0]
[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.static_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.static_idx_chan]
except ValueError:
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.
Elevation angle is in degree.
"""
def __init__(self, info_src):
"""**Constructor**
:param info_src: single source information {'elev','dir','intensity','name'}.
Elev and dir in deg.
:type info_src: astropy.table.Row or dict
"""
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):
"""This function must be renamed as 'get_coef_hors_axis'.
Verbose argument is unused.
"""
return self.sin_elev
[docs]class ListModelPointSrcFromCatalog(object):
"""Manage list of point source
"""
def __init__(self, o_dpix, use_irf=True):
"""**Constructor**
o_dpix used to define dpix channel
"""
assert isinstance(o_dpix, ECLAIRsDetectorEffect)
self._dpix = o_dpix
self.use_irf = use_irf
# declare logger
self.logger = logging.getLogger(__name__)
self.logger.info('creating an instance of ListModelPointSrcFromCatalog')
[docs] def set_use_irf(self, use_irf):
assert use_irf in [True, False]
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)
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 = {}
self.logger.info(f'Number of sources in FOV: {self.cat_fov.get_nb_element()}')
for isrc in range(self.cat_fov.get_nb_element()):
info_src = self.cat_fov._catalog[isrc]
mdl = self._get_model_point_src(info_src, verbose)
d_src[mdl.name] = mdl
if verbose:
self.logger.info(f'adding {info_src["name"]} source')
if verbose:
self.logger.info(f"Selected {len(d_src)} sources in FOV with SWIFT/BAT catalog")
return d_src
[docs] def init_with_cat_fov(self, cat_fov, verbose=False):
"""Initialize catalog with sources in FOV.
"""
assert isinstance(cat_fov, cat.CatalogFovEclairs)
d_src = {}
for isrc in range(cat_fov.get_nb_element()):
info_src = cat_fov._catalog[isrc]
if verbose:
self.logger.info(f'adding {info_src["name"]} source')
mdl = self._get_model_point_src(info_src, verbose)
d_src[mdl.name] = mdl
if verbose:
self.logger.info(f"Selected {len(d_src)} sources in FOV with user catalog")
return d_src