Source code for ixpeobssim.irfgen.astar

#!/usr/bin/env python
# Copyright (C) 2020, 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
# 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.

from __future__ import print_function, division

"""Small convenience module for interfacing with the astar database.

import numpy
import os

from ixpeobssim import IXPEOBSSIM_IRFGEN
from ixpeobssim.utils.logging_ import logger
from ixpeobssim.core.spline import xInterpolatedUnivariateSpline
from ixpeobssim.utils.matplotlib_ import plt
from ixpeobssim.irfgen.constants import H_MASS, C_MASS, O_MASS, DME_MASS
from ixpeobssim.irfgen.constants import dme_density

IXPEOBSSIM_ASTAR_DATA = os.path.join(IXPEOBSSIM_IRFGEN, 'data', 'astar')

[docs] class xAlphaStoppingPowerTable: """Basic interface to the output files of the ASTAR database. """ def __init__(self, identifier, density=None): """Constructor. There has been a lot of back and forth on this class, essentially because we want to code things in a sensible way for elements and at the same time we want to make it practical to handle compounds (e.g., DME). A few noticeable things: * here we are using MeV and cm (instead of keV and mm) throughout, as these are the "natural" units in this business; * we support an optional density argument in the constructor, that allows to handle solids in a natural fashion; this implies some bookkeeping with the units, as the numbers have a different meaning depending of whether they are normalized to the density or not. """ self.identifier = identifier self.density = density # Read the underlying data file. file_name = '%s.txt' % identifier.lower() file_path = os.path.join(IXPEOBSSIM_ASTAR_DATA, file_name)'Parsing ASTAR data file %s...' % file_path) self._energy, _, _, self._stopping_power, self._csda_range, _, _ = \ numpy.loadtxt(file_path, unpack=True) self.process_data()
[docs] def process_data(self): """Do all the post-processing of the underlying data. """ # Some book-keeping. energy_label = '$\\alpha$ kinetic energy [MeV]' stopping_power_label = 'Stopping power [%s]' range_label = 'CSDA range [%s]' if self.density is None: stopping_power_label = stopping_power_label % 'MeV cm$^2$ g$^{-1}$' range_label = range_label % 'g cm$^{-2}$' else: stopping_power_label = stopping_power_label % 'MeV cm$^{-1}$' range_label = range_label % 'cm' # If we define the density of the material, we need to rescale the # relevant quantities. self._stopping_power *= self.density self._csda_range /= self.density # Build a spline with the stopping power as a function of the energy. args = self._energy, self._stopping_power fmt = dict(xlabel=energy_label, ylabel=stopping_power_label) self.stopping_power = xInterpolatedUnivariateSpline(*args, **fmt) # Calculate the inverse stopping power spline, which will be handy for # later use. args = self._energy, 1. / self._stopping_power fmt = dict(xlabel=energy_label, ylabel='Inverse stopping power') self.inverse_stopping_power = xInterpolatedUnivariateSpline(*args, **fmt) # Now to the CSDA range spline. args = self._energy, self._csda_range fmt = dict(xlabel=energy_label, ylabel=range_label) self.csda_range = xInterpolatedUnivariateSpline(*args, **fmt)
[docs] def bragg_curve(self, energy=10., energy_step=0.01): """Return a spline representing the Bragg curve for the element or compound. """ depth = 0. x = [] y = [] while energy > 0.: stopping_power = self.stopping_power(energy) x.append(depth) if energy < energy_step: y.append(0.) break y.append(stopping_power) energy -= energy_step depth += energy_step / stopping_power x = numpy.array(x) y = numpy.array(y) fmt = dict(xlabel='Path Length [cm]', ylabel=self.stopping_power.ylabel) return xInterpolatedUnivariateSpline(x, y, **fmt)
def _energy_profile(self, energy=10., energy_step=0.01): """Return a spline with the longitudinal energy profile of the alpha at a given energy as a function of the depth in the material. Note that for this to make sense the density has to be defined. """ assert self.density is not None depth = 0. range_ = self.csda_range(energy) x = [depth] y = [energy] while energy > 0.: stopping_power = self.stopping_power(energy) energy -= energy_step depth += energy_step / stopping_power if energy < energy_step: # Add the last step. x.append(range_) y.append(0.) break x.append(depth) y.append(energy) x = numpy.array(x) y = numpy.array(y) fmt = dict(xlabel='Path Length [cm]', ylabel=self.stopping_power.xlabel) # Note that we set the ext spline parameter to 'const' so that the # energy profile returns 0. beyond the particle range. return xInterpolatedUnivariateSpline(x, y, ext='const', **fmt)
[docs] def energy_loss(self, max_energy=10., energy_step=0.01): """ """ # Build the energy profile. _profile = self._energy_profile(max_energy, energy_step) # Invert the energy profile (mind you have to flip the arrays, as # y is decreasing). # Note we pass the axis=0 argument explicitly in order for this to # work with numpy versions prior to 1.15 _x = numpy.flip(_profile.y, 0) _y = numpy.flip(_profile.x, 0) _inverse_profile = xInterpolatedUnivariateSpline(_x, _y, ext='const') def _loss(energy, path_length): """ """ d = _inverse_profile(energy) return _profile(d) - _profile(d + path_length) return _loss
[docs] class xDMEAlphaStoppingPowerTable(xAlphaStoppingPowerTable): """Specialized subclass for the DME compound. """ def __init__(self, temperature, pressure): """Overloaded constructor. """ self.identifier = 'DME' self.density = dme_density(temperature, pressure) h = load_alpha_stopping_power_data('H') c = load_alpha_stopping_power_data('C') o = load_alpha_stopping_power_data('O') self._energy = h._energy self._stopping_power = h._stopping_power * (6. * H_MASS / DME_MASS) +\ c._stopping_power * (2. * C_MASS / DME_MASS) +\ o._stopping_power * (1. * O_MASS / DME_MASS) self.__calculate_csda_range() self.process_data() def __calculate_csda_range(self): """Integrate the inverse of the stopping power, yielding the range as a function of energy in CSDA approximation. """ s = xInterpolatedUnivariateSpline(self._energy, 1. / self._stopping_power) self._csda_range = numpy.zeros(self._energy.shape, 'd') emin = self._energy.min() for i, ebar in enumerate(self._energy): self._csda_range[i] = s.integral(emin, ebar) self._csda_range += emin * s(0.)
[docs] def load_alpha_stopping_power_data(identifier, density=None): """Load the stopping power data for a given element or compound. """ return xAlphaStoppingPowerTable(identifier, density)
[docs] def load_dme_alpha_stopping_power_data(temperature=20., pressure=800.): """Load the stopping power data for a given element or compound. """ return xDMEAlphaStoppingPowerTable(temperature, pressure)