ctools

import numpy as np

from astropy.io import fits

from scipy.interpolate import interp1d

import matplotlib.pyplot as plt

import gammalib

import ctools

import cscripts

sim_cr = ctools.ctmodel()

sim_cr['inobs'] = 'obs_CTA_GC_survey_ROI_5deg.xml'

sim_cr['emin'] = 0.03

sim_cr['emax'] = 100.0

sim_cr['inmodel'] = 'cr.xml'

sim_cr['incube'] = 'bin_generic_GC_survey_ROI_12x12deg_100EBINS_0p01deg.fits.gz'

sim_cr['outcube'] = 'ctmodel_CR_GC_survey_ROI_12x12deg_100Ebins_0p01deg.fits.gz'

sim_cr.execute()

sim_expcube = ctools.ctexpcube()

sim_expcube['inobs'] = 'obs_CTA_GC_survey_ROI_5deg.xml'

sim_expcube['emin'] = 0.03

sim_expcube['emax'] = 100.0

sim_expcube['incube'] = 'bin_generic_GC_survey_ROI_12x12deg_100EBINS_0p01deg.fits.gz'

sim_expcube['outcube'] = 'ctexpcube_GC_survey_ROI_12x12deg_100Ebins_0p01deg.fits.gz'

sim_expcube.execute()

expcube = fits.open('ctexpcube_GC_survey_ROI_12x12deg_100Ebins_0p01deg.fits.gz')

map_expcube = np.array(expcube[0].data).astype('float')

E_expcube = np.array(expcube[1].data).astype('float')

f_expsoure = interp1d(E_expcube*1e-03, map_expcube, axis = 0)

cr = fits.open('ctmodel_CR_GC_survey_ROI_12x12deg_100Ebins_0p01deg.fits.gz')

map_cr = np.array(cr[0].data).astype('float') 

E_cr = np.array([np.sqrt(cr[2].data[i][0] * cr[2].data[i][1]) * 1e-06 for i in range(len(cr[2].data))]) ###calculate residual cosmic-ray flux at an energy bin's central energy (w.r.t. logarithmic-spacing)

all_cr_fluxes = []

for i, E in enumerate(E_cr):

    exposure_shape = np.array([1.0 if x != 0. else 0. for x in f_expsoure(E).flatten()]) ### Both maps might have zeros at different positions, I try to lift them on the same footing by imposing the structure of the exposure map on the CR background map

    cr_flux = map_cr[i].flatten() * exposure_shape / np.array([x if x != 0. else 1.0 for x in f_expsoure(E).flatten()])  ### To avoid divisony by zero, since both maps have their zeros at the same position.s

    cr_flux = cr_flux.sum()

    all_cr_fluxes.append(cr_flux)

plt.plot(E_cr, all_cr_fluxes)

plt.yscale('log')

plt.xscale('log')

plt.xlabel('Energy [GeV]')

plt.ylabel('Residual cosmic-ray flux [arbitrary units]')

plt.show()