Fermi Gamma-ray Space Telescope

pyLikelihood Usage Notes

Analysis object creation using the hoops/ape interface

The UnbinnedAnalysis and BinnedAnalysis modules contain functions that use the hoops/pil/ape interface to take advantage of the gtlike.par file for specifying inputs. Usage of this interface may be more convenient than creating the UnbinnedObs, UnbinnedAnalysis, BinnedObs, and BinnedAnalysis objects directly. Here are some examples of their use:

  • In this example, the call to the unbinnedAnalysis function is made without providing any arguments. The gtlike.par file is read from the user's PFILES path, and the various parameters are prompted for just as when running gtlike (except that the order is a bit different). If one is using the ipython interface with readline enabled, then tab-completion works as well.
>>> from UnbinnedAnalysis import *
>>> like = unbinnedAnalysis()
Response functions to use[P6_V3_DIFFUSE]
Spacecraft file[test_scData_0000.fits] test_scData_0000.fits
Event file[filtered.fits] filtered.fits
Unbinned exposure map[none]
Exposure hypercube file[expCube.fits]
Source model file[anticenter_model.xml]
Optimizer (DRMNFB|NEWMINUIT|MINUIT|DRMNGB|LBFGS) [MINUIT]
>>> like.model
Crab
Spectrum: PowerLaw2
0 Integral: 1.540e+01 0.000e+00 1.000e-05 1.000e+03 ( 1.000e-06)
1 Index: -2.190e+00 0.000e+00 -5.000e+00 0.000e+00 ( 1.000e+00)
2 LowerLimit: 2.000e+01 0.000e+00 2.000e+01 3.000e+05 ( 1.000e+00) fixed
3 UpperLimit: 2.000e+05 0.000e+00 2.000e+01 3.000e+05 ( 1.000e+00) fixed

Geminga
Spectrum: PowerLaw2
4 Integral: 1.020e+01 0.000e+00 1.000e-05 1.000e+03 ( 1.000e-06)
5 Index: -1.660e+00 0.000e+00 -5.000e+00 0.000e+00 ( 1.000e+00)
6 LowerLimit: 2.000e+01 0.000e+00 2.000e+01 3.000e+05 ( 1.000e+00) fixed
7 UpperLimit: 2.000e+05 0.000e+00 2.000e+01 3.000e+05 ( 1.000e+00) fixed

PKS 0528+134
Spectrum: PowerLaw2
8 Integral: 9.802e+00 0.000e+00 1.000e-05 1.000e+03 ( 1.000e-06)
9 Index: -2.460e+00 0.000e+00 -5.000e+00 0.000e+00 ( 1.000e+00)
10 LowerLimit: 2.000e+01 0.000e+00 2.000e+01 3.000e+05 ( 1.000e+00) fixed
11 UpperLimit: 2.000e+05 0.000e+00 2.000e+01 3.000e+05 ( 1.000e+00) fixed
  • Alternatively, one can give all of the parameters explicitly. This is useful for running in scripts.
>>> like2 = unbinnedAnalysis(evfile='filtered.fits', scfile='test_scData_0000.fits', irfs='P6_V3_DIFFUSE', expcube='expCube.fits', srcmdl='anticenter_model.xml', optimizer='minuit', expmap='none')
  • The mode='h' option is available as well. In this case, one can set a specific parameter, leaving the remaining ones to be read silently from the gtlike.par file.
>>> like3 = unbinnedAnalysis(evfile='filtered.fits', mode='h')
  • A BinnedAnalysis object may be created in a similar fashion using the binnedAnalysis function:
>>> from BinnedAnalysis import *
>>> like = binnedAnalysis()
Response functions to use[P6_V3_DIFFUSE::FRONT]
Counts map file[smaps_0_20_inc_front.fits]
Binned exposure map[binned_expmap_0_20_inc_front.fits]
Exposure hypercube file[expCube_0_20_inc.fits]
Source model file[Vela_model_0_20_inc_front.xml]
Optimizer (DRMNFB|NEWMINUIT|MINUIT|DRMNGB|LBFGS) [MINUIT]

» Back to Top

Photon and Energy Flux Calculations

A facility has been added for calculating photon (ph/cm^2/s) and energy (MeV/cm^2/s) fluxes over a selectable energy range. Here are some usage examples:

>>> like.model
Extragalactic Diffuse
Spectrum: PowerLaw
0 Prefactor: 7.527e-02 7.807e-01 1.000e-05 1.000e+02 ( 1.000e-07)
1 Index: -2.421e+00 1.968e+00 -3.500e+00 -1.000e+00 ( 1.000e+00)
2 Scale: 1.000e+02 0.000e+00 5.000e+01 2.000e+02 ( 1.000e+00) fixed

