Transient Background Modeling in CDS
This notebook demonstrates a background-modeling method for transients. Traditionally, one models the background by fitting the light curve before and after the burst. However, that approach only uses the count information after projecting to the Time axis. cosipy operates in Compton Data Space (CDS), so the background must be modeled in CDS as well. We therefore need a different approach that preserves the multi-dimensional structure of the data.
Overview
The general idea is the same:
Define background windows before and after the burst.
Sum the data in those background windows and scale the result to the burst interval.
We provide two scaling methods: duration and fitting.
Scaling Methods
duration
Convert the summed background model into a unit-rate background (i.e., a 1 s background).
Multiply by the burst duration so the expected total counts scale to the burst interval.
fitting
This will be supported in the future.
Fit the burst background level in a light curve using the background windows.
Use the fitted result to estimate the number of background counts in the burst interval.
Scale the summed background model to match the fitted counts.
[1]:
from pathlib import Path
from histpy import Histogram
from cosipy import TransientBackgroundEstimation
import matplotlib.pyplot as plt
/home/sheng2/Applications/conda_envs/cosipy_yong_dev/lib/python3.10/site-packages/threeML/io/package_data.py:4: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
import pkg_resources
22:21:08 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.
22:21:09 WARNING The ebltable package is not available. Models that depend on it will not be absorption.py:33 available
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
22:21:09 WARNING ROOT minimizer not available minimization.py:1345
WARNING Multinest minimizer not available minimization.py:1357
WARNING PyGMO is not available minimization.py:1369
22:21:10 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?
22:21:11 WARNING No fermitools installed lat_transient_builder.py:44
22:21:11 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
Load and examine data
[2]:
data_dir = Path("") # Current directory by default. Modify if you want a different path
[4]:
# Read the data
file_name = "GRB_bn110605183_with_bkg.hdf5"
data = Histogram.open(data_dir / file_name)
[5]:
# initialize background estimator
estimator = TransientBackgroundEstimation(data)
[6]:
# view the full light curve
%matplotlib inline
estimator.plot_lightcurve()
[7]:
# You can zoom in the light curve by using `plot_limts` parameters
estimator.plot_lightcurve(plot_limits = [1.841896e9 + 400, 1.841896e9+450])
Define Burst Window
You can estimate the burst window by the zoom-in light curve.
The burst windows can be recognized by Bayesian blocks. It will be implemented in future updates.
[8]:
# add a burst window
estimator.add_burst_windows([1.8418964e9+4, 1.8418964e9+48])
[9]:
# plot the light curve with the burst window
estimator.plot_lightcurve(burst_windows=True)
Define Background Windows
You can define the background windows like you did for the burst window.
[10]:
estimator.add_bkg_windows([1.841896e9+200, 1.841896e9+300], [1.841896e9+500, 1.841896e9+600])
[11]:
estimator.plot_lightcurve(burst_windows=True, bkg_windows=True)
Make the background model
[12]:
# It only supports scaling by duratiom now.
bkg_model = estimator.make_background_model(scaling = "duration", save_path = data_dir / "bkg_model.hdf5")
[13]:
bkg_model.project("PsiChi").plot()
plt.title("Background model (PsiChi)")
Could not determine default ax_kw['coord']. WCSAxes do not support intricsic coordinates <SpacecraftFrame Frame (attitude=None, obstime=None, location=None)>. Fallling back to plotting onto ICRS coodinates
Failed to transform from 'spacecraftframe' to 'icrs'. Rasterizing in 'spacecraftframe' frame. ERROR: Spacecraft coordinates need attitude to transform from ICRS
[13]:
Text(0.5, 1.0, 'Background model (PsiChi)')
Compare with the true background
[14]:
bkg_true = Histogram.open("True_GRB_bn110605183_background.hdf5")
bkg_true.project("PsiChi").plot()
plt.title("Real background (PsiChi)")
Could not determine default ax_kw['coord']. WCSAxes do not support intricsic coordinates <SpacecraftFrame Frame (attitude=None, obstime=None, location=None)>. Fallling back to plotting onto ICRS coodinates
Failed to transform from 'spacecraftframe' to 'icrs'. Rasterizing in 'spacecraftframe' frame. ERROR: Spacecraft coordinates need attitude to transform from ICRS
[14]:
Text(0.5, 1.0, 'Real background (PsiChi)')
[15]:
bkg_diff = bkg_true.project("PsiChi") - bkg_model.project("PsiChi")
bkg_diff.plot()
plt.title("Real background (PsiChi)")
Could not determine default ax_kw['coord']. WCSAxes do not support intricsic coordinates <SpacecraftFrame Frame (attitude=None, obstime=None, location=None)>. Fallling back to plotting onto ICRS coodinates
Failed to transform from 'spacecraftframe' to 'icrs'. Rasterizing in 'spacecraftframe' frame. ERROR: Spacecraft coordinates need attitude to transform from ICRS
[15]:
Text(0.5, 1.0, 'Real background (PsiChi)')
The difference map is centered around \(0\) with variations up to \(\pm5\)
[ ]: