import numpy as np
from NuRadioReco.utilities import units
from scipy.interpolate import interp1d
from scipy import constants
import logging
logger = logging.getLogger("cross sections")
[docs]def param(energy, inttype='cc', parameterization='ctw'):
"""
Parameterization and constants as used in
get_nu_cross_section()
See documentation there for details
"""
if(np.any(energy < 1e4 * units.GeV)):
logger.warning(f"CTW / BGR neutrino nucleon cross sections not valid for energies below 1e4 GeV, ({energy/units.GeV}GeV was requested)")
if(hasattr(energy, "__len__")):
return np.nan * np.ones_like(energy)
else:
return np.nan
if parameterization == 'ctw':
"""
Phys.Rev.D83:113009,2011 Amy Connolly, Robert S. Thorne, David Waters
"""
if inttype == 'cc':
c = (-1.826, -17.31, -6.406, 1.431, -17.91) # nu, CC
elif inttype == 'nc':
c = (-1.826, -17.31, -6.448, 1.431, -18.61) # nu, NC
elif inttype == 'cc_bar':
c = (-1.033, -15.95, -7.247, 1.569, -17.72) # nu_bar, CC
elif inttype == 'nc_bar':
c = (-1.033, -15.95, -7.296, 1.569, -18.30) # nu_bar, NC
elif inttype == 'nc_up':
c = (-1.456, 32.23, -32.32, 5.881, -49.41) # nu, NC
elif inttype == 'cc_up':
c = (-1.456, 33.47, -33.02, 6.026, -49.41) # nu, CC
elif inttype == 'nc_bar_up':
c = (-2.945, 143.2, -76.70, 11.75, -142.8) # nu_bar, NC
elif inttype == 'cc_bar_up':
c = (-2.945, 144.5, -77.44, 11.9, -142.8) # nu_bar, CC
elif inttype == 'nc_down':
c = (-15.35, 16.16, 37.71, -8.801, -253.1) # nu, NC
elif inttype == 'cc_down':
c = (-15.35, 13.86, 39.84, -9.205, -253.1) # nu, CC
elif inttype == 'nc_bar_down':
c = (-13.08, 15.17, 31.19, -7.757, -216.1) # nu_bar, NC
elif inttype == 'cc_bar_down':
c = (-13.08, 12.48, 33.52, -8.191, -216.1) # nu_bar, CC
else:
logger.error("Type {0} of interaction not defined for 'ctw'".format(inttype))
raise NotImplementedError
elif parameterization == 'hedis_bgr18':
"""
Parameterization as above fitted to GENIE HEDIS module (with BGR18) cross sections
as in arXiv:2004.04756v2 (prepared for JCAP)
Precalculated xsec tables for 'nu_mu(_bar)_O16/tot_cc(nc)'/16 for isoscalar target
obtained from GHE19_00a_00_000.root in https://github.com/pochoarus/GENIE-HEDIS/tree/nupropearth/genie_xsec
Fitted in the energy range above 1 TeV; do not use below
"""
if inttype == 'cc':
c = (-1.6049779136562436, -17.7480299104706, -6.748861524562085, 1.5569481852252935, -16.545379184836094) # nu, CC
elif inttype == 'nc':
c = (-1.9625311094497564, -17.576550328008224, -6.444583672267122, 1.4702739736023922, -18.6167800243672) # nu, NC
elif inttype == 'cc_bar':
c = (-2.28879962998228, -15.725804320703244, -5.273935123272873, 1.0314821502761589, -23.15773837113476) # nu_bar, CC
elif inttype == 'nc_bar':
c = (-2.582585867636026, -15.742658435090945, -5.075692336968196, 0.9963850387362603, -24.870843546539973) # nu_bar, NC
else:
logger.error("Type {0} of interaction not defined for 'hedis_bgr18'".format(inttype))
raise NotImplementedError
else:
logger.error("Parameterization {0} of interaction cross section not defined".format(parameterization))
raise NotImplementedError
epsilon = np.log10(energy / units.GeV)
l_eps = np.log(epsilon - c[0])
crscn = c[1] + c[2] * l_eps + c[3] * l_eps ** 2 + c[4] / l_eps
crscn = np.power(10, crscn) * units.cm ** 2
return crscn
[docs]def csms(energy, inttype, flavors):
"""
Neutrino cross sections according to
Amanda Cooper-Sarkar, Philipp Mertsch, Subir Sarkar
JHEP 08 (2011) 042
"""
if type(inttype) == str:
inttype = np.array([inttype] * energy.shape[0])
if isinstance(flavors, (int, np.integer)):
flavors = np.array([flavors] * energy.shape[0])
neutrino = np.array((
[50, 0.32, 0.10],
[100, 0.65, 0.20],
[200, 1.3, 0.41],
[500, 3.2, 1.0],
[1000, 6.2, 2.0],
[2000, 12., 3.8],
[5000, 27., 8.6],
[10000, 47., 15.],
[20000, 77., 26.],
[50000, 140., 49.],
[100000, 210., 75.],
[200000, 310., 110.],
[500000, 490., 180.],
[1e6, 690., 260.],
[2e6, 950., 360.],
[5e6, 1400., 540.],
[1e7, 1900., 730.],
[2e7, 2600., 980.],
[5e7, 3700., 1400.],
[1e8, 4800., 1900.],
[2e8, 6200., 2400.],
[5e8, 8700., 3400.],
[1e9, 11000., 4400.],
[2e9, 14000., 5600.],
[5e9, 19000., 7600.],
[1e10, 24000., 9600.],
[2e10, 30000., 12000.],
[5e10, 39000., 16000.],
[1e11, 48000., 20000.],
[2e11, 59000., 24000.],
[5e11, 75000., 31000.]
))
neutrino[:, 0] *= units.GeV
neutrino[:, 1] *= units.picobarn # CC
neutrino[:, 2] *= units.picobarn # NC
neutrino_cc = interp1d(neutrino[:, 0], neutrino[:, 1], bounds_error=True)
neutrino_nc = interp1d(neutrino[:, 0], neutrino[:, 2], bounds_error=True)
antineutrino = np.array((
[50, 0.15, 0.05],
[100, 0.33, 0.12],
[200, 0.69, 0.24],
[500, 1.8, 0.61],
[1000, 3.6, 1.20],
[2000, 7., 2.4],
[5000, 17., 5.8],
[10000, 31., 11.],
[20000, 55., 19.],
[50000, 110., 39.],
[100000, 180., 64.],
[200000, 270., 99.],
[500000, 460., 170.],
[1e6, 660., 240.],
[2e6, 920., 350.],
[5e6, 1400., 530.],
[1e7, 1900., 730.],
[2e7, 2500., 980.],
[5e7, 3700., 1400.],
[1e8, 4800., 1900.],
[2e8, 6200., 2400.],
[5e8, 8700., 3400.],
[1e9, 11000., 4400.],
[2e9, 14000., 5600.],
[5e9, 19000., 7600.],
[1e10, 24000., 9600.],
[2e10, 30000., 12000.],
[5e10, 39000., 16000.],
[1e11, 48000., 20000.],
[2e11, 59000., 24000.],
[5e11, 75000., 31000.]
))
antineutrino[:, 0] *= units.GeV
antineutrino[:, 1] *= units.picobarn # CC
antineutrino[:, 2] *= units.picobarn # NC
antineutrino_cc = interp1d(antineutrino[:, 0], antineutrino[:, 1], bounds_error=True)
antineutrino_nc = interp1d(antineutrino[:, 0], antineutrino[:, 2], bounds_error=True)
crscn = np.zeros_like(energy)
particles_cc = np.where((flavors >= 0) & (inttype == 'cc'))
particles_nc = np.where((flavors >= 0) & (inttype == 'nc'))
antiparticles_cc = np.where((flavors < 0) & (inttype == 'cc'))
antiparticles_nc = np.where((flavors < 0) & (inttype == 'nc'))
#
crscn[particles_cc] = neutrino_cc(energy[particles_cc])
crscn[particles_nc] = neutrino_nc(energy[particles_nc])
crscn[antiparticles_cc] = antineutrino_cc(energy[antiparticles_cc])
crscn[antiparticles_nc] = antineutrino_nc(energy[antiparticles_nc])
return crscn
[docs]def get_nu_cross_section(energy, flavors, inttype='total', cross_section_type='ctw'):
"""
return neutrino cross-section
Parameters
----------
energy: float/ array of floats
neutrino energies/momenta in standard units
flavors: float / array of floats
neutrino flavor (integer) encoded as using PDG numbering scheme,
particles have positive sign, anti-particles have negative sign, relevant are:
* 12: electron neutrino
* 14: muon neutrino
* 16: tau neutrino
inttype: str, array of str
interaction type. Options:
* nc : neutral current
* cc : charged current
* total: total (for non-array type)
cross_section_type: {'ctw', 'ghandi', 'csms'}, default 'ctw'
defines model of cross-section. Options:
* ctw : A. Connolly, R. S. Thorne, and D. Waters, Phys. Rev.D 83, 113009 (2011).
cross-sections for all interaction types and flavors
* ghandi : according to Ghandi et al. Phys.Rev.D58:093009,1998
only one cross-section for all interactions and flavors
* csms : A. Cooper-Sarkar, P. Mertsch, S. Sarkar, JHEP 08 (2011) 042
"""
if cross_section_type == 'ghandi':
crscn = 7.84e-36 * units.cm ** 2 * np.power(energy / units.GeV, 0.363)
elif cross_section_type == 'ctw' or cross_section_type == 'hedis_bgr18':
crscn = np.zeros_like(energy)
if type(inttype) == str:
if inttype == 'total':
if isinstance(flavors, (int, np.integer)):
if flavors >= 0:
crscn = param(energy, 'nc', parameterization=cross_section_type) + param(energy, 'cc', parameterization=cross_section_type)
else:
crscn = param(energy, 'nc_bar', parameterization=cross_section_type) + param(energy, 'cc_bar', parameterization=cross_section_type)
else:
antiparticles = np.where(flavors < 0)
particles = np.where(flavors >= 0)
crscn[particles] = param(energy[particles], 'nc', parameterization=cross_section_type) + param(energy[particles], 'cc', parameterization=cross_section_type)
crscn[antiparticles] = param(energy[antiparticles], 'nc_bar', parameterization=cross_section_type) + param(energy[antiparticles], 'cc_bar', parameterization=cross_section_type)
elif inttype == 'total_up':
if isinstance(flavors, (int, np.integer)):
if flavors >= 0:
crscn = param(energy, 'nc_up') + param(energy, 'cc_up', parameterization=cross_section_type)
else:
crscn = param(energy, 'nc_bar_up') + param(energy, 'cc_bar_up', parameterization=cross_section_type)
else:
antiparticles = np.where(flavors < 0)
particles = np.where(flavors >= 0)
crscn[particles] = param(energy[particles], 'nc_up') + param(energy[particles], 'cc_up', parameterization=cross_section_type)
crscn[antiparticles] = param(energy[antiparticles], 'nc_bar_up') + param(energy[antiparticles], 'cc_bar_up', parameterization=cross_section_type)
elif inttype == 'total_down':
if isinstance(flavors, (int, np.integer)):
if flavors >= 0:
crscn = param(energy, 'nc_down') + param(energy, 'cc_down', parameterization=cross_section_type)
else:
crscn = param(energy, 'nc_bar_down') + param(energy, 'cc_bar_down', parameterization=cross_section_type)
else:
antiparticles = np.where(flavors < 0)
particles = np.where(flavors >= 0)
crscn[particles] = param(energy[particles], 'nc_down') + param(energy[particles], 'cc_down', parameterization=cross_section_type)
crscn[antiparticles] = param(energy[antiparticles], 'nc_bar_down') + param(energy[antiparticles], 'cc_bar_down', parameterization=cross_section_type)
else:
if isinstance(flavors, (int, np.integer)):
crscn = param(energy, inttype, parameterization=cross_section_type)
else:
antiparticles = np.where(flavors < 0)
particles = np.where(flavors >= 0)
crscn[particles] = param(energy[particles], inttype, parameterization=cross_section_type)
crscn[antiparticles] = param(energy[antiparticles], inttype, parameterization=cross_section_type)
else:
if isinstance(flavors, (int, np.integer)):
particles_cc = np.where(inttype == 'cc')
particles_nc = np.where(inttype == 'nc')
if flavors >= 0:
crscn[particles_cc] = param(energy[particles_cc], 'cc', parameterization=cross_section_type)
crscn[particles_nc] = param(energy[particles_nc], 'nc', parameterization=cross_section_type)
else:
crscn[particles_cc] = param(energy[particles_cc], 'cc_bar', parameterization=cross_section_type)
crscn[particles_nc] = param(energy[particles_nc], 'nc_bar', parameterization=cross_section_type)
else:
particles_cc = np.where((flavors >= 0) & (inttype == 'cc'))
particles_nc = np.where((flavors >= 0) & (inttype == 'nc'))
antiparticles_cc = np.where((flavors < 0) & (inttype == 'cc'))
antiparticles_nc = np.where((flavors < 0) & (inttype == 'nc'))
crscn[particles_cc] = param(energy[particles_cc], 'cc', parameterization=cross_section_type)
crscn[particles_nc] = param(energy[particles_nc], 'nc', parameterization=cross_section_type)
crscn[antiparticles_cc] = param(energy[antiparticles_cc], 'cc_bar', parameterization=cross_section_type)
crscn[antiparticles_nc] = param(energy[antiparticles_nc], 'nc_bar', parameterization=cross_section_type)
elif cross_section_type == 'csms':
crscn = csms(energy, inttype, flavors)
else:
logger.error("Cross-section {} not defined".format(cross_section_type))
raise NotImplementedError
return crscn
[docs]def get_interaction_length(Enu, density=.917 * units.g / units.cm ** 3, flavor=12, inttype='total',
cross_section_type='ctw'):
"""
calculates interaction length from cross section
Parameters
----------
Enu: float
neutrino energy
density: float (optional)
density of the medium, default density of ice = 0.917 g/cm**3
flavors: float / array of floats
neutrino flavor (integer) encoded as using PDG numbering scheme,
particles have positive sign, anti-particles have negative sign, relevant are:
* 12: electron neutrino
* 14: muon neutrino
* 16: tau neutrino
inttype: str, array of str
interaction type. Options:
* nc : neutral current
* cc : charged current
* total: total (for non-array type)
cross_section_type: {'ctw', 'ghandi', 'csms'}, default 'ctw'
defines model of cross-section. Options:
* ctw: A. Connolly, R. S. Thorne, and D. Waters, Phys. Rev.D 83, 113009 (2011).
cross-sections for all interaction types and flavors
* ghandi: according to Ghandi et al. Phys.Rev.D58:093009,1998
only one cross-section for all interactions and flavors
* csms: A. Cooper-Sarkar, P. Mertsch, S. Sarkar, JHEP 08 (2011) 042
Returns
-------
L_int: float
interaction length
"""
m_n = constants.m_p * units.kg # nucleon mass, assuming proton mass
L_int = m_n / get_nu_cross_section(Enu, flavors=flavor, inttype=inttype, cross_section_type=cross_section_type) / density
return L_int
if __name__ == "__main__": # this part of the code gets only executed it the script is directly called
n_points = 100
inttype = np.array(['cc'] * n_points)
# inttype = 'cc'
flavors = np.array([14] * n_points)
# flavors = 14
energy = np.logspace(13, 20, n_points) * units.eV
cc = get_nu_cross_section(energy, flavors, inttype=inttype) / units.picobarn
cc_new = get_nu_cross_section(energy, flavors, inttype=inttype, cross_section_type='csms') / units.picobarn
cc_hedis_bgr18 = get_nu_cross_section(energy, flavors, inttype=inttype, cross_section_type='hedis_bgr18') / units.picobarn
import matplotlib.pyplot as plt
plt.figure()
plt.loglog(energy / units.GeV, cc, label='CTW')
plt.loglog(energy / units.GeV, cc_new, label='CSMS')
plt.loglog(energy / units.GeV, cc_hedis_bgr18, label='HEDIS-BGR')
plt.xlabel("Energy [GeV]")
plt.ylabel("CC [pb]")
plt.legend()
plt.show()