r''' A module to convert between the flux and count data.
Some classes supporting this algorithm.
Differential flux conversion to count rate
- :class:`irfpy.cena.cena_flux.Flux2Count_Simple`.
- :class:`irfpy.cena.cena_flux.Flux2Count_Matrix`.
Count rate to differential flux.
- :class:`irfpy.cena.cena_flux.Count2Flux`.
.. |cm2| replace:: cm\ :sup:`2`
'''
import logging
import numpy as np
from scipy.integrate import quad
from irfpy.cena import gfactor
from irfpy.cena import cena_time
from irfpy.cena import energy
[docs]class Count2Flux:
''' Implementation of count to flux conversion.
According to the :file:`cena_calibration_issue_1_revision_1.pdf`,
count to flux conversion is implemented.
It is essential to understand the g-factor (a sort of conversion factor).
The g-factor, however, has not been solidly calculated.
This means that the values of g-factor and its affecting factors
my be changed in the near future.
The :class:`Count2Flux` class implementation is based on the calibration
issue 1 revision 1. In the near future, API would also be changed.
**Usage**
This class is a factory class.
Seveal G-factor implementaions can be set through g-factor class.
See :mod:`irfpy.cena.gfactor` for details.
To instance the conversion class, just call the constructor.
>>> c2f = Count2Flux()
The default is to use the g-factor class of
:class:`irfpy.cena.gfactor.GFactorH_Component`.
If one wants to change the gfactor class, for example to
:class:`irfpy.cena.gfactor.GFactorH_Table`,
you can specify the instance in the constructor.
The instance that has method :meth:`geometricFactor` can be used.
>>> from irfpy.cena.gfactor import GFactorH_Table
>>> gf = GFactorH_Table()
>>> c2f_tbl = Count2Flux(g_factor=gf)
**Todo**
- Support of *sigma*, the ambiguity.
'''
def __init__(self, g_factor=None):
'''
:keyword g_factor: G factor matrix.
By default (``None``), the g-factor is obtained from
:meth:`irfpy.cena.gfactor.GFactorH_Component`.
If you want to specify the g-factor, it should be a class instance
that have a method named ``geometricFactor`` that returns a ``np.masked_array``
with a shape of ``(16, 7)``. See :mod:`irfpy.cena.gfactor` module for details.
:type g_factor: ``None`` or ``np.masked_array`` with a shape of ``(16, 7)``.
'''
if g_factor is None:
g_factor = gfactor.GFactorH_Component()
try:
# Substitute gfactor from an instance.
# self.gfactor is [16][7] masked array.
self.gfactor = g_factor.geometricFactor()
except AttributeError as e:
logging.error('The g_factor keyword should be used with ' +
'an instance that has geometricFactor method')
raise e
self.ene = energy.getEnergyE16()
self.dt = cena_time.deltaT()
[docs] def getFlux(self, raw_count_rate, energy_step, sector_nr):
r''' Return the flux from the count rate.
Return the flux.
:param raw_count_rate: Raw count rate.
:type raw_count_rate: ``int``
:param sector_nr: Sector number, ranges *0..6*.
:type sector_nr: ``int``
:param energy_step: General energy step,
i.e in terms of E16, ranges *0..15*.
:type energy_step: ``int``
:returns: The particle flux, in the unit of #/|cm2| sr eV s.
:rtype: ``Float``.
The particle flux is calculated as in (3.19):
.. math::
J_{n, i, m} = \frac{C_{n, i, m}}{G_{n, i, m}
{\cdot}E_i{\cdot}\Delta t}
where *n*, *i* and *m* is the sector number, energy step,
and the mass index, respectively.
'''
return (raw_count_rate / self.gfactor[energy_step][sector_nr]
/ self.ene[energy_step] / self.dt)
[docs]class Flux2Count_Simple:
''' Class converting flux to count. Simplest implementation.
This provides inverse functions of :class:`Count2Flux`.
The flux to count conversion is not simple, because the flux is a
3-D continuous function.
This class provides a simple conversion.
It assumes 1-D energy spectra for the direction of your interest.
The spectra should **not** be a continuous function, but the
sampled value at the central energy of the CENA sensor.
The energy values are obtained from :meth:`irfpy.cena.energy.getEnergyE16` method.
'''
def __init__(self, g_factor=None):
'''
'''
if g_factor is None:
g_factor = gfactor.GFactorH_Component()
try:
self.gfactor = g_factor.geometricFactor()
except AttributeError as e:
self.logger.error('The g-factor keyword should be used with ' +
'an instance that has geometricFactor() method.')
raise
self.ene = energy.getEnergyE16()
self.dt = cena_time.deltaT()
[docs] def get_counts(self, flux_array, sector_nr):
''' Obtain the expected counts corresponding to the given flux.
:param flux_array: Differential flux in the unit of [#/cm2 sr eV s].
This parameter should be the ``np.array`` with the shape of (16,)
:type flux_array: ``np.array``
:param sector_nr: Sector number, 0 to 6.
:type sector_nr: ``int``
This is a array version of :meth:`getCount`.
Internally just call getCount 16 times.
>>> j = np.ones(16) * 1e5 # A constant flux over energy
>>> f2c = Flux2Count_Simple()
>>> c = f2c.get_counts(j, 3) # Get an array of counts
>>> c_8 = f2c.getCount(1e5, 8, 3) # Obtain count rate by each energy
>>> c[8] == c_8
True
'''
return np.array([self.getCount(flux_array[i], i, sector_nr) for i in range(16)])
[docs] def getCount(self, flux, energy_step, sector_nr):
r''' Return the expected count rate for the given flux and energy and sector.
:param flux: Flux in #/|cm2| sr eV s.
:type flux: ``float``
:param energy_step: The energy step, 0 to 15.
:type energy_step: ``int``
:param sector_nr: The sector number, 0 to 6.
:type sector_nr: ``int``
:returns: The expected count rate. #/sample.
:rtype: ``float``
The count rate will be estimated from (3.19):
.. math::
C_{n, i, m} = J_{n, i, m} {\cdot} G_{n,i,m} {\cdot} E_i {\cdot} \Delta t
where *n*, *i*, *m* is the sector number, energy step, and the mass index, respectively.
A sample for understanding.
The energy step for CENA can be found by:
>>> etbl = energy.getEnergyE16()
>>> print('%.1f' % etbl[8])
193.0
We take the central channel (Ch=3).
Then, assuming the flux at the direction 3 at energy of 193 eV as
10\ :sup:`5` #/|cm2| sr eV s, the expected count rate is calculated by the following.
>>> f2c = Flux2Count_Simple()
>>> f = 1e5
>>> c = f2c.getCount(f, 8, 3)
>>> print('%.3g' % c)
6.16
This class is the inversed function of the :class:`Count2Flux`. Thus, the value
can be reverted.
>>> c2f = Count2Flux()
>>> finv = c2f.getFlux(c, 8, 3)
>>> print('%.3e' % finv)
1.000e+05
'''
raw_cnt = flux * self.gfactor[energy_step][sector_nr] * self.ene[energy_step] * self.dt
return raw_cnt
[docs]class Flux2Count_Matrix:
''' Class converting flux to count. Matrix implementation.
See also :ref:`gfactor_matrix`.
'''
def __init__(self, energy_response_type="simulated"):
self.gfactor = gfactor.GFactorH_Component().geometricFactor() # (16, 7)
# This is the classical g-factor (integral)
self.dt = cena_time.deltaT()
self.dgfactor_matrix = self.get_differential_gfactor_matrix(
self.gfactor, energy_response_type=energy_response_type) # (16, 16, 7)
'''Differential g-factor matrix'''
[docs] @staticmethod
def get_differential_gfactor_matrix(g_factor, energy_response_type="simulated"):
r''' Return the differential gfactor matrix.
:param g_factor: G-factor (cm2 sr eV/eV)
array in ``np.array`` with the shape of (16, 7).
Simplest way of getting is just
``irfpy.cena.gfactor.GFactorH_Component().geometricFactor()``
:keyword energy_response_type: Energy response functions.
``simulated`` uses :mod:`irfpy.cena.energy.response_simulation`.
``simple_triangle`` uses :mod:`irfpy.cena.energy.response_simple_triangle`.
:returns: Differential gfactor matrix in ``np.array``
with the shape of (ie=16, je=16, 7).
Here *je* defines the beam (particle) energy and *ie* defines
the sensor energy.
The unit of the returned array is in cm:sup:`^2` sr eV.
:math:`\Delta E` has already be multiplied.
'''
# G-factor, Gi is to be ge_i, by multiplying Ei.
etbl = energy.getEnergyE16() # Shape (16)
delta_e = energy.getDeltaE16() # Shape 16
ge_i = (g_factor.T * etbl).T # Shape (16,7) -- ge,i (cm2 sr eV)
full_energy_response = np.zeros([16, 16])
gfactor_matrix = np.zeros([16, 16, 7])
for ic in range(7):
for je in range(16): # ie defines the beam (incoming) energy
if energy_response_type == "simulated":
energy_response = energy.response_simulation(
etbl[je], etbl) # Shape = (16,)
elif energy_response_type == "simple_triangle":
energy_response = energy.response_simple_triangle(
etbl[je], etbl) # Shape (16,)
else:
raise ValueError('Unknown response type (%s)' % energy_response_type)
k_ie_je = energy_response # Shape (16,), relative
full_energy_response[:, je] = k_ie_je * delta_e # Still relative but integral values
# Normalize. Use ge_i as conservation along j.
for ie in range(16):
full_energy_response_normalization = full_energy_response[ie, :].sum()
full_energy_response[ie, :] = full_energy_response[ie, :] / full_energy_response_normalization * ge_i[ie, ic]
gfactor_matrix[:, :, ic] = full_energy_response
# print gfactor_matrix[8, :, 3]
return gfactor_matrix
[docs] def get_counts(self, flux_array, sector_nr):
''' Obtain the expected counts corresponding to the given flux.
:param flux_array: Differential flux in the unit of [#/cm2 sr eV s].
This parameter should be the ``np.array`` with the shape of (16,)
:type flux_array: ``np.array``
:param sector_nr: Sector number, 0 to 6.
:type sector_nr: ``int``
'''
return (self.dgfactor_matrix[:, :, sector_nr]).dot(np.array(flux_array)) * self.dt
[docs]class Flux2Count_Integral:
''' Class converting flux to count. Integral implementation.
See also :ref:`gfactor_matrix`.
Integration is depending on the ``quad`` in ``scipy.integrate``, this class is heavy.
For my iMac (2007 late), it takes about 1.5 minutes to instance.
'''
def __init__(self, energy_response_type='simulated'):
# Classical g-factor.
self.gfactor = gfactor.GFactorH_Component().geometricFactor() # (16, 7)
self.dt = cena_time.deltaT()
self.energy_response_type = energy_response_type
self.__differential_gfactor(self.gfactor, energy_response_type=energy_response_type)
def __differential_gfactor(self, g_factor, energy_response_type='simulated'):
''' Return the differentil gfactor function
:param g_factor: Classical g-factor in cm2 sr eV/eV.
``np.array`` with (16, 7) shape. Usually you can get
:meth:`irfpy.cena.gfactor.GFactorH_Component.geometricFactor`.
:keyword energy_response_type: Energy response function.
``'simulates'`` uses :mod:`irfpy.cena.energy.response_simulation`.
``'simple_triangle'`` uses :mod:`irfpy.cena.energy.response_simple_triangle`.
:returns: Array of differential gfactor functions.
Array is 2-D array of functions, with length of 16x7.
Each function has an argument of energy in ``eV`` to return the differential flux in ``/cm2 sr eV s``.
'''
self.etbl = energy.getEnergyE16() # (16,)
delta_e = energy.getDeltaE16() # (16,)
self.ge_i = (g_factor.T * self.etbl).T # (16, 7) ge,i (cm2 sr eV)
self.diff_gf = [[None for j in range(7)] for i in range(16)]
class __GFactorInteg:
''' A helper class of gfactor integration.
'''
def __init__(self, ie, ic, ene, ge_i, energy_response_type):
'''
:param ie: Energy step
:param ic: Channel
:param ene: Energy in eV.
:param ge_i: Gfactor in cm2 sr eV s
:param energy_response_type: 'simulated' or 'simple_triangle'
'''
# from pudb import set_trace; set_trace()
self.ie = ie
self.ic = ic
self.ene = ene
self.ge_i = ge_i
if energy_response_type == 'simulated':
self.func = lambda eev: energy.response_simulation(self.ene, eev)
elif energy_response_type == 'simple_triangle':
self.func = lambda eev: energy.response_simple_triangle(self.ene, eev)
else:
raise ValueError('Unknown response type (%s)' % energy_response_type)
self.integ = quad(self.func, 0, 20 * self.ene)[0]
print(self.integ)
def differential_gfactor(self, eev):
return self.func(eev) / self.integ * self.ge_i
def integrate(self, flux_function):
tobeinteg = lambda eev: flux_function(eev) * self.differential_gfactor(eev)
return quad(tobeinteg, 0, 20 * self.ene)[0]
def integrate2(self, flux_function):
esplit = np.logspace(0, 4, num=1000) # 1 eV to 10000 eV.
esplit_delta = esplit[1:] - esplit[:-1]
esplit_center = (esplit[1:] + esplit[:-1]) / 2.
return (flux_function(esplit_center) * self.differential_gfactor(esplit_center) * esplit_delta).sum()
[docs] def calculate_diff_gf(self, ie, ic):
# The g-factor function, gfun, should have this two conditions.
# 1) gfun(E) = energy_response_func(E) * something
# 2) integral(gfun(E)) = ge_i[ie, ic]
logger = logging.getLogger('calculate_diff_gf')
logger.setLevel(logging.DEBUG)
logger.debug('Calculate ie=%d ic=%d' % (ie, ic))
self.diff_gf[ie][ic] = self.__GFactorInteg(ie, ic, self.etbl[ie], self.ge_i[ie, ic], self.energy_response_type)
[docs] def getCount(self, flux_function, energy_step, channel_step):
''' Return the count rate for the given function.
:param flux_function: A function which the differential flux follows.
It is a function with an argument of energy in ``eV`` returning the differential flux ``/cm2 sr eV s``.
:type flux_function: ``callable``
:param energy_step: Energy index
:param channel_step: Channel index.
'''
if self.diff_gf[energy_step][channel_step] is None:
self.calculate_diff_gf(energy_step, channel_step)
return self.diff_gf[energy_step][channel_step].integrate(flux_function) * self.dt
[docs] def getCount2(self, flux_function, energy_step, channel_step):
''' Same as :meth:getCount:, but faster. Precision would be ok.
Note that all the function should take np.array.
'''
if self.diff_gf[energy_step][channel_step] is None:
self.calculate_diff_gf(energy_step, channel_step)
return self.diff_gf[energy_step][channel_step].integrate2(flux_function) * self.dt