''' Calculate mass ring number from energy and mass.
The way of calculating the mass ring number from energy and mass
is described in the MEXImaTables.pdf file.
This module implements that.
.. autosummary::
rm
dm
massline
Users may use :func:`rm` function to calculate the mass ring number
from energy and mass together with PACC index.
The mass ring number looks to have 0-based index.
It may say that the range of mass ring number could be (-0.5, 31.5).
>>> print('{:.2f}'.format(rm(100, 32, 4))) # Mass ring number for 100 eV/q, 32 amu/q particle with PACC=4.
10.34
:func:`dm` will give the width, although it is almost always 1.2.
Users also want to know the iso-mass line. This will be done by
:func:`massline` function.
>>> energies, masses = massline(32, 7, enestep=np.logspace(1, 4, 96))
>>> print('{e:.2f} {m:.2f}'.format(e=energies[15], m=masses[15]))
29.76 7.84
See also the example in :mod:`snippet_mima.mima_masscontour`.
*Details*: The :attr:`table` contains the table needed for calculation.
It is 8-element list, but only used for 0, 4, and 7th element.
Each corresponds to the table for each PACC.
The table is (5, 12) shaped. Contents is in :func:`_load_table6` function.
>>> print(table[0].shape)
(5, 12)
>>> print(table[3])
None
'''
import numpy as np
from scipy.interpolate import InterpolatedUnivariateSpline
def _load_table6():
''' Load the table 6. This is constant over the mission.
>>> tbl = _load_table6()
>>> print(tbl.shape)
(15, 12)
:returns: A table. (15, 12 shaped).
[0:5, :] is for low pacc, [5:10, :] is for mid pac,
and [10:15, :] is for high pac.
For column, refer to the table below.
Column
::
0: PaccIndex
1: Pacc
2: Elimit
3: #
4: GfitP0
5: GfitP1
6: GfitD0
7: GfitD1
8: Kpacc
9: Kmass
10: Dmass
11: Omass
'''
tbl = '''0 300. 200. 0 1.913E+01 -2.389E+00 1.200E+00 1.200E+00 2.737E-01 1.000E+00 -5.000E+01 1.000E+00
0 300. 200. 1 2.479E-01 2.392E+00 5.416E-16 -1.203E-16 1.517E+00 2.000E+00 -5.000E+01 2.000E+00
0 300. 200. 2 -1.246E-02 -2.307E-02 -2.604E-17 -4.879E-19 -7.909E-01 1.700E+01 -4.900E+00 1.600E+01
0 300. 200. 3 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 2.800E+01 0.000E+00 3.200E+01
0 300. 200. 4 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 3.900E+01 1.500E+00 5.000E+01
1 2433. 340. 0 1.945E+01 -3.707E+00 1.200E+00 1.200E+00 1.067E-01 1.000E+00 -4.000E+01 1.000E+00
1 2433. 340. 1 -1.134E+00 2.117E+00 -2.996E-16 8.844E-17 1.957E+00 2.000E+00 -2.500E+01 2.000E+00
1 2433. 340. 2 8.625E-03 1.202E-02 1.765E-17 -1.240E-17 -9.828E-01 1.700E+01 -3.500E+00 1.600E+01
1 2433. 340. 3 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 2.600E+01 0.000E+00 3.200E+01
1 2433. 340. 4 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 3.300E+01 1.500E+00 5.000E+01
2 4216. 499. 0 1.132E+01 -1.500E+00 1.200E+00 1.200E+00 3.594E-02 1.000E+00 -3.000E+01 1.000E+00
2 4216. 499. 1 -4.321E-01 1.647E+00 -1.451E-16 4.208E-17 3.741E-01 1.800E+00 -2.200E+01 2.000E+00
2 4216. 499. 2 1.041E-02 -1.065E-02 2.290E-17 1.632E-17 -1.731E-01 1.700E+01 -3.600E+00 1.600E+01
2 4216. 499. 3 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 3.200E+01 0.000E+00 3.200E+01
2 4216. 499. 4 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 4.500E+01 1.500E+00 5.000E+01'''
tbl = np.fromstring(tbl, sep=' ').reshape(15, 12)
return tbl
def _load_table6_pacc(pi):
''' Return the table 6 contents according to PACC.
:param pi: Pacc index. 0, 4 or 7.
:returns: (5, 12) shaped array. See :func:`load_table6` for contents.
'''
tbl = _load_table6()
if pi == 0:
tbl = tbl[:5, :]
elif pi == 4:
tbl = tbl[5:10, :]
elif pi == 7:
tbl = tbl[10:, :]
else:
raise ValueError('PI={} not used.'.format(pi))
return tbl
table = [_load_table6_pacc(0), None, None, None,
_load_table6_pacc(4), None, None,
_load_table6_pacc(7)]
''' Table for Mass Ring calculation.
An list of numpy array with a shape of (5,12).
This is eight elements, each corresponds to PACC index.
>>> print(table[0].shape)
(5, 12)
>>> print(table[3])
None
'''
from irfpy.mima import energy as _energy
def _meff(m_per_q, pi):
''' Calculate Meff.
Meff = Lin(M/Q, Kmass_i(PI), Omass_i(PI)
Here Lin(X, Ytab, Xtab) is the linear interpolation of X on the tabulated Xtab and Ytab.
:param m_per_q: Mass per charge. Float or array.
:param pi: Pacc index. 0, 4 or 7.
:returns: Meff value
:rtype: np.array (that's why one need float() for the later examples)
>>> print(_meff(1, 0), _meff(1, 4), _meff(1, 7))
1.0 1.0 1.0
>>> print(_meff(16, 0), _meff(16, 4), _meff(16, 7))
17.0 17.0 17.0
>>> print(_meff(32, 0), _meff(32, 4), _meff(32, 7))
28.0 26.0 32.0
>>> print('{:.3f} {:.3f} {:.3f}'.format(float(_meff(44, 0)), float(_meff(44, 4)), float(_meff(44, 7))))
35.333 30.667 40.667
>>> print('{:.3f}'.format(float(_meff(128, 4))))
63.333
'''
tbl = table[pi]
kmass = tbl[:, 9]
omass = tbl[:, 11]
f = InterpolatedUnivariateSpline(omass, kmass, k=1) # k=1 means linear
return f(m_per_q)
def _pacceff(m_per_q, pi):
r''' Calculate Pacceff.
.. math::
Pacc_{eff} = Pacc(PI) \cdot (Kpacc0(PI) + Kpacc1(PI) / Meff + Kpacc2(PI) / Meff ^ 2)
>>> peff = _pacceff(12.8, 0) # M/q = 12.8 for PACC=0
>>> print('{:.2f}'.format(float(peff)))
114.36
>>> peff = _pacceff(37.6, 4) # M/q = 37.6 for PACC=4
>>> print('{:.2f}'.format(float(peff)))
425.57
>>> peff = _pacceff(3.75, 7) # M/q = 3.75 for PACC=7
>>> print('{:.2f}'.format(float(peff)))
524.49
'''
tbl = table[pi]
pacc = tbl[0, 1] # PACC in the table (tbl[1, :] is all same)
kpacc0 = tbl[0, 8]
kpacc1 = tbl[1, 8]
kpacc2 = tbl[2, 8]
Meff = _meff(m_per_q, pi)
coeffic = kpacc0 + kpacc1 / Meff + kpacc2 / Meff / Meff
return pacc * coeffic
def _geff(e_per_q, m_per_q, pi):
r''' Geff calculation.
.. math::
G_{eff} = 10^3 / \sqrt{(E/Q + Pacc_{eff})\cdot M_{eff}}
:param e_per_q: Given by the user.
Technical note shows it as the function of energy index,
it is complicated to get the energy from the index,
so that this is specified by the user.
Pacceff and Meff is obtained from the other functions (
:func:`_pacceff` and :func:`_meff`).
>>> Geff = _geff(1200, 12.8, 0) # M/q = 12.8 for PACC=0
>>> print('{:.4f}'.format(float(Geff)))
7.4874
>>> Geff = _geff(2400, 37.6, 4) # M/q = 37.6 for PACC=4
>>> print('{:.4f}'.format(float(Geff)))
3.5440
>>> Geff = _geff(750, 3.75, 7) # M/q = 3.75 for PACC=7
>>> print('{:.4f}'.format(float(Geff)))
14.5624
'''
PaccEff = _pacceff(m_per_q, pi)
Meff = _meff(m_per_q, pi)
geff = 1e3 / np.sqrt((e_per_q + PaccEff) * Meff)
return geff
def _glim(m_per_q, pi):
r''' Glim
.. math::
G_{lim} = 10^3 / \sqrt{(Elimit(PI) + Pacc_{eff})\cdot(M_{eff})}
>>> Glim = _glim(12.8, 0) # M/q = 12.8 for PACC=0
>>> print('{:.3f}'.format(float(Glim)))
15.310
>>> Glim = _glim(37.6, 4) # M/q = 37.6 for PACC=4
>>> print('{:.3f}'.format(float(Glim)))
6.809
>>> Glim = _glim(3.75, 7) # M/q = 3.75 for PACC=7
>>> print('{:.3f}'.format(float(Glim)))
16.250
'''
## Elimit is obtained from table [0, 2]
tbl = table[pi]
elimit = tbl[0, 2]
Meff = _meff(m_per_q, pi)
PaccEff = _pacceff(m_per_q, pi)
glim = 1e3 / np.sqrt((elimit + PaccEff) * Meff)
return glim
[docs]def rm(e_per_q, m_per_q, pi):
''' Return Rm.
>>> Rm = rm(1200, 12.8, 0) # M/q = 12.8 for PACC=0
:param e_per_q: Energy per charge (in eV/q)
:param m_per_q: Mass per charge (in amu/q)
:param pi: Pacc index, 0, 4 or 7.
:returns: Mass ring number.
>>> print('{:.2f}'.format(float(Rm)))
14.23
>>> Rm = rm(2400, 37.6, 4) # M/q = 37.6 for PACC=4
>>> print('{:.2f}'.format(float(Rm)))
3.95
>>> Rm = rm(750, 3.75, 7) # M/q = 3.75 for PACC=7
>>> print('{:.2f}'.format(float(Rm)))
20.23
>>> Rm = rm(100, 12.8, 0) # M/q = 12.8 for PACC=0
>>> print('{:.2f}'.format(float(Rm)))
28.26
>>> Rm = rm(250, 37.6, 4) # M/q = 37.6 for PACC=4
>>> print('{:.2f}'.format(float(Rm)))
10.82
>>> Rm = rm(350, 3.75, 7) # M/q = 3.75 for PACC=7
>>> print('{:.2f}'.format(float(Rm)))
22.35
'''
tbl = table[pi]
elimit = tbl[0, 2]
if e_per_q < elimit:
return _rm_low(e_per_q, m_per_q, pi)
else:
return _rm_high(e_per_q, m_per_q, pi)
[docs]def massline(m_per_q, pi, enestep=None):
""" Return the iso-mass contour.
For the given mass and post acceleration,
the mass ring numbers as a function of energies are
returned.
For normal data analysis, use :func:`central_mass_channel` function.
:param m_per_q: Mass per charge (amu/q)
:param pi: Post acceleration index. 0, 4, or 7.
:keyword enestep: Energy step. Default is logspace(0, 4, 96).
:return: (enestep, massring) with each element np.array
with shape of given enestep.
"""
if enestep is None:
enestep = np.logspace(0, 4, 96)
ml = np.zeros_like(enestep)
for i in range(len(enestep)):
ml[i] = rm(enestep[i], m_per_q, pi)
return np.array(enestep), np.array(ml)
__central_mass_channel = {}
[docs]def central_mass_channel(m_per_q, promsection, paccindex, fill_nan_invalid_energy=True):
""" Return the central mass channel, Rm, for given PROM section and PACC index.
:param m_per_q: Mass per charge in amu
:param promsection: PROM section, 0 to 16
:param paccindex: PACC index, 0 to 7
:keyword fill_nan: If *True*, the step for invalid energy step is filled by nan.
:return: Central mass channel for given m_per_q. 96 steps corresponding to the energy table.
:rtype: Tuple
>>> m1curve = central_mass_channel(1, 15, 4) # Proton for prom=15 and pacc=4
>>> print(f'{m1curve[0]:.2f}')
12.92
>>> m16curve = central_mass_channel(32, 16, 4) # OO+ for prom-16 and pacc=4
>>> print(f'{m16curve[30]:.2f}')
7.18
>>> print(m16curve[90]) # The non-valid energy step corresponds to nan
nan
If you do not want to fill nan for all the energy step,
>>> m16curve_float = central_mass_channel(32, 16, 4, fill_nan_invalid_energy=False) # Dont fill nan for any energy step
>>> print(f'{m16curve_float[90]:.2f}') # Value shall be used with care!
9.44
"""
prms = (m_per_q, promsection, paccindex, fill_nan_invalid_energy)
if prms in __central_mass_channel:
rm = __central_mass_channel[prms]
else:
etable = _energy.get_default_table(promsection, keep_negative=fill_nan_invalid_energy) # Nagative to be nan.
etable = np.where(etable >= 0, etable, np.nan)
rm = tuple(massline(m_per_q, paccindex, enestep=etable)[1])
__central_mass_channel[prms] = rm
return rm
[docs]def dm0(*args):
""" Simple version of :func:`dm`. It just returns 1.2.
:param args: Whatever.
:return: 1.2
"""
return 1.2
[docs]def dm(e_per_q, m_per_q, pi):
''' Return Dm.
:param e_per_q: Energy per charge (in eV/q)
:param m_per_q: Mass per charge (in amu/q)
:param pi: Pacc index, 0, 4 or 7.
:returns: Mass ring number.
Dm is almost always 1.2 (as fitting parameters on Geff dependence is 1e-16 or less..).
I wonder if there are some errors in the calibration document, I could not find them.
So, a simple way to use dm is just use :func:`dm0`, which always returns ``1.2``.
>>> Dm = dm(1200, 12.8, 0) # M/q = 12.8 for PACC=0
>>> print('{:.2f}'.format(float(Dm)))
1.20
>>> Dm = dm(2400, 37.6, 4) # M/q = 37.6 for PACC=4
>>> print('{:.2f}'.format(float(Dm)))
1.20
>>> Dm = dm(750, 3.75, 7) # M/q = 3.75 for PACC=7
>>> print('{:.2f}'.format(float(Dm)))
1.20
>>> Dm = dm(100, 12.8, 0) # M/q = 12.8 for PACC=0
>>> print('{:.2f}'.format(float(Dm)))
1.20
>>> Dm = dm(250, 37.6, 4) # M/q = 37.6 for PACC=4
>>> print('{:.2f}'.format(float(Dm)))
1.20
>>> Dm = dm(350, 3.75, 7) # M/q = 3.75 for PACC=7
>>> print('{:.2f}'.format(float(Dm)))
1.20
'''
tbl = table[pi]
elimit = tbl[0, 2]
if e_per_q < elimit:
return _dm_low(e_per_q, m_per_q, pi)
else:
return _dm_high(e_per_q, m_per_q, pi)
def _rm_high(e_per_q, m_per_q, pi):
r''' Return Rm when the energy is higher than limit.
.. math::
Rm = GfitP1_0(PI) + GfitP1_i(PI)\cdot G_{eff} + GfitP1_2(PI)\cdot G_{eff}^2
'''
tbl = table[pi] # (5, 12) shape
gfp10 = tbl[0, 5]
gfp11 = tbl[1, 5]
gfp12 = tbl[2, 5]
Geff = _geff(e_per_q, m_per_q, pi)
Rm = gfp10 + gfp11 * Geff + gfp12 * Geff * Geff
return Rm
def _rm_low(e_per_q, m_per_q, pi):
r''' Return Rm for E<Elimit
.. math::
dM = GfitP0_0(PI) + GfitP0_1(PI)\cdot G_{lim} + GfitP0_2(PI)\cdot G_{lim}^2 - (GfitP1_0(PI) + GfitP1_1(PI)\cdot G_{lim} + GfitP1_2(PI)\cdot G_{lim}^2 ) \\
Rm = GfitP0_0(PI) + GfitP0_1(PI)\cdot G_{eff} + GfitP0_2(PI)\cdot G_{eff}^2 - dM
'''
tbl = _load_table6_pacc(pi) # (5, 12) array
gfp00 = tbl[0, 4]
gfp01 = tbl[1, 4]
gfp02 = tbl[2, 4]
gfp10 = tbl[0, 5]
gfp11 = tbl[1, 5]
gfp12 = tbl[2, 5]
Glim = _glim(m_per_q, pi)
Geff = _geff(e_per_q, m_per_q, pi)
dM = gfp00 + gfp01 * Glim + gfp02 * Glim * Glim - (gfp10 + gfp11 * Glim + gfp12 * Glim * Glim)
Rm = gfp00 + gfp01 * Geff + gfp02 * Geff * Geff - dM
return Rm
def _dm_high(e_per_q, m_per_q, pi):
r''' Dm for high energy part.
.. math::
Dm = GfitD1_0(PI) + GfitD1_1(PI)\cdot G_{eff} + GfitD1_2(PI)\cdot G_{eff}^2
'''
tbl = table[pi]
gfd10 = tbl[0, 7]
gfd11 = tbl[1, 7]
gfd12 = tbl[2, 7]
Geff = _geff(e_per_q, m_per_q, pi)
dm = gfd10 + gfd11 * Geff + gfd12 * Geff * Geff
return dm
def _dm_low(e_per_q, m_per_q, pi):
r''' Dm for low energy part.
.. math::
Dm = GfitD0_0(PI) + GfitD0_1(PI)\cdot G_{eff} + GfitD0_2(PI)\cdot G_{eff}^2
'''
tbl = table[pi]
gfd10 = tbl[0, 6]
gfd11 = tbl[1, 6]
gfd12 = tbl[2, 6]
Geff = _geff(e_per_q, m_per_q, pi)
dm = gfd10 + gfd11 * Geff + gfd12 * Geff * Geff
return dm