Generating DRMs for GBM detectors¶
The main purpose for gbm_drm_gen is a fast and easy way to generate detector response matrices for the Fermi-GBM detectors in a pythonic way.
When we want to generate DRMS for GBM detectors, we simply need to obtain the proper data
[1]:
import matplotlib.pyplot as plt
import astropy.io.fits as fits
import numpy as np
%matplotlib inline
from jupyterthemes import jtplot
jtplot.style(context="talk", grid=False)
from gbm_drm_gen import DRMGen
from gbm_drm_gen.utils.package_data import get_path_of_data_file
WARNING: AstropyDeprecationWarning: The private astropy._erfa module has been made into its own package, pyerfa, which is a dependency of astropy and can be imported directly using "import erfa" [astropy._erfa]
[WARNING ] The naima package is not available. Models that depend on it will not be available
[WARNING ] The GSL library or the pygsl wrapper cannot be loaded. Models that depend on it will not be available.
[WARNING ] The ebltable package is not available. Models that depend on it will not be available
[INFO ] Starting 3ML!
[WARNING ] no display variable set. using backend for graphics without display (agg)
[WARNING ] ROOT minimizer not available
[WARNING ] Multinest minimizer not available
[WARNING ] PyGMO is not available
[WARNING ] The cthreeML package is not installed. You will not be able to use plugins which require the C/C++ interface (currently HAWC)
[WARNING ] Could not import plugin FermiLATLike.py. Do you have the relative instrument software installed and configured?
[WARNING ] Could not import plugin HAWCLike.py. Do you have the relative instrument software installed and configured?
[WARNING ] Env. variable OMP_NUM_THREADS is not set. Please set it to 1 for optimal performances in 3ML
[WARNING ] Env. variable MKL_NUM_THREADS is not set. Please set it to 1 for optimal performances in 3ML
[WARNING ] Env. variable NUMEXPR_NUM_THREADS is not set. Please set it to 1 for optimal performances in 3ML
quick start¶
To create a DRM generator for TTE data, we need the TTE, CSPEC, and the TRIGDAT data files. * The CSPEC data contains the output side of the DRM’s energy bounds. * The TRIGDAT data contains the spacecraft orientation data.
[2]:
trigdat_file = get_path_of_data_file('example_data/glg_trigdat_all_bn110721200_v01.fit')
cspec_file = get_path_of_data_file('example_data/glg_cspec_n6_bn110721200_v00.pha')
# create the generator
gbm_n6_generator = DRMGen.from_128_bin_data(det_name= "n6",
time=0, # time relative to T0 or trigger time.
trigdat = trigdat_file,
mat_type = 2, # direct response + atmospheric scattering
cspecfile = cspec_file)
We can set the location of the source directly. The first run can be a bit slow as numba is used in the background to be very fast.
[3]:
gbm_n6_generator.set_location(ra=329,dec = -38.2)
We can now checkout the matrix object created in the background:
[4]:
gbm_n6_generator.matrix
[4]:
array([[1.21529588e-02, 2.66685346e-02, 3.35702294e-02, ...,
1.71696052e-02, 1.79967322e-02, 0.00000000e+00],
[2.40656348e-03, 9.02243768e-03, 2.44969585e-02, ...,
1.76004235e-02, 1.84410737e-02, 0.00000000e+00],
[0.00000000e+00, 6.64336828e-04, 4.79652366e-03, ...,
1.75084305e-02, 1.83343065e-02, 0.00000000e+00],
...,
[0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
1.18137462e-01, 1.18589681e-01, 0.00000000e+00],
[0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
1.19105016e-01, 1.19567934e-01, 0.00000000e+00],
[0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
2.49905465e+00, 2.50882594e+00, 0.00000000e+00]])
Or we can input RA and DEC to create a 3ML style OGIP response directly:
[5]:
response = gbm_n6_generator.to_3ML_response(ra=329,dec = -38.2)
[6]:
fig = response.plot_matrix()
WARNING MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. This has been deprecated since 3.3 and in 3.6, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = mpl.cm.get_cmap("ocean").copy()
To see how the effective area varies with location, we can loop through various angles.
[7]:
fig, ax = plt.subplots()
bounds = np.vstack((gbm_n6_generator.monte_carlo_energies[:-1],gbm_n6_generator.monte_carlo_energies[1:])).T
de = np.diff(bounds)
ene = np.mean(bounds,axis=1)
for ra in np.linspace(260, 350, 10):
gbm_n6_generator.set_location(ra=ra,dec = -38.2)
ax.loglog(ene,gbm_n6_generator.matrix.sum(axis=0),label=r'%d$^{\circ}$'%ra)
ax.set_ylim(1)
ax.legend()
ax.set_xlabel(r'Effective Area (cm$^2$)');
ax.set_xlabel('MC Energies');
Into the details¶
Ok, now let’s go through the various specifics of the DRM generator constructor.
First, the generator needs to know:
The current location and oreintation of GBM for the data of interest.
The channel to PHA reconstructed energy from the CSPEC files.
For the first, one can get the current spacecraft from the either the triggers trigdat file, or for non-trggered data, one can obtain the position history file. These are available at the NASA database. If using a position history file. You need to specify T0=<Fermi MET> in the constructor so that the time coordinate will be relative to this MET.
Internally, the class uses gbmgeometry to convert RA,Dec to the approriate spacecraft coordinates. However, one can also create DRMs in spacecraft coordinates directly.
custom energy binning¶
Maybe you are a curious person and want to investigate a response with finer input energies to model line features in solar flares?
It is possible to add a custom array of input energies. To do this we need to import the NaITTEEdges and BgoTTEEdges classes.
[8]:
from gbm_drm_gen import NaiTTEEdges, BgoTTEEdges
These objects allow you to specify input (monte carlo) energies either from an array or in log spaced binning.
Note: The number of energies must be off and include the low and high end points which are specific to the NaI [5.0, 50000.0] and BGO [100., 200000.0] detectors.
Log-spaced energies¶
First we will try with VERY fine log spaced binning above the typical 140 input energies
[9]:
custom_edges = NaiTTEEdges.from_log_bins(n_bins=531)
gbm_n6_generator = DRMGen.from_128_bin_data(det_name= "n6",
trigdat = trigdat_file,
mat_type = 2, # direct response + atmospheric scattering
cspecfile = cspec_file,
# pass the custom edges
custom_input_edges=custom_edges
)
response = gbm_n6_generator.to_3ML_response(ra=329,dec = -38.2)
fig = response.plot_matrix()
WARNING MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. This has been deprecated since 3.3 and in 3.6, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = mpl.cm.get_cmap("ocean").copy()
And now with a much coarser input binning:
[10]:
custom_edges = NaiTTEEdges.from_log_bins(n_bins=73)
gbm_n6_generator = DRMGen.from_128_bin_data(det_name= "n6",
trigdat = trigdat_file,
mat_type = 2, # direct response + atmospheric scattering
cspecfile = cspec_file,
# pass the custom edges
custom_input_edges=custom_edges
)
response = gbm_n6_generator.to_3ML_response(ra=329,dec = -38.2)
fig = response.plot_matrix()
It is easier to see the difference with simple matrix plotting:
[11]:
custom_edges = NaiTTEEdges.from_log_bins(n_bins=91)
gbm_n6_generator = DRMGen.from_128_bin_data(det_name= "n6",
trigdat = trigdat_file,
mat_type = 2, # direct response + atmospheric scattering
cspecfile = cspec_file,
# pass the custom edges
custom_input_edges=custom_edges
)
gbm_n6_generator.set_location(ra=329,dec = -38.2)
plt.matshow(gbm_n6_generator.matrix.T)
custom_edges = NaiTTEEdges.from_log_bins(n_bins=141)
gbm_n6_generator = DRMGen.from_128_bin_data(det_name= "n6",
trigdat = trigdat_file,
mat_type = 2, # direct response + atmospheric scattering
cspecfile = cspec_file,
# pass the custom edges
custom_input_edges=custom_edges
)
gbm_n6_generator.set_location(ra=329,dec = -38.2)
plt.matshow(gbm_n6_generator.matrix.T)
custom_edges = NaiTTEEdges.from_log_bins(n_bins=1541)
gbm_n6_generator = DRMGen.from_128_bin_data(det_name= "n6",
trigdat = trigdat_file,
mat_type = 2, # direct response + atmospheric scattering
cspecfile = cspec_file,
# pass the custom edges
custom_input_edges=custom_edges
)
gbm_n6_generator.set_location(ra=329,dec = -38.2)
plt.matshow(gbm_n6_generator.matrix.T)
[11]:
<matplotlib.image.AxesImage at 0x7f1d1443b820>
WARNING RuntimeWarning: overflow encountered in double_scalars
WARNING RuntimeWarning: overflow encountered in double_scalars
a custom array¶
And we can even supply and entirely custom array of energies:
[12]:
edges = np.linspace(5., 50000., 1001)
custom_edges = NaiTTEEdges.from_custom_array(edges)
gbm_n6_generator = DRMGen.from_128_bin_data(det_name= "n6",
trigdat = trigdat_file,
mat_type = 2, # direct response + atmospheric scattering
cspecfile = cspec_file,
# pass the custom edges
custom_input_edges=custom_edges
)
response = gbm_n6_generator.to_3ML_response(ra=329,dec = -38.2)
fig = response.plot_matrix()
Creating RSP2 files¶
[13]:
from gbm_drm_gen import create_rsp2
[14]:
# create the generator
gbm_n6_generator = DRMGen.from_128_bin_data(det_name= "n6",
time=0, # time relative to T0 or trigger time.
trigdat = trigdat_file,
mat_type = 2, # direct response + atmospheric scattering
cspecfile = cspec_file)
[15]:
output_file_name = "my_new_rsp.rsp2" # you must call it an RSP2 file!
[16]:
create_rsp2(output_file_name,
response_generator=gbm_n6_generator,
ra=0,
dec=0,
tstart=0,
tstop=10,
delta_time=2,
overwrite=True
)
[17]:
with fits.open(output_file_name) as f:
f.info()
print(f[0].header['DRM_NUM'])
print(f[2].header['TSTART'])
print(f[2].header['TSTOP'])
Filename: my_new_rsp.rsp2
No. Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 5 ()
1 EBOUNDS 1 BinTableHDU 28 128R x 3C [K, E, E]
2 SPECRESP MATRIX 1 BinTableHDU 42 140R x 6C [E, E, I, I, I, 128D]
3 SPECRESP MATRIX 2 BinTableHDU 42 140R x 6C [E, E, I, I, I, 128D]
4 SPECRESP MATRIX 3 BinTableHDU 42 140R x 6C [E, E, I, I, I, 128D]
5 SPECRESP MATRIX 4 BinTableHDU 42 140R x 6C [E, E, I, I, I, 128D]
6 SPECRESP MATRIX 5 BinTableHDU 42 140R x 6C [E, E, I, I, I, 128D]
5
332916465.760476
332916467.760476
[18]:
from threeML import TimeSeriesBuilder
[19]:
tsb = TimeSeriesBuilder.from_gbm_cspec_or_ctime("cspec", cspec_or_ctime_file=cspec_file, rsp_file="my_new_rsp.rsp2")
tsb.view_lightcurve();
[WARNING ] The default choice for MATRIX extension failed:KeyError("Extension ('MATRIX', 1) not found.")available: None 'EBOUNDS' 'SPECRESP MATRIX' 'SPECRESP MATRIX' 'SPECRESP MATRIX' 'SPECRESP MATRIX' 'SPECRESP MATRIX'
WARNING FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
WARNING FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
[WARNING ] The default choice for MATRIX extension failed:KeyError("Extension ('MATRIX', 1) not found.")available: None 'EBOUNDS' 'SPECRESP MATRIX' 'SPECRESP MATRIX' 'SPECRESP MATRIX' 'SPECRESP MATRIX' 'SPECRESP MATRIX'
[WARNING ] The default choice for MATRIX extension failed:KeyError("Extension ('MATRIX', 2) not found.")available: None 'EBOUNDS' 'SPECRESP MATRIX' 'SPECRESP MATRIX' 'SPECRESP MATRIX' 'SPECRESP MATRIX' 'SPECRESP MATRIX'
[WARNING ] The default choice for MATRIX extension failed:KeyError("Extension ('MATRIX', 3) not found.")available: None 'EBOUNDS' 'SPECRESP MATRIX' 'SPECRESP MATRIX' 'SPECRESP MATRIX' 'SPECRESP MATRIX' 'SPECRESP MATRIX'
[WARNING ] The default choice for MATRIX extension failed:KeyError("Extension ('MATRIX', 4) not found.")available: None 'EBOUNDS' 'SPECRESP MATRIX' 'SPECRESP MATRIX' 'SPECRESP MATRIX' 'SPECRESP MATRIX' 'SPECRESP MATRIX'
[WARNING ] The default choice for MATRIX extension failed:KeyError("Extension ('MATRIX', 5) not found.")available: None 'EBOUNDS' 'SPECRESP MATRIX' 'SPECRESP MATRIX' 'SPECRESP MATRIX' 'SPECRESP MATRIX' 'SPECRESP MATRIX'
[20]:
tsb.set_active_time_interval('1-5')
[INFO ] Interval set to 0.7039839625358582-4.800044000148773 for cspec