Continuum Background Estimation

This example shows how to use cosipy for estimating the continuum background. In short, the method is based on the traditional on-off analysis, where the background for a source region on the sky is estimated by performing some kind of interpolation of the nearby surrounding region. The main difference with a Compton telescope is that we are now performing the on-off analysis in the Compton data space.

In this particular example, we want to estimate the background for the Crab. We start with the full dataset. This contains the full background, which includes instrumental + astrophysical, where the latter consists of all astrophysical sources other than the Crab. Then, for each bin of Em and Phi, we mask the Crab in the PsiChi plane based on the point source response. Finally, we interpolate over the masked region using an inpainting method.

The accuracy of this method depends strongly on two things:

1. Choosing the proper source region to mask

2. Making an accurate interpolation

The current source code takes a very simple appoach for both of these, but eventually they need to be improved. For the first point, ultimately we should use the ARM measurement, and we’ll need to figure out the optimal percentage of counts to mask. A major challenge here is that point sources have really long tails in the ARM distribution, extending over the entire sky. So we probably can’t use a typical confidence level of 95% for the masking. For the second point, we are currently using the simplest possible inpainting algorithm. More sophisticated methods are needed. The best alrorithm will likely come from deep convolution neural networks. An alterantive to inpainting methods is to use background templates, which can be scaled outside of the masked region. The code also needs to be developed in a way that will make this easy to do. Finally, the code is super slow and needs to be vectorized.

from cosipy.background_estimation import ContinuumEstimation
from cosipy.spacecraftfile import SpacecraftFile
from cosipy.util import fetch_wasabi_file
import os
import logging
import astropy.units as u
from astropy.coordinates import SkyCoord
logging.basicConfig()
logging.getLogger().setLevel(logging.INFO)
%matplotlib inline

15:24:51 WARNING   The naima package is not available. Models that depend on it will not be         functions.py:48
                  available                                                                                        

         WARNING   The GSL library or the pygsl wrapper cannot be loaded. Models that depend on it  functions.py:69
                  will not be available.                                                                           

15:24:52 WARNING   The ebltable package is not available. Models that depend on it will not be     absorption.py:33
                  available                                                                                        

15:24:53 INFO      Starting 3ML!                                                                     __init__.py:39

         WARNING   WARNINGs here are NOT errors                                                      __init__.py:40

         WARNING   but are inform you about optional packages that can be installed                  __init__.py:41

         WARNING    to disable these messages, turn off start_warning in your config file            __init__.py:44

         WARNING   no display variable set. using backend for graphics without display (agg)         __init__.py:50

15:24:54 WARNING   ROOT minimizer not available                                                minimization.py:1345

         WARNING   Multinest minimizer not available                                           minimization.py:1357

         WARNING   PyGMO is not available                                                      minimization.py:1369

15:24:54 WARNING   The cthreeML package is not installed. You will not be able to use plugins which  __init__.py:94
                  require the C/C++ interface (currently HAWC)                                                     

         WARNING   Could not import plugin FermiLATLike.py. Do you have the relative instrument     __init__.py:144
                  software installed and configured?                                                               

         WARNING   Could not import plugin HAWCLike.py. Do you have the relative instrument         __init__.py:144
                  software installed and configured?                                                               

15:24:57 WARNING   No fermitools installed                                              lat_transient_builder.py:44

15:24:57 WARNING   Env. variable MKL_NUM_THREADS is not set. Please set it to 1 for optimal         __init__.py:387
                  performances in 3ML                                                                              

         WARNING   Env. variable NUMEXPR_NUM_THREADS is not set. Please set it to 1 for optimal     __init__.py:387
                  performances in 3ML                                                                              

The notebook requires the following files:

1. DC3_final_530km_3_month_with_slew_15sbins_GalacticEarth.ori

2. SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.earthocc.h5

3. crab_bkg_binned_data_galactic.hdf5

4. inputs_crab.yaml

They can be downloaded using the cells below.

# 20280301_3_month_modifies.ori
fetch_wasabi_file('COSI-SMEX/DC3/Data/Orientation/DC3_final_530km_3_month_with_slew_15sbins_GalacticEarth.ori')

# Detector response file
# Make sure to unzip file after downloading!
fetch_wasabi_file('COSI-SMEX/DC2/Responses/SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.earthocc.h5.zip')

# crab_bkg_binned_data_galactic.hdf5
fetch_wasabi_file('COSI-SMEX/cosipy_tutorials/background_estimation/crab_bkg_binned_data_galactic.hdf5')

# inputs_crab.yaml
fetch_wasabi_file('COSI-SMEX/cosipy_tutorials/background_estimation/inputs_crab.yaml')

Define instance of class:

instance = ContinuumEstimation()

In order to estimate the background, we need the point source response. If you don’t already have this, you can calculate it, as shown below. Note that the coordinates of the Crab need to be passed as a tuple, giving Galactic longitude and latitude in degrees.

data_path = "your/path"

# Orientatin file:
ori_file = os.path.join(".","DC3_final_530km_3_month_with_slew_15sbins_GalacticEarth.ori")

# Spacecraft orientation:
sc_orientation = SpacecraftFile.parse_from_file(ori_file)

# Detector response:
dr = os.path.join(data_path,\
        "SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.earthocc.h5")

crab = SkyCoord(l=184.56*u.deg,b=-5.78*u.deg,frame="galactic")
psr = instance.calc_psr(sc_orientation, dr, crab)

# Alternatively, you can load a PSR from file with:
#psr = instance.load_psr_from_file("psr_file_name.h5")

Now let’s calculate the estimated background. To make a short example, we’ll only consider 1 Em bin and 2 Phi bins, as specified by the optional keywords e_loop and s_loop, respectively. We’ll also make plots here for demonstrational purposes.

Note that the current code has not yet been optimized for speed, as it uses a simple nested for loop. The time required to generate the estimated background using all bins is roughyly 4 hours. The option to use a subset of the Em bins and/or Phi bins may be useful for analyses that also use a given subset, but at this point the main motivation for this option is for demonstrational purposes, and when using this option, nothing is done with the other bins.

data_file = "crab_bkg_binned_data_galactic.hdf5"
data_yaml = "inputs_crab.yaml"
estimated_bg = instance.continuum_bg_estimation(data_file, data_yaml, psr, make_plots=True, e_loop=(2,3), s_loop=(4,6))

INFO:yayc.configurator:Using configuration file at inputs_crab.yaml
  0%|          | 0/2 [00:00<?, ?it/s]INFO:cosipy.background_estimation.ContinuumEstimation:Bin 2 4

100%|██████████| 2/2 [00:42<00:00, 21.49s/it]INFO:cosipy.background_estimation.ContinuumEstimation:Bin 2 5

100%|██████████| 2/2 [01:25<00:00, 42.71s/it]

Finaly, let’s save the estimated background to file:

estimated_bg.write("estimated_background",overwrite=True)