''' An implementation of CENA g-factor.
G-Factor implementation according to the section 3.5 of the calibration report.
The implementation is done in :class:`GFactorH` class.
But this is somehow a test purpose, so the class will be re-implemented some day.
There are two implementations.
:class:`GFactorH_Table` starts with a tabulared G_i in Table 3.5 in the calibration report.
:class:`GFactorH_Component` starts with a calculation from component-based parameters.
One can start from
.. code-block:: py
from irfpy.cena.gfactor import GFactorH_Table as GFactorH
or
.. code-block:: py
from irfpy.cena.gfactor import GFactorH_Component as GFactorH
.. warning::
Apart from the updated definition in 3.19 in calibration report i1r2,
the dead time is in the g-factor.
It concerns the class:`GFactorH_Component`.
Also, to calculate the flux, you must use dt=0.5 sec.
'''
import logging
import math
import numpy
from irfpy.cena import fov
from irfpy.cena import cena_time
from irfpy.cena import energy
[docs]class GFactorH_Base():
r''' An base class of GFactor classes.
This is an abstract class. Derived classes should implement
:meth:`energyDependentGfactor` (:math:`G_{n,i}`\ ).
'''
def __init__(self):
''' Constructor.
'''
self.prm = GFactorParams()
[docs] def G_n_i(self, sigma=0):
''' Alias to :meth:`geometricFactor`.
'''
return self.geometricFactor(sigma=sigma)
[docs] def G_i(self, sigma=0):
''' Alias to :meth:`energyDependentGfactor`.
'''
return self.energyDependentGfactor(sigma=sigma)
[docs] def energyDependentGfactor(self, sigma=0):
''' Return the energy dependent g-factor. (Abstract class)
Override this method in the subclass.
'''
raise NotImplementedError("Override this method.")
[docs] def geometricFactor(self, sigma=0):
r''' Return the g-factor.
The g-factor is a kind of conversion factor between the differential flux, J,
and the counts, C.
Here, we denote :math:`J_{n,i}^{H}` [#/cm^2 sr eV s] is the differential number
flux seen by each sector n and energy step i for hydrogen ENA.
The corresponding count rate, :math:`C_{n,i}^{H}` [#] is calculated by
.. math::
C_{n,i}^{H} = J_{n,i}^{H} \cdot G_{n,i}^{H} \cdot E_i \cdot \Delta t
as in (3.15).
:param sigma: You can give an error parameter :math:`\sigma`.
:type sigma: ``float``
:returns: geometric factor as a ``numpy.ma.masked_array``.
The dimension is [i=16][n=7]
:rtype: Instance of ``numpy.ma.masked_array`` with dimensions of [i=16][n=7].
**Usage**
>>> gfH = GFactorH_Table()
>>> gf = gfH.geometricFactor()
>>> gf.shape
(16, 7)
.. code-block:: py
> print '%.1e' % gf[8][3]
5.8e-07
# While this method returns 6.6e-07 ...
**Note**
Indeed, the calculated value is different by about 10% from the calibration report.
Last column in Table 3.5 shows the gf[3][:], but the values are different....
.. todo::
The g-factor related classes is an objective of reimplementation.
The inconsist g-factors in the calibration report.
cena_flux class also should be re-implemented.
'''
geom_i = self.G_i(sigma=sigma).__array__() # dim is [i=16]
sensitivity_n = self.prm.s_n() # dim is [n=7]
solid_angle_n = numpy.array([fov.solid_angle(n) for n in range(7)]) # dim is [n=7]
sensXsolang_n = sensitivity_n * solid_angle_n # dim is [n=7]
gf = numpy.zeros([16, 7])
for nnn in range(7):
for iii in range(16):
gf[iii, nnn] = geom_i[iii] * sensXsolang_n[nnn]
gf = numpy.ma.masked_where(gf<0, gf)
return gf
[docs]class GFactorH_Table(GFactorH_Base):
''' An implementation of G-factor for hydrogen ENA.
This is based on section 3.5.2 in the calibration report.
**For developer**
The class is low-level class. Thus, in the near future, more high-level methods
may be needed. However, so far I used this low-level class also for ebugging purpose.
Thus, the class itself would also needed to be re-implemented.
'''
def __init__(self):
''' Constructor.
'''
GFactorH_Base.__init__(self)
[docs] def energyDependentGfactor(self, sigma=0):
r''' Return the energy dependent g-factor.
The energy dependent g-factor for :math:`H^0`, :math:`G_i^H`, is listed in ``Table 3.5``.
There are uncertainty of factor 3, corresponding to :math:`1\sigma`.
:param sigma: You can give a error parameter :math:`\sigma`.
For default value (``None``), central value of g-factor will be returned.
If you give +1, the g-factor becomes 3 times larger than the default
according to the uncertainty.
If you give -1, the g-factor becomes 3 times smaller than the default.
:type sigma: float
:returns: Energy depenent g-factor, as an array of 16 elements.
Each index corresponds to the generic energy steps.
(See also :mod:`irfpy.cena.energy`)
However, the data is only available for 4 to 11.
Thus, the unavailable element is a masked numpy array.
:rtype: ``numpy.ma.masked_array``
>>> gf = GFactorH_Table()
>>> g_i_H = gf.energyDependentGfactor()
>>> numpy.ma.is_masked(g_i_H[0])
True
>>> numpy.ma.is_masked(g_i_H[5])
False
>>> print((g_i_H[8]))
4.5e-05
>>> g_i_H_lower = gf.energyDependentGfactor(sigma=-1)
>>> print('%.2f' % (g_i_H_lower[5]/g_i_H[5]))
0.33
>>> g_i_H_upper = gf.energyDependentGfactor(sigma=1)
>>> print('%.2f' % (g_i_H_upper[5]/g_i_H[5]))
3.00
'''
arr = numpy.array([-9e30, -9e30, -9e30, -9e30,
2.1e-5, 2.4e-5, 2.8e-5, 3.5e-5,
4.5e-5, 6.0e-5, 8.2e-5, 1.2e-4,
-9e30, -9e30, -9e30, -9e30])
factor = math.pow(3, sigma)
return numpy.ma.masked_where(arr<0, arr) * factor
[docs]class GFactorH_Component(GFactorH_Base):
''' Another implementation of G-factor for hydrogen ENA.
This is based on section 3.5.2 in the calibration report.
The energy dependent geometric factors are calculated from scratch.
**For developer**
The class is low-level class. Thus, in the near future, more high-level methods
may be needed. However, so far I used this low-level class also for ebugging purpose.
Thus, the class itself would also needed to be re-implemented.
'''
def __init__(self):
GFactorH_Base.__init__(self)
[docs] def energyDependentGfactor(self, sigma=0):
r''' Return the energy dependent g-factor.
The energy dependent g-factor for :math:`H^0`, :math:`G_i^H`, is listed in ``Table 3.5``.
However, the values are inconsistent with those calculate from the equations.
Apaart from the inconsistency, There are uncertainty of factor 3, corresponding
to :math:`1\sigma`.
:param sigma: You can give a error parameter :math:`\sigma`.
For default value (``None``), central value of g-factor will be returned.
If you give +1, the g-factor becomes 3 times larger than the default
according to the uncertainty.
If you give -1, the g-factor becomes 3 times smaller than the default.
:type sigma: float
:returns: Energy depenent g-factor, as an array of 16 elements.
Each index corresponds to the generic energy steps.
(See also :mod:`irfpy.cena.energy`)
However, the data is only available for 4 to 11.
Thus, the unavailable element is a masked numpy array.
:rtype: ``numpy.ma.masked_array``
>>> gf = GFactorH_Component()
>>> g_i_H = gf.energyDependentGfactor()
>>> numpy.ma.is_masked(g_i_H[0])
False
>>> numpy.ma.is_masked(g_i_H[5])
False
>>> print(('%.2e' % g_i_H[8]))
4.38e-05
>>> g_i_H_lower = gf.energyDependentGfactor(sigma=-1)
>>> print('%.2f' % (g_i_H_lower[5]/g_i_H[5]))
0.33
>>> g_i_H_upper = gf.energyDependentGfactor(sigma=1)
>>> print('%.2f' % (g_i_H_upper[5]/g_i_H[5]))
3.00
**Definition of energy dependent g-factor**
From equation 3.19 in calibration report, :math:`G_i^H` is defined as
.. math::
G_i^H = \sigma{\cdot}A{\cdot}e_{i,H}{\cdot}\frac{\Delta E_i}{E_i}
\sigma = T_1{\cdot}T_2{\cdot}r_{CS}{\cdot}\eta_{START}{\cdot}
\eta_{MCP}{\cdot}\eta_{TOF}{\cdot}(1-\frac{t_D}{t_{SLOT}})
A = 1.09 \textrm{cm}^2
'''
p = self.prm
T1 = p.T1()
T2 = p.T2()
rcs = p.rcs()
eta_start = p.eta_start()
eta_mcp = p.eta_mcp()
eta_tof = p.eta_tof()
td = cena_time.dead_time()
tslot = cena_time.slot_time()
sigma_ = T1 * T2 * rcs * eta_start * eta_mcp * eta_tof * ( 1 - td / tslot)
logging.debug('sigma = %e' % sigma_)
A = p.A()
energy_in_eV = energy.getEnergyE16() # [16] array
delta_e_in_eV = energy.getFwhmRangeE16() # [2][16] array
delta_e_in_eV = delta_e_in_eV[1] - delta_e_in_eV[0]
e_i_H = numpy.array([p.eta_i_H_plus(energy_in_eV[i]) for i in range(16)])
g_i_H = sigma_ * A * e_i_H * delta_e_in_eV / energy_in_eV
factor = math.pow(3, sigma)
return numpy.ma.masked_where(g_i_H<0, g_i_H) * factor
[docs]class GFactorParams:
# def s_n(self):
# ''' Alias to :meth:`sectorSensitivity`.
# '''
# return self.getSectorSensitivity()
[docs] def getSectorSensitivity(self):
''' Returns the sector sensitivities.
Sensitivity depends on the sector of interest.
This method returns the sensitivities, :math:`s_n`.
Returned is 7-element ``numpy.array``, which is listed in ``Table 3.4``
in the calibration report.
'''
return numpy.array([0.525, 0.188, 0.443, 0.165, 0.613, 0.270, 0.285])
s_n = getSectorSensitivity
''' Alias to :meth:`getSectorSensitivity`.'''
[docs] def centralPixelsSensitivity(self):
''' Return the sensitivity of the sector 2-3-4.
It is sometimes worth to consider the sectors 2-3-4 is a single sector.
This method returns the sensitivity of this aggregated sector.
See also :meth:`sectorSensitivity`.
'''
return 1.0
s_234 = centralPixelsSensitivity
''' Alias to :meth:`centralPixelsSensitivity` '''
[docs] def getTransparencyGrid1(self):
''' Returns the first grid transparency, T1
This is geometrically defined.
T1 = 0.88 (see calibration report)
'''
return 0.88
T1 = getTransparencyGrid1
''' Alias to :meth:`getTransparencyGrid1` '''
[docs] def getTransparencyGrid2(self):
''' Returns the second grid transparency, T2
This is geometrically defined.
T2 = 0.88 (see calibration report)
'''
return 0.88
T2 = getTransparencyGrid2
''' Alias to :meth:`getTransparencyGrid2` '''
[docs] def getCollectionEfficiencyCS(self):
''' Returns the collection efficiency of wave system from CS.
Collection efficiency of wave system from conversion surface, r_CS is 0.7
'''
return 0.7
rcs = getCollectionEfficiencyCS
''' Alias to :meth:`getCollectionEfficiencyCS`. '''
[docs] def getPositiveIonYieldCS(self, energy_in_eV):
''' Returns the positive ion yield at CS.
Positive ion yield at the conversion surface is taken from Wieser et al. 2002.
The value depends on the energy of the impinging particle.
:param energy_in_ev: Energy of particle in eV
:type energy_in_ev: float
:returns: The yield
:rtype: float
'''
energy_ref = 1000.
return 0.005 + 0.05 * energy_in_eV / energy_ref
eta_i_H_plus = getPositiveIonYieldCS
''' Alias to :meth:`getPositiveIonYieldCS`. '''
def getStartMcpEfficiency(self):
r''' Start surface electron yield, i.e. :math:`\eta_{START}=1.0`.
It depends on MCP bias.
This is for post acceleration of ions with 2.4 kV.
'''
return 1.0
eta_start = getStartMcpEfficiency
''' Alias to :meth:`getStartMcpEfficiency`. '''
[docs] def getStartMcpEfficiency(self):
''' TOF efficiency.
Defining factors are found correted rate and start rate fro inflight data.
Using 2.1 kV start MCP bias, the quantiy is defined to be 0.15
'''
return 0.15
eta_mcp = getStartMcpEfficiency
''' Alias to :meth:`getStartMcpEfficiency` '''
[docs] def getTofEfficiency(self):
''' Sector sensitivity
The quantity is obtained from flight data.
'''
return 0.045
eta_tof = getTofEfficiency
''' Alias to :meth:`getTofEfficiency`. '''
[docs] def getOpenAperture(self):
''' Return the aperture in cm^2, based on (3.9).
'''
return 1.09
A = getOpenAperture
'''Alias to :meth:`getOpenAperture`'''