Source code for ixpeobssim.irf.rmf

# Copyright (C) 2015--2022, the ixpeobssim team.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

"""Definition of the response matrix.
"""

from __future__ import print_function, division

import numpy

from ixpeobssim.core.fitting import fit
from ixpeobssim.core.modeling import xGaussian
from ixpeobssim.core.rand import xUnivariateAuxGenerator
from ixpeobssim.core.spline import xInterpolatedUnivariateSplineLinear
from ixpeobssim.irf.base import xResponseBase
from ixpeobssim.utils.matplotlib_ import plt, last_line_color, setup_gca

# pylint: disable=invalid-name, too-many-ancestors




[docs] class xEnergyDispersionMatrix(xUnivariateAuxGenerator): """Class encapsulating the energy dispersion matrix, as stored in the MATRIX extension of a .rmf file. We don't downsample the matrix to avoid the failure of the fit procedure in XSPEC. Arguments --------- hdu : FITS hdu The MATRIX hdu in the .rmf FITS file. """ def __init__(self, hdu): """Constructor. """ data = hdu.data _x = numpy.arange(0, len(data['MATRIX'][0]), 1) _y = 0.5 * (data['ENERG_LO'] + data['ENERG_HI']) _z = data['MATRIX'].T fmt = dict(xlabel='Channel', ylabel='Energy [keV]') xUnivariateAuxGenerator.__init__(self, _x, _y, _z, **fmt)
[docs] @staticmethod def digitize(value, dtype=numpy.int16): """Digitize an energy value. This is essentially rounding an array to the nearest integer and cast the result to the specified type. """ return numpy.ndarray.astype(numpy.rint(value), dtype)
[docs] def rvs(self, aux): """Overloaded method. This is the original rvs() overloade method that was initially implemented in the energy dispersion code, and it is digitizing the measured energy before returning it, just like the detector would do in real life. We discovered along the way that it is handy to keep track of the measured (smneared) energy *before* the digitization (i.e., with all the significant digits), as this allows to apply corrections (e.g., for the GEM charging) without running into rounding problems. This is why we added the rvs2() class method, which is returning both. """ val = xUnivariateAuxGenerator.rvs(self, aux) return self.digitize(val)
[docs] def rvs2(self, aux): """Alterntative rvs() implementation where we return the measured energy both *before* and *after* the digitization. Note that the energy before the digitization is still in the pha space, and needs to be explicitely converted to physical units (keV), if that it what one is interested. When that's the case, the operation is deferred to the xEnergyDispersion class, who is aware of both the energy dispersion matric and the energy bounds. """ val = xUnivariateAuxGenerator.rvs(self, aux) return val, self.digitize(val)
[docs] class xEnergyDispersionBounds(xInterpolatedUnivariateSplineLinear): """Class encapsulating the bounds for the energy dispersion matrix, as stored in the EBOUNDS extension of a .rmf file. This is essentially a univariate spline with the channel number on the x axis and the average energy of the corresponding interval on the y axis. At the same time we do store the bounds of the intervals themselves. """ def __init__(self, hdu): """Constructor. """ _bounds = hdu.data self.chans = _bounds['CHANNEL'] self.emin = _bounds['E_MIN'] self.emax = _bounds['E_MAX'] _x = self.chans _y = 0.5 * (self.emin + self.emax) fmt = dict(xlabel='Channel', ylabel='Average energy [keV]') xInterpolatedUnivariateSplineLinear.__init__(self, _x, _y, **fmt)
[docs] def num_channels(self): """Return the number of channels. """ return len(self.chans)
[docs] def min(self): """Return the minimum energy of the first channel. """ return self.emin[0]
[docs] def max(self): """Return the maximum energy of the last channel. """ return self.emax[-1]
[docs] def channel_to_energy(self, channel): """Convert a channel number to energy. """ return self(channel)
[docs] def energy_to_channel(self, energy): """Convert a physical energy (in keV) into a channel number. This is using a simple binary search on the underlying energy bounds. """ if not isinstance(energy, numpy.ndarray): energy = numpy.array(energy) ch = numpy.searchsorted(self.emax, energy, side='left') if numpy.isscalar(ch): return numpy.int16(ch) return numpy.ndarray.astype(ch, numpy.int16)
[docs] class xEnergyDispersion(xResponseBase): """Class representing the energy dispersion. Arguments --------- file_path : str The path to the .rmf FITS file containing the energy dispersion tables. """ def __init__(self, file_path): """Constructor. """ xResponseBase.__init__(self, file_path, 'rmf') self.matrix = xEnergyDispersionMatrix(self.hdu_list['MATRIX']) self.ebounds = xEnergyDispersionBounds(self.hdu_list['EBOUNDS'])
[docs] def channel_to_energy(self, channel): """Convenient proxy to the underlying xEnergyDispersionBounds method. """ return self.ebounds.channel_to_energy(channel)
[docs] def energy_to_channel(self, energy): """Convenient proxy to the underlying xEnergyDispersionBounds method. """ return self.ebounds.energy_to_channel(energy)
[docs] @staticmethod def pha_to_pi(pha): """Simple placeholder method to convert a PHA value to pulse invariant. Note ---- Since we're not handling the gain corrections, yet, this is simply converting the values from integer to floating point numbers. """ return pha.astype(numpy.double)
[docs] def pha_analysis(self, energy): """Perform a pulse-height analysis on a given array of energies. """ pha = self.energy_to_channel(energy) pi = self.pha_to_pi(pha) return pha, pi
[docs] @staticmethod def scale_energy(energy): """Hook to allow to apply distorsions to the energy response of the detector. In this default implementation this is just returning the energy, untouched, but the method can be overridden to simulate a variety of systematic effects. """ return energy
[docs] def convolve_energy(self, mc_energy): """Convolve an array of mc_energy with the energy dispersion. This function is less trivial than the name might suggest, and it is performing several different operations at the same time, namely: * call xEnergyDispersionMatrix.rvs2() to smear the true energy and get the measured energy (in channel space) both before and after the digitization; * convert the measured energy before digitization to keV; * convert the pha to pulse invariant. """ # Smear the Monte Carlo energy with the energy dispersion, and retrieve # the measured energy (in channel space) both before and after the # digitization. energy, pha = self.matrix.rvs2(self.scale_energy(mc_energy)) # Convert the measured energy from channel to keV. energy = self.ebounds(energy) # Calculate the naive pulse invariant. pi = self.pha_to_pi(pha) return energy, pha, pi
[docs] def plot(self): """Plot the energy dispersion. """ plt.figure('%s energy dispersion' % self.base_name) self.matrix.plot() plt.figure('%s energy resolution' % self.base_name) for energy in [2.7, 5.9, 10.]: # Retrieve the energy dispersion 1-d spline at a given energy pdf = self.matrix.slice(energy) # Initial fit with a Gaussian... ch = self.energy_to_channel(energy) model = fit(xGaussian(), pdf.x, pdf.y, p0=(1., ch, 20.)) # ... refine the fit around the peak. xmin = model.peak - 1.5 * model.sigma xmax = model.peak + 2.5 * model.sigma model = fit(model, pdf.x, pdf.y, xmin=xmin, xmax=xmax) # Retrieve the energy resolution. eres = 100. * model.resolution() pdf.plot() label = '$\\sigma_E$/E = %.1f%% FWHM @ %.2f keV' % (eres, energy) x = model.peak y = pdf(x) plt.text(x - 28, y, label, va='bottom', color=last_line_color()) _, ymax = plt.ylim() setup_gca(ymin=0., ymax=1.15 * ymax)