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()

#

# print('done with cube computation')

#### Christopher's way

# 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)

#### my way

# extract data from FITS

map_exp = fits.getdata('ctexpcube_GC_survey_ROI_12x12deg_100Ebins_0p01deg.fits.gz',0)

E_exp = np.array(fits.getdata('ctexpcube_GC_survey_ROI_12x12deg_100Ebins_0p01deg.fits.gz',1))['Energy']

map_cr =fits.getdata('ctmodel_CR_GC_survey_ROI_12x12deg_100Ebins_0p01deg.fits.gz',0)

# convert energies from keV to GeV

E_exp /= 1.e6

# exposure is calculated on bin edges, determine bin centers using geometric mean

E_cr = np.sqrt(E_exp[1:] * E_exp[:-1])

# set mask, all bins where exposure at lower edge is 0

mask = np.ones([np.shape(map_exp)[0] - 1,np.shape(map_exp)[1],np.shape(map_exp)[2]])

mask[map_exp[:-1,:,:] == 0.] = 0.

#### all bins

# # derive exposure at bin center as PL interp of edges values

# # index

# idx = (np.log(map_exp[1:]) - np.log(map_exp[:-1])) / (np.log(E_exp[1:,np.newaxis,np.newaxis]) - np.log(E_exp[:-1,np.newaxis,np.newaxis]))

# # set NaNs to 0, this occurs for bins with exposure = 0 that are going to be taken care of later

# idx[np.isnan(idx)] = 0.

# # exposure at bin center

# exp = map_exp[:-1] * np.power(E_cr[:,np.newaxis,np.newaxis]/E_exp[:-1,np.newaxis,np.newaxis], idx)

#

# # divide cr counts by exposure to get "flux"

# cr_fluxes = map_cr / exp

#

# # # apply mask to flux map

# cr_fluxes = cr_fluxes * mask

#

# # sum over spatial bins

# all_cr_fluxes = np.nansum(cr_fluxes,axis=(1,2))

#### just central bin

ipix = 600

# derive exposure at bin center as PL interp of edges values

# index

idx = (np.log(map_exp[1:,ipix,ipix]) - np.log(map_exp[:-1,ipix,ipix])) / (np.log(E_exp[1:]) - np.log(E_exp[:-1]))

# set NaNs to 0, this occurs for bins with exposure = 0 that are going to be taken care of later

idx[np.isnan(idx)] = 0.

# exposure at bin center

exp = map_exp[:-1,ipix,ipix] * np.power(E_cr/E_exp[:-1], idx)

# divide cr counts by exposure to get "flux"

cr_fluxes = map_cr[:,ipix,ipix] / exp

all_cr_fluxes = cr_fluxes

# plot

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.savefig('bin_600.png')