* [EBL Attenuation models (20 June 2009)|#EBL Attenuation models (20 June 2009)]
* [SummedLikelihood (11 June 2009)|#SummedLikelihood. Unbinned and binned likelihood analysis objects can now be combined in joint analyses. (11 June 2009)]
* [SmoothBrokenPowerLaw (11 June 2009)|#SmoothBrokenPowerLaw (11 June 2009)]
* [Flux, energy flux, and upper limit calculations for diffuse sources (22 Feb 2009)|#Flux, energy flux, and upper limit calculations for diffuse sources (22 Feb 2009)]
* [Diffuse response calculations and FITS image templates (22 Feb 2009)|#Diffuse response calculations and FITS image templates (22 Feb 2009)]
* [Analysis object creation using the hoops/ape interface|#Analysis object creation using the hoops/ape interface]
* [Photon and Energy Flux Calculations|#Photon and Energy Flux Calculations]
h2. EBL Attenuation models (20 June 2009)
I've added access to the EBL attenuation models that are available to *gtobssim* (and Gleam) simulations via the [celestialSources/eblAtten|http://www-glast.stanford.edu/cgi-bin/viewcvs/celestialSources/eblAtten/] package. Here is an example xml definition for the PowerLaw2 spectral model:
{code}
<spectrum type="EblAtten::PowerLaw2">
<parameter free="1" max="1000.0" min="1e-05" name="Integral" scale="1e-06" value="2.0"/>
<parameter free="1" max="-1.0" min="-5.0" name="Index" scale="1.0" value="-2.0"/>
<parameter free="0" max="200000.0" min="20.0" name="LowerLimit" scale="1.0" value="20.0"/>
<parameter free="0" max="200000.0" min="20.0" name="UpperLimit" scale="1.0" value="2e5"/>
<parameter free="1" max="10" min="0" name="tau_norm" scale="1.0" value="1"/>
<parameter free="0" max="10" min="0" name="redshift" scale="1.0" value="0.5"/>
<parameter free="0" max="8" min="0" name="ebl_model" scale="1.0" value="0"/>
</spectrum>
{code}
As seen in this example, three parameters are added to the underlying spectral model: tau_norm, redshift, and ebl_model. The models with EBL attenuation that are available are
| EblAtten::PowerLaw2 | EblAtten::BrokenPowerLaw2 | EblAtten::LogParabola |
| EblAtten::BandFunction | EblAtten::SmoothBrokenPowerLaw | EblAtten::FileFunction |
| EblAtten::ExpCutoff | EblAtten::BPLExpCutoff | EblAtten::PLSuperExpCutoff |
The current list of EBL models are
|| model id || model name ||
| 0 | Kneiske |
| 1 | Primack05 |
| 2 | Kneiske_HighUV |
| 3 | Stecker05 |
| 4 | Franceschini |
| 5 | Finke |
| 6 | Gilmore |
The definitive list can be found in the enums defined in the [EblAtten|http://www-glast.stanford.edu/cgi-bin/viewcvs/celestialSources/eblAtten/eblAtten/EblAtten.h?view=markup] class, and references for each model are given in the header comments for the [implementation of the optical depth functions|http://www-glast.stanford.edu/cgi-bin/viewcvs/celestialSources/eblAtten/src/IRB_routines.cxx?view=markup].
celestialSources/eblAtten v0r7p1, optimizers v2r17, Likelihood v15r2
h2. SummedLikelihood. Unbinned and binned likelihood analysis objects can now be combined in joint analyses. (11 June 2009)
Here is a script that does a joint analysis of unbinned likelihood with binned front and binned back event likelihoods:
{code}
from UnbinnedAnalysis import *
from BinnedAnalysis import *
from SummedLikelihood import *
from UpperLimits import *
scfile = 'test_scData_0000.fits'
ltcube = 'expCube.fits'
irfs = 'P6_V3_DIFFUSE'
optimizer = 'MINUIT'
srcmdl = 'anticenter_model.xml'
like = unbinnedAnalysis(evfile='filtered.fits',
scfile=scfile, expmap='expMap.fits',
expcube=ltcube, irfs=irfs,
optimizer=optimizer, srcmdl=srcmdl)
like_u = unbinnedAnalysis(evfile='filtered_gt_3GeV.fits',
scfile=scfile, expmap='expMap.fits',
expcube=ltcube, irfs=irfs,
optimizer=optimizer, srcmdl=srcmdl)
like_f = binnedAnalysis(irfs=irfs, expcube=ltcube, srcmdl=srcmdl,
optimizer=optimizer,
cmap='smaps_front_3GeV.fits',
bexpmap='binned_exp_front_lt_3GeV.fits')
like_b = binnedAnalysis(irfs=irfs, expcube=ltcube, srcmdl=srcmdl,
optimizer=optimizer,
cmap='smaps_back_3GeV.fits',
bexpmap='binned_exp_back_lt_3GeV.fits')
summed_like = SummedLikelihood()
summed_like.addComponent(like_u)
summed_like.addComponent(like_f)
summed_like.addComponent(like_b)
like.fit(0)
print "Unbinned results:"
print like.model
summed_like.fit(0)
print "Combined unbinned and binned results:"
print summed_like.model
ul = UpperLimits(like)
ul_summed = UpperLimits(summed_like)
print "Unbinned upper limit example:"
print ul['PKS 0528+134'].compute()
print "Combined unbinned and binned upper limit example:"
print ul_summed['PKS 0528+134'].compute()
print like.Ts('Crab')
print summed_like.Ts('Crab')
{code}
Here is the output:
{noformat}
Unbinned results:
Crab
Spectrum: PowerLaw2
0 Integral: 1.955e+01 3.198e+00 1.000e-05 1.000e+03 ( 1.000e-07)
1 Index: -2.039e+00 1.091e-01 -5.000e+00 0.000e+00 ( 1.000e+00)
2 LowerLimit: 1.000e+02 0.000e+00 2.000e+01 3.000e+05 ( 1.000e+00) fixed
3 UpperLimit: 5.000e+05 0.000e+00 2.000e+01 5.000e+05 ( 1.000e+00) fixed
Extragalactic Diffuse
Spectrum: PowerLaw
4 Prefactor: 1.049e+00 5.228e-01 1.000e-05 1.000e+02 ( 1.000e-07)
5 Index: -2.136e+00 1.772e-01 -3.500e+00 -1.000e+00 ( 1.000e+00)
6 Scale: 1.000e+02 0.000e+00 5.000e+01 2.000e+02 ( 1.000e+00) fixed
GalProp Diffuse
Spectrum: ConstantValue
7 Value: 1.229e+00 6.405e-02 0.000e+00 1.000e+01 ( 1.000e+00)
Geminga
Spectrum: PowerLaw2
8 Integral: 3.463e+01 3.295e+00 1.000e-05 1.000e+03 ( 1.000e-07)
9 Index: -1.700e+00 5.405e-02 -5.000e+00 0.000e+00 ( 1.000e+00)
10 LowerLimit: 1.000e+02 0.000e+00 2.000e+01 3.000e+05 ( 1.000e+00) fixed
11 UpperLimit: 5.000e+05 0.000e+00 2.000e+01 5.000e+05 ( 1.000e+00) fixed
PKS 0528+134
Spectrum: PowerLaw2
12 Integral: 1.214e+01 3.081e+00 1.000e-05 1.000e+03 ( 1.000e-07)
13 Index: -2.555e+00 2.340e-01 -5.000e+00 0.000e+00 ( 1.000e+00)
14 LowerLimit: 1.000e+02 0.000e+00 2.000e+01 3.000e+05 ( 1.000e+00) fixed
15 UpperLimit: 5.000e+05 0.000e+00 2.000e+01 5.000e+05 ( 1.000e+00) fixed
Combined unbinned and binned results:
Crab
Spectrum: PowerLaw2
0 Integral: 1.927e+01 3.175e+00 1.000e-05 1.000e+03 ( 1.000e-07)
1 Index: -2.026e+00 1.083e-01 -5.000e+00 0.000e+00 ( 1.000e+00)
2 LowerLimit: 1.000e+02 0.000e+00 2.000e+01 3.000e+05 ( 1.000e+00) fixed
3 UpperLimit: 5.000e+05 0.000e+00 2.000e+01 5.000e+05 ( 1.000e+00) fixed
Extragalactic Diffuse
Spectrum: PowerLaw
4 Prefactor: 7.895e-01 6.264e-01 1.000e-05 1.000e+02 ( 1.000e-07)
5 Index: -2.088e+00 2.110e-01 -3.500e+00 -1.000e+00 ( 1.000e+00)
6 Scale: 1.000e+02 0.000e+00 5.000e+01 2.000e+02 ( 1.000e+00) fixed
GalProp Diffuse
Spectrum: ConstantValue
7 Value: 1.269e+00 7.695e-02 0.000e+00 1.000e+01 ( 1.000e+00)
Geminga
Spectrum: PowerLaw2
8 Integral: 3.494e+01 3.327e+00 1.000e-05 1.000e+03 ( 1.000e-07)
9 Index: -1.701e+00 5.386e-02 -5.000e+00 0.000e+00 ( 1.000e+00)
10 LowerLimit: 1.000e+02 0.000e+00 2.000e+01 3.000e+05 ( 1.000e+00) fixed
11 UpperLimit: 5.000e+05 0.000e+00 2.000e+01 5.000e+05 ( 1.000e+00) fixed
PKS 0528+134
Spectrum: PowerLaw2
12 Integral: 1.176e+01 3.008e+00 1.000e-05 1.000e+03 ( 1.000e-07)
13 Index: -2.534e+00 2.311e-01 -5.000e+00 0.000e+00 ( 1.000e+00)
14 LowerLimit: 1.000e+02 0.000e+00 2.000e+01 3.000e+05 ( 1.000e+00) fixed
15 UpperLimit: 5.000e+05 0.000e+00 2.000e+01 5.000e+05 ( 1.000e+00) fixed
Unbinned upper limit example:
0 12.1361746852 -0.000307310165226 1.21960706796e-06
1 13.368492501 0.0833717813548 1.34368135193e-06
2 14.6008103168 0.305775613851 1.46778865044e-06
3 15.8331281326 0.646238077657 1.59193260068e-06
4 17.0654459483 1.08809170817 1.71611295547e-06
5 17.8801837016 1.42911471373 1.7982332395e-06
(1.7803859938766175e-06, 17.703116306269838)
Combined unbinned and binned upper limit example:
0 11.7635006963 -0.000623093163085 1.18207829753e-06
1 12.9668988172 0.0780380520519 1.30323312058e-06
2 14.1702969382 0.291922864584 1.42440474102e-06
3 15.3736950591 0.621846241655 1.54561984399e-06
4 16.57709318 1.05134415647 1.66686579099e-06
5 17.7804913009 1.56725279011 1.78814776158e-06
(1.7382504833201744e-06, 17.285394704120577)
282.100179467
271.963534683
{noformat}
It would be wise to ensure that the datasets used in the summed likelihoods do not intersect.
Likelihood v15r0p3, pyLikelihood v1r14p1
h2. SmoothBrokenPowerLaw (11 June 2009)
Benoit implemented and added a new spectral model. Documentation can be found in the workbook [here|http://glast-ground.slac.stanford.edu/workbook/pages/sciTools_likelihoodTutorial/likelihoodTutorial.htm#createSourceModel].
Likelihood v15r1p2.
h2. Flux, energy flux, and upper limit calculations for diffuse sources (22 Feb 2009)
pyLikelihood v1r10p1; Likelihood v14r4
Flux, energy flux and upper limit calculations can now be made for diffuse sources in the python interface:
{code}
>>> like.model
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)
{code}
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.
h2. Diffuse response calculations and FITS image templates (22 Feb 2009)
Prior to Likelihood v14r3, there had been a long-standing problem with computing diffuse responses (with *gtdiffrsp*) for diffuse sources that have FITS image templates for discrete sources, such as SNRs or molecular clouds. In summary, *gtdiffrsp* performs the following integral for each diffuse source component {latex}i{latex} in the xml model definition:
{latex}
\newcommand{\e}{{\varepsilon}}
\newcommand{\ep}{{\varepsilon^\prime}}
\newcommand{\phat}{{\hat{p}}}
\newcommand{\phatp}{{\hat{p}^\prime}}
\begin{equation}
\int d\phat S_i(\phat, \ep_j) A(\ep_j, \phat) P(\phatp_j; \ep_j, \phat)
\end{equation}
{latex}
where {latex}$\hat{p}_j${latex} and {latex}$\varepsilon^\prime_j${latex} are the measured direction and measured energy of detected photon {latex}$j${latex}, and the integral is performed over (true) directions on the sky {latex}$\hat{p}${latex}. {latex}$A(...)${latex} and {latex}$P(...)${latex} are the effective area and point spread function, and I have made the approximation equating the true photon energy with the measured value, {latex}$\varepsilon = \varepsilon^\prime${latex}.
The problem arises for a discrete diffuse source when its spatial distribution {latex}$S_i(\hat{p})${latex} is only significantly different from zero far from the location of event {latex}$j${latex}. Operationally, the integral is evaluated using an adaptive Romberg integrator that samples the integrand at theta and phi values that are referenced to {latex}$\hat{p}_j${latex}. For very compact sources, the integrator will often miss the source entirely and evaluate the integral to zero; and unless the measured photon direction lies directly on a bright part of the extended source, the integral will usually not be very accurate, even if non-zero.
In Likelihood v14r3, the code has been modified to restrict the range of theta and phi values over which the integration is performed based on the boundaries of the input map. As long as the input map is not mostly filled with zero or near zero-valued entries, restricting the angular range this way seems to produce accurate diffuse response values.
This means that care should be taken in preparing the FITS image templates used for defining diffuse sources in Likelihood. The map should be as small as possible and should exclude regions where the emission is negligible. Here are examples of bad (left) and good (right) template files for the same source:
!snr_images.png|thumbnail!
The FITS file on the left is 1000x1000 pixels most of which are zero. This is almost optimally bad. Except for events lying right on the source, the diffuse response integrals are almost all zero. The image on the right is the same data, but cropped to the 17x23 pixel region that actually contains relevant data. Performing an unbinned likelihood fit using this smaller map produces accurate results even in the presence of significant Galactic diffuse emission.
h2. 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.
{code}
>>> from UnbinnedAnalysis import *
>>> like = unbinnedAnalysis()
Response functions to use[P6_V1_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
{code}
* Alternatively, one can give all of the parameters explicitly. This is useful for running in scripts.
{code}
>>> like2 = unbinnedAnalysis(evfile='filtered.fits', scfile='test_scData_0000.fits', irfs='P6_V1_DIFFUSE', expcube='expCube.fits', srcmdl='anticenter_model.xml', optimizer='minuit', expmap='none')
{code}
* 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.
{code}
>>> like3 = unbinnedAnalysis(evfile='filtered.fits', mode='h')
{code}
* A BinnedAnalysis object may be created in a similar fashion using the binnedAnalysis function:
{code}
>>> from BinnedAnalysis import *
>>> like = binnedAnalysis()
Response functions to use[P6_V1_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]
{code}
h2. Photon and Energy Flux Calculations
With Likelihood v13r18 and pyLikelihood v1r6 (ST v9r8), a facility has been added for calculating photon (ph/cm^2/s) and energy (MeV/cm^2/s) fluxes over a selectable energy range. The errors on these quantities are computed using the procedure described in [this presentation|Likelihood Usage Notes^Likelihood_dev.pdf] made at the [September 2008 Collaboration meeting|SCIGRPS:Issues on high-level analyses]. Here are some usage examples:
{code}
>>> 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
{code}
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:
{code}
>>> like.fit(covar=True)
{code} |