'''script to rebinned the mask in order to be used by the reconstruction algorithm
Created on 19 mars 2018
@author: Catalano Camille
'''
import os
import time
import matplotlib.pyplot as plt
import numpy as np
import common.sky.catalog as catX
import simu.lib.InstruX as simX
from common import add_path_data_ref_eclairs, add_path_tests_eclairs
from common.io import fits_tools
from common.instru.ref_geom import ECLAIRsGeom, add_kw_geom_eclairs_to_mask
from common.instru.array_convention import yzegp_to_ijdet_array
[docs]def mask_rebinning_v1():
'''build the rebinned mask by ray tracing simulation
the rebinned mask is constructed by simulation :
- projection of the mask on a larger detector 120x120pix
- source in the middle of the fov
- detector of 120x120 pix and of the same size of the mask (54cm)
rebinned mask saved as maskECLAIRs_rebinned.fits with same keywords as maskECLAIRs.fits
'''
#source in the middle of the fov with a large number of photon
ocat = catX.CatalogFovBasic()
ocat.add_src([ 90, 0, 1000, 'src_00']) # 1000 ph et 37s = 4min
#modification of the ECLAIRs' instrument for the mask projection
#pix_size = 4.5mm as before
#pix_nb = 120 in order to match the size of the mask (54cm)
simu_obs_x = simX.SimuECLAIRsRayTracing(ocat, p_nb_pixel=120)
#simulation (=projection of the mask on the larger detector)
simu_obs_x.simu_catalog(37)
simu_obs_x.inst.detect.plot("rebinned mask")
#normalization and binning to discrete values
rebinned_mask = simu_obs_x.inst.detect.ar_val.copy()
#additional simulation (first has reach memory limit)
simu_obs_x = simX.SimuECLAIRsRayTracing(ocat, p_nb_pixel=120)
simu_obs_x.simu_catalog(37)
rebinned_mask = rebinned_mask + simu_obs_x.inst.detect.ar_val
simu_obs_x = simX.SimuECLAIRsRayTracing(ocat, p_nb_pixel=120)
simu_obs_x.simu_catalog(37)
rebinned_mask = rebinned_mask + simu_obs_x.inst.detect.ar_val
simu_obs_x = simX.SimuECLAIRsRayTracing(ocat, p_nb_pixel=120)
simu_obs_x.simu_catalog(37)
rebinned_mask = rebinned_mask + simu_obs_x.inst.detect.ar_val
#custom histogram method
#zeros_value = rebinned_mask.max()/4
#one_value = rebinned_mask.max() - rebinned_mask.max()/4
nb_bin = 12
zeros_value = rebinned_mask.max()/nb_bin
one_value = rebinned_mask.max() - rebinned_mask.max()/nb_bin
step_bin = (one_value - zeros_value)/(nb_bin-2)
bins = np.concatenate(([0,zeros_value],(np.arange(nb_bin-2-1)+1)*step_bin+zeros_value,[one_value,rebinned_mask.max()+1]))
print(bins)
bins_indices = np.digitize(rebinned_mask,bins)
simu_obs_x.inst.detect.ar_val = (bins_indices-1)/(len(bins)-2)
#save the rebinned mask
geom = ECLAIRsGeom()
geom.save_with_kw_geom_eclairs(yzegp_to_ijdet_array(simu_obs_x.inst.detect.ar_val), add_path_data_ref_eclairs('maskECLAIRs_rebinned.fits'), "IJDET")
mask = fits_tools.read_fits(add_path_data_ref_eclairs('maskECLAIRs.fits'), ugts_standard=1)
geom.save_with_kw_geom_eclairs(mask, add_path_data_ref_eclairs('maskECLAIRs.fits'), "IJDET")
simu_obs_x.inst.detect.plot("rebinned mask")
print(len(np.where(simu_obs_x.inst.detect.ar_val == 0.)[0])/14400)
print(len(np.where(simu_obs_x.inst.detect.ar_val == 1.)[0])/14400)
[docs]def pixel_in_hole(p_ar, p_i, p_j):
if (p_ar[p_i+1, p_j-1]*p_ar[p_i+1, p_j]*p_ar[p_i+1, p_j+1] == 0):
return False
if (p_ar[p_i, p_j-1]*p_ar[p_i, p_j+1] == 0):
return False
if (p_ar[p_i-1, p_j-1]*p_ar[p_i-1, p_j]*p_ar[p_i-1, p_j+1] == 0):
return False
return True
[docs]def force_to_1_pixel_in_hole(p_ar):
for li in range(1, p_ar.shape[0]-1):
for lj in range(1, p_ar.shape[1]-1):
if pixel_in_hole(p_ar, li, lj):
p_ar[li, lj] = 1
[docs]def cannot_be_superius_to_1(p_ar, p_correct=True):
li, lj = np.where(p_ar > 1)
print("%d values superius to 1"%li.size)
if p_correct:
p_ar[li, lj] = 1
else:
print("pas de correction")
[docs]def mask_rebinning_v2(p_nb_ph_pix=100, p_nb_ite=4, p_path_out=None, p_prefix_out='mask'):
'''
build the rebinned mask by ray tracing simulation
the rebinned mask is constructed by simulation :
- projection of the mask on a larger detector 120x120pix
- source in the middle of the fov
- detector of 120x120 pix and of the same size of the mask (54cm)
rebinned mask saved as maskECLAIRs_rebinned.fits with same keywords as maskECLAIRs.fits
:param p_nb_ph_pix:
:type p_nb_ph_pix:
:param p_nb_ite:
:type p_nb_ite:
:param p_path_out:
:type p_path_out:
:param p_prefix_out:
:type p_prefix_out:
'''
#source in the middle of the fov with a large number of photon
ocat = catX.CatalogFovBasic()
intensity_src = 1 # ph.cm-2.s-1
ocat.add_src([ 90, 0, intensity_src, 'src_00'])
o_sim = simX.SimuECLAIRsRayTracing(ocat, p_nb_pixel=120)
time_exp_s = p_nb_ph_pix/(o_sim.inst.detect.surf_cell*intensity_src)
print("%.1fs de pose pour %d ph/pix"%(time_exp_s, p_nb_ph_pix))
a_cumul = np.zeros_like(o_sim.inst.detect.ar_val)
o_sim.b_photon_noise = False
for l_i in range(1, p_nb_ite+1):
print('itération %d:'%l_i)
o_sim.raz_nb_photon()
o_sim.simu_catalog(time_exp_s)
#o_sim.inst.detect.plot("rebinned mask")
if p_path_out:
file_name = p_prefix_out+"_%d_%d.fits"%(p_nb_ph_pix, l_i)
file_out = os.path.join(p_path_out,file_name)
o_sim.save_stack_evt_image(file_out)
a_cumul += o_sim.inst.detect.ar_val
a_cumul_nor = a_cumul/(p_nb_ite*p_nb_ph_pix)
plt.figure()
plt.title("Raw resampling mask for %d ph/pix"%(p_nb_ite*p_nb_ph_pix))
plt.imshow(yzegp_to_ijdet_array(a_cumul_nor))
plt.colorbar()
# correction
cannot_be_superius_to_1(a_cumul_nor, False)
force_to_1_pixel_in_hole(a_cumul_nor)
cannot_be_superius_to_1(a_cumul_nor)
plt.figure()
plt.title("Corrected resampling mask for %d ph/pix"%(p_nb_ite*p_nb_ph_pix))
plt.imshow(yzegp_to_ijdet_array(a_cumul_nor))
plt.colorbar()
# surface estimation
est_surf = a_cumul_nor.sum()*o_sim.inst.detect.surf_cell
print("surface estimation %.2f cm2"%est_surf)
print("surface théorique %.2f cm2"%o_sim.inst.get_open_surface_mask())
if p_path_out:
file_name = p_prefix_out+"_%d.fits"%(p_nb_ite*p_nb_ph_pix)
file_out = os.path.join(p_path_out,file_name)
o_sim.inst.detect.ar_val = a_cumul_nor
o_sim.save_stack_evt_image(file_out)
add_kw_geom_eclairs_to_mask(file_out, file_out)
[docs]def mask_proj_maskproj(p_path_out=None, p_prefix_out='mask'):
'''
build the rebinned mask by mask projection simulation
the rebinned mask is constructed by simulation :
- projection of the mask on a larger detector 120x120pix
- source in the middle of the fov
- detector of 120x120 pix and of the same size of the mask (54cm)
rebinned mask saved as maskECLAIRs_rebinned.fits with same keywords as maskECLAIRs.fits
:param p_path_out:
:type p_path_out:
:param p_prefix_out:
:type p_prefix_out:
'''
#source in the middle of the fov with a large number of photon
ocat = catX.CatalogFovBasic()
intensity_src = 1 # ph.cm-2.s-1
ocat.add_src([ 90, 0, intensity_src, 'src_00'])
o_sim = simX.SimuECLAIRsMaskProjection(p_nb_pixel=120)
o_sim.set_catalog(ocat)
o_sim.raz_nb_photon()
o_sim.b_photon_noise = False
o_sim.simu_catalog(1)
print("Durée calcul: %.3fs"%(time.time()-time_start))
plt.figure()
plt.title("resampling mask with projection method")
plt.imshow(yzegp_to_ijdet_array(o_sim.det_tmp.ar_val))
plt.colorbar()
# surface estimation
a_cumul_nor = o_sim.det_tmp.ar_val
est_surf = o_sim.det_tmp.ar_val.sum()*o_sim.inst.detect.surf_cell
print("surface estimation %.2f cm2"%est_surf)
print("surface théorique %.2f cm2"%o_sim.inst.get_open_surface_mask())
if p_path_out:
file_name = p_prefix_out
file_out = os.path.join(p_path_out,file_name+".fits")
o_sim.inst.detect.ar_val = a_cumul_nor
o_sim.save_stack_evt_image(file_out)
add_kw_geom_eclairs_to_mask(file_out, file_out)
if __name__ == '__main__':
time_start = time.time()
#mask_rebinning_v2(1000, 10, add_path_tests_eclairs(""))
mask_proj_maskproj(add_path_tests_eclairs(""), 'proj')
print("Durée fonction: %.3fs"%(time.time()-time_start))
plt.show()
#plt.show()