GalProp Diffuse
Spectrum: ConstantValue
3 Value: 1.186e+00 5.273e-02 0.000e+00 1.000e+01 ( 1.000e+00)

Vela
Spectrum: BrokenPowerLaw2
4 Integral: 9.176e-02 3.730e-03 1.000e-03 1.000e+03 ( 1.000e-04)
5 Index1: -1.683e+00 5.131e-02 -5.000e+00 -1.000e+00 ( 1.000e+00)
6 Index2: -3.077e+00 2.266e-01 -5.000e+00 -1.000e+00 ( 1.000e+00)
7 BreakValue: 1.716e+03 2.250e+02 3.000e+01 1.000e+04 ( 1.000e+00)
8 LowerLimit: 1.000e+02 0.000e+00 2.000e+01 2.000e+05 ( 1.000e+00) fixed
9 UpperLimit: 3.000e+05 0.000e+00 2.000e+01 5.000e+05 ( 1.000e+00) fixed

>>> like.flux('Vela', emin=100, emax=3e5)
9.2005203669330761e-06

>>> like.flux('Vela')
9.2005203669330761e-06

>>> like.fluxError('Vela')
3.73060211564e-07

>>> like.energyFlux('Vela')
0.0047870395747710718

>>> like.energyFluxError('Vela')
0.000260743911378

The default energy range for the flux and energy flux calculations is (emin, emax) = (100, 3e5) MeV. Either or both of these may be set as keyword arguments to the function call. The errors are available as separate function calls and require that the covariance matrix has been computed using "covar=True" keyword option to the fit function:

>>> like.fit(covar=True)

Flux, energy flux, and upper limit calculations for diffuse sources

Flux, energy flux and upper limit calculations can now be made for diffuse sources in the python interface:

Extragalactic Diffuse
Spectrum: PowerLaw
0 Prefactor: 1.450e+00 4.286e-01 1.000e-05 1.000e+02 ( 1.000e-07)
1 Index: -2.054e+00 1.366e-01 -3.500e+00 -1.000e+00 ( 1.000e+00)
2 Scale: 1.000e+02 0.000e+00 5.000e+01 2.000e+02 ( 1.000e+00) fixed

SNR_map
Spectrum: PowerLaw2
3 Integral: 6.953e+03 2.566e+03 0.000e+00 1.000e+10 ( 1.000e-06)
4 Index: -2.039e+00 1.332e-01 -5.000e+00 -1.000e+00 ( 1.000e+00)
5 LowerLimit: 2.000e+01 0.000e+00 2.000e+01 2.000e+05 ( 1.000e+00) fixed
6 UpperLimit: 5.000e+05 0.000e+00 2.000e+01 5.000e+05 ( 1.000e+00) fixed


>>> like.flux('Extragalactic Diffuse', 100, 5e5)
0.00017279621316112551

>>> like.fluxError('Extragalactic Diffuse', 100, 5e5)

3.32420171199e-05

>>> like.flux('SNR_map', 100, 5e5)
1.3482399819603797e-06

>>> like.fluxError('SNR_map', 100, 5e5)
2.54429802134e-07

>>> from UpperLimits import UpperLimits
>>> ul = UpperLimits(like)
>>> ul['Extragalactic Diffuse'].compute()
0 1.44990964752 0.000405768508472 0.000172805047298
1 1.62134142735 0.073917905358 0.000185125831974
2 1.79277320719 0.262761447131 0.000197078187151
3 1.96420498702 0.53963174062 0.00020864241284
4 2.13563676685 0.882155849953 0.000219869719057
5 2.33677514646 1.33948412506 0.000232753706737
6 2.37253076663 1.42681888532 0.000234966127989
(0.00023314676502589186, 2.34312748234)

>>> ul['SNR_map'].compute()
0 6953.47703963 7.03513949247e-05 1.34891151045e-06
1 7979.85337599 0.072031366891 1.43605234286e-06
2 9006.22971235 0.246899738264 1.51625454475e-06
3 10032.6060487 0.487519676687 1.59037705975e-06
4 11058.9823851 0.770755618154 1.65973696844e-06
5 12629.8902716 1.25421175769 1.75816622567e-06
6 13039.134979 1.38671547284 1.77783091654e-06
(1.7731240676197149e-06, 12941.1800697)

The compute() command will use profile likelihood to compute the upper limit and will thus scan in normalization parameter value. The screen output comprises the scan values with columns index, parameter value, delta(log-likelihood), flux. The values that are returned are the total flux upper limit (i.e., integrated over all angles) and the corresponding normalization parameter value.

Last updated by: Elizabeth Ferrara - 02/03/2010