This tutorial demonstrates how to perform an analysis of the Draco dSph galaxy. This tutorial assumes that you have first gone through the PG 1553 analysis tutorial. In this example we will use the following data selection which is chosen to match the selection used in the 6-year LAT Dwarf Analysis.
The P8 dSph paper used a multi-component analysis that used the four PSF event types in a joint likelihood. In this example we will perform a single component analysis using all SOURCE-class events (evtype=3).
%matplotlib inline
import os
import numpy as np
from fermipy.gtanalysis import GTAnalysis
from fermipy.plotting import ROIPlotter, SEDPlotter
import matplotlib.pyplot as plt
import matplotlib
WARNING: AstropyDeprecationWarning: astropy.extern.six will be removed in 4.0, use the six module directly if it is still needed [astropy.extern.six]
In this thread we will use a pregenerated data set which is contained in a tar archive in the data directory of the fermipy-extra repository.
if os.path.isfile('../data/draco.tar.gz'):
!tar xzf ../data/draco.tar.gz
else:
!curl -OL https://raw.githubusercontent.com/fermiPy/fermipy-extras/master/data/draco.tar.gz
!tar xzf draco.tar.gz
We will begin by looking at the contents of the configuration file. The configuration is similar to our PG1553 example except for the addition of a 'draco' component in the ROI model. We also set the ROI coordinates explicitly since the ROI center isn't at the position of a 3FGL source.
!cat draco/config.yaml
logging : verbosity : 3 data: evfile: draco_ft1.fits scfile: null ltcube : ltcube_239557414_428903014_z100_r180_gti.fits binning: # Binning roiwidth : 10.0 npix : null binsz : 0.1 # spatial bin size in deg binsperdec : 8 # nb energy bins per decade coordsys : 'GAL' selection: # Data selections emin : 500 emax : 500000 zmax : 100 evclass : 128 evtype : 3 tmin : 239557414 tmax : 428903014 # 6 years filter : null # Set the ROI center to these coordinates ra: 260.05167 dec: 57.91528 gtlike: # IRFs edisp : True irfs : 'P8R3_SOURCE_V3' # Settings for ROI model model: # Include catalog sources within this distance from the ROI center src_radius : null # Include catalog sources within a box of width roisrc. src_roiwidth : 15.0 galdiff : '$FERMI_DIFFUSE_DIR/gll_iem_v07.fits' isodiff : '$FERMI_DIFFUSE_DIR/iso_P8R3_SOURCE_V3_v1.txt' # List of catalogs to be used in the model. catalogs : - '4FGL' sources : - { name : 'draco', ra : 260.05167, dec : 57.91528, SpectrumType : 'PowerLaw', Index : 2.0, Prefactor : { 'value' : 0.0, 'scale' : 1e-13 } } components: null
Note that the setup for a joint analysis is identical to the above except for the modification to the components section. The following example shows the components configuration one would use to define a joint analysis with the four PSF event types:
components:
- { selection : { evtype : 4 } }
- { selection : { evtype : 8 } }
- { selection : { evtype : 16 } }
- { selection : { evtype : 32 } }
To get started we will first instantiate a GTAnalysis instance using the config file in the draco directory and the run the setup() method. This will prepare all the ancillary files and create the pylikelihood instance for binned analysis. Note that in this example these files have already been generated so the routines that would normally be executed to create these files will be skipped.
gta = GTAnalysis('draco/config.yaml')
matplotlib.interactive(True)
gta.setup()
gta.write_roi('fit0')
2021-06-18 15:57:57 INFO GTAnalysis.__init__(): -------------------------------------------------------------------------------- fermipy version v1.0.1 ScienceTools version 2.0.18 2021-06-18 15:57:58 INFO GTAnalysis.setup(): Running setup. 2021-06-18 15:57:58 INFO GTBinnedAnalysis.setup(): Running setup for component 00 2021-06-18 15:57:58 INFO GTBinnedAnalysis._select_data(): Skipping data selection. 2021-06-18 15:57:58 INFO GTBinnedAnalysis.setup(): Using external LT cube. /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/fermipy/irfs.py:51: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result. log_ratio = np.log(x[xs1] / x[xs0]) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/fermipy/irfs.py:52: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result. return 0.5 * (y[ys0] * x[xs0] + y[ys1] * x[xs1]) * log_ratio 2021-06-18 15:57:58 INFO GTBinnedAnalysis._create_expcube(): Skipping gtexpcube. 2021-06-18 15:57:58 INFO GTBinnedAnalysis._create_srcmaps(): Skipping gtsrcmaps. 2021-06-18 15:57:58 INFO GTBinnedAnalysis.setup(): Finished setup for component 00 2021-06-18 15:57:58 INFO GTBinnedAnalysis._create_binned_analysis(): Creating BinnedAnalysis for component 00. 2021-06-18 15:58:07 INFO GTAnalysis.setup(): Initializing source properties 2021-06-18 15:58:12 INFO GTAnalysis.setup(): Finished setup. 2021-06-18 15:58:12 INFO GTBinnedAnalysis.write_xml(): Writing /Users/manuelmeyer/Python/fermipy-extra/notebooks/draco/fit0_00.xml... 2021-06-18 15:58:12 INFO GTAnalysis.write_fits(): Writing /Users/manuelmeyer/Python/fermipy-extra/notebooks/draco/fit0.fits... WARNING: Format %s cannot be mapped to the accepted TDISPn keyword values. Format will not be moved into TDISPn keyword. [astropy.io.fits.column] WARNING: Format %f cannot be mapped to the accepted TDISPn keyword values. Format will not be moved into TDISPn keyword. [astropy.io.fits.column] 2021-06-18 15:58:22 INFO GTAnalysis.write_roi(): Writing /Users/manuelmeyer/Python/fermipy-extra/notebooks/draco/fit0.npy...
We can print the ROI object to see a list of sources in the model along with their distance from the ROI center (offset), TS, and number of predicted counts (Npred). Since we haven't yet fit any sources, the ts of all sources will initially be assigned as nan.
gta.print_roi()
2021-06-18 15:58:22 INFO GTAnalysis.print_roi(): name SpatialModel SpectrumType offset ts npred -------------------------------------------------------------------------------- draco PointSource PowerLaw 0.000 nan 0.0 4FGL J1725.5+5851 PointSource PowerLaw 1.174 nan 352.2 4FGL J1741.2+5739 PointSource PowerLaw 2.821 nan 103.5 4FGL J1722.6+6104 PointSource PowerLaw 3.173 nan 220.8 4FGL J1742.5+5944 PointSource PowerLaw 3.415 nan 268.3 4FGL J1657.0+6010 PointSource PowerLaw 3.734 nan 189.2 4FGL J1705.4+5436 PointSource PowerLaw 3.894 nan 124.1 4FGL J1733.4+5428 PointSource PowerLaw 3.899 nan 98.9 4FGL J1658.4+6150 PointSource PowerLaw 4.782 nan 101.4 4FGL J1740.6+5346 PointSource PowerLaw 5.037 nan 136.6 4FGL J1647.0+6040 PointSource PowerLaw 5.044 nan 67.8 4FGL J1756.3+5522 PointSource PowerLaw 5.568 nan 119.7 4FGL J1638.1+5721 PointSource PowerLaw 5.656 nan 32.5 4FGL J1647.1+6149 PointSource PowerLaw 5.695 nan 55.8 4FGL J1632.4+5800 PointSource PowerLaw 6.327 nan 1.6 4FGL J1740.5+5211 PointSource PowerLaw 6.417 nan 110.1 4FGL J1649.4+5235 PointSource PowerLaw 6.888 nan 256.6 4FGL J1645.6+6329 PointSource PowerLaw 6.990 nan 17.8 4FGL J1815.1+5917 PointSource LogParabola 7.280 nan 4.9 4FGL J1623.6+5743 PointSource PowerLaw 7.515 nan 1.2 4FGL J1806.3+5345 PointSource PowerLaw 7.680 nan 4.7 4FGL J1626.5+6257 PointSource PowerLaw 8.292 nan 1.4 4FGL J1810.7+5335 PointSource PowerLaw 8.297 nan 1.9 4FGL J1807.2+6429 PointSource PowerLaw 8.657 nan 1.3 4FGL J1647.5+4950 PointSource PowerLaw 9.390 nan 4.4 isodiff ConstantValue FileFunction ----- nan 10541.7 galdiff MapCubeFunctio PowerLaw ----- nan 16069.6
We can assess the quality of our pre-fit model by running the residmap method. This will generate four maps
resid = gta.residmap('draco_prefit',model={'SpatialModel' : 'PointSource', 'Index' : 2.0})
fig = plt.figure(figsize=(14,6))
ROIPlotter(resid['data'],roi=gta.roi).plot(vmin=400,vmax=1000,subplot=121,cmap='magma')
plt.gca().set_title('Data')
ROIPlotter(resid['model'],roi=gta.roi).plot(vmin=400,vmax=1000,subplot=122,cmap='magma')
plt.gca().set_title('Model')
fig = plt.figure(figsize=(14,6))
ROIPlotter(resid['sigma'],roi=gta.roi).plot(vmin=-5,vmax=5,levels=[-5,-3,3,5],subplot=121,cmap='RdBu_r')
plt.gca().set_title('Significance')
ROIPlotter(resid['excess'],roi=gta.roi).plot(vmin=-100,vmax=100,subplot=122,cmap='RdBu_r')
plt.gca().set_title('Excess')
2021-06-18 15:58:22 INFO GTAnalysis.residmap(): Generating residual maps 2021-06-18 15:58:22 INFO GTAnalysis.add_source(): Adding source residmap_testsource 2021-06-18 15:58:24 INFO GTAnalysis.delete_source(): Deleting source residmap_testsource 2021-06-18 15:58:24 INFO GTAnalysis.residmap(): Finished residual maps 2021-06-18 15:58:30 WARNING GTAnalysis.residmap(): Saving maps in .npy files is disabled b/c of incompatibilities in python3, remove the maps from the /Users/manuelmeyer/Python/fermipy-extra/notebooks/draco/draco_prefit_pointsource_powerlaw_2.00_residmap.npy 2021-06-18 15:58:30 INFO GTAnalysis.residmap(): Execution time: 7.59 s /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/fermipy/plotting.py:146: MatplotlibDeprecationWarning: Passing raw data via parameters data and lut to register_cmap() is deprecated since 3.3 and will become an error two minor releases later. Instead use: register_cmap(cmap=LinearSegmentedColormap(name, data, lut)) plt.register_cmap(name='ds9_b', data=ds9_b) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/fermipy/plotting.py:301: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. In future versions, 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 = copy.copy(mpl.cm.get_cmap("magma")) colormap.set_under(colormap(0)) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:202: MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it. return super().imshow(X, *args, origin=origin, **kwargs) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/fermipy/plotting.py:301: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. In future versions, 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 = copy.copy(mpl.cm.get_cmap("RdBu_r")) colormap.set_under(colormap(0))
Text(0.5, 1.0, 'Excess')
/Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe)
/Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe)
Now we will run the optimize method. This method will iteratively optimize the parameters of all components in the ROI in several stages:
Running optimize gives us a baseline model that we can use as a starting point for subsequent stages of the analysis. We will also save the results of the analysis with write_roi. By saving the analysis state we can restore the analysis to this point at any time with the load_roi method. (NOTE: This step is computationally intensive and can take up to 5-10 minutes)
gta.optimize()
gta.write_roi('fit1')
2021-06-18 15:58:33 INFO GTAnalysis.optimize(): Starting
Joint fit ['galdiff', 'isodiff', '4FGL J1725.5+5851', '4FGL J1742.5+5944', '4FGL J1649.4+5235']
/Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/scipy/interpolate/fitpack2.py:253: UserWarning: The maximal number of iterations maxit (set to 20 by the program) allowed for finding a smoothing spline with fp=s has been reached: s too small. There is an approximation returned but the corresponding weighted sum of squared residuals does not satisfy the condition abs(fp-s)/s < tol. warnings.warn(message)
Fitting shape galdiff TS: 585.797 Fitting shape isodiff TS: 544.666 Fitting shape 4FGL J1649.4+5235 TS: 413.410 Fitting shape 4FGL J1725.5+5851 TS: 310.493 Fitting shape 4FGL J1742.5+5944 TS: 190.030 Fitting shape 4FGL J1756.3+5522 TS: 103.030 Fitting shape 4FGL J1740.6+5346 TS: 77.101 Fitting shape 4FGL J1657.0+6010 TS: 67.301 Fitting shape 4FGL J1705.4+5436 TS: 52.557 Fitting shape 4FGL J1722.6+6104 TS: 49.037 Fitting shape 4FGL J1658.4+6150 TS: 43.143 Fitting shape 4FGL J1733.4+5428 TS: 25.890
2021-06-18 15:58:39 INFO GTAnalysis.optimize(): Finished 2021-06-18 15:58:39 INFO GTAnalysis.optimize(): LogLike: nan Delta-LogLike: nan 2021-06-18 15:58:39 INFO GTAnalysis.optimize(): Execution time: 6.45 s 2021-06-18 15:58:39 INFO GTBinnedAnalysis.write_xml(): Writing /Users/manuelmeyer/Python/fermipy-extra/notebooks/draco/fit1_00.xml... 2021-06-18 15:58:39 INFO GTAnalysis.write_fits(): Writing /Users/manuelmeyer/Python/fermipy-extra/notebooks/draco/fit1.fits...
Fitting shape 4FGL J1647.1+6149 TS: 25.573
WARNING: Format %s cannot be mapped to the accepted TDISPn keyword values. Format will not be moved into TDISPn keyword. [astropy.io.fits.column] WARNING: Format %f cannot be mapped to the accepted TDISPn keyword values. Format will not be moved into TDISPn keyword. [astropy.io.fits.column] 2021-06-18 15:58:49 INFO GTAnalysis.write_roi(): Writing /Users/manuelmeyer/Python/fermipy-extra/notebooks/draco/fit1.npy...
After running optimize we can rerun print_roi to see a summary of the updated model. All sources that were fit in this step now have ts values and an Npred value the reflects the optimized normalization of that source. Note that model components that were not fit during the optimize step still have ts=nan.
gta.print_roi()
2021-06-18 15:58:49 INFO GTAnalysis.print_roi(): name SpatialModel SpectrumType offset ts npred -------------------------------------------------------------------------------- draco PointSource PowerLaw 0.000 nan 0.0 4FGL J1725.5+5851 PointSource PowerLaw 1.174 313.90 387.3 4FGL J1741.2+5739 PointSource PowerLaw 2.821 16.34 80.3 4FGL J1722.6+6104 PointSource PowerLaw 3.173 51.56 211.4 4FGL J1742.5+5944 PointSource PowerLaw 3.415 189.03 284.4 4FGL J1657.0+6010 PointSource PowerLaw 3.734 69.52 196.3 4FGL J1705.4+5436 PointSource PowerLaw 3.894 52.66 81.2 4FGL J1733.4+5428 PointSource PowerLaw 3.899 27.34 154.6 4FGL J1658.4+6150 PointSource PowerLaw 4.782 44.35 129.8 4FGL J1740.6+5346 PointSource PowerLaw 5.037 74.47 141.1 4FGL J1647.0+6040 PointSource PowerLaw 5.044 4.78 31.6 4FGL J1756.3+5522 PointSource PowerLaw 5.568 99.93 113.2 4FGL J1638.1+5721 PointSource PowerLaw 5.656 -0.00 0.0 4FGL J1647.1+6149 PointSource PowerLaw 5.695 26.53 81.3 4FGL J1632.4+5800 PointSource PowerLaw 6.327 -0.00 0.0 4FGL J1740.5+5211 PointSource PowerLaw 6.417 6.80 101.6 4FGL J1649.4+5235 PointSource PowerLaw 6.888 414.99 273.3 4FGL J1645.6+6329 PointSource PowerLaw 6.990 0.00 2.1 4FGL J1815.1+5917 PointSource LogParabola 7.280 1.00 53.7 4FGL J1623.6+5743 PointSource PowerLaw 7.515 -0.00 0.0 4FGL J1806.3+5345 PointSource PowerLaw 7.680 0.53 30.3 4FGL J1626.5+6257 PointSource PowerLaw 8.292 -0.00 0.0 4FGL J1810.7+5335 PointSource PowerLaw 8.297 -0.00 0.0 4FGL J1807.2+6429 PointSource PowerLaw 8.657 0.01 3.8 4FGL J1647.5+4950 PointSource PowerLaw 9.390 -0.00 0.3 isodiff ConstantValue FileFunction ----- 10037.09 14832.3 galdiff MapCubeFunctio PowerLaw ----- 13181.25 16586.4
gta.print_params()
2021-06-18 15:58:49 INFO GTAnalysis.print_params(): idx parname value error min max scale free -------------------------------------------------------------------------------- draco 73 Prefactor 0 0 0 1e+03 1e-13 * 74 Index 2 0 0 5 -1 75 Scale 1 0 0.001 1e+03 1e+03
To evaluate the quality of the optimized model we can rerun the residmap method. In the updated residual map the positive residuals from before are strongly reduced.
resid = gta.residmap('draco_postfit',model={'SpatialModel' : 'PointSource', 'Index' : 2.0})
fig = plt.figure(figsize=(14,6))
ROIPlotter(resid['sigma'],roi=gta.roi).plot(vmin=-5,vmax=5,levels=[-5,-3,3,5],subplot=121,cmap='RdBu_r')
plt.gca().set_title('Significance')
ROIPlotter(resid['excess'],roi=gta.roi).plot(vmin=-100,vmax=100,subplot=122,cmap='RdBu_r')
plt.gca().set_title('Excess')
2021-06-18 15:58:49 INFO GTAnalysis.residmap(): Generating residual maps 2021-06-18 15:58:49 INFO GTAnalysis.add_source(): Adding source residmap_testsource 2021-06-18 15:58:50 INFO GTAnalysis.delete_source(): Deleting source residmap_testsource 2021-06-18 15:58:51 INFO GTAnalysis.residmap(): Finished residual maps 2021-06-18 15:58:57 WARNING GTAnalysis.residmap(): Saving maps in .npy files is disabled b/c of incompatibilities in python3, remove the maps from the /Users/manuelmeyer/Python/fermipy-extra/notebooks/draco/draco_postfit_pointsource_powerlaw_2.00_residmap.npy 2021-06-18 15:58:57 INFO GTAnalysis.residmap(): Execution time: 7.67 s /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/fermipy/plotting.py:146: MatplotlibDeprecationWarning: Passing raw data via parameters data and lut to register_cmap() is deprecated since 3.3 and will become an error two minor releases later. Instead use: register_cmap(cmap=LinearSegmentedColormap(name, data, lut)) plt.register_cmap(name='ds9_b', data=ds9_b) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/fermipy/plotting.py:301: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. In future versions, 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 = copy.copy(mpl.cm.get_cmap("RdBu_r")) colormap.set_under(colormap(0)) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:202: MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it. return super().imshow(X, *args, origin=origin, **kwargs)
Text(0.5, 1.0, 'Excess')
/Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe)
Another diagnostic for the quality of the ROI model is the TS map. The tsmap method can be used to generate a TS map of the ROI with a given test source model. Here we use the same source model we did for the residual map -- a point source with a power-law index of 2.
tsmap_postfit = gta.tsmap('draco_postfit',model={'SpatialModel' : 'PointSource', 'Index' : 2.0})
2021-06-18 15:58:58 INFO GTAnalysis.tsmap(): Generating TS map 2021-06-18 15:58:58 INFO GTAnalysis._make_tsmap_fast(): Fitting test source. 2021-06-18 15:59:08 INFO GTAnalysis.tsmap(): Finished TS map 2021-06-18 15:59:16 WARNING GTAnalysis.tsmap(): Saving TS maps in .npy files is disabled b/c of incompatibilities in python3, remove the maps from the /Users/manuelmeyer/Python/fermipy-extra/notebooks/draco/draco_postfit_pointsource_powerlaw_2.00_tsmap.npy 2021-06-18 15:59:16 INFO GTAnalysis.tsmap(): Execution time: 17.81 s
Here we see that the same excess in the south west part of the ROI that was also visible in the residual map. This excess is detected as a new point source with TS close to 25:
tsmap_postfit['sqrt_ts'].data.max()**2.
21.66518724742127
fig = plt.figure(figsize=(14,6))
ROIPlotter(tsmap_postfit['sqrt_ts'],roi=gta.roi).plot(levels=[0,3,5,7],vmin=0,vmax=5,subplot=121,cmap='magma')
plt.gca().set_title('Sqrt(TS)')
ROIPlotter(tsmap_postfit['npred'],roi=gta.roi).plot(vmin=0,vmax=100,subplot=122,cmap='magma')
plt.gca().set_title('NPred')
/Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/fermipy/plotting.py:146: MatplotlibDeprecationWarning: Passing raw data via parameters data and lut to register_cmap() is deprecated since 3.3 and will become an error two minor releases later. Instead use: register_cmap(cmap=LinearSegmentedColormap(name, data, lut)) plt.register_cmap(name='ds9_b', data=ds9_b) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/fermipy/plotting.py:301: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. In future versions, 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 = copy.copy(mpl.cm.get_cmap("magma")) colormap.set_under(colormap(0)) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:202: MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it. return super().imshow(X, *args, origin=origin, **kwargs)
Text(0.5, 1.0, 'NPred')
/Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe)
We can add this source into the model by running the find_sources
method. This method will generate a TS map of the region and add a point-source at the location of each peak with TS > 16.
src = gta.find_sources(sqrt_ts_threshold=4.0)
2021-06-18 15:59:18 INFO GTAnalysis.find_sources(): Starting. 2021-06-18 15:59:18 INFO GTAnalysis.tsmap(): Generating TS map 2021-06-18 15:59:18 INFO GTAnalysis._make_tsmap_fast(): Fitting test source. 2021-06-18 15:59:27 INFO GTAnalysis.tsmap(): Finished TS map 2021-06-18 15:59:34 WARNING GTAnalysis.tsmap(): Saving TS maps in .npy files is disabled b/c of incompatibilities in python3, remove the maps from the /Users/manuelmeyer/Python/fermipy-extra/notebooks/draco/sourcefind_00_pointsource_powerlaw_2.00_tsmap.npy 2021-06-18 15:59:34 INFO GTAnalysis.tsmap(): Execution time: 16.64 s 2021-06-18 15:59:34 INFO GTAnalysis._build_src_dicts_from_peaks(): Found source name: PS J1726.8+5954 ts: 21.665187 2021-06-18 15:59:34 INFO GTAnalysis.add_source(): Adding source PS J1726.8+5954 2021-06-18 15:59:36 INFO GTAnalysis.free_source(): Fixing parameters for PS J1726.8+5954 : ['Prefactor'] 2021-06-18 15:59:36 INFO GTAnalysis._find_sources_iterate(): Performing spectral fit for PS J1726.8+5954. 2021-06-18 15:59:36 INFO GTAnalysis.free_source(): Freeing parameters for PS J1726.8+5954 : ['Prefactor', 'Index'] 2021-06-18 15:59:36 INFO GTAnalysis.fit(): Starting fit. 2021-06-18 15:59:36 INFO GTAnalysis.fit(): Fit returned successfully. Quality: 3 Status: 0 2021-06-18 15:59:36 INFO GTAnalysis.fit(): LogLike: -60755.226 DeltaLogLike: nan 2021-06-18 15:59:36 INFO GTAnalysis._find_sources_iterate(): {'Index': {'error': 0.2631627116784879, 'value': -2.008301461402123}, 'Prefactor': {'error': 5.3406068480684565e-14, 'value': 9.841199185890776e-14}, 'Scale': {'error': nan, 'value': 1000.0}} 2021-06-18 15:59:36 INFO GTAnalysis.free_source(): Fixing parameters for PS J1726.8+5954 : ['Prefactor', 'Index'] 2021-06-18 15:59:36 INFO GTAnalysis.find_sources(): Found 1 sources in iteration 0. 2021-06-18 15:59:36 INFO GTAnalysis.tsmap(): Generating TS map 2021-06-18 15:59:37 INFO GTAnalysis._make_tsmap_fast(): Fitting test source. 2021-06-18 15:59:47 INFO GTAnalysis.tsmap(): Finished TS map 2021-06-18 15:59:55 WARNING GTAnalysis.tsmap(): Saving TS maps in .npy files is disabled b/c of incompatibilities in python3, remove the maps from the /Users/manuelmeyer/Python/fermipy-extra/notebooks/draco/sourcefind_01_pointsource_powerlaw_2.00_tsmap.npy 2021-06-18 15:59:55 INFO GTAnalysis.tsmap(): Execution time: 18.55 s 2021-06-18 15:59:55 INFO GTAnalysis.find_sources(): Found 0 sources in iteration 1. 2021-06-18 15:59:55 INFO GTAnalysis.find_sources(): Done. 2021-06-18 15:59:55 INFO GTAnalysis.find_sources(): Execution time: 37.15 s
From the log we can see that a new source PS J1726.8+5954 was found. Rerunning the tsmap
method we can see that this source is located at the peak of the TS Map that was previously present in the northeast corner of the ROI. Some residual emission is still present, though.
print(gta.roi['PS J1726.8+5954'])
tsmap_newsrcs = gta.tsmap('draco_newsrcs',model={'SpatialModel' : 'PointSource', 'Index' : 2.0})
fig = plt.figure(figsize=(14,6))
ROIPlotter(tsmap_newsrcs['sqrt_ts'],roi=gta.roi).plot(levels=[0,3,5,7],vmin=0,vmax=5,subplot=121,cmap='magma')
plt.gca().set_title('Sqrt(TS)')
ROIPlotter(tsmap_newsrcs['npred'],roi=gta.roi).plot(vmin=0,vmax=100,subplot=122,cmap='magma')
plt.gca().set_title('NPred')
2021-06-18 15:59:55 INFO GTAnalysis.tsmap(): Generating TS map
Name : PS J1726.8+5954 Associations : ['PS J1726.8+5954'] RA/DEC : 261.718/ 59.905 GLON/GLAT : 88.707/ 33.743 TS : 21.63 Npred : 58.27 Flux : 1.962e-10 +/- 9.51e-11 EnergyFlux : 6.645e-07 +/- 2.96e-07 SpatialModel : PointSource SpectrumType : PowerLaw Spectral Parameters b'Prefactor' : 9.841e-14 +/- 5.341e-14 b'Index' : -2.008 +/- 0.2632 b'Scale' : 1000 +/- nan
2021-06-18 15:59:56 INFO GTAnalysis._make_tsmap_fast(): Fitting test source. 2021-06-18 16:00:05 INFO GTAnalysis.tsmap(): Finished TS map 2021-06-18 16:00:13 WARNING GTAnalysis.tsmap(): Saving TS maps in .npy files is disabled b/c of incompatibilities in python3, remove the maps from the /Users/manuelmeyer/Python/fermipy-extra/notebooks/draco/draco_newsrcs_pointsource_powerlaw_2.00_tsmap.npy 2021-06-18 16:00:13 INFO GTAnalysis.tsmap(): Execution time: 18.01 s /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/fermipy/plotting.py:146: MatplotlibDeprecationWarning: Passing raw data via parameters data and lut to register_cmap() is deprecated since 3.3 and will become an error two minor releases later. Instead use: register_cmap(cmap=LinearSegmentedColormap(name, data, lut)) plt.register_cmap(name='ds9_b', data=ds9_b) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/fermipy/plotting.py:301: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. In future versions, 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 = copy.copy(mpl.cm.get_cmap("magma")) colormap.set_under(colormap(0)) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:202: MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it. return super().imshow(X, *args, origin=origin, **kwargs)
Text(0.5, 1.0, 'NPred')
/Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe) /Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/astropy/visualization/wcsaxes/core.py:458: MatplotlibDeprecationWarning: The 'inframe' parameter of draw() was deprecated in Matplotlib 3.3 and will be removed two minor releases later. Use Axes.redraw_in_frame() instead. If any parameter follows 'inframe', they should be passed as keyword, not positionally. super().draw(renderer, inframe=inframe)
After optimizing the ROI model we are ready to perform our analysis of the source of interest (draco). We will begin by freeing draco along with all other point sources within 3 deg of the ROI center and refitting their normalizations.
gta.free_sources(distance=3.0,pars='norm')
gta.free_sources(distance=3.0,pars='shape',minmax_ts=[100.,None])
fit_results = gta.fit()
2021-06-18 16:00:15 INFO GTAnalysis.free_source(): Freeing parameters for 4FGL J1725.5+5851 : ['Prefactor'] 2021-06-18 16:00:15 INFO GTAnalysis.free_source(): Freeing parameters for PS J1726.8+5954 : ['Prefactor'] 2021-06-18 16:00:15 INFO GTAnalysis.free_source(): Freeing parameters for 4FGL J1741.2+5739 : ['Prefactor'] 2021-06-18 16:00:15 INFO GTAnalysis.free_source(): Freeing parameters for isodiff : ['Normalization'] 2021-06-18 16:00:15 INFO GTAnalysis.free_source(): Freeing parameters for galdiff : ['Prefactor'] 2021-06-18 16:00:15 INFO GTAnalysis.free_source(): Freeing parameters for 4FGL J1725.5+5851 : ['Index'] 2021-06-18 16:00:15 INFO GTAnalysis.free_source(): Freeing parameters for galdiff : ['Index'] 2021-06-18 16:00:15 INFO GTAnalysis.fit(): Starting fit. 2021-06-18 16:00:16 INFO GTAnalysis.fit(): Fit returned successfully. Quality: 3 Status: 0 2021-06-18 16:00:16 INFO GTAnalysis.fit(): LogLike: -60754.928 DeltaLogLike: 0.297
gta.print_roi()
2021-06-18 16:00:16 INFO GTAnalysis.print_roi(): name SpatialModel SpectrumType offset ts npred -------------------------------------------------------------------------------- draco PointSource PowerLaw 0.000 -0.00 0.0 4FGL J1725.5+5851 PointSource PowerLaw 1.174 304.12 379.3 PS J1726.8+5954 PointSource PowerLaw 2.168 21.97 59.7 4FGL J1741.2+5739 PointSource PowerLaw 2.821 15.97 80.8 4FGL J1722.6+6104 PointSource PowerLaw 3.173 51.56 211.4 4FGL J1742.5+5944 PointSource PowerLaw 3.415 189.03 284.4 4FGL J1657.0+6010 PointSource PowerLaw 3.734 69.52 196.3 4FGL J1705.4+5436 PointSource PowerLaw 3.894 52.66 81.2 4FGL J1733.4+5428 PointSource PowerLaw 3.899 27.34 154.6 4FGL J1658.4+6150 PointSource PowerLaw 4.782 44.35 129.8 4FGL J1740.6+5346 PointSource PowerLaw 5.037 74.47 141.1 4FGL J1647.0+6040 PointSource PowerLaw 5.044 4.78 31.6 4FGL J1756.3+5522 PointSource PowerLaw 5.568 99.93 113.2 4FGL J1638.1+5721 PointSource PowerLaw 5.656 -0.00 0.0 4FGL J1647.1+6149 PointSource PowerLaw 5.695 26.53 81.3 4FGL J1632.4+5800 PointSource PowerLaw 6.327 -0.00 0.0 4FGL J1740.5+5211 PointSource PowerLaw 6.417 6.80 101.6 4FGL J1649.4+5235 PointSource PowerLaw 6.888 414.99 273.3 4FGL J1645.6+6329 PointSource PowerLaw 6.990 0.00 2.1 4FGL J1815.1+5917 PointSource LogParabola 7.280 1.00 53.7 4FGL J1623.6+5743 PointSource PowerLaw 7.515 -0.00 0.0 4FGL J1806.3+5345 PointSource PowerLaw 7.680 0.53 30.3 4FGL J1626.5+6257 PointSource PowerLaw 8.292 -0.00 0.0 4FGL J1810.7+5335 PointSource PowerLaw 8.297 -0.00 0.0 4FGL J1807.2+6429 PointSource PowerLaw 8.657 0.01 3.8 4FGL J1647.5+4950 PointSource PowerLaw 9.390 -0.00 0.3 isodiff ConstantValue FileFunction ----- 399.13 14913.9 galdiff MapCubeFunctio PowerLaw ----- 452.96 16385.0
After running the fit completes we can execute the spectral analysis of draco with the sed method.
sed_draco = gta.sed('draco')
gta.write_roi('fit_sed')
2021-06-18 16:00:16 INFO GTAnalysis.sed(): Computing SED for draco 2021-06-18 16:00:17 INFO GTAnalysis._make_sed(): Fitting SED 2021-06-18 16:00:17 INFO GTAnalysis.free_source(): Fixing parameters for 4FGL J1725.5+5851 : ['Index'] 2021-06-18 16:00:17 INFO GTAnalysis.free_source(): Fixing parameters for galdiff : ['Index'] 2021-06-18 16:00:22 INFO GTAnalysis.sed(): Finished SED 2021-06-18 16:00:29 INFO GTAnalysis.sed(): Execution time: 12.10 s 2021-06-18 16:00:29 INFO GTBinnedAnalysis.write_xml(): Writing /Users/manuelmeyer/Python/fermipy-extra/notebooks/draco/fit_sed_00.xml... 2021-06-18 16:00:29 INFO GTAnalysis.write_fits(): Writing /Users/manuelmeyer/Python/fermipy-extra/notebooks/draco/fit_sed.fits... WARNING: Format %s cannot be mapped to the accepted TDISPn keyword values. Format will not be moved into TDISPn keyword. [astropy.io.fits.column] WARNING: Format %f cannot be mapped to the accepted TDISPn keyword values. Format will not be moved into TDISPn keyword. [astropy.io.fits.column] 2021-06-18 16:00:39 INFO GTAnalysis.write_roi(): Writing /Users/manuelmeyer/Python/fermipy-extra/notebooks/draco/fit_sed.npy...
We can now inspect the fit results by looking at the elements of the output dictionary. By default the sed
method will perform a likelihood scan in each energy bin which is saved in the dloglike_scan
array. In the following example we plot the likelihood profile in the first energy bin and overplot the flux upper limit in that bin (vertical black line). fermiPy uses the delta-log-likelihood method to evaluate ULs and we can see that the 95% CL flux upper limit intersects with the point at which the log-likelihood has decreased by 2.71/2 from its maximum value (horizontal red line).
# E^2 x Differential flux ULs in each bin in units of MeV cm^{-2} s^{-1}
print(sed_draco['e2dnde_ul95'])
e2dnde_scan = sed_draco['norm_scan']*sed_draco['ref_e2dnde'][:,None]
plt.figure()
plt.plot(e2dnde_scan[0],
sed_draco['dloglike_scan'][0]-np.max(sed_draco['dloglike_scan'][0]))
#plt.gca().set_ylim(-5,1)
plt.gca().axvline(sed_draco['e2dnde_ul95'][0],color='k')
plt.gca().axhline(-2.71/2.,color='r')
[2.40285095e-25 7.40765731e-26 3.41231872e-25 2.19682560e-25 3.74476648e-25 6.95933941e-26 1.03795947e-25 1.50025889e-25 1.15250358e-25 1.40418183e-25 1.51014898e-25 1.39455554e-24 1.84020521e-25 2.04472491e-25 2.07466174e-25 2.25251861e-25 2.19321609e-25 2.24051095e-25 2.27586992e-25 2.30241745e-25 2.32198521e-25 2.33659895e-25 2.34711785e-25 2.35533187e-25]
<matplotlib.lines.Line2D at 0x16d88a450>
We can also visualize the results of the scan with the SEDPlotter class. This class accepts a source object as its argument and creates a visualization of the SED as a sequence of points with errors. Setting showlnl=True
overplots the likelihood function in each bin as a color gradient (the so-called castro plot).
fig = plt.figure(figsize=(6,4))
ylim=[1E-8,1E-5]
fig.add_subplot(111)
SEDPlotter(sed_draco).plot()
plt.gca().set_ylim(ylim)
fig = plt.figure(figsize=(6,4))
fig.add_subplot(111)
SEDPlotter(sed_draco).plot(showlnl=True,ylim=ylim)
plt.gca().set_ylim(ylim)
/Users/manuelmeyer/anaconda3/envs/fermipy-1.0.1-test/lib/python3.7/site-packages/fermipy/plotting.py:719: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3. Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading']. This will become an error two minor releases later. linewidth=0)
(1e-08, 1e-05)
Now that we have run a spectral analysis we can use the bin-by-bin likelihoods to gamma-ray flux from DM annihilations in Draco. In the following sample code we demonstrate how to calculate the UL on the DM cross section for a given DM spectral model.
import pyLikelihood
# Load the sed data structure
data = sed_draco
# Instantiate a DM Fit Function for a DM particle spectrum given the following parameters
# Mass = 100 GeV
# Cross-Section: 3 x 10^{-26} cm^{3} s^{-1}
# J-Factor: 10^19 GeV^2 cm^{-5}
# Channel: b-bbar
dmf = pyLikelihood.DMFitFunction()
dmf.readFunction(os.path.expandvars('$FERMIPY_ROOT/data/gammamc_dif.dat'))
dmf.setParam('norm',1E19)
dmf.setParam('sigmav',3E-26)
dmf.setParam('mass',100.0)
dmf.setParam('bratio',1.0)
dmf.setParam('channel0',4)
def integrate_eflux(fn,ebins,nstep=10):
"""Compute energy flux within a sequence of energy bins."""
loge = np.linspace(ebins[0],ebins[-1],100)
dfde = [fn(pyLikelihood.dArg(10**x)) for x in loge]
dfde = np.array(dfde)
x = ebins
dx = (x[1:] - x[:-1])
yedge = x[1:,np.newaxis] + np.linspace(0,1,nstep)[np.newaxis,:]*dx[:,np.newaxis]
dy = 10**yedge[:,1:]-10**yedge[:,:-1]
y = 0.5*(yedge[:,1:]+yedge[:,:-1])
eflux = np.interp(np.ravel(y),loge,dfde)
eflux = np.sum(eflux.reshape(y.shape)*10**y*dy,axis=1)
return eflux
class SEDLike(object):
def __init__(self,sed):
self._sed = sed
self._eflux_scan = sed['norm_scan']*sed['ref_eflux'][:,None]
def __call__(self,eflux):
lnl = np.zeros(eflux.shape)
for i, ectr in enumerate(self._sed['e_ctr']):
v = np.interp(eflux[i],
self._eflux_scan[i],
self._sed['dloglike_scan'][i])
lnl[i] += v
return np.sum(lnl,axis=0)
ebins = np.log10(np.array(list(data['e_min']) + list([data['e_max'][-1]])))
eflux = integrate_eflux(dmf,ebins)
sigmav = 3.E-26*np.logspace(-3.,1.,101)
eflux = eflux[:,np.newaxis]*np.logspace(-3,1,101)[np.newaxis,:]
slike = SEDLike(data)
lnl = slike(eflux)
lnl -= np.max(lnl)
# Plot log-likelihood profile
plt.figure()
plt.plot(sigmav,lnl)
plt.gca().set_xscale('log')
plt.gca().axhline(-2.71/2.,color='k')
plt.gca().set_ylim(-4,1)
plt.gca().set_xlabel('Cross Section [cm$^{3}$ s$^{-1}$]')
sigmav_ul = float(np.interp(2.71/2.,-lnl,sigmav))
print ('Sigma-V Upper Limit: ', sigmav_ul)
Sigma-V Upper Limit: nan