"""
Functions to read an astrophysical catalogue (e.g. Swift/BAT), compute count rate through the ARF
and save infos in a fits file that will be read by the simulator.
"""
import logging
import numpy as np
from astropy.io import fits
import astropy.units as u
import ecpi.common.sky.catalog as cata
from ecpi.common.instru.model_effect import ECLAIRsDetectorEffect
logger = logging.getLogger(__name__)
[docs]def swift_catalog_to_eclairs_catalog(cat_filename, dpix):
"""return a catalog with mean flux in channel
:param cat_filename: path to fits catalog
:type cat_filename: string
:param dpix: object calibration ECLAIRs
:type dpix:
:param verbose:
:type verbose:
"""
assert isinstance(dpix, ECLAIRsDetectorEffect)
f_fits = fits.open(cat_filename)
sources_cat = f_fits[1].data
f_fits.close()
# in keV. Energy band of the Swift/BAT catalogue. Assumed to be always the same.
band_bat = np.array([14., 195.])
spectral_index = sources_cat['GAMMA']
# (given in 1e-12 erg/s/cm^2 in the Swift/BAT catalogue)
flux_integrated_kev = ((sources_cat['FLUX'] * 1e-12) * u.erg).to('keV').value
catalog_of_sources = cata.CatalogAstroWithEnergySpecSampling(dpix.chan_center)
for idx_src in range(len(sources_cat)):
spec_index = spectral_index[idx_src]
src_name = sources_cat['COUNTERPART_NAME'][idx_src]
# logger.debug(f"process '{src_name}' spec_idx={spec_index:.2f}")
fact_indice = (2 - spec_index)
if np.isclose(spec_index, 2, atol=1e-3):
phi0 = flux_integrated_kev[idx_src] / np.log(band_bat[1] / band_bat[0])
else:
band_power = np.power(band_bat, 2 - spec_index)
phi0 = fact_indice * flux_integrated_kev[idx_src] / (band_power[1] - band_power[0])
# spec_index
fact_indice = (1 - spec_index)
if np.isclose(spec_index, 1, atol=1e-3):
boundary_log = np.log(dpix.chan_boundary)
flux_in_channel = phi0 * (boundary_log[1:] - boundary_log[:-1])
logger.info(f'{src_name} closed to 1 {spec_index}')
else:
boundary_pow = np.power(dpix.chan_boundary, 1 - spec_index)
flux_in_channel = phi0 * (boundary_pow[1:] - boundary_pow[:-1]) / fact_indice
# logger.debug(f'current source name: {src_name}')
try:
assert (flux_in_channel > 0).all()
except AssertionError:
logger.error(f'source name: {src_name} / spectral index: {spec_index}')
raise AssertionError
my_src = (sources_cat['RA'][idx_src], sources_cat['DEC'][idx_src],
flux_in_channel, src_name)
catalog_of_sources.add_src(my_src)
return catalog_of_sources