Source code for gammapy.spectrum.extract
# Licensed under a 3-clause BSD style license - see LICENSE.rst
from __future__ import absolute_import, division, print_function, unicode_literals
import logging
import numpy as np
import astropy.units as u
from regions import CircleSkyRegion
from ..utils.scripts import make_path
from ..irf import PSF3D
from .core import PHACountsSpectrum
from .observation import SpectrumObservation, SpectrumObservationList
__all__ = ["SpectrumExtraction"]
log = logging.getLogger(__name__)
[docs]class SpectrumExtraction(object):
"""Creating input data to 1D spectrum fitting.
This class is responsible for extracting a
`~gammapy.spectrum.SpectrumObservation` from a
`~gammapy.data.DataStoreObservation`. The background estimation is done
beforehand, using e.g. the
`~gammapy.background.ReflectedRegionsBackgroundEstimator`. For point
sources analyzed with 'full containment' IRFs, a correction for PSF
leakage out of the circular ON region can be applied.
For more info see :ref:`spectral_fitting`.
For a usage example see :gp-extra-notebook:`spectrum_analysis`
Parameters
----------
observations : `~gammapy.data.Observations`
Observations to process
bkg_estimate : `~gammapy.background.BackgroundEstimate`
Background estimate, e.g. of
`~gammapy.background.ReflectedRegionsBackgroundEstimator`
e_reco : `~astropy.units.Quantity`, optional
Reconstructed energy binning
e_true : `~astropy.units.Quantity`, optional
True energy binning
containment_correction : bool
Apply containment correction for point sources and circular on regions.
max_alpha : float
Maximum alpha value to accept, if the background was estimated using
reflected regions this is 1 / minimum number of regions.
use_recommended_erange : bool
Extract spectrum only within the recommended valid energy range of the
effective area table (default is True).
"""
DEFAULT_TRUE_ENERGY = np.logspace(-2, 2.5, 109) * u.TeV
"""True energy axis to be used if not specified otherwise"""
DEFAULT_RECO_ENERGY = np.logspace(-2, 2, 73) * u.TeV
"""Reconstruced energy axis to be used if not specified otherwise"""
def __init__(
self,
observations,
bkg_estimate,
e_reco=None,
e_true=None,
containment_correction=False,
max_alpha=1,
use_recommended_erange=True,
):
self.observations = observations
self.bkg_estimate = bkg_estimate
self.e_reco = e_reco if e_reco is not None else self.DEFAULT_RECO_ENERGY
self.e_true = e_true if e_true is not None else self.DEFAULT_TRUE_ENERGY
self.containment_correction = containment_correction
self.max_alpha = max_alpha
self.use_recommended_erange = use_recommended_erange
self.spectrum_observations = SpectrumObservationList()
self.containment = None
self._on_vector = None
self._off_vector = None
self._aeff = None
self._edisp = None
[docs] def run(self):
"""Run all steps.
"""
log.info("Running {}".format(self))
for obs, bkg in zip(self.observations, self.bkg_estimate):
if not self._alpha_ok(bkg):
continue
self.spectrum_observations.append(self.process(obs, bkg))
def _alpha_ok(self, bkg):
"""Check if observation fulfills alpha criterion"""
condition = bkg.a_off == 0 or bkg.a_on / bkg.a_off > self.max_alpha
if condition:
msg = "Skipping because {} / {} > {}"
log.info(msg.format(bkg.a_on, bkg.a_off, self.max_alpha))
return False
else:
return True
[docs] def process(self, observation, bkg):
"""Process one observation.
Parameters
----------
observation : `~gammapy.data.DataStoreObservation`
Observation
bkg : `~gammapy.background.BackgroundEstimate`
Background estimate
Returns
-------
spectrum_observation : `~gammapy.spectrum.SpectrumObservation`
Spectrum observation
"""
log.info("Process observation\n {}".format(observation))
self.make_empty_vectors(observation, bkg)
self.extract_counts(bkg)
self.extract_irfs(observation)
if self.containment_correction:
self.apply_containment_correction(observation, bkg)
else:
self.containment = np.ones(self._aeff.energy.nbins)
spectrum_observation = SpectrumObservation(
on_vector=self._on_vector,
aeff=self._aeff,
off_vector=self._off_vector,
edisp=self._edisp,
)
if self.use_recommended_erange:
try:
spectrum_observation.hi_threshold = observation.aeff.high_threshold
spectrum_observation.lo_threshold = observation.aeff.low_threshold
except KeyError:
log.warning("No thresholds defined for obs {}".format(observation))
return spectrum_observation
[docs] def make_empty_vectors(self, observation, bkg):
"""Create empty vectors.
This method copies over all meta info and sets up the energy binning.
Parameters
----------
observation : `~gammapy.data.DataStoreObservation`
Observation
bkg : `~gammapy.background.BackgroundEstimate`
Background estimate
"""
log.info("Update observation meta info")
offset = observation.pointing_radec.separation(bkg.on_region.center)
log.info("Offset : {}\n".format(offset))
self._on_vector = PHACountsSpectrum(
energy_lo=self.e_reco[:-1],
energy_hi=self.e_reco[1:],
backscal=bkg.a_on,
offset=offset,
livetime=observation.observation_live_time_duration,
obs_id=observation.obs_id,
)
self._off_vector = self._on_vector.copy()
self._off_vector.is_bkg = True
self._off_vector.backscal = bkg.a_off
[docs] def extract_counts(self, bkg):
"""Fill on and off vector for one observation.
Parameters
----------
bkg : `~gammapy.background.BackgroundEstimate`
Background estimate
"""
log.info("Fill events")
self._on_vector.fill(bkg.on_events)
self._off_vector.fill(bkg.off_events)
[docs] def extract_irfs(self, observation):
"""Extract IRFs.
Parameters
----------
observation : `~gammapy.data.DataStoreObservation`
Observation
"""
log.info("Extract IRFs")
offset = self._on_vector.offset
self._aeff = observation.aeff.to_effective_area_table(
offset, energy=self.e_true
)
self._edisp = observation.edisp.to_energy_dispersion(
offset, e_reco=self.e_reco, e_true=self.e_true
)
[docs] def apply_containment_correction(self, observation, bkg):
"""Apply PSF containment correction.
Parameters
----------
observation : `~gammapy.data.DataStoreObservation`
observation
bkg : `~gammapy.background.BackgroundEstimate`
background esimate
"""
# TODO: This should be split out into a separate class
if not isinstance(bkg.on_region, CircleSkyRegion):
raise TypeError(
"Incorrect region type for containment correction."
" Should be CircleSkyRegion."
)
log.info("Apply containment correction")
# First need psf
angles = np.linspace(0.0, 1.5, 150) * u.deg
offset = self._on_vector.offset
if isinstance(observation.psf, PSF3D):
psf = observation.psf.to_energy_dependent_table_psf(theta=offset)
else:
psf = observation.psf.to_energy_dependent_table_psf(offset, angles)
center_energies = self._aeff.energy.nodes
containment = []
for index, energy in enumerate(center_energies):
try:
cont_ = psf.integral(energy, 0.0 * u.deg, bkg.on_region.radius)
except:
msg = "Containment correction failed for bin {}, energy {}."
log.warning(msg.format(index, energy))
cont_ = 1
finally:
containment.append(cont_)
self.containment = np.array(containment)
self._aeff.data.data *= self.containment
[docs] def compute_energy_threshold(self, **kwargs):
"""Compute and set the safe energy threshold for all observations.
See `SpectrumObservation.compute_energy_threshold` for full
documentation about the options.
"""
for obs in self.spectrum_observations:
obs.compute_energy_threshold(**kwargs)
[docs] def write(self, outdir, ogipdir="ogip_data", use_sherpa=False, overwrite=False):
"""Write results to disk.
Parameters
----------
outdir : `~gammapy.extern.pathlib.Path`
Output folder
ogipdir : str, optional
Folder name for OGIP data, default: 'ogip_data'
use_sherpa : bool, optional
Write Sherpa compliant files?
overwrite : bool
Overwrite existing files?
"""
outdir = make_path(outdir)
log.info("Writing OGIP files to {}".format(outdir / ogipdir))
outdir.mkdir(exist_ok=True, parents=True)
self.spectrum_observations.write(
outdir / ogipdir, use_sherpa=use_sherpa, overwrite=overwrite
)
# TODO : add more debug plots etc. here