This tutorial shows how to configure the model and how to run Fermipy-LAT g-tools with Fermipy. Many parts of this tutorial are taken directly from the documentation page of Fermipy: fermipy.readthedocs. I suggest to visit it to find further informations.
The model can be configured with a configuration file or directly running Fermipy with your script.
The configuration file defines the data selection and analysis parameters. The configuration file has a hierarchical structure that groups parameters into dictionaries that are keyed to a section name (data, binning, etc.).
When creating a class instance, the configuration is initialized by passing either a configuration dictionary or configuration file path to the class constructor. Keyword arguments can be passed to the constructor to override configuration parameters in the input dictionary. In the following example the config dictionary defines values for the parameters emin and emax. By passing a dictionary for the selection keyword argument, the value of emax in the keyword argument (10000) overrides the value of emax in the input dictionary.
The configuration file can be created also with a yaml file. Below I report a sample of configuration applied for an analysis of the SMC:
The configuration file has the same structure as the configuration dictionary such that one can read/write configurations using the load/dump methods of the yaml module:
DATA
The data section defines the input data files for the analysis (FT1, FT2, and livetime cube). evfile and scfile can either be individual files or group of files. The optional ltcube option can be used to choose a pre-generated livetime cube. If ltcube is null a livetime cube will be generated at runtime with gtltcube.
Below an example of the data part:
The options for the data component are the following:
BINNING
Options in the binning section control the spatial and spectral binning of the data.
The roiwidth is the ROI width, npix specifies the number of pixels, binsz indicates the pixel size, binsperdec is the number of energy bins per decade and projtype is the type of projection and can be WCS and HEALPIX.
Other options are:
COMPONENTS
The components section can be used to define analysis configurations for independent subselections of the data. Each subselection will have its own binned likelihood instance that is combined in a global likelihood function for the ROI (implemented with the SummedLikelihood class in pyLikelihood). The components section is optional and when set to null (the default) only a single likelihood component will be created with the parameters of the root analysis configuration.
The component section is defined as a list of dictionaries where each element sets analysis parameters for a different subcomponent of the analysis. The component configurations follow the same structure and accept the same parameters as the root analysis configuration. Parameters not defined in a given element will default to the values set in the root analysis configuration.
The following example illustrates how to define a Front/Back analysis with two components. Files associated to each component will be given a suffix according to their order in the list (e.g. file_00.fits, file_01.fits, etc.).
The following example illustrates how to define an analysis divided in four components. Each of them has a different zmax and uses a different evtype and isotropic template.
EXTENSION
The options in extension control the default behavior of the extension method. For more information about using this method see the Extension Fitting page. The main options are the following:
There is a notebook dedicated to this analysis.
FILEIO
The fileio section collects options related to file bookkeeping. The outdir option sets the root directory of the analysis instance where all output files will be written. If outdir is null then the output directory will be automatically set to the directory in which the configuration file is located. Enabling the usescratch option will stage all output data files to a temporary scratch directory created under scratchdir.
Sample fileio Configuration
The options for the fileio are the following:
GTLIKE
Options in the gtlike section control the setup of the likelihood analysis include the IRF name (irfs).
A subsample of options for the gtlike section is listed below:
The Likelihood wights map is created and used in order to account for example for the uncertainty in the interstellar emission.
MODEL
The model section collects options that control the inclusion of point-source and diffuse components in the model. galdiff and isodiff set the templates for the Galactic IEM and isotropic diffuse respectively. catalogs defines a list of catalogs that will be merged to form a master analysis catalog from which sources will be drawn. Valid entries in this list can be FITS files or XML model files. sources can be used to insert additional point-source or extended components beyond those defined in the master catalog. src_radius and src_roiwidth set the maximum distance from the ROI center at which sources in the master catalog will be included in the ROI model.
Sample for the Model section:
In the sample above the official IEM and isotropic template are used together with the 3FGL catalog given with a fits file and by an additional catalog given with an xml file following the notation of the Fermi Science Tools fermi.fssc.spectralmodel. Then two sources named as SourceA and SourceB are added.
To define two or more galactic diffuse components you can optionally define the galdiff and isodiff parameters as lists. A separate component will be generated for each element in the list with the name galdiffXX or isodiffXX where XX is an integer position in the list.
To explicitly set the name of a component you can define any element as a dictionary containing name and file fields:
The list of sources for inclusion in the ROI model is set by defining a list of catalogs with the catalogs parameter. Catalog files can be in either XML or FITS format. Sources from the catalogs in this list that satisfy either the src_roiwidth or src_radius selections are added to the ROI model. If a source is defined in multiple catalogs the source definition from the last file in the catalogs list takes precedence.
Fermipy supports four spatial models which are defined with the SpatialModel property:
The spatial extension of RadialDisk and RadialGaussian can be controlled with the SpatialWidth parameter which sets the 68% containment radius in degrees. Note for ST releases prior to 11-01-01, RadialDisk and RadialGaussian sources will be represented with the SpatialMap type.
OPTIMIZER
This section provides informations for the setup used for the optimization. The list of options are the following:
PLOTTING
This section selects the options to be used to make the plots:
SELECTION
The selection section collects parameters related to the data selection and target definition. The majority of the parameters in this section are arguments to gtselect and gtmktime. The ROI center can be set with the target parameter by providing the name of a source defined in one of the input catalogs (defined in the model section). Alternatively the ROI center can be defined by giving explicit sky coordinates with ra and dec or glon and glat.
A sample of options for this section are the following:
First of all you need to load the configuration file, create the object gta and run the tool gta.setup that implements the ST gtselect, gtmktime, gtbin, gtexpcube, gtsrcmap tools
With the commands below we import the packages needed to run the script.
In [1]:
%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
from IPython.display import Image
Then we untar the file ../data/SMC_data.tar.gz. This will copy the config.yaml and ft1 file in the notebook directory.
In [2]:
if os.path.isfile('../data/SMC_data.tar.gz'):
!tar xzf ../data/SMC_data.tar.gz
else:
!curl -OL https://raw.githubusercontent.com/fermiPy/fermipy-extras/master/data/SMC_data.tar.gz
!tar xzf SMC_data.tar.gz
In [3]:
gta = GTAnalysis('config.yaml')
matplotlib.interactive(True)
gta.setup()
The setup() method performs the data preparation and response calculations needed for the analysis (selecting the data, creating counts and exposure maps, etc.). Depending on the data selection and binning of the analysis this will often be the slowest step in the analysis sequence. The output of setup() is cached in the analysis working directory so subsequent calls to setup() will run much faster.
The g-tools performed with gta.setup() are the following:
In order to have a summary of the source model you can use the command:
In [4]:
gta.print_model()
The table above contains:
First of all, let's make a fit using gta.fit tool. Source fitting with fermipy is generally performed with the optimize and fit methods. fit is a wrapper on the pyLikelihood fit method and performs a likelihood fit of all free parameters of the model. This method can be used to manually optimize of the model by calling it after freeing one or more source parameters.
In [5]:
gta.free_sources()
fitresult=gta.fit()
In [6]:
gta.print_model()
Now you see that the ts column is filled and the free column has all asterix because we have freed the sources.
Fermipy gives the possibility to access many informations on the sources in the model. Below a few examples.
NAME AND CHARACTERISTICS OF THE SOURCE
In [8]:
print gta.roi.sources[0]['name'] #name of the source
print gta.roi.sources[0]['Source_Name'] #name of the source
print gta.roi.sources[0]['SpatialModel'] #spatial model
print gta.roi.sources[0]['SpatialWidth'] #spatial size parameter
print gta.roi.sources[0]['SpatialType'] #spatial size parameter
print gta.roi.sources[0]['SourceType'] #Source type
print gta.roi.sources[0]['SpectrumType'] #Spectrum type string
print gta.roi.sources[0]['Spatial_Filename'] #Path to spatial template
print gta.roi.sources[0]['Spectrum_Filename'] #Path to the SED source template
print gta.roi.sources[0]['correlation'] #Dictionary of correlation coefficients.
print gta.roi.sources[0]['model_counts'] #Vector of predicted counts for this source
CHARACTERISTICS OF THE POSITION OF THE SOURCE
In [9]:
print gta.roi.sources[0]['ra'] #ra
print gta.roi.sources[0]['dec'] #dec
print gta.roi.sources[0]['glon'] #glon
print gta.roi.sources[0]['glat'] #glat
print gta.roi.sources[0]['ra_err'] #error for ra
print gta.roi.sources[0]['dec_err'] #error for dec
print gta.roi.sources[0]['glon_err'] #error for glon
print gta.roi.sources[0]['glat_err'] #error for glat
print gta.roi.sources[0]['pos_err'] #error for the position in deg
print gta.roi.sources[0]['pos_r68'] #68% CL error for the position
print gta.roi.sources[0]['pos_r95'] #95% CL error for the position
print gta.roi.sources[0]['pos_r99'] #99% CL error for the position
print gta.roi.sources[0]['pos_err_semimajor'] #1-sigma uncertainty (deg) along major axis of uncertainty ellipse.
print gta.roi.sources[0]['pos_err_semiminor'] #1-sigma uncertainty (deg) along minor axis of uncertainty ellipse.
print gta.roi.sources[0]['offset_ra'] #Right ascension offset from ROI center in local celestial projection (deg).
print gta.roi.sources[0]['offset_dec'] #Declination offset from ROI center in local celestial projection (deg).
print gta.roi.sources[0]['offset_glon'] #Galactic longitude offset from ROI center in local galactic projection (deg).
print gta.roi.sources[0]['offset_glat'] #Galactic latitude offset from ROI center in local galactic projection (deg).
print gta.roi.sources[0]['offset'] #Angular offset from ROI center (deg).
The error for the position is nan because we have not calculated yet the relocalization. Let's run the localization for a source with gta.localize tool.
In [10]:
gta.free_sources(free=True)
gta.print_model()
gta.free_sources(skydir=gta.roi[gta.roi.sources[2].name].skydir,distance=[3.0],free=True)
gta.print_model()
localsmc = gta.localize(gta.roi.sources[2].name, update=True, make_plots=True)
gta.print_model()
Using the option make_plots=True a few control plots are created like 3fgl_j0023.9-7203_localize_peak.png with the result for the likelihood analysis for the besy fit position and error for the position.
In [13]:
Image(filename='3fgl_j0023.9-7203_localize_peak.png')
Out[13]:
In [14]:
print gta.roi.sources[2]['ra'] #ra
print gta.roi.sources[2]['dec'] #dec
print gta.roi.sources[2]['glon'] #glon
print gta.roi.sources[2]['glat'] #glat
print gta.roi.sources[2]['ra_err'] #error for ra
print gta.roi.sources[2]['dec_err'] #error for dec
print gta.roi.sources[2]['glon_err'] #error for glon
print gta.roi.sources[2]['glat_err'] #error for glat
print gta.roi.sources[2]['pos_err'] #error for the position in deg
print gta.roi.sources[2]['pos_r68'] #68% CL error for the position
print gta.roi.sources[2]['pos_r95'] #95% CL error for the position
print gta.roi.sources[2]['pos_r99'] #99% CL error for the position
print gta.roi.sources[2]['pos_err_semimajor'] #1-sigma uncertainty (deg) along major axis of uncertainty ellipse.
print gta.roi.sources[2]['pos_err_semiminor'] #1-sigma uncertainty (deg) along minor axis of uncertainty ellipse.
print gta.roi.sources[2]['offset_ra'] #Right ascension offset from ROI center in local celestial projection (deg).
print gta.roi.sources[2]['offset_dec'] #Declination offset from ROI center in local celestial projection (deg).
print gta.roi.sources[2]['offset_glon'] #Galactic longitude offset from ROI center in local galactic projection (deg).
print gta.roi.sources[2]['offset_glat'] #Galactic latitude offset from ROI center in local galactic projection (deg).
print gta.roi.sources[2]['offset'] #Angular offset from ROI center (deg).
SED PARAMETERS FLUX AND SCAN OF FLUX
In [15]:
print gta.roi.sources[0]['param_names'] #Names of spectral parameters.
print gta.roi.sources[0]['param_values'] #Spectral parameter values.
print gta.roi.sources[0]['param_errors'] #Spectral parameters errors.
print gta.roi.sources[0]['ts'] #Source test statistic.
print gta.roi.sources[0]['loglike'] #Log-likelihood of the model evaluated at the best-fit normalization of the source.
print gta.roi.sources[0]['flux_scan'] #Flux values for scan of source normalization.
print gta.roi.sources[0]['norm_scan'] #Normalization parameters values for scan of source normalization.
print gta.roi.sources[0]['npred'] #Number of predicted counts from this source integrated over the analysis energy range.
print gta.roi.sources[0]['flux'] #Photon flux (cm−2 s−1) integrated over analysis energy range
print gta.roi.sources[0]['flux_err'] #Photon flux uncertainty (cm−2 s−1) integrated over analysis energy range
print gta.roi.sources[0]['flux_ul95'] #95% CL upper limit on the photon Differential photon flux (cm−2 s−1 MeV−1tegrated over analysis energy range
print gta.roi.sources[0]['dnde'] #Differential photon flux (cm−2 s−1 MeV−1) evaluated at the pivot energy.
In [16]:
gta.write_roi('model_test',make_plots=True,save_model_map=True)
As we saw before you can customize your model directly in the config.yaml file. However, sometimes you want to change something in your model directly in your Python script.
You can free or fix sources using gta.free_sources. First let's fix the SED parameter of all the sources.
In [17]:
gta.free_sources(free=False)
gta.print_model()
Now we free all parameters:
In [18]:
gta.free_sources(free=True)
gta.print_model()
Now we free all the SED parameters of sources within 3 degrees from 3FGL J0059.0-7242e:
In [19]:
gta.free_sources(free=False)
gta.free_sources(skydir=gta.roi['3FGL J0059.0-7242e'].skydir,distance=3.0)
gta.print_model()
Now we want to delete the source 3FGL J0021.6-6835:
In [20]:
gta.delete_source('3FGL J0021.6-6835')
gta.print_model()
You can also delete all the sources as a function of the npred/ts or position using the options minmax_npred/minmax_ts/skydir options. Finally, you can delete sources in a given list using the option names.
In [21]:
help(gta.delete_sources)
In the example below we delete all the sources that have an npred=[0,500]
In [22]:
gta.delete_sources(minmax_npred=[0,500])
gta.print_model()
We can aldo add a new source in the model using the gta.add_source function.
In [23]:
help(gta.add_source)
In the example below we add a Source called Source_PL with a PLSuperExpCutoff and pointlike and a source called Source_Gauss that is spatial extended with a PowerLaw SED and with a RadialGaussian template with an extension of 1 deg.
In [24]:
gta.add_source('Source_PL',{ 'glon' : 300., 'glat' : -46.,'SpectrumType' : 'PLSuperExpCutoff', 'Index1':-1.5, 'Index2' : 1.0,'Scale' : 1000,'Prefactor':1e-9,'SpatialModel' : 'PointSource' })
gta.add_source('Source_Gauss',{ 'glon' : 302., 'glat' : -45.,'SpectrumType' : 'PowerLaw', 'Index':2.0,'Scale' : 1000,'Prefactor':1e-9,'SpatialModel' : 'RadialGaussian', 'SpatialWidth': 1.0 })
gta.print_model()
Now I load the model saved into model_test to have back the intial model.
In [25]:
gta.load_roi('model_test')
gta.print_model()
It is also possible to modify the SED parameters using the functions gta.set_norm, gta.set_parameter, gta.set_parameter_bounds, gta.set_parameter_error .
In [26]:
help(gta.set_norm)
help(gta.set_parameter)
In [27]:
print gta.roi['3FGL J0059.0-7242e']['param_names']
print gta.roi['3FGL J0059.0-7242e']['param_values']
print gta.roi['3FGL J0059.0-7242e']['param_errors']
In the example below we fix the normalization to 1e-11 and the slope to 2.0.
In [28]:
gta.set_norm('3FGL J0059.0-7242e',value=1.)
gta.set_parameter('3FGL J0059.0-7242e',par='Index',value=2.0)
print gta.roi['3FGL J0059.0-7242e']['param_names']
print gta.roi['3FGL J0059.0-7242e']['param_values']
print gta.roi['3FGL J0059.0-7242e']['param_errors']
It is also possible to set the SED shape of a source using gta.set_source_spectrum tool.
In [30]:
help(gta.set_source_spectrum)
In [33]:
gta.load_roi('model_test')
gta.print_model()
gta.roi['3FGL J0059.0-7242e']['SpectrumType']
Out[33]:
In [34]:
gta.set_source_spectrum('3FGL J0059.0-7242e',spectrum_type='LogParabola')
gta.roi['3FGL J0059.0-7242e']['SpectrumType']
Out[34]:
In [35]:
gta.fit()
Out[35]:
In [36]:
print gta.roi['3FGL J0059.0-7242e']['param_names']
print gta.roi['3FGL J0059.0-7242e']['param_values']
print gta.roi['3FGL J0059.0-7242e']['param_errors']
As you can see above the tool gta.set_source_spectrum has modified the SED shape from PowerLaw to LogParabola.
It is also possible to set the spatial morphology of a source using gta.set_source_morphology
In [37]:
help(gta.set_source_morphology)
In [38]:
gta.set_source_morphology(name='3FGL J0029.1-7045',spatial_model='RadialGaussian',spatial_pars={'SpatialWidth': 1.0} )
In [39]:
gta.print_model()
print gta.roi['3FGL J0029.1-7045']['SpatialType']
print gta.roi['3FGL J0029.1-7045']['SpatialWidth']