"""Geometric description of the X-ray instrument
"""
import logging
import os.path as osp
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from matplotlib.patches import Rectangle
from matplotlib.collections import PatchCollection
from ecpi.common.io import fits_tools
from ecpi.common.instru.array_convention import yzegp_to_ijdet_array
from ecpi.common.num.array_square_cell import ArraySquareCell
from ecpi.common import add_path_data_ref_eclairs, add_path_caldb_eclairs
s_logger = logging.getLogger(__name__)
[docs]class InstruXbase(object):
"""
Base class for X-ray instrument
mask convention : 1 for hole, 0 stop photon
"""
def __init__(self, p_mask: ArraySquareCell, p_detec: ArraySquareCell, p_dist_mask_detec):
'''**Constructor**
'''
self.mask = ArraySquareCell()
self.mask = p_mask
self.mask.set_origin_center()
self.detect = ArraySquareCell()
self.detect = p_detec
assert isinstance(self.mask, ArraySquareCell)
assert isinstance(self.detect, ArraySquareCell)
self.detect.set_origin_center()
self.dist_mask_detect = p_dist_mask_detec # cm
# declare logger
s_logger.info('creating an instance of InstruXbase')
[docs] def set_pos_source_fov(self, elevation, direction):
"""
set definition of the translation vector linked to the source
vec_trans_source is also position on mask where photon passed to go on origin
"""
self.vec_trans_source = np.array([0, 0], dtype=np.float32)
dd = self.dist_mask_detect / np.tan(elevation)
self.vec_trans_source[0] = np.cos(direction)
self.vec_trans_source[1] = np.sin(direction)
self.vec_trans_source *= dd
[docs] def get_surface_rect_proj(self):
"""
compute intersection detector and shadow mask then compute its surface
:return: intersection surface in cm2
:rtype: float
"""
self.rect_proj = self.detect.intersection(self.mask, -self.vec_trans_source)
# why -self.vec_trans_source ?
# we want intersection in detector referentiel so we translate position mask to detector
inter_rec = self.rect_proj
surface = (inter_rec[1, 0] - inter_rec[0, 0]) * (inter_rec[1, 1] - inter_rec[0, 1])
return surface
[docs] def get_full_surface_mask(self):
"""in cm2"""
return self.mask.get_surface()
[docs] def get_open_surface_mask(self):
"""in cm2"""
return self.mask.ar_val.sum() * self.mask.surf_cell
[docs] def get_percent_open_mask(self):
return self.get_open_surface_mask() / self.get_full_surface_mask()
[docs] def print_total_count(self):
return self.detect.ar_val.sum()
[docs] def get_elevation_limit_fov_side(self):
"""for special case where mask and detector are square and with the same center
return deg
"""
half_size_mask = (self.get_full_surface_mask() ** .5) / 2
half_size_detect = (self.detect.get_surface() ** .5) / 2
angle = np.arctan(self.dist_mask_detect / (half_size_mask + half_size_detect))
return np.rad2deg(angle)
[docs] def get_elevation_limit_fov_corner(self):
"""for special case where mask and detector are square and with the same center
return deg
"""
half_size_mask = (self.get_full_surface_mask() ** .5) / (2 ** .5)
half_size_detect = (self.detect.get_surface() ** .5) / (2 ** .5)
angle = np.arctan(self.dist_mask_detect / (half_size_mask + half_size_detect))
return np.rad2deg(angle)
[docs] def skypix_to_elevdir(self, pix_y, pix_z):
"""convert pixel coordinates into elevation and direction
(0,0) is at the center of the sky array.
elevation from 0 to 90.
direction from -180 to +180
:param pix_y: y coordinate of the source, in pixels.
:type pix_y: float
:param pix_z: z coordinates of the source, in pixels.
:type pix_z: float
:return: elevation and direction of the source in deg
:rtype: float, float
"""
direction = np.rad2deg(np.arctan2(pix_z, pix_y))
norm_vec_trans_source = np.sqrt(pix_y ** 2 + pix_z ** 2) * self.detect.size_cm
elev = np.rad2deg(np.arctan2(self.dist_mask_detect, norm_vec_trans_source))
return elev, direction
[docs] def elevdir_to_skypix(self, elev, direction):
"""convert elevation/direction coordinates into sky pixels coordinates
(0,0) is at the center of the sky array
:param elev: elevation in degrees
:type elev: float or array(float)
:param direction: direction in degrees
:type direction: float or array(float)
:return: pixels y,z coordinates
:rtype: float, float or array2D(float)
"""
dd = self.dist_mask_detect / np.tan(np.radians(elev)) / self.detect.size_cm
pix_y = dd * np.cos(np.radians(direction))
pix_z = dd * np.sin(np.radians(direction))
return pix_y, pix_z
[docs]class InstruXMultiMask(InstruXbase):
"""To manage multiple mask for X-ray instrument
"""
def __init__(self, p_detec: ArraySquareCell, p_dist_mask_detec):
"""**Constructor**
"""
self.detect = p_detec
self.detect.set_origin_center()
self.dist_mask_detect = p_dist_mask_detec
self._l_mask = []
self._d_mask = {}
self.mask = ArraySquareCell()
[docs] def add_mask(self, p_mask: ArraySquareCell, p_name=None):
p_mask.perct_open = p_mask.ar_val.sum() * 1.0 / p_mask.ar_val.size
if p_name == None:
p_mask.name = str(len(self._l_mask) + 1)
else:
p_mask.name = p_name
self._d_mask[p_mask.name] = p_mask
self._l_mask.append(p_mask)
[docs] def set_mask(self, p_mask):
self.mask = p_mask
[docs] def select_mask(self, p_idx):
self.mask = self._l_mask[p_idx]
[docs] def get_mask_name(self, p_name):
return self._d_mask[p_name]
[docs] def plot_submasks(self): # pragma: no cover
for mask in self._l_mask:
mask.plot()
plt.figure()
s_logger.debug(f'lenght of _l_mask: {len(self._l_mask)}')
plt.title('Pos mask')
for mask in self._l_mask:
plt.plot(mask.to_origin_cm[0], mask.to_origin_cm[1], 'o')
[docs] def plot_mask(self, mes="mask"): # pragma: no cover
# plt.figure()
_, ax = plt.subplots(1)
plt.title(mes)
ax.patch.set_facecolor('red')
amax = 27.01
ax.set_xlim([-amax, amax])
ax.set_ylim([-amax, amax])
l_holes = []
# rect = Rectangle((0, 0), 1, 1)
# l_holes.append(rect)
# rect = Rectangle((3, 3), 1, 1)
# l_holes.append(rect)
for mask in self._l_mask:
assert isinstance(mask, ArraySquareCell)
for idx_x in range(mask.get_nb_cell_x()):
for idx_y in range(mask.get_nb_cell_y()):
idx = np.array([[idx_x, idx_y]])
if mask.value(idx[0]) == 1:
pos = mask.idx2pos_inf(idx)[0]
# print(pos[0], pos[1])
rect = Rectangle((pos[0], pos[1]), mask.size_cm, mask.size_cm)
l_holes.append(rect)
pc = PatchCollection(l_holes, facecolor='white', alpha=1)
# Add collection to axes
ax.add_collection(pc)
plt.xlabel("Y ECLLos, cm [RADIATOR]*")
plt.ylabel("Z ECLLos, cm")
[docs] def get_full_surface_mask(self):
raise
[docs] def get_open_surface_mask(self):
"""in cm2"""
s_surf = 0
for r_mask in self._l_mask:
s_surf += r_mask.ar_val.sum() * r_mask.surf_cell
return s_surf
[docs] def get_surface_rect_proj_allmasks(self):
"""compute the surface on the detector of the projection of the masks
(the coded detector surface). It is source dependent
:return: coded detector surface in cm2
:rtype: float
"""
inter_rec = np.array([[99999, 99999], [-99999, -99999]])
for mask in self._l_mask:
self.set_mask(mask)
self.rect_proj = self.detect.intersection(self.mask, -self.vec_trans_source)
inter_rec[0, 0] = min(inter_rec[0, 0], self.rect_proj[0, 0])
inter_rec[0, 1] = min(inter_rec[0, 1], self.rect_proj[0, 1])
inter_rec[1, 0] = max(inter_rec[1, 0], self.rect_proj[1, 0])
inter_rec[1, 1] = max(inter_rec[1, 1], self.rect_proj[1, 1])
coded_surface = (inter_rec[1, 0] - inter_rec[0, 0]) * (inter_rec[1, 1] - inter_rec[0, 1])
self.select_mask(0)
return coded_surface
[docs]class InstruECLAIRs(InstruXMultiMask):
"""Implementation of the ECLAIRs mask geometry in ECL_Los referential
"""
def __init__(self, p_file_mask=None, p_nb_pix=None):
"""**Constructor**
"""
if p_file_mask is None:
p_file_mask = add_path_data_ref_eclairs(osp.join("instru", "maskECLAIRs.fits"))
# referential detector is center of detector
mask, _, header = fits_tools.set_array_and_header_from_file(p_file_mask,
ugts_standard=True)
self._raw_mask = mask
# get the characteristics of the instrument from keywords
if p_nb_pix is None:
# number of detector pixel
self.p_nb_pix = header['DETPIX']
else:
self.p_nb_pix = p_nb_pix
# thickness of the mask cross in cm
self.half_cross = header['MASNERV'] / 2
# size of a detector pixel in cm
self.det_pix_size = header['DETPSIZE']
# pitch between detector pixels in cm
self.pitch_pix_size = header['PIXPITCH']
# distance between mask and detector in cm
self.dist_mask_detect = header['DETMAS']
# size of a mask pixel in cm
self.mask_pix_size = header['MASPSIZE']
# number of mask pixels in a row
self.mask_nb_pix = header['MASPIX']
# size of the mask in cm
self.mask_total_size = header['MASSIZE']
# average open fraction of the mask
self.mask_aperture = header['ALPHA']
# compute half FOV
size_mask_dpix = self.mask_total_size + self.p_nb_pix * self.pitch_pix_size
arg = np.sqrt(2) * size_mask_dpix / (2 * self.dist_mask_detect)
# TODO: test value self.half_fov
self.half_fov = np.rad2deg(np.arctan(arg))
detect = ArraySquareCell(self.p_nb_pix, self.p_nb_pix, self.pitch_pix_size)
super().__init__(detect, self.dist_mask_detect)
# Now origin is the center of detector
pm23_hc = self.mask_pix_size * self.mask_nb_pix / 2 + self.half_cross
# left bottom
m1 = ArraySquareCell(self.mask_nb_pix // 2, self.mask_nb_pix // 2, self.mask_pix_size)
m1.set_array(mask[:self.mask_nb_pix // 2, :self.mask_nb_pix // 2])
pos_corner_lb = np.array([-pm23_hc, -pm23_hc], dtype=np.float32)
m1.set_pos_corner_lb(pos_corner_lb)
self.add_mask(m1, "lb")
# right bottom
m2 = ArraySquareCell(self.mask_nb_pix // 2, self.mask_nb_pix // 2, self.mask_pix_size)
m2.set_array(mask[self.mask_nb_pix // 2:, :self.mask_nb_pix // 2])
pos_corner_lb = np.array([self.half_cross, -pm23_hc], dtype=np.float32)
m2.set_pos_corner_lb(pos_corner_lb)
self.add_mask(m2, "rb")
# left up
m3 = ArraySquareCell(self.mask_nb_pix // 2, self.mask_nb_pix // 2, self.mask_pix_size)
m3.set_array(mask[:self.mask_nb_pix // 2, self.mask_nb_pix // 2:])
pos_corner_lb = np.array([-pm23_hc, self.half_cross], dtype=np.float32)
m3.set_pos_corner_lb(pos_corner_lb)
self.add_mask(m3, "lu")
# right up
m4 = ArraySquareCell(self.mask_nb_pix // 2, self.mask_nb_pix // 2, self.mask_pix_size)
m4.set_array(mask[self.mask_nb_pix // 2:, self.mask_nb_pix // 2:])
pos_corner_lb = np.array([self.half_cross, self.half_cross], dtype=np.float32)
m4.set_pos_corner_lb(pos_corner_lb)
self.add_mask(m4, "ru")
self.select_mask(0)
def __repr__(self):
desc = """\n
InstruECLAIRs object
=====================
Mask aperture: {}
Half-cross size (cm): {}
Detector pixel size (cm): {}
Detector pixel pitch size (cm): {}
Mask-detector distance (cm): {}
Mask pixel size (cm): {}
Number of mask pixels: {}
Mask total size (cm): {}
Left-bottom submask
-------------------
y range: [{}, {}] | z range: [{}, {}]
Right-bottom submask
-------------------
y range: [{}, {}] | z range: [{}, {}]
Left-up submask
-------------------
y range: [{}, {}] | z range: [{}, {}]
Right-up submask
-------------------
y range: [{}, {}] | z range: [{}, {}]
"""
lb_y_range, lb_z_range = self.get_mask_name('lb').pos_center_cells()
rb_y_range, rb_z_range = self.get_mask_name('rb').pos_center_cells()
lu_y_range, lu_z_range = self.get_mask_name('lu').pos_center_cells()
ru_y_range, ru_z_range = self.get_mask_name('ru').pos_center_cells()
return desc.format(self.mask_aperture,
self.half_cross,
self.det_pix_size,
self.pitch_pix_size,
self.dist_mask_detect,
self.mask_pix_size,
self.mask_nb_pix,
self.mask_total_size,
lb_y_range[0], lb_y_range[-1], lb_z_range[0], lb_z_range[-1],
rb_y_range[0], rb_y_range[-1], rb_z_range[0], rb_z_range[-1],
lu_y_range[0], lu_y_range[-1], lu_z_range[0], lu_z_range[-1],
ru_y_range[0], ru_y_range[-1], ru_z_range[0], ru_z_range[-1],)
[docs] def get_full_surface_mask(self):
s_surf = 4 * self.mask.get_surface()
# add cross without center
s_surf += (2 * self.half_cross) * (4 * self.mask.get_sizecm()[0])
# add center cross
s_surf += (2 * self.half_cross) ** 2
return s_surf
[docs] def get_half_fov(self):
"""return in deg diagonal the half of field of view of ECLAIRs instrument
"""
return self.half_fov
[docs] def plot_eclairs_on_axis(self, a_pos, p_title=""): # pragma: no cover
_, ax = plt.subplots(1)
lb_corner = tuple(self.detect.pos_corner_lb())
for mask in self._l_mask:
lb_corner = tuple(mask.pos_corner_lb())
size_mask = mask.get_sizecm()
rect = patches.Rectangle(lb_corner, size_mask[0], size_mask[1],\
linewidth=1, edgecolor='r', facecolor='none')
# Add the patch to the Axes
ax.add_patch(rect)
plt.plot(a_pos[:, 0], a_pos[:, 1], '.')
plt.grid()
plt.title(p_title)
[docs] def plot_raw_mask(self): # pragma: no cover
plt.figure()
plt.title("raw fits mask")
plt.imshow(yzegp_to_ijdet_array(self._raw_mask))