Energy Dispersion in Binned Likelihood

Detector Response Matrix (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:

where are the counts in measured energy bin and are the counts in true energy bin . The DRM calculation in Likelihood follows that performed in gtrspgen:

Here is the energy dispersion function, is the effective area, and is the integrated livetime as a function of detector coordinates associated with the specified sky position. is the logarithmic center of the th true energy bin. The integral over measured energy is taken over the width of the th 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:

1. 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,
   
   where is the observed count in pixel and 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.
2. 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 , the convolution is an operation. For , 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 , 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

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

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

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

  <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 this python script, using the PLSuperExpCutoff model that is available from pyLikelihood. This same model was used in the spectral fitting to ensure consistency.

  Results

• P7SOURCE_V6
  
  

  edisp-off	edisp handling turned off
  edisp-on	edisp handling turned on
  MC	edisp handling turned off, but using MCENERGY values

  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.
• 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.
  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:
  
  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.
• The xml model that was input to gtobssim

  <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

  <?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.
• Input model to gtobssim:

  <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

  <?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.

    Error rendering macro 'deck'

    java.lang.NullPointerException

Usage

To enable the energy dispersion handling (in ST-09-26-00 and later), set the USE_BL_EDISP environment variable:

% setenv USE_BL_EDISP

To disable it, unset the env var:

% unsetenv USE_BL_EDISP

This functionality can also be turned on and off in python at runtime:

>>> import os
>>> os.environ['USE_BL_EDISP'] = "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:

>>> like = binnedAnalysis(...)

>>> like.logLike.set_edisp_flag(True)    # Turn on energy dispersion handling
>>> print like.logLike.use_edisp()       # Check if edisp handling is turned 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.

User Performance Studies

• A closer look at the energy dispersion - Stefan
• Energy dispersion studies using the smoothly broken powerlaw model - Jim