''' Implementation of CENA energy.
**Energy table**
CENA energy table based on the calibration report Issue 1 Rev 1.
Energy bin is generalized. There are 16 (TBC) energy settings in the sensor.
Three energy tables are used.
============ ========= ========= =========
Energy bin Table 1 Table 2 Table 3
============ ========= ========= =========
0 193 8 652 11 3300 15
1 129 7 435 10 2200 14
2 86 6 290 9 1467 13
3 57 5 193 8 978 12
4 38 4 129 7 652 11
5 25 3 86 6 435 10
6 17 2 57 5 290 9
7 11 1 38 4 193 8
============ ========= ========= =========
The definition of the energy ranges are illustrated :ref:`here <energy_definition>`.
One definition is the range between the boundaries of neighboring energy steps.
This is important for analysis to calculate the integral for example.
Another definition is defined by the range of the energy channel response, which
is also important for calibration, flux calculation, or error estimation.
In this module, the response range is distinguished by using ``Fwhm`` in the method name.
**Energy response**
While we do not exactly know the energy response, the function is implemented
based on the Kazama 2006 report.
This is approximated by a triangle, whose shape changed according to the incoming
beam energy. See :func:`response_simulation` for details.
'''
import logging
import numpy as np
[docs]def getEnergyE16():
''' Returns the list of the central energy in eV for each energy bin.
16-element energy in eV is returned as numpy.array instance.
>>> etbl = getEnergyE16()
>>> print(len(etbl))
16
>>> print(etbl[5])
57.0
'''
return np.array([7, 11, 17, 25, 38, 57, 86, 129, 193, 290,
435, 652, 978, 1467, 2200, 3300], dtype=float)
[docs]def getEnergy(nr_table):
''' Returns the list of the central energy in eV for each energy bin.
Corresponding energy bins (8) for the given table.
>>> etbl = getEnergy(1)
>>> print(len(etbl))
8
>>> print(etbl[6])
17.0
>>> etbl = getEnergy(3)
>>> print(etbl[1])
2200.0
'''
etbl = getEnergyE16()
if nr_table == 1:
etbl = etbl[8:0:-1]
elif nr_table == 2:
etbl = etbl[11:3:-1]
elif nr_table == 3:
etbl = etbl[15:7:-1]
else:
raise IndexError('Energy table for CENA should be 1, 2 or 3')
return etbl
[docs]def getRangeE16():
''' Returns the range of the energy.
Bounds are the average (in multiplication domain) of the neighboring bins.
Returned is 2x16 array.
>>> erng = getRangeE16()
>>> print(erng.shape)
(2, 16)
'''
ec = list(getEnergyE16())
ec.append(ec[15] * 1.5)
ec.insert(0, ec[0] / 1.5)
ec = np.array(ec)
elow = np.sqrt(ec[0:16] * ec[1:17])
ehigh = np.sqrt(ec[1:17] * ec[2:18])
return np.array([elow, ehigh])
[docs]def getDeltaE16():
''' Return the dE of the energy. Using average of boundaries as dE.
:returns: The ndarray of the dE.
:rtype: ``ndarray (shape=(16,))``
See also :meth:`getFwhmDeltaE16()`.
'''
erng = getRangeE16()
return erng[1] - erng[0]
[docs]def getFwhmDeltaE16():
''' Return the dE of the energy. Using FWHM as dE.
:returns: The ndarray of the dE.
:rtype: ``ndarray (shape=(16,))``
See also :meth:`getDeltaE16()`.
'''
erng = getFwhmRangeE16()
return erng[1] - erng[0]
[docs]def getBoundaryE16():
''' Return the boundary of the energy steps. np.array with 17 element.
'''
elow, ehigh = getRangeE16()
e = [e for e in elow]
e.append(ehigh[-1])
return np.array(e)
[docs]def getRange(nr_table):
''' Returns the list of the energy range for the given table.
Bounds are calculated by the average of the neighboring bins.
See also getRangeE16() funciton.
Returned is 2x8 array.
>>> erng = getRange(2)
>>> print(erng.shape)
(2, 8)
>>> print('%.1f' % erng[0, 3]) # Lower bound for the energy step 3
157.8
>>> print('%.1f' % erng[1, 3]) # Lower bound for the energy step 3
236.6
'''
er = getRangeE16()
if nr_table == 1:
etbl = er[:, 8:0:-1]
elif nr_table == 2:
etbl = er[:, 11:3:-1]
elif nr_table == 3:
etbl = er[:, 15:7:-1]
else:
raise IndexError('Energy table for CENA should be 1, 2 or 3')
return etbl
[docs]def getFwhmRangeE16():
''' Returns the list of the effective energy range of the given table.
Bounds are calculated from the FWHM (or effective) energy range.
According to the calibration report, (0.6-1.6) x E.
Returned is 2x16 array
'''
ec = getEnergyE16()
el = ec * 0.6
eh = ec * 1.6
return np.array([el, eh])
[docs]def getFwhmRange(nr_table):
''' Returns the list of the energy range for the given table.
Bounds are calculated by the average of the neighboring bins.
See also getRangeE16() funciton.
Returned is 2x8 array
>>> erng = getFwhmRange(2)
>>> print(erng.shape)
(2, 8)
>>> print('%.1f' % erng[0, 3]) # Lower bound for the energy step 3
115.8
>>> print('%.1f' % erng[1, 3]) # Lower bound for the energy step 3
308.8
'''
er = getFwhmRangeE16()
if nr_table == 1:
etbl = er[:, 8:0:-1]
elif nr_table == 2:
etbl = er[:, 11:3:-1]
elif nr_table == 3:
etbl = er[:, 15:7:-1]
else:
raise IndexError('Energy table for CENA should be 1, 2 or 3')
return etbl
[docs]def step_function(x, n):
''' A step function
:param x: Scalar or np.array.
:param n: Step function parameter.
:returns: 0 if x<n whereas 1 for x>=n.
'''
x2 = np.array(x)
mask = (x2 >= n)
x2 = np.ma.zeros(x2.shape)
x2.mask = mask
return x2.filled(1) * 1.
[docs]def response_simulation(incoming_energy, detecting_energies, normalize="peak"):
r''' Energy response function derived from Kazama's simulation.
According to Kazama [2006], the energy response may be
well fitted by triangle, while there may be a high energy tail.
:param incoming_energy: Incoming particle energy in ``eV``
:param detecting_energies: Energies to detect.
:type detecting_energies: ``iterable`` or ``float``
:keyword normalize: Way of normalize. "peak" or "integral".
:returns: The energy response. Normalized (thus this is relative).
The maximum of the response is set to 1 (\``normalize``=="peak")
or the integral over the energy (in eV) is set to 1 (\``normalize``=="integral").
See also the script in :mod:`cena_energy_resolution_kazama`.
For example, beam energy 25 eV is examined.
The response is maximum (=1) at 25 eV detecting energy
>>> print(response_simulation(25, 25))
1.0
The response falls down if it goes down to 12.5 eV or up to 63 eV.
>>> print(response_simulation(25, [12.5, 63.]))
[0. 0.]
The response is about half at 18.75 eV and 44. eV.
>>> print('%.2f' % response_simulation(25, 18.75))
0.50
>>> print('%.2f' % response_simulation(25, 44.))
0.49
If you normalize by the ``integral``, the integral response in energy in eV
should be 1.
>>> e = np.array([12.5, 18.75, 25., 44., 63.])
>>> r = response_simulation(25, e, normalize='integral')
>>> print('%.4f %.4f %.4f %.4f %.4f' % (r[0], r[1], r[2], r[3], r[4]))
0.0000 0.0200 0.0400 0.0197 0.0000
Numerical integral to be roughly 1.
>>> print('%.2f' % (
... r[1] * (e[2]-e[0])/2. + r[2] * (e[3]-e[1])/2. + r[3] * (e[4]-e[2])/2.))
1.00
'''
# Relative to incoming_energy
peak = 1.0
limit0 = 0.5
limit1 = 2.510 * np.exp(-1.5597e-4 * incoming_energy)
logger = logging.getLogger('response_simulation')
logger.setLevel(logging.WARN)
logger.debug('%g %g %g' % (limit0, peak, limit1))
# Relative to incoming_energy
e2 = np.ma.array(detecting_energies, dtype=np.float64) / incoming_energy
# Mask the zero efficiency
e2 = np.ma.masked_less(e2, limit0)
e2 = np.ma.masked_greater(e2, limit1)
ke = np.ma.array((((e2 - peak) / (peak - limit0) + 1) * (1. - step_function(e2, peak))
+ ((e2 - peak) / (peak - limit1) + 1) * step_function(e2, peak)))
if normalize == "integral":
ke_integral = (limit1 - limit0) / 2.0 * incoming_energy
ke /= ke_integral
return ke.filled(0) * 1. # * 1. is for scalar input to return as scalar
[docs]def response_simple_triangle(incoming_energy, detecting_energies, normalize="peak"):
peak = 1.0
limit0 = 0.2
limit1 = 2.2
# Relative to incoming_energy
e2 = np.ma.array(detecting_energies, dtype=np.float64) / incoming_energy
# Mask the zero efficiency
e2 = np.ma.masked_less(e2, limit0)
e2 = np.ma.masked_greater(e2, limit1)
ke = np.ma.array((((e2 - peak) / (peak - limit0) + 1) * (1. - step_function(e2, peak))
+ ((e2 - peak) / (peak - limit1) + 1) * step_function(e2, peak)))
if normalize == "integral":
ke_integral = (limit1 - limit0) / 2.0 * incoming_energy
ke /= ke_integral
return ke.filled(0) * 1. # * 1. is for scalar input to return as scalar