"""
Created on 29 nov. 2017
@author: user
"""
import logging
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
[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
self.logger = logging.getLogger(__name__)
self.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 get_total_photon(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
: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()
self.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): # pragma: no cover
_, ax = plt.subplots(1)
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.show()
[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.
"""
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('instru/maskECLAIRs.fits')
# referential detector is center of detector
mask, _, header = fits_tools.set_array_and_header_from_file(p_file_mask,
ugts_standard=1)
self._raw_mask = mask
# get the characteristics of the instrument from keywords
if p_nb_pix is None:
self.p_nb_pix = header['DETPIX']
else:
self.p_nb_pix = p_nb_pix
self.half_cross = header['MASNERV'] / 2
self.det_pix_size = header['DETPSIZE']
self.pitch_pix_size = header['PIXPITCH']
self.dist_mask_detect = header['DETMAS']
self.mask_pix_size = header['MASPSIZE']
self.mask_nb_pix = header['MASPIX']
self.mask_total_size = header['MASSIZE']
self.mask_aperture = header['ALPHA']
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
=====================
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.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 plot_eclairs_on_axis(self, a_pos, p_title=""):
fig, 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):
plt.figure()
plt.title("raw fits mask")
plt.imshow(yzegp_to_ijdet_array(self._raw_mask))