...
...
...
...
(DRM) Implementation
#Convolution Approximations
#Performance
#Usage
#User Performance Studies
Here is a description of the binned likelihood calculation.
Detector Response Matrix (DRM) Implementation
The DRM is the matrix that transforms a binned counts spectrum from true energies to measured energies:
Latex |
---|
Implementation|#Detector Response Matrix (DRM) Implementation] [#Convolution Approximations] [#Performance] Here is a [description of the binned likelihood calculation|https://confluence.slac.stanford.edu/download/attachments/3342362/binned.pdf?version=1&modificationDate=1254102725000]. h3. Detector Response Matrix (DRM) Implementation The DRM is the matrix that transforms a binned counts spectrum from true energies to measured energies: {latex} \begin{equation} n_{k^\prime} = \sum_k D_{kk^\prime} n_k \nonumber \end{equation} {latex} where {latex} |
where
Latex |
---|
$n_{k^\prime}$ |
...
...
the
...
counts
...
in
...
measured
...
energy
...
bin
...
Latex |
---|
...
$k^\prime$ |
Latex |
---|
$n^k$ |
Latex |
---|
$k$ |
Latex |
---|
{latex} and {latex}$n^k${latex} are the counts in true energy bin {latex}$k${latex}. The [DRM calculation in Likelihood|http://glast.stanford.edu/cgi-bin/viewcvs/Likelihood/src/Drm.cxx?view=markup] follows that performed in [gtrspgen|http://glast.stanford.edu/cgi-bin/viewcvs/rspgen/src/PointResponse.cxx?view=markup]: {latex} \begin{equation} D_{kk^\prime} = \frac{\int d\theta d\phi \left[\int_{\Delta E_{k^\prime}} dE^\prime D(E^\prime; E_k, \theta, \phi)\right] A(E_k, \theta, \phi) lt(\theta, \phi)}{\int d\theta d\phi A(E_k, \theta, \phi) lt(\theta, \phi)} \nonumber \end{equation} {latex} Here {latex} |
Here
Latex |
---|
$D(E^\prime; E, \theta, \phi$) |
...
...
the
...
energy
...
dispersion
...
function,
...
Latex |
---|
...
$A(E, \theta, \phi)$ |
...
...
the
...
effective
...
area,
...
and
...
Latex |
---|
...
$lt(\theta, \phi)$ |
...
...
the
...
integrated
...
livetime
...
as
...
a
...
function
...
of
...
detector
...
coordinates
...
associated
...
with
...
the
...
specified
...
sky
...
position.
...
Latex |
---|
...
$E_k$ |
...
...
the
...
logarithmic
...
center
...
of
...
the
...
Latex |
---|
...
$k$ |
...
...
true
...
energy
...
bin.
...
The
...
integral
...
over
...
measured
...
energy
...
is
...
taken
...
over
...
the
...
width
...
of
...
the
...
Latex |
---|
...
$k^\prime$ |
...
...
bin.
...
In
...
principle,
...
the
...
DRM
...
should
...
be
...
evaluated
...
at
...
each
...
sky
...
pixel
...
position
...
in
...
the
...
binned
...
counts
...
map,
...
but
...
we
...
make
...
the
...
approximation
...
that
...
the
...
DRM
...
does
...
not
...
change
...
much
...
over
...
the
...
counts
...
map
...
region
...
and
...
just
...
evaluate
...
it
...
at
...
the
...
map
...
center
...
and
...
assume
...
it
...
applies
...
everywhere
...
on
...
the
...
map.
...
This
...
is
...
supported
...
by
...
these
...
plots
...
of
...
the
...
energy
...
dispersion
...
evaluated
...
at
...
various
...
points
...
on
...
the
...
sky
...
for
...
a
...
one-day
...
survey
...
mode
...
integration:
...
Convolution Approximations
The model counts cubes used by the binned analysis are 2D spatial arrays of sky pixels with a counts spectrum at each pixel location forming the third dimension. When energy dispersion is neglected, the model counts spectrum is a function of true energy, but with the other aspects of the instrument response applied (i.e.,
...
effective
...
area
...
and
...
PSF).
...
To
...
complete
...
the
...
forward-folding
...
of
...
the
...
model
...
counts,
...
we
...
simply
...
need
...
to
...
apply
...
the
...
DRM
...
to
...
the
...
true
...
counts
...
spectrum
...
at
...
each
...
pixel.
...
Unfortunately,
...
there
...
are
...
two
...
constraints
...
that
...
prevent
...
us
...
from
...
doing
...
this:
...
- To
...
- save
...
- execution
...
- time,
...
- the
...
- log-likelihood
...
- is
...
- evaluated
...
- using
...
- only
...
- the
...
- pixels
...
- that
...
- have
...
- non-zero
...
- data
...
- counts
...
- (as
...
- opposed
...
- to
...
- model
...
- counts).
...
- This
...
- can
...
- be
...
- seen
...
- in
...
- the
...
- expression
...
- for
...
- the
...
- Poisson
...
- log-likelihood,
...
Latex \begin{equation} \log{\cal L} = \sum_j \left[n_j\log\theta_j - \theta_j\right] \nonumber \end{equation}
...
whereLatex $n_j$
...
is
...
- the
...
- observed
...
- count
...
- in
...
- pixel
...
Latex
...
and$j$
Latex $\theta_j$
...
is
...
- the
...
- model
...
- counts
...
- for
...
- that
...
- pixel.
...
- In
...
- the
...
- current
...
- implementation,
...
- this
...
- means
...
- that
...
- model
...
- counts
...
- are
...
- only
...
- computed
...
- for
...
- those
...
- non-zero
...
- pixels.
...
- For
...
- sparse
...
- datasets
...
- (e.g.,
...
- short
...
- integration
...
- times
...
- or
...
- at
...
- higher
...
- energies),
...
- this
...
- speeds
...
- up
...
- the
...
- calculation
...
- substantially:
...
- at
...
- least
...
- an
...
- order
...
- of
...
- magnitude
...
- for
...
- 1
...
- day
...
- integrations.
...
- Applying
...
- the
...
- DRM
...
- to
...
- each
...
- pixel
...
- counts
...
- spectrum
...
- would
...
- require
...
- the
...
- model
...
- calcuation
...
- to
...
- be
...
- made
...
- for
...
- every
...
- pixel
...
- in
...
- the
...
- counts
...
- cube,
...
- even
...
- if
...
- it
...
- is
...
- empty.
...
- However,
...
- assuming
...
- binned
...
- analysis
...
- is
...
- generally
...
- used
...
- for
...
- longer
...
- observations,
...
- and
...
- most
...
- pixels
...
- are
...
- occupied,
...
- this
...
- may
...
- not
...
- be
...
- a
...
- real
...
- limitation.
...
- A
...
- more
...
- serious
...
- problem
...
- is
...
- applying
...
- the
...
- DRM
...
- convolution
...
- to
...
- the
...
- true
...
- counts
...
- spectrum
...
- at
...
- each
...
- pixel
...
- location.
...
- If
...
- the
...
- number
...
- of
...
- energy
...
- bins
...
- is
...
Latex
...
$n_e$
...
,
...
- the
...
- convolution
...
- is
...
- an
...
Latex
...
${\cal O}(n_e^2)$
...
operation.
...
- For
...
Latex
...
$n_e = 30$
...
,
...
- this
...
- results
...
- in
...
- about
...
- a
...
- factor
...
- of
...
- 900
...
- slow
...
- down
...
- in
...
- computing
...
- the
...
- expected
...
- pixel
...
- counts
...
- that
...
- go
...
- into
...
- the
...
- log-likelihood.
...
Since
...
we
...
need
...
to
...
compute
...
the
...
total
...
predicted
...
counts
...
for
...
each
...
component,
...
we
...
have
...
the
...
total
...
counts
...
spectrum
...
in
...
true
...
energy
...
space
...
(This
...
is
...
easily
...
obtained
...
for
...
each
...
source
...
by
...
integrating
...
each
...
energy
...
plane
...
over
...
angle
...
in
...
the
...
associated
...
source
...
map
...
and
...
multiplying
...
by
...
the
...
spectral
...
function.)
...
We
...
can
...
then
...
convolve
...
the
...
spatially
...
integrated
...
true
...
energy
...
counts
...
spectrum
...
with
...
the
...
DRM
...
to
...
obtain
...
the
...
overall
...
measured
...
energy
...
counts
...
spectrum.
...
We
...
form
...
the
...
ratio
...
of
...
the
...
convolved
...
model
...
counts
...
to
...
unconvolved
...
model
...
counts
...
in
...
each
...
energy
...
bin.
...
When
...
computing
...
the
...
contribution
...
to
...
the
...
log-likelihood
...
from
...
each
...
pixel
...
Latex |
---|
...
$j$ |
...
...
we
...
multiply
...
the
...
unconvolved
...
model
...
counts
...
for
...
that
...
pixel
...
by
...
the
...
ratio
...
of
...
convolved
...
to
...
unconvolved
...
model
...
counts
...
from
...
the
...
spatially
...
integrated
...
spectrum.
...
This
...
procedure
...
does
...
not
...
really
...
account
...
for
...
any
...
differences
...
in
...
counts
...
spectra
...
that
...
arise
...
from
...
spatial
...
variation
...
of
...
the
...
source
...
spectrum
...
(for
...
diffuse
...
sources),
...
exposure
...
variations
...
across
...
the
...
map
...
region,
...
and
...
the
...
effect
...
of
...
the
...
energy-dependence
...
of
...
the
...
PSF
...
on
...
the
...
counts
...
spectrum
...
(the
...
observed
...
counts
...
spectrum
...
for
...
a
...
point
...
source
...
should
...
be
...
softer
...
near
...
the
...
source
...
location).
...
All
...
of
...
these
...
neglected
...
effects
...
(they
...
may
...
be
...
others
...
I
...
am
...
missing)
...
are
...
ameliorated
...
somewhat
...
by
...
the
...
fact
...
that
...
this
...
procedure
...
has
...
the
...
right
...
local
...
effect,
...
i.e.,
...
energy
...
bins
...
on
...
the
...
falling
...
part
...
of
...
the
...
counts
...
spectrum
...
tend
...
to
...
have
...
their
...
count
...
increased
...
while
...
bins
...
on
...
the
...
rising
...
part
...
have
...
their
...
count
...
decreased.
...
This
...
procedure
...
will
...
probably
...
not
...
work
...
very
...
well
...
for
...
sharp
...
spectral
...
features.
Performance
Composition Setup |
---|
Deck of Cards | ||||||||
---|---|---|---|---|---|---|---|---|
| =
|
| }
{||||||
Card | ||
---|---|---|
| ||
I've simulated a single point source with indices 1.25, 1.50, 2.00, 2.50, 3.00, generating 70-100k events for each case for a week long observation with idealized +/-50 deg rocking and fit those data with and without energy dispersion handling enabled:
The blue model curves and residual points correspond to the fit with energy dispersion turned off, while the red curves and residual points correspond to energy dispersion turned on. The corresponding data files have columns of
Emin | Emax | obs. counts | model counts w/ edisp | model counts w/out edisp |
The input spectra to gtobssim covered 20 to 1e6 MeV and the energy range selection for the binned analysis was 30 to 3e5 MeV, with 40 logarithmically spaced bins.
Photon index best-fit values
MC index | edisp handling | no edisp handling |
---|---|---|
1.25 | -1.253 +/- 1.86e-03 | -1.260 +/- 1.86e-03 |
1.50 | -1.505 +/- 2.19e-03 | -1.516 +/- 2.19e-03 |
2.00 | -2.007 +/- 3.20e-03 | -2.031 +/- 3.20e-03 |
2.50 | -2.525 +/- 4.17e-03 | -2.557 +/- 4.03e-03 |
3.00 | -3.036 +/- 4.97e-03 | -3.063 +/- 4.74e-03 |
Card |
---|
Card | ||
---|---|---|
| ||
For these tests, a point source with an exponentially cut-off power-law spectrum is modeled and fit for different values of the photon index. The cutoff energy is fixed at 1 GeV, and values of photon index = 1.25, 1.5, 2.0, 2.5, 3.0 are used. These tests are similar to those shown in fig. 68 of v1r0 of the LAT performance paper.
Simulations are performed for P7SOURCE_V6 and DC1A IRFs and these simulations are analyzed with and without the energy dispersion handling turned on for the standard gtobssim output and without the energy dispersion handling but using the MC values for the ENERGY column which is attained simply by swapping the labels on the ENERGY and MCENERGY columns in the FT1 file. Note that this swap was done prior to making an energy range selection.
- Input XML model to gtobssim
Code Block <source_library title="cutoff_pl_source"> <source name="cutoff_pl_source_3"> <spectrum escale="MeV" flux="20" particle_name="gamma"> <SpectrumClass name="FileSpectrum" params = "flux=1, specFile=Cutoff_PowerLaw_-3.00.txt"/> <celestial_dir ra="266.4" dec="-28.9"/> </spectrum> </source> <source name="cutoff_pl_source_2.5"> <spectrum escale="MeV" flux="7.8" particle_name="gamma"> <SpectrumClass name="FileSpectrum" params = "flux=1, specFile=Cutoff_PowerLaw_-2.50.txt"/> <celestial_dir ra="266.4" dec="-28.9"/> </spectrum> </source> <source name="cutoff_pl_source_2"> <spectrum escale="MeV" flux="3" particle_name="gamma"> <SpectrumClass name="FileSpectrum" params = "flux=1, specFile=Cutoff_PowerLaw_-2.00.txt"/> <celestial_dir ra="266.4" dec="-28.9"/> </spectrum> </source> <source name="cutoff_pl_source_1.5"> <spectrum escale="MeV" flux="1.14" particle_name="gamma"> <SpectrumClass name="FileSpectrum" params = "flux=1, specFile=Cutoff_PowerLaw_-1.50.txt"/> <celestial_dir ra="266.4" dec="-28.9"/> </spectrum> </source> <source name="cutoff_pl_source_1.25"> <spectrum escale="MeV" flux="0.73" particle_name="gamma"> <SpectrumClass name="FileSpectrum" params = "flux=1, specFile=Cutoff_PowerLaw_-1.25.txt"/> <celestial_dir ra="266.4" dec="-28.9"/> </spectrum> </source> </source_library>
...
The
...
- flux
...
- for
...
- each
...
- source
...
- was
...
- adjusted
...
- so
...
- that
...
- 100000
...
- events
...
- are
...
- generated
...
- for
...
- a
...
- week
...
- long
...
- idealized
...
- survey
...
- mode
...
- observation
...
- using
...
- P7SOURCE_V6.
...
- The
...
- template
...
- spectra
...
- were
...
- generated
...
- with
...
...
...
...
- ,
...
- using
...
- the
...
- PLSuperExpCutoff
...
- model
...
- that
...
- is
...
- available
...
- from
...
- pyLikelihood.
...
- This
...
- same
...
- model
...
- was
...
- used
...
- in
...
- the
...
- spectral
...
- fitting
...
- to
...
- ensure
...
- consistency.
...
Results
- P7SOURCE_V6
The red and green points substantially agree, as one would hope, and the trend is towards a harder measured spectra, as expected. However, there seems to be a residual bias in the edisp-on and MC results towards harder spectra.edisp-off
edisp handling turned off
edisp-on
edisp handling turned on
MC
edisp handling turned off, but using MCENERGY values
- DC1A (source fluxes were adjust to produce 200000 events)
Here the differences between the three cases is much smaller, but edisp-on and MC still agree better. The same residual bias persists.Card
Figure 68 of the instrument performance paper (v1r0) shows the "fractional" error on the photon index for a power-law and cutoff power-law model for different input indices:Card label Comparison with Instrument Performance Paper Figure
Using the same data that went into those simulations, here are the comparisons using energy dispersion handling:In each plot, the black points are obtained from refitting the data without energy dispersion handling (consistent with the paper figure), and the red points are the results obtained with energy dispersion turned on.Card Card label Soft power-law source, isotropic diffuse background - The xml model that was input to gtobssim
Code Block <source_library title="soft_source
...
"> <source
...
flux="1" name="soft_source"> <spectrum
...
escale="
...
MeV">
...
<particle
...
name="
...
gamma"> <power_law emax="1000000.0"
...
emin="
...
20.0"
...
gamma="
...
2.5"/>
...
</particle> <celestial_dir dec="-28.9" ra="266.4"/> </spectrum> </source>
...
<source
...
name="
...
test_isotropic"> <spectrum escale="MeV"> <SpectrumClass name="
...
Isotropic"
...
params="
...
100, 2.1, 20., 2e5"/>
...
<use_spectrum
...
frame="
...
galaxy"/> </spectrum> </source> </source_library>
- The model used in the binned analysis
Code Block <?xml version="1.0" ?> <source_library title="source library"> <source name="test_source" type="PointSource"> <!-- point source units are cm^-2 s^-1 MeV^-1 -->
...
<spectrum type="
...
PowerLaw2"> <parameter free="
...
1" max="
...
1000.0" min="1e-
...
05" name="
...
Integral" scale="
...
1e-06" value="
...
1.
...
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"/> </spectrum> <spatialModel type="SkyDirFunction"> <parameter free="0" max="
...
360.0" min="
...
-360.0" name="
...
RA" scale="1.0" value="
...
266.
...
4"/> <parameter free="0" max="
...
90.0" min="
...
-90.0" name="
...
DEC" scale="1.0" value="
...
-28.
...
9"/> </spatialModel> </
...
source> <source name="Extragalactic Diffuse" type="DiffuseSource">
...
<spectrum type="
...
PowerLaw2"> <parameter free="
...
1" max="
...
1000.0" min="0.0" name="
...
Integral" scale="
...
1e-06" value="1.0"/> <parameter free="
...
1" max="
...
-1.0" min="
...
-5.0" name="
...
Index" scale="1.0" value="-2.1
...
"/>
...
...
<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="200000.0"/> </spectrum> <spatialModel type="ConstantValue"> <parameter free="0" max="10.0" min="0.0" name="Value" scale="1.0" value="1.0"/> <parameter free="0" max="10.0" min="0.0" name="Value" scale="1.0" value="1.0"/> </spatialModel> </source> </source_library>
- Distributions from simulations. The leftmost plot in each row is the parameter distribution obtained from the standard analysis with the energy dispersion handling turned off. The middle plot has energy dispersion handling turned on. The rightmost plot uses the MC energy values in the analysis and has energy dispersion handling turned off.
- Photon Index of the soft point source.
- Photon Index of the isotropic component.
Card Card label Exponential cut-off with isotropic diffuse
- Photon Index of the soft point source.
- Input model to gtobssim:
Code Block <source_library title="cutoff_pl_source"> <source name="cutoff_pl_source"> <spectrum escale="MeV" flux="1" particle_name="gamma"> <SpectrumClass name="FileSpectrum" params = "flux=1e-2, specFile=$(rootdir)/Cutoff_PowerLaw.txt"/> <celestial_dir ra="266.4" dec="-28.9"/> </spectrum> </source> <source name="test_isotropic"> <spectrum escale="MeV"> <SpectrumClass name="Isotropic" params="100, 2.1, 20., 2e5"/> <use_spectrum frame="galaxy"/> </spectrum> </source> </source_library>
...
- Likelihood model
Code Block <?xml version="1.0" ?> <source_library title="source library"> <source name="test_source" type="PointSource"> <!-- point source units are cm^-2 s^-1 MeV^-1 --> <spectrum type="PLSuperExpCutoff"> <parameter free="1" max="1e3" min="1e-5" name="Prefactor" scale="1e-07" value="1.0"/> <parameter free="1" max="0.0" min="-5.0" name="Index1" scale="1.0" value="-1.7"/> <parameter free="0" max="1000.0" min="50.0" name="Scale" scale="1.0" value="200.0"/> <parameter free="1" max="30000.0" min="500.0" name="Cutoff" scale="1.0" value="3000.0"/> <parameter free="0" max="5.0" min="0.0" name="Index2" scale="1.0" value="1.0"/> </spectrum> <spatialModel type="SkyDirFunction"> <parameter free="0" max="360.0" min="-360.0" name="RA" scale="1.0" value="266.4"/> <parameter free="0" max="90.0" min="-90.0" name="DEC" scale="1.0" value="-28.9"/> </spatialModel> </source> <source name="Extragalactic Diffuse" type="DiffuseSource"> <spectrum type="PowerLaw2"> <parameter free="1" max="1e3" min="1e-5" name="Integral" scale="1e-06" value="1.0"/> <parameter free="1" max="-1.0" min="-5.0" name="Index" scale="1.0" value="-2.1"/> <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"/> </spectrum> <spatialModel type="ConstantValue"> <parameter free="0" max="10.0" min="0.0" name="Value" scale="1.0" value="1.0"/> </spatialModel> </source> </source_library>
...
- Distributions
- Photon index (Index1)
...
- of
...
- the
...
- exponentially
...
- cutoff
...
- power-law
...
- of
...
- the
...
- point
...
- source.
- Cutoff energy of the point source.
- Photon index of the isotropic component.
Card Deck of Cards
- source.
Usage
To enable the energy dispersion handling (in ST-09-26-00
...
and later), set the USE_BL_EDISP environment variable:
Code Block |
---|
% setenv USE_BL_EDISP
|
To disable it, unset the env var:
Code Block |
---|
% unsetenv USE_BL_EDISP
|
This functionality can also be turned on and off in python at runtime:
Code Block |
---|
>>> import os >>> os.environ[' later), set the USE_BL_EDISP environment variable: {code} % setenv USE_BL_EDISP {code} To disable it, unset the env var: {code} % unsetenv USE_BL_EDISP {code} This functionality can also be turned on and off in python at runtime: {code} >>> import os >>> os.environ['USE_BL_EDISP'] = "1" '] = "1" # The precise value doesn't matter >>> del os.environ['USE_BL_EDISP'] # In a script, this should probably be enclosed in a try/except block |
In ST-09-28-00 and later, energy dispersion handling may be turned on in python for individual binned analysis objects:
Code Block |
---|
>>> like = binnedAnalysis(...) >>> like.logLike.set_edisp_flag(True) # TheTurn preciseon valueenergy doesn'tdispersion matterhandling >>> delprint os.environ['USE_BL_EDISP']like.logLike.use_edisp() # InCheck if aedisp script,handling thisis shouldturned probably be enclosed in a try/except block {code}on >>> like.logLike.set_edisp_flag(False) # Turn off |
This will allow binned analysis objects within a single python script to have energy dispersion handling turned on or off independently, whereas the environment variable switch will affect all analysis objects within the script the same way.