| Wiki Markup |
|---|
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}$n_{k^\prime}${latex} are the counts in measured energy bin {latex}$k^\prime${latex} and {latex}$n^k${latex} are the counts in true energy bin {latex}$k${latex}. The DRM calcuation 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}$D(E^\prime; E, \theta, \phi$){latex} is the energy dispersion function, {latex}$A(E, \theta, \phi)${latex} is the effective area, and {latex}$lt(\theta, \phi)${latex} is the integrated livetime as a function of detector coordinates associated with the specified sky position. {latex}$E_k${latex} is the logarithmic center of the {latex}$k${latex}th true energy bin. The integral over measured energy is taken over the width of the {latex}$k^\prime${latex}th bin.
In princple, 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:
!edisp_plots.png|thumbnail!
h3. 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}
{latex}
where {latex}$n_j${latex} is the observed count in pixel
{latex}$j${latex} and {latex}$\theta_j${latex} 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${latex}, the convolution is an {latex}${\cal
O}(n_e^2)${latex} operation. For {latex}$n_e = 30${latex}, this
results in about a factor of 900 slow down in computing 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${latex}, 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 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.
h3. Performance
Two simulation tests are applied.
h4. Soft power-law source, isotropic diffuse background.
* The xml model that was input to gtobssim
{code}
<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>
{code}
* The fitted model used in the binned analysis
{code}
<?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="0" 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>
{code}
* Distributions of photon index for the point source
|
Page History
Overview
Content Tools