"""ECLAIRs ray tracing simulation method
"""
import logging
import numpy as np
from numba import njit, int8, int16, int32, void
import ecpi.common.instru.model_geom as iX
from ecpi.simu.lib.instru_x import SimuInstruXBase
from ecpi.common.instru.model_geom import InstruXbase
from ecpi.common.sky.catalog import CatalogFovBasic
logger = logging.getLogger(__name__)
[docs]class SimuInstruXRayTracing(SimuInstruXBase):
"""
old class
..warning
simulation instru X with
* catalog FOV => stable pointing
* evt in pixel => no energy
* CXB model uniforme noise on detector
"""
def __init__(self, instru_x: InstruXbase, catalog_x: CatalogFovBasic = None):
"""**Constructor**
"""
self.inst = instru_x
if catalog_x:
self.cat_fov = catalog_x
super().__init__()
self.det_pixel_surface = self.inst.detect.size_cm ** 2
self._use_numba = True
[docs] def simu_catalog(self, time_expose, verbose=None):
self.time_expose_s = time_expose
# loop on source and call ray_tracing_for_one_src()
src_idx = 0
# init individual source shadowgram and the
# intersection coordinates pixels on
# the detector only for a raytrac object
if not isinstance(self, SimuECLAIRsRayTracing):
arr_nb_cell = [self.inst.detect.get_nb_cell_x() - 1,
self.inst.detect.get_nb_cell_y() - 1]
arr_inter_src = [arr_nb_cell, [0, 0]]
self.shadows_sources = np.zeros((len(self.cat_fov.get_catalog()),
self.inst.detect.get_nb_cell_x(),
self.inst.detect.get_nb_cell_y()))
self.intersect_sources = np.zeros((len(self.cat_fov.get_catalog()), 2, 2))\
+arr_inter_src
for src in self.cat_fov.get_catalog():
# compute rectangle mask projection in detector
el_rad = np.deg2rad(src['elev'])
dr_rad = np.deg2rad(src['dir'])
self.inst.set_pos_source_fov(el_rad, dr_rad)
surface = self.inst.get_surface_rect_proj()
if verbose is True:
self.logger.info(f'surface: {surface}')
if np.isclose(0, surface, 1e-7):
# src out of detector
continue
# compute nb photon with
# - time
# - surface
# - angle between source and detector
# - intensity of source
# - noise photon of source
surface_open = surface * np.sin(el_rad)
nb_photon_theo = src['intensity'] * surface_open * time_expose
# add photon noise
if self.b_photon_noise:
nb_photon = np.random.poisson(nb_photon_theo)
else:
nb_photon = nb_photon_theo
self.nb_photon_src = nb_photon
apos = self.inst.detect.random_pos_rect(int(round(nb_photon)), self.inst.rect_proj)
self.ray_tracing_for_one_src(apos)
min_inter_src = np.minimum(self.intersect_sources[src_idx, 0],
self.intersect_det_subsrc[0])
max_inter_src = np.maximum(self.intersect_sources[src_idx, 1],
self.intersect_det_subsrc[1])
self.intersect_sources[src_idx] = [min_inter_src, max_inter_src]
self.shadows_sources[src_idx] += self.shadow_src
src_idx += 1
if verbose is True:
self.logger.info(f'process simulation {src["name"]} with {nb_photon} photons')
# self.inst.detect.plot()
[docs] def ray_tracing_for_one_src(self, apos):
"""
result: select by ray tracing photon can be arrived on detector by hole mask
"""
if self._use_numba:
# if False:
det = self.inst.detect
self.shadow_src = np.zeros(det.ar_val.shape, dtype=np.int32)
shadow = self.shadow_src
mask = self.inst.mask.ar_val
vec_trans_source = self.inst.vec_trans_source
size_p_mask = self.inst.mask.size_cm
size_p_det = det.size_cm
trans_mask = self.inst.mask.to_origin_cm
trans_det = det.to_origin_cm
idx_det = ray_tracing_for_one_src_numba(shadow, mask, apos, vec_trans_source,
size_p_mask, size_p_det,
trans_mask, trans_det)
det.ar_val += shadow
# intersection on the detector pix for this source and this mask
self.intersect_det_subsrc = [[idx_det[:, 0].min(), idx_det[:, 1].min()],
[idx_det[:, 0].max(), idx_det[:, 1].max()]]
else:
self.ray_tracing_for_one_src_numpy(apos)
[docs] def ray_tracing_for_one_src_numpy(self, apos):
"""
result: select by ray tracing photon can be arrived on detector by hole mask
"""
apos_mask = apos + self.inst.vec_trans_source
idx_det = self.inst.detect.pos2idx_inf(apos)
idx_mask = self.inst.mask.pos2idx_inf(apos_mask)
self.logger.debug(f'index detector: {idx_det}')
self.logger.debug(f'index mask: {idx_mask}')
det = self.inst.detect
mask = self.inst.mask
# apply mask on photon in field of view
if idx_mask.size:
aval_mask = (mask.value_array(idx_mask) + 0.5).astype(np.int8)
self.logger.debug(f'shape of aval_mask: {aval_mask.shape}')
self.logger.debug(f'sum of aval_mask: {aval_mask.sum()}')
self.logger.debug(f'normalized sum of aval_mask: {aval_mask.sum()/aval_mask.shape[0]}')
self.logger.debug(f'detect.ar_val.sum(): {self.inst.detect.ar_val.sum()}')
self.shadow_src = np.zeros(det.ar_val.shape, dtype=np.int32)
if self._use_numba:
self.logger.info('using numba')
summation_loop_for_ray_tracing(idx_det[:, 0], idx_det[:, 1],
aval_mask, self.shadow_src)
else:
it = np.nditer([idx_det[:, 0], idx_det[:, 1], aval_mask], [],
[['readonly'], ['readonly'], ['readonly']])
for il, ic, val in it:
# save individual source shadowgram for this source and this mask
self.shadow_src[il, ic] += val
det.ar_val += self.shadow_src
# intersection on the detector pix for this source and this mask
self.intersect_det_subsrc = [
[idx_det[:, 0].min(), idx_det[:, 1].min()],
[idx_det[:, 0].max(), idx_det[:, 1].max()]
]
self.logger.debug(f'sum of detect.ar_val: {self.inst.detect.ar_val.sum()}')
self.logger.debug(f'max of detect.ar_val: {self.inst.detect.ar_val.max()}')
[docs]@njit(cache=True)
def ray_tracing_for_one_src_numba(shadow, mask, apos, vec_trans_source,
size_p_mask, size_p_det,
lbpos_mask, lbpos_det): # pragma: no cover
"""compute shadowgram count with ray tracing method with photon position on detector
full computation with numba, ie compute index position instead call array_square_cell method
:param shadow: {out} shadowgram with count
:type shadow: 2D array
:param mask:
:type mask: 2D array
:param apos: photon position on detector
:type apos: [cm] 2D array (n,2)
:param vec_trans_source: translation vector of the source
:type vec_trans_source: 1D array (2,)
:param size_p_mask: size of pixel mask
:type size_p_mask: [cm] real
:param size_p_det: size of pixel detector
:type size_p_det: [cm] real
:param lbpos_mask: position of left bottom corner position of mask
:type lbpos_mask: [cm] 1D array (2,)
:param lbpos_det: position of left bottom corner position of detector
:type lbpos_det: [cm] 1D array (2,)
:return: index of the photon on detector
:rtype: [int] 2D array (n,2)
"""
# index of the photon on detector
idx_det = np.floor((apos + lbpos_det) / size_p_det).astype(np.int16)
# index of the photon on mask
idx_mask = np.floor((apos + vec_trans_source + lbpos_mask) / size_p_mask).astype(np.int16)
if idx_mask.size:
size_ar = apos.shape[0]
for idx in range(size_ar):
shadow[idx_det[idx, 0], idx_det[idx, 1]] += mask[idx_mask[idx, 0], idx_mask[idx, 1]]
return idx_det
[docs]@njit
def ray_tracing_photon_list_numba(on_detect, mask, apos, vec_trans_source,
size_p_mask, lbpos_mask): # pragma: no cover
"""compute shadowgram count with ray tracing method with photon position on detector
full computation with numba, ie compute index position instead call array_square_cell method
:param shadow: {out} shadowgram with count
:type shadow: 2D array
:param mask:
:type mask: 2D array
:param apos: photon position on detector
:type apos: [cm] 2D array (n,2)
:param vec_trans_source: translation vector of the source
:type vec_trans_source: 1D array (2,)
:param size_p_mask: size of pixel mask
:type size_p_mask: [cm] real
:param lbpos_mask: position of left bottom corner position of mask
:type lbpos_mask: [cm] 1D array (2,)
:return: index of the photon on detector
:rtype: [int] 2D array (n,2)
"""
# index of the photon on mask
idx_mask = np.floor((apos + vec_trans_source + lbpos_mask) / size_p_mask).astype(np.int16)
if idx_mask.size:
size_ar = apos.shape[0]
for idx in range(size_ar):
on_detect[idx] = mask[idx_mask[idx, 0], idx_mask[idx, 1]]
[docs]@njit(void(int16[:], int16[:], int8[:], int32[:, :]))
def summation_loop_for_ray_tracing(row, col, aval_mask, shadow): # pragma: no cover
"""compute shadowgram count with ray tracing method with index position of photon on mask
summation of the photons that have passed through the mask
:param row: row index of the mask where the photon arrives
:type row: [int] 1D array (n,)
:param col: row index of the mask where the photon arrives
:type col: [int] 1D array (n,)
:param aval_mask: pixel mask value (1 is open)
:type aval_mask: [0,1] 2D array
:param shadow: photon count on pixel detector
:type shadow: [int] 2D array
"""
size_ar = row.shape[0]
# print(f'loop size: {size_ar}')
for idx in range(size_ar):
shadow[row[idx], col[idx]] += aval_mask[idx]
[docs]class SimuECLAIRsRayTracing(SimuInstruXRayTracing):
""" Ray tracing simulation for ECLAIRs
Essentially loop on sub-mask
"""
def __init__(self, catalog_x: CatalogFovBasic = None, nb_pixel=None):
eclairs = iX.InstruECLAIRs(p_nb_pix=nb_pixel)
self.inst = eclairs
if catalog_x:
self.cat_fov = catalog_x
super().__init__(eclairs)
[docs] def simu_catalog(self, time_expose, verbose=None, show=None):
self.time_expose_s = time_expose
# init individual source shadowgram and the intersection
# coordinates pixels on the detector
self.shadows_sources = np.zeros((len(self.cat_fov.get_catalog()),
self.inst.p_nb_pix, self.inst.p_nb_pix))
arr = [[self.inst.p_nb_pix - 1, self.inst.p_nb_pix - 1], [0, 0]]
self.intersect_sources = np.zeros((len(self.cat_fov.get_catalog()), 2, 2)) + arr
for mask in self.inst._l_mask:
if verbose is True:
self.logger.info(f'mask name: {mask.name}')
self.inst.set_mask(mask)
super().simu_catalog(time_expose, verbose=verbose)
if show is True:
self.inst.detect.plot()