''' Response function. See also :ref:`geometricfactor`.
'''
import numpy as np
import matplotlib.pyplot as plt
[docs]def param_sigma(param_fwhm, param_p):
r''' Return the sigma for general gaussian.
As in :ref:`geometricfactor`, it is derived as
.. math::
\sigma = \frac{\mathrm{FWHM}}{2(2\ln 2)^\frac{1}{2p}}
>>> print('%.2f' % param_sigma(22.5, 10))
11.07
'''
x = 2 * np.log(2)
x = np.power(x, 0.5 / param_p)
x = x * 2
sigma = param_fwhm / x
return sigma
[docs]def genfunc_general_gaussian(param_sigma, param_p):
r''' Return the general gaussian function
It is in the form of:
.. math::
f(x) = \exp\Bigl(-\frac{1}{2}\Bigl|\frac{x}{\sigma}\Bigr|^{2p}\Bigr)
:param param_sigma: :math:`\sigma`.
:param param_p: :math:`p`.
:returns: A function, f(x).
'''
def gg(x):
return np.exp(-0.5 * np.power(np.abs(x / param_sigma), 2 * param_p))
return gg
[docs]def resp_phi1(angle_degrees):
''' Return the response function
This is outdated. (Too unrealistic response for MIMA)
'''
phi = angle_degrees
fwhm = 22.5
p = 1
## Calculate parameter sigma from fwhm and p.
sigma = param_sigma(fwhm, p)
print(sigma)
respfunc = genfunc_general_gaussian(sigma, p)
return respfunc(angle_degrees)
[docs]class MimaAzimResp:
''' Mex/IMA azimuthal response.
The data is taken from Fig 30 in Barabash et al. 2007.
>>> aresp = MimaAzimResp.get_instance()
>>> print('%.4f' % aresp(0))
0.9985
'''
def __str__(self):
return 'MimaAzimResp:Fig30/SB07/G3DATA/LINEXPL'
g3data = '''
-180 0
-17.9893238434 0
-16.0676156584 0
-15 0.0632480957563
-13.9857651246 0.133405488714
-13.024911032 0.223150840888
-12.0106761566 0.3259524123
-10.9964412811 0.423857356945
-9.98220640569 0.525026719434
-8.96797153025 0.611506201619
-8.00711743772 0.707780389484
-6.99288256228 0.792627662747
-5.97864768683 0.86604947355
-5.01779359431 0.928047274037
-4.05693950178 0.975355194219
-2.98932384342 1.0144991655
-1.97508896797 1.04711575324
-1.01423487544 1.06667612173
0 1.07970620239
1.06761565836 1.08130936845
2.02846975089 1.07148997634
2.98932384342 1.04045186823
4.00355871886 0.963710458141
5.07117437722 0.820047030077
5.97864768683 0.6159962283
6.99288256228 0.418471357928
8.00711743772 0.224210905402
9.02135231317 0.093606600862
10.0355871886 0.0446127424595
12.0106761566 0
180 0
'''
singleton_instance = None
def __init__(self):
self.dset = np.array([float(s) for s in self.g3data.split()]).reshape(31, 2)
from scipy import interpolate
self.dset[:, 1] /= self.dset[:, 1].max()
self.lin_intp = interpolate.interp1d(self.dset[:, 0], self.dset[:, 1])
[docs] @classmethod
def get_instance(cls):
if cls.singleton_instance is None:
cls.singleton_instance = MimaAzimResp()
return cls.singleton_instance
def __call__(self, angle_degrees):
return self.get_response(angle_degrees)
[docs] def get_response(self, angle_degrees):
return self.lin_intp(angle_degrees)
[docs]class MimaElevRespTrapezoid:
''' MEX/IMA elevation response approximated by trapezoid
The original function comes from :class:`MimaElevResp`,
No dependence on PACC.
>>> lresp = MimaElevRespTrapezoid.get_instance(-2000)
>>> print(lresp(np.array([-5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5]))) # doctest: +NORMALIZE_WHITESPACE
[0. 0.33333333 0.66666667 1. 1. 0.41176471
0. 0. 0. 0. 0. ]
'''
def __str__(self):
return 'MimaElevResp:Fig30/SB07/SIMPLE/TRAPEZOID/PACC%d' % self.pacc
singleton_instance = None
def __init__(self, pacc):
self.a = -5.
self.b = -2.
self.c = -1.
self.d = 0.7
self.pacc = pacc
[docs] @classmethod
def get_instance(cls, pacc):
if cls.singleton_instance is None:
cls.singleton_instance = MimaElevRespTrapezoid(pacc)
return cls.singleton_instance
def __call__(self, angle_degrees):
return self.get_response(angle_degrees)
[docs] def get_response(self, x):
a, b, c, d = self.a, self.b, self.c, self.d
return np.clip(np.min([(d-x)/(d-c), (x-a) / (b-a)], axis=0), 0, 1)
[docs]class MimaElevResp:
''' Mex/IMA elevation response.
The data is taken from Fig 30 in Barabash et al. 2007.
>>> elres0 = MimaElevResp(0)
>>> elres4 = MimaElevResp(-2000)
>>> print(elres4.get_response(-10))
0.0
>>> elres7 = MimaElevResp(-4000)
>>> print(elres7.get_response(-10))
0.0
'''
def __str__(self):
return 'MimaElevResp:Fig30/SB07/G3DATA/LINEXPL/PAC%d' % self.pacc
g3data = {0: np.array([-180,0,
-4.64,0,
-3.89074346485, 2.40146211151e-05,
-3.01121680669, 5.20958846334e-05,
-1.99878831465, 0.000172907114925,
-1.01671297658, 0.000234880863541,
-0.0013032155427, 3.44418182846e-05,
0.100, 0,
180,0]).reshape(9, 2),
-2000: np.array([-180,0,
-4.3,0,
-4.29434801085, 2.56694588368e-05,
-4.00506031186, 9.64229762167e-05,
-3.01001383849, 0.000374671105825,
-2.02092990313, 0.000745867860924,
-1.01179649788, 0.000610641996466,
-0.00737035947392, 0.000181583423095,
0.369489937127, 2.58643103805e-05,
0.5,0,
180,0]).reshape(11, 2),
-4000: np.array([-180,0,
-5.32,0,
-5.2945376091, 2.70382702455e-05,
-5.00537195036, 0.000105018428167,
-4.00218364879, 0.000438386959566,
-2.99639763329, 0.000898093975469,
-2.00604099245, 0.00126148335041,
-1.00602573742, 0.00125622136158,
-0.00280256720386, 0.000519408340055,
0.999130463208, 3.05828743002e-05,
1.05033506587, 2.46678093405e-05,
1.1,0,
180,0]).reshape(13, 2)
} # Note that the range should be, inpriniple, -90 and 90, but
# considering the future use as -180 and 180.
# This is because in the future, the offset may be added,
# and the maximum could be 180 degrees to be called.
singleton_instance = {0: None, -2000: None, -4000: None}
def __init__(self, pacc):
self.dset = self.g3data[pacc]
self.dset[:, 1] /= self.dset[:, 1].max()
self.pacc = pacc
from scipy import interpolate
self.lin_intp = interpolate.interp1d(self.dset[:, 0], self.dset[:, 1])
[docs] @classmethod
def get_instance(cls, pacc):
if cls.singleton_instance[pacc] is None:
cls.singleton_instance[pacc] = MimaElevResp(pacc)
return cls.singleton_instance[pacc]
def __call__(self, angle_degrees):
return self.get_response(angle_degrees)
[docs] def get_response(self, angle_degrees):
return self.lin_intp(angle_degrees)
[docs]class MimaEnerResp:
''' Energy response for a single channel.
I do not have good reference on it. So that my guess work is here.
Generalized Maxwellian with p = 2 and dE=7%.
See also :ref:`geometricfactor`.
'''
def __str__(self):
return 'MimaEnerResp:GENMAX/p=.2f/dE=%.3f' % (self.p, self.resol)
def __init__(self, Ec):
self.Ec = Ec
self.resol = 0.07
self.p = 2
self.sigma = param_sigma(Ec * self.resol, self.p)
self.respfunc = genfunc_general_gaussian(self.sigma, self.p)
[docs] def get_response(self, energy):
return self.respfunc(energy - self.Ec)
def __call__(self, energy):
return self.get_response(energy)
[docs]class MimaResponse:
''' IMA response for a "single" channel.
'''
def __init__(self, energy0, elev0, azim0, pacc=-4000):
''' Instance IMA response
:param energy0: The center energy in eV.
:param elev0: The elevation angle in degrees.
:param azim0: The azimuth angle in degrees.
:keyword pacc: Post acceleration in V. -4000, -2000, and 0 is supported
>>> ene = 1000.0
>>> elev0 = 15.0
>>> azim0 = -45.0
>>> respfunc = MimaResponse(ene, elev0, azim0)
>>> print('%.3f' % respfunc(ene, elev0 - 1.5, azim0)) # Maximum elevation hist at -1.5
0.999
Note, it was 0.996 before (using :class:`MimaElevResp`. Now using :class:`MimaElevRespTrapezoid` as default)
>>> res = respfunc([999, 1000, 1001], [13, 14, 15], [-46, -45, -44])
>>> print(res.shape)
(3,)
.. warning::
I have not yet considered plus and minus in the definition of elevatio/azimuth.
'''
self.eresp = MimaEnerResp(energy0) # Reference is energy0
# self.tresp = MimaElevResp.get_instance(pacc) # Reference is 0 degrees.
self.tresp = MimaElevRespTrapezoid.get_instance(pacc) # Reference is 0 degrees.
self.presp = MimaAzimResp.get_instance() # Reference is 0 degrees.
self.pacc = pacc
self.energy0 = energy0
self.elev0 = elev0
self.azim0 = azim0
[docs] def get_eresp(self, ener, normalized2keV=False):
''' Return the energy response (1D)
If normalized2keV is given, the energy efficiency is normalized by 2 keV g-factor values
(see :mod:`gfactor0` for the relative g-factor).
This means that the maximul total efficiency other than 2keV becomes not 1.
'''
if normalized2keV:
from . import gfactor0
if self.pacc == 0:
effic = gfactor0.GL_CR4.Proton.PAC0()
elif self.pacc == -2000:
effic = gfactor0.GL_CR4.Proton.PAC1()
elif self.pacc == -4000:
effic = gfactor0.GL_CR4.Proton.PAC2()
else:
raise ValueError('No such pacc %g' % self.pacc)
effic2000 = effic(2000.)
effic_value = effic(ener) / effic2000
else:
effic_value = 1.
return self.eresp(np.array(ener)) * effic_value
[docs] def get_tresp(self, theta, cos_corr=False):
''' Return the elevation response (1D)
'''
if cos_corr:
cospart = np.cos(np.deg2rad(theta))
else:
cospart = 1
return self.tresp(np.array(theta) - self.elev0) * cospart
[docs] def get_presp(self, phi):
''' Return the azimuth (phi) response (1D)
'''
phi2 = np.array(phi) - self.azim0
while phi2.min() < -180.:
phi2 = np.where(phi2 < -180., phi2 + 360, phi2)
while phi2.max() > 180.:
phi2 = np.where(phi2 > 180., phi2 - 360, phi2)
return self.presp(phi2)
[docs] def get_response(self, ener, theta, phi, normalized2keV=False, cos_corr=False):
''' Get the responce (3D).
If ``cos_corr`` is given, the theta efficiency is multiplied by cos(theta).
Usually, the final integral is over the "omgea", so that domega=cos(theta)dtheta dphi.
If ``cos_corr`` is set, the integral may be done over theta domain.
Use this option only if you understand it.
'''
# print self.get_eresp(ener)
# print self.get_tresp(theta)
# print self.get_presp(phi)
return self.get_eresp(ener, normalized2keV=normalized2keV) * self.get_tresp(theta, cos_corr=cos_corr) * self.get_presp(phi)
[docs] def get_meshed_response(self, ene1d, theta1d, phi1d, normalized2keV=False, cos_corr=False):
''' Get the response (3D).
This provides the same interface to :meth:`get_response`, but faster.
'''
r_ene = self.get_eresp(ene1d, normalized2keV=normalized2keV)
r_theta = self.get_tresp(theta1d, cos_corr=cos_corr)
r_phi = self.get_presp(phi1d)
r_ene, r_theta, r_phi = np.meshgrid(r_ene, r_theta, r_phi, copy=False) # Shape is (ntheta, nene, nphi)
resp = r_ene * r_theta * r_phi
return resp.transpose((1, 0, 2))
def __call__(self, *args, **kwds):
return self.get_response(*args, **kwds)
def __str__(self):
return 'MEX/IMA efficiency\n' + '\n'.join([str(self.eresp), str(self.tresp), str(self.presp)])
import unittest
import doctest
[docs]def doctests():
return unittest.TestSuite((
doctest.DocTestSuite(),
))
if __name__ == '__main__':
unittest.main(defaultTest='doctests')