# Source code for irfpy.vels.bible.gfactor

```''' Module for g-factor.

G-factor, well, not simple.

But for simplicity, use :func:`gfactor`.

>>> print('%.2e' % gfactor(3021, 5))   # Return the g-factor of 3021 eV for CH-5
7.40e-06
>>> print('%.2e' % gfactor(3021, 10))   # Return the g-factor of 3021 eV for CH-5
8.09e-06

The unit is ``cm2 sr eV/eV``.

The module refers to table 6 in :ref:`vexelsbible` (see :func:`lab_gfactor`).
They are fitted using a linear function, and the fitted function will be returned from :func:`gfactor_func`.
The returned function is the function to the energy.
:func:`gfactor` is the simple way of getting the g-factor for the specified energy and the anode.
This is the way that the :ref:`vexelsbible` recommends.
'''

import numpy as np
from scipy.stats import linregress

[docs]def lab_gfactor():
''' Return 4 energy and 4 x 16 element np.array g-factor cm2 sr eV/eV.  Table 6 in :ref:`vexelsbible`.
'''

g199 = np.array([8.85756e-6, 9.63588e-6, 9.47093e-6, 8.70769e-6, 8.25067e-6, 8.07105e-6, 7.23690e-6, 8.70944e-6,
8.51592e-6, 9.39149e-6, 8.33294e-6, 1.03845e-5, 1.04415e-5, 1.14662e-5, 1.17197e-5, 1.37654e-5,])

g970 = np.array([1.08817e-5, 1.07470e-5, 9.99465e-6, 9.52572e-6, 9.26740e-6, 8.73019e-6, 8.07841e-6, 9.17226e-6,
9.72213e-6, 1.06782e-5, 9.62592e-6, 1.12857e-5, 1.29351e-5, 1.42535e-5, 1.44310e-5, 1.49384e-5,])

g3021 = np.array([8.04262e-6, 8.32344e-6, 8.43561e-6, 8.44829e-6, 8.10313e-6, 7.78427e-6, 7.01710e-6, 7.42745e-6,
8.29093e-6, 9.33263e-6, 8.62126e-6, 9.49013e-6, 1.00476e-5, 1.07045e-5, 1.06524e-5, 1.04071e-5,])

g5994 = np.array([6.94579e-6, 7.52185e-6, 7.54510e-6, 7.37566e-6, 6.92202e-6, 6.40499e-6, 5.74396e-6, 5.87616e-6,
6.61914e-6, 6.99206e-6, 6.89954e-6, 6.46349e-6, 9.05845e-6, 9.34059e-6, 8.89006e-6, 8.82440e-6,])

g = np.array([g199, g970, g3021, g5994])

e = np.array([199.990, 970., 3021., 5994.])

return (e, g)

[docs]def lab_loggfactor():
e, g = lab_gfactor()
return (np.log10(e), np.log10(g))

[docs]def gfactor_fit_func(x, y):
fit1 = linregress(x, y)
p0 = fit1[1]
p1 = fit1[0]
return p0, p1

[docs]def gfactor_func(anode, above=True):
''' Return the g-factor as a function of energy.

You can get g-factor from this function.
The function returns the g-factor of the specified anode.

:param anode: The anode number
:returns: A function of the g-factor in the unit of cm2 sr eV/eV
:rtype: ``function``

>>> g0 = gfactor_func(0)   # Return the function of g-factor of CH-0
>>> print('%.1e' % g0(1000))   # 1 keV g-factor.
1.1e-05

According to the :ref:`vexelsbible`, the g-factor should be get
via the fitting of the laboratory data (:func:`lab_gfactor`).

The method will return the function (energy dependence) of the g-factor.
'''

e, g = lab_loggfactor()
loge = e
logg = g[:, anode]

# Two functions fitted for above and below 970 eV (according to g-factor used by Rudy Frahm)
if not above:
p0, p1 = gfactor_fit_func(loge[:2], logg[:2])
fitfunc = lambda x: np.power(10, p0 + np.log10(x) * p1)
return fitfunc
if above:
p0, p1 = gfactor_fit_func(loge[1:], logg[1:])
fitfunc = lambda x: np.power(10, p0 + np.log10(x) * p1)
return fitfunc

__gfactor_funcs_above970 = None
__gfactor_funcs_below970 = None

[docs]def gfactor(energy, anode):
''' Simplest g-factor accessor. Return the g-factor in the unit of ``cm2 sr eV/eV``.

:param energy: Energy in eV.
:param anode: Anode id, 0 to 15.
:returns: G-factor in the unit of ``cm2 sr eV/eV``

For example, you can ge the g-factor for 10, 100 and 1000 eV for anode 3 as follows.

>>> print('%.3e' % gfactor(10, 3))
7.344e-06
>>> print('%.3e' % gfactor(100, 3))
8.371e-06
>>> print('%.3e' % gfactor(1000, 3))
9.579e-06
'''
# For performance.
global __gfactor_funcs_above970
global __gfactor_funcs_below970

if energy > 970:
if __gfactor_funcs_above970 is None:
__gfactor_funcs_above970 = [gfactor_func(ian, above=True) for ian in range(16)]
return __gfactor_funcs_above970[anode](energy)
else:
if __gfactor_funcs_below970 is None:
__gfactor_funcs_below970 = [gfactor_func(ian, above=False) for ian in range(16)]
return __gfactor_funcs_below970[anode](energy)

[docs]def gfactor1d(energies, anode):
""" Simplest g-factor accessor, with multiple dimensional version, in cm2 sr eV/eV.

:param energies: Energies in eV. Any shape of array
:param anode: Anodes. 0 to 15 (inclusive), integer.
:returns: Gfactor.

Gfactor returned also have the same shape as energies.

>>> g = gfactor1d([10., 100., 1000.], 3)    # Gfactor of 10, 100 and 1000 eV for anode 3
>>> print('%.3e' % g[0])
7.344e-06
>>> print('%.3e' % g[1])
8.371e-06
>>> print('%.3e' % g[2])
9.579e-06
"""
energies = np.array(energies)

# For performance.
global __gfactor_funcs_above970
global __gfactor_funcs_below970

if __gfactor_funcs_above970 is None:
__gfactor_funcs_above970 = [gfactor_func(ian, above=True) for ian in range(16)]
if __gfactor_funcs_below970 is None:
__gfactor_funcs_below970 = [gfactor_func(ian, above=False) for ian in range(16)]

g_a = __gfactor_funcs_above970[anode](energies)
g_b = __gfactor_funcs_below970[anode](energies)

g = np.where(energies > 970, g_a, g_b)

return g
```