Galactic Diffuse Continuum
This tutorial demonstrates how to do a spectral fit of the Galactic diffuse continuum emission using an input GALPROP model in healpix format.
[8]:
import logging
import sys
import os
logger = logging.getLogger('cosipy')
logger.setLevel(logging.INFO)
logger.addHandler(logging.StreamHandler(sys.stdout))
os.environ["OMP_NUM_THREADS"] = "1"
os.environ["MKL_NUM_THREADS"] = "1"
os.environ["NUMEXPR_NUM_THREADS"] = "1"
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[8], line 5
2 import sys
3 import os
----> 5 sys.path.remove('/uni-mainz.de/homes/sgallego/software/COSItools/cosipy_israel')
8 logger = logging.getLogger('cosipy')
9 logger.setLevel(logging.INFO)
ValueError: list.remove(x): x not in list
[5]:
# imports:
from cosipy import BinnedData
from cosipy.statistics import PoissonLikelihood
from cosipy.background_estimation import FreeNormBinnedBackground
from cosipy.interfaces import ThreeMLPluginInterface
from cosipy.response import BinnedThreeMLModelFolding, BinnedInstrumentResponse, BinnedThreeMLExtendedSourceResponse
from cosipy.data_io import EmCDSBinnedData
from cosipy.spacecraftfile import SpacecraftHistory
from cosipy.response.FullDetectorResponse import FullDetectorResponse
from cosipy.response.ExtendedSourceResponse import ExtendedSourceResponse
from cosipy.threeml.custom_functions import GalpropHealpixModel
from cosipy.util import fetch_wasabi_file
from threeML import PointSource, Model, JointLikelihood, DataList, update_logging_level
from threeML.analysis_results import *
from astromodels import *
from astromodels.functions import GalPropTemplate_3D
import numpy as np
import matplotlib.pyplot as plt
17:24:41 WARNING The naima package is not available. Models that depend on it will not be functions.py:43 available
WARNING The GSL library or the pygsl wrapper cannot be loaded. Models that depend on it functions.py:65 will not be available.
WARNING The ebltable package is not available. Models that depend on it will not be absorption.py:33 available
17:24:42 INFO Starting 3ML! __init__.py:44
WARNING WARNINGs here are NOT errors __init__.py:45
WARNING but are inform you about optional packages that can be installed __init__.py:46
WARNING to disable these messages, turn off start_warning in your config file __init__.py:47
WARNING no display variable set. using backend for graphics without display (agg) __init__.py:53
WARNING ROOT minimizer not available minimization.py:1208
WARNING Multinest minimizer not available minimization.py:1218
WARNING PyGMO is not available minimization.py:1228
WARNING The cthreeML package is not installed. You will not be able to use plugins which __init__.py:95 require the C/C++ interface (currently HAWC)
WARNING Could not import plugin FermiLATLike.py. Do you have the relative instrument __init__.py:136 software installed and configured?
WARNING Could not import plugin HAWCLike.py. Do you have the relative instrument __init__.py:136 software installed and configured?
17:24:43 WARNING No fermitools installed lat_transient_builder.py:44
Get the data
[6]:
#path to download your data
data_path = Path("") # /path/to/files. Current dir by default
[7]:
# ori file
fetch_wasabi_file('COSI-SMEX/DC3/Data/Orientation/DC3_final_530km_3_month_with_slew_15sbins_GalacticEarth_SAA.ori',
output = data_path / "DC3_final_530km_3_month_with_slew_15sbins_GalacticEarth_SAA.ori" , checksum = 'e5e71e3528e39b855b0e4f74a1a2eebe')
A file named /localscratch/sgallego/linkToXauron/COSIpyData/DC3_data/DC3_final_530km_3_month_with_slew_15sbins_GalacticEarth_SAA.ori already exists with the specified checksum (e5e71e3528e39b855b0e4f74a1a2eebe). Skipping.
[ ]:
# background file
fetch_wasabi_file('COSI-SMEX/DC3/Data/Backgrounds/Ge/AlbedoPhotons_3months_unbinned_data_filtered_with_SAAcut.fits.gz',
output = data_path / "AlbedoPhotons_3months_unbinned_data_filtered_with_SAAcut.fits.gz" ,checksum = '191a451ee597fd2e4b1cf237fc72e6e2')
[ ]:
# source file
fetch_wasabi_file('COSI-SMEX/DC3/Data/Sources/GalTotal_SA100_F98_3months_unbinned_data_filtered_with_SAAcut.fits.gz',
output = data_path / "GalTotal_SA100_F98_3months_unbinned_data_filtered_with_SAAcut.fits.gz" ,checksum = '9fda5a7b15a90358abc2b886979f9fef')
[18]:
# precomputed point source response
fetch_wasabi_file('COSI-SMEX/DC3/Data/Responses/extended_source_response/extended_source_response_continuum_merged.h5.gz',
output = data_path / "extended_source_response_continuum_merged.h5.gz",unzip = True, checksum = '92ed7e22b1dafce6b57611d5cdb6cf70')
A file named /localscratch/sgallego/linkToXauron/COSIpyData/DC3_data/extended_source_response_continuum_merged.h5.gz with the same ETag (e1a6ee94cb0e557df15eb9cc25869ae1-72) as the requested file already exists. Skipping.
[5]:
# GALPROP input model
fetch_wasabi_file('COSI-SMEX/cosipy_tutorials/galactic_diffuse_continuum/total_healpix_57_SA100_F98_example.gz',
output = data_path / "total_healpix_57_SA100_F98_example.gz" ,checksum = '82cbeb9a86d86637f19f31c762f379fc')
Downloading cosi-pipeline-public/COSI-SMEX/cosipy_tutorials/galactic_diffuse_continuum/total_healpix_57_SA100_F98_example.gz
Input files:
[11]:
ori_file = data_path / "DC3_final_530km_3_month_with_slew_15sbins_GalacticEarth_SAA.ori"
BG_file = data_path / "AlbedoPhotons_3months_unbinned_data_filtered_with_SAAcut.fits.gz"
src_file = data_path / "GalTotal_SA100_F98_3months_unbinned_data_filtered_with_SAAcut.fits.gz"
psr_file = data_path / "extended_source_response_continuum_merged.h5"
galprop_model_file = data_path / "total_healpix_57_SA100_F98_example.gz"
Make the dataset and bin
This step only needs to be run once. Afterwards, the files can be loaded directly using the cell below.
[12]:
# Bin galdiff:
galdiff = BinnedData("galdiff.yaml")
galdiff.get_binned_data(unbinned_data=src_file, output_name=data_path/"galdiff_binned_data")
# Bin background:
bg_tot = BinnedData("galdiff.yaml")
bg_tot.get_binned_data(unbinned_data=BG_file, output_name=data_path/ "albedo_photons_binned_data")
binning data...
Time unit: s
Em unit: keV
Phi unit: deg
PsiChi unit: None
binning data...
Time unit: s
Em unit: keV
Phi unit: deg
PsiChi unit: None
Load binned files:
[12]:
#load signal
galdiff = BinnedData("galdiff.yaml")
galdiff.load_binned_data_from_hdf5(binned_data= data_path / "galdiff_binned_data.hdf5")
# Load background:
bg_tot = BinnedData("galdiff.yaml")
bg_tot.load_binned_data_from_hdf5(binned_data= data_path / "albedo_photons_binned_data.hdf5")
Define GALPROP model
Below is how to define the custom GALPROP model. We will save the model to a yaml file so that it can be directly uploaded in the future (as shown at the bottom).
[13]:
# defining the model:
galprop_model = GalpropHealpixModel()
galprop_model.load_file(galprop_model_file)
# The spectrum is defined in the data cube,
# and so we use a dummy model for defining an extended source in astromodels.
# NB: This has no impact on the results - just make sure the parameter is fixed!
spectrum = Constant()
spectrum.k.value = 0.0
spectrum.k.free = False
src = ExtendedSource("galprop_source", spatial_shape=galprop_model, spectral_shape=spectrum)
model = Model(src)
model.save("galprop_model.yaml", overwrite=True)
# uncomment below to load saved model:
#model = load_model('galprop_model.yaml')
loading GALPROP model: /localscratch/sgallego/linkToXauron/COSIpyData/DC3_data/total_healpix_57_SA100_F98_example.gz
Setup and perform fit
Setup background parameter, extended source , orientation file and the LH fit
[14]:
#open the orientation file
sc_orientation = SpacecraftHistory.open(ori_file)
#open the extended source response
dr = ExtendedSourceResponse.open(psr_file)
bkg_dist = {"background":bg_tot.binned_data.project('Em', 'Phi', 'PsiChi')}
# Workaround to avoid inf values. Out bkg should be smooth, but currently it's not.
# Reproduces results before refactoring. It's not _exactly_ the same, since this fudge value was 1e-12, and
# it was added to the expectation, not the normalized bkg
for bckfile in bkg_dist.keys() :
bkg_dist[bckfile] += sys.float_info.min
#combine the data + the bck like we would get for real data
data = EmCDSBinnedData(galdiff.binned_data.project('Em', 'Phi', 'PsiChi')
+ bg_tot.binned_data.project('Em', 'Phi', 'PsiChi')
)
bkg = FreeNormBinnedBackground(bkg_dist,
sc_history=sc_orientation,
copy = False)
# Currently using the same NnuLambda, Ei and Pol axes as the underlying FullDetectorResponse,
# matching the behavior of v0.3. This is all the current BinnedInstrumentResponse can do.
# In principle, this can be decoupled, and a BinnedInstrumentResponseInterface implementation
# can provide the response for an arbitrary directions, Ei and Pol values.
esr = BinnedThreeMLExtendedSourceResponse(data = data,
precomputed_psr = dr,
polarization_axis = dr.axes['Pol'] if 'Pol' in dr.axes.labels else None,
)
##====
response = BinnedThreeMLModelFolding(data = data, extended_source_response = esr)
like_fun = PoissonLikelihood(data, response, bkg)
cosi = ThreeMLPluginInterface('cosi',
like_fun,
response,
bkg)
cosi.bkg_parameter['background'] = Parameter("background", # background parameter
1, # initial value of parameter
unit = u.Hz,
min_value=0, # minimum value of parameter
max_value=10, # maximum value of parameter
delta=0.05, # initial step used by fitting engine
)
WARNING RuntimeWarning: divide by zero encountered in divide
Instantiate the COSI 3ML plugin and perform the fit
[15]:
# Optional: if you want to call get_log_like manually, then you also need to set the model manually
# 3ML does this internally during the fit though
cosi.set_model(model)
# Gather all plugins and combine with the model in a JointLikelihood object, then perform maximum likelihood fit
plugins = DataList(cosi) # If we had multiple instruments, we would do e.g. DataList(cosi, lat, hawc, ...)
like = JointLikelihood(model, plugins, verbose = False)
like.fit()
10:47:35 INFO set the minimizer to minuit joint_likelihood.py:1017
... Reading Extended source response ...
--> done (source name : galprop_source)
Interpolating GALPROP map...
Integrating intensity over energy bins...
Best fit values:
| result | unit | |
|---|---|---|
| parameter | ||
| galprop_source.GalpropHealpixModel.K | (7.721 +/- 0.008) x 10^-1 | |
| background | 4.5842 +/- 0.0013 | Hz |
Correlation matrix:
| 1.00 | -0.72 |
| -0.72 | 1.00 |
Values of -log(likelihood) at the minimum:
| -log(likelihood) | |
|---|---|
| cosi | -171655715.58496818 |
| total | -171655715.58496818 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | -343311427.16988426 |
| BIC | -343311406.47479194 |
[15]:
( value negative_error \
galprop_source.GalpropHealpixModel.K 0.772051 -0.000787
background 4.584189 -0.001347
positive_error error unit
galprop_source.GalpropHealpixModel.K 0.000823 0.000805
background 0.001292 0.001320 Hz ,
-log(likelihood)
cosi -171655715.58496818
total -171655715.58496818)
[16]:
results = like.results
print(results.display())
parameters = {par.name:results.get_variates(par.path)
for par in results.optimized_model["galprop_source"].parameters.values()
if par.free}
results_err = results.propagate(results.optimized_model["galprop_source"].spectrum.main.shape.evaluate_at, **parameters)
print(results.optimized_model["galprop_source"])
Best fit values:
| result | unit | |
|---|---|---|
| parameter | ||
| galprop_source.GalpropHealpixModel.K | (7.721 +/- 0.008) x 10^-1 | |
| background | 4.5842 +/- 0.0013 | Hz |
Correlation matrix:
| 1.00 | -0.72 |
| -0.72 | 1.00 |
Values of -log(likelihood) at the minimum:
| -log(likelihood) | |
|---|---|
| cosi | -171655715.58496818 |
| total | -171655715.58496818 |
Values of statistical measures:
| statistical measures | |
|---|---|
| AIC | -343311427.16988426 |
| BIC | -343311406.47479194 |
None
* galprop_source (extended source):
* shape:
* K:
* value: 0.7720511842515105
* desc: Normalization factor (unitless)
* min_value: 0.0
* max_value: 1000.0
* unit: ''
* is_normalization: true
* spectrum:
* main:
* Constant:
* k:
* value: 0.0
* desc: Constant value
* min_value: null
* max_value: null
* unit: keV-1 s-1 cm-2
* is_normalization: false
* polarization: {}
Make plots
Compare best-fit to injected source:
[11]:
%matplotlib inline
# Get expected counts from likelihood scan (i.e. best-fit convolved with response):
total_expectation = response.expectation(data)
print()
print("galdiff expected counts:")
print(total_expectation.project('Em').todense().contents)
print()
# Plot:
fig,ax = plt.subplots()
binned_energy_edges = galdiff.binned_data.axes['Em'].edges.value
binned_energy = galdiff.binned_data.axes['Em'].centers.value
ax.stairs(total_expectation.project('Em').todense().contents, binned_energy_edges, color='purple', label = "Best-fit convolved with response")
ax.errorbar(binned_energy, total_expectation.project('Em').todense().contents, yerr=np.sqrt(total_expectation.project('Em').todense().contents), color='purple', linewidth=0, elinewidth=1)
ax.stairs(galdiff.binned_data.project('Em').todense().contents, binned_energy_edges, color = 'black', ls = ":", label = "Simulated source counts")
ax.errorbar(binned_energy, galdiff.binned_data.project('Em').todense().contents, yerr=np.sqrt(galdiff.binned_data.project('Em').todense().contents), color='black', linewidth=0, elinewidth=1)
ax.set_xlabel("Energy (keV)")
ax.set_ylabel("Counts")
plt.yscale('log')
plt.xscale('log')
ax.legend()
plt.savefig("injected_model_comparison.pdf")
plt.show()
plt.close()
galdiff expected counts:
[7.75333154e+05 1.89566158e+06 1.76323392e+06 1.10662500e+06
6.56070357e+05 3.87738124e+05 2.35773666e+05 9.77338158e+04
2.13603691e+04 1.03920517e+03]
percent difference:
[12]:
diff = (galdiff.binned_data.project('Em').todense().contents - total_expectation.project('Em').todense().contents)/total_expectation.project('Em').todense().contents
plt.semilogx(binned_energy,diff,ls="--",marker="o")
plt.xlabel("Energy [keV]")
plt.ylabel("(data - model) / model")
plt.savefig("percent_diff.pdf")
plt.show()
plt.close()
Compare best-fit to injected for total counts:
[17]:
# Plot:
fig,ax = plt.subplots()
expectation_bkg = bkg.expectation(copy = True)
data_combined = galdiff.binned_data.project('Em').todense().contents+bg_tot.binned_data.project('Em').todense().contents
ax.stairs(total_expectation.project('Em').todense().contents \
+(expectation_bkg.project('Em').todense().contents) \
, binned_energy_edges, color='purple', label = "Best fit convolved with response plus BG")
ax.errorbar(binned_energy, total_expectation.project('Em').todense().contents \
+(expectation_bkg.project('Em').todense().contents) \
, yerr=np.sqrt(total_expectation.project('Em').todense().contents \
+(expectation_bkg.project('Em').todense().contents)), \
color='purple', linewidth=0, elinewidth=1)
ax.stairs(data_combined
, binned_energy_edges, \
color = 'black', ls = ":", label = "Simulated total counts")
ax.errorbar(binned_energy,data_combined, \
yerr=np.sqrt(data_combined), color='black', linewidth=0, elinewidth=1)
ax.set_xlabel("Energy (keV)")
ax.set_ylabel("Counts")
ax.legend()
plt.yscale('log')
plt.xscale('log')
plt.savefig("injected_total_comparison.pdf")
plt.show()
plt.close()
[18]:
mod_tot = total_expectation.project('Em').todense().contents \
+(expectation_bkg.project('Em').todense().contents)
diff = (data_combined - mod_tot)/mod_tot
plt.semilogx(binned_energy,diff,ls="--",marker="o")
plt.xlabel("Energy [keV]")
plt.ylabel("(data - model) / model")
plt.savefig("percent_diff.pdf")
plt.show()
plt.close()
Plot average intensity (averaged over full sky):
[58]:
# Get parameter error manually:
K = results.optimized_model["galprop_source"].parameters["galprop_source.GalpropHealpixModel.K"].value
Kerr = K - results.get_variates('galprop_source.GalpropHealpixModel.K').equal_tail_interval()[0]
# We will pass the response nside in order to use the same spatial sampling as the fit:
nside = 8
# We will also use the same energy values as was used in the fit:
binned_energy = galdiff.binned_data.axes['Em'].centers.to(u.MeV).value
binned_energy_edges = galdiff.binned_data.axes['Em'].edges.to(u.MeV).value
energy_err = np.diff(binned_energy_edges)/2.0
# Below we will pass avg_int=True in order to get the average intensity. Otherwise, the function returns the total intensity by default.
intensity = results.optimized_model["galprop_source"].spatial_shape.get_total_spatial_integral(binned_energy, avg_int=True, nside=nside)
intensity = intensity.value
yerr = (Kerr/K)*(intensity)
yerr *= (binned_energy**2)
print("intensity error:")
print(yerr)
intensity *= (binned_energy**2)
fig,ax = plt.subplots()
ax.loglog(binned_energy, intensity, ls="", marker="o", color="black", label = "Best-fit")
ax.errorbar(binned_energy, intensity, xerr=energy_err, yerr=yerr, ls="", marker="o", color="black", label = "_nolabel_")
# Plot model specturm with galpy:
# This is optional and requires galpy package:
# https://github.com/ckarwin/galpy
#from galpy import GalMapsHeal
#instance = GalMapsHeal()
#instance.read_healpix_file("GALPROP_DC3/total_healpix_57_SA100_F98_example.gz")
#instance.make_spectrum()
#gal_energy = instance.energy
#gal_spec = instance.spectra_list
#ax.loglog(gal_energy, gal_spec, ls="-", marker="", color="red", label = "GALPROP model")
plt.ylabel("dN/dE [$\mathrm{MeV \ cm^{-2} \ s^{-1} \ sr^{-1}}$]")
plt.xlabel("Energy [MeV]")
plt.title("Average Intensity")
ax.legend()
plt.xlim(5e-2,20)
plt.ylim(1e-4,5e-3)
plt.savefig("intensity.pdf")
plt.show()
plt.close()
using nside=8 from user input in evaluate method
loading GALPROP model: /localscratch/sgallego/linkToXauron/COSIpyData/DC3_data/total_healpix_57_SA100_F98_example.gz
Interpolating GALPROP map...
intensity error:
[6.33011004e-06 6.69254441e-06 6.98040174e-06 7.37073109e-06
7.79093418e-06 8.50441800e-06 9.24401661e-06 1.03413466e-05
1.13966358e-05 1.28929686e-05]
Below we plot the best-fit spectrum just for demonstration. Again, this is just a dummy model since the spectrum is contained in the 3D data cube. This has no impact on the fit.
[23]:
energy = np.linspace(100.,10000.,201)*u.keV
flux = results.optimized_model["galprop_source"].spectrum.main.shape(energy)
fig,ax = plt.subplots()
ax.plot(energy, flux, label = "Best fit")
plt.ylabel("dN/dE [$\mathrm{ph \ cm^{-2} \ s^{-1} \ keV^{-1}}$]", fontsize=14)
plt.xlabel("Energy [keV]", fontsize=14)
ax.legend()
plt.savefig("best_fit_model.pdf")
plt.show()
plt.close()