Index: trunk/RELEASE_NOTES
===================================================================
--- trunk/RELEASE_NOTES (revision 16)
+++ trunk/RELEASE_NOTES (revision 17)
@@ -1,65 +1,74 @@
Release notes for nuCraft
+2014-06-16, rev17, Marius Wallraff
+-------------------------------------------------------------------------------
+
+Minor changes:
+- Added or improved various docstrings and brought them in line with Python's
+ recommended docstring conventions, PEP 257.
+
+
+
2014-05-20, rev16, Marius Wallraff
-------------------------------------------------------------------------------
Minor changes:
- Cleaned up the example script.
- Removed useless output of one line of code for every warning.
- Made code compatible with Python 3 (tested Python 3.3.2 and scipy 0.12.0).
2014-02-06, rev14, Marius Wallraff
-------------------------------------------------------------------------------
New feature:
- Added optional keyword argument numPrec to CalcWeights() to govern the
numerical precision with which the ZVODE ODE solver solves the Schroedinger
equation. The actual meaning of the parameter is the allowable upper limit
for the deviation of the sum of the oscillation probabilities from one, i.e.,
from the unitarity condition. An additional cross-check was implemented to
ensure the precision is met; if not, a warning is issued.
The new default value for numPrec is 5e-4 and gives slightly lower precision
then the old hard-coded default value 5e-5.
2013-11-05, rev13, Marius Wallraff
-------------------------------------------------------------------------------
Major change:
- Switched from Geant-style particle IDs to the PDG Monte Carlo particle
numbering scheme (DOI:10.1103/PhysRevD.86.010001, chapter 39) for known
flavors; sterile neutrinos are using the numbers +/- 81+ (reserved for MC
internal use), because PDG only reserved numbers for one additional
generation of fermions (7,8 for quarks, 17,18 for charged lepton and
neutrino).
This breaks backwards compatibility! Errors due the use of the old
convention will be raised as exceptions.
Conversion table for your convenience:
name: old --> new
nu_e 66 12
nu_e_bar 67 -12
nu_mu 68 14
nu_mu_bar 69 -14
nu_tau 133 16
nu_tau_bar 134 -16
nu_ster_i -2*i 80+i
nu_ster_i_bar -2*i-1 -80-i
(i=1,2,3...)
2013-xx-xx, Marius Wallraff
-------------------------------------------------------------------------------
older changes to be documented
Index: trunk/NuCraft.py
===================================================================
--- trunk/NuCraft.py (revision 16)
+++ trunk/NuCraft.py (revision 17)
@@ -1,746 +1,786 @@
#!/usr/bin/env python
-from __future__ import print_function, division
-
"""
+Main and only module of the nuCraft package that computes neutrino oscillation
+probabilities for atmospheric neutrinos by directly solving the Schroedinger
+equation.
+
+It contains the main class NuCraft with its main method CalcWeights, as well as
+the auxiliary class EarthModel.
+
+For default usage, only CalcWeights and the NuCraft constructor are needed:
+
+nC = NuCraft([1., 7.50e-5, 7.50e-5+2.32e-3], [(1,2,33.89),(1,3,9.12),(2,3,45.00)])
+nC.CalcWeights([(type1, energy1, zenith1),(type2, energy2, zenith2),...])
+
+For more information please see the docstrings of the respective classes and
+methods, or check out the example script.
+
+################################################################################
+
Copyright (c) 2013, Marius Wallraff (mwallraff#physik.rwth-aachen.de)
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name of the RWTH Aachen University nor the
names of its contributors may be used to endorse or promote products
derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL MARIUS WALLRAFF BE LIABLE FOR ANY
DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+################################################################################
+
+NuCraft heavily relies on NumPy and SciPy. If you want to use nuCraft, you might
+want to cite scipy (see www.scipy.org/Citing_SciPy) and the publication related
+to the ODE solver ZVODE, http://dx.doi.org/10.1137/0910062.
+
+NuCraft is shipped with a table of data points sampled from the Preliminary
+Earth Reference Model PREM, see http://dx.doi.org/10.1016/0031-9201(81)90046-7.
+
+Also, please cite our publication; a link and bibTeX entry will be added to this
+file as soon as it is available.
"""
+from __future__ import print_function, division
+
from math import sqrt as ssqrt
from math import cos as scos
from math import fabs as fabs # more restrictive than abs() (in a good way)
# from math import pow as spower # slower than **
from numpy import *
from scipy import interpolate, integrate
# from scipy.linalg.matfuncs import expm2 as expm2
from scipy.stats import lognorm
import warnings
def CustomWarning(message, category, filename, lineno):
print("%s:%s: %s:%s" % (filename, lineno, category.__name__, message))
warnings.showwarning = CustomWarning
from os import path
set_printoptions(precision=5, linewidth=150)
-"""
-NuCraft heavily relies on NumPy and SciPy. If you want to use nuCraft, you might
-want to cite scipy (see www.scipy.org/Citing_SciPy) and the publication related to
-the ODE solver ZVODE, http://dx.doi.org/10.1137/0910062.
-NuCraft is shipped with a table of data points sampled from the Preliminary Earth
-Model PREM, see http://dx.doi.org/10.1016/0031-9201(81)90046-7.
-
-Also, please cite our publication; a link and bibTeX entry will be added to this
-file as soon as it is available.
-"""
class EarthModel:
"""
- Auxiliary class with informations regarding the matter properties at a given
- distance from the center of the Earth, for use with nuCraft.
+ Auxiliary class with informations regarding Earth matter properties
+ at a given distance from the center of the Earth, for use with nuCraft.
The class can be constructed from an entry of the included dictionary of models,
in which case the parameter 'name' should just be the name of the model, i.e.,
the key value of the dictionary entry; or it can be constructed from a file, in
which case 'name' should be the path to the profile file. An example profile
file is provided as EarthModelPrem.txt. Names in the dictionary have priority.
In both cases, values for the (relative) electron numbers y and for the radii
of mantle (including crust), outer and inner core that are given by the model
can be overwritten with keyword arguments:
y: tuple of electron numbers, (y_mantle, y_oCore, y_iCore), 0 <= y <= 1
rE: radius of the Earth in km
rOCore: radius of the outer core in km
rICore: radius of the inner core in km, 0 <= rICore <= rOCore <= rE
"""
# models is a dictionary: name:(r values, rho values, y, rE, rOCore, rICore),
# where r are radii in km, rho the corresponding matter density values in kg/cm^3,
# y the tuple of electron numbers for mantle, outer and inner core, and the last
# three entries the corresponding radii (in km)
models = {"prem":
([0., 200., 400., 600., 800., 1000., 1200., 1221.5, 1221.5005, 1400.,
1600., 1800., 2000., 2200., 2400., 2600., 2800., 3000., 3200., 3400.,
3480., 3480.0005, 3600., 3630., 3630.0005, 3800., 4000., 4200., 4400.,
4600., 4800., 5000., 5200., 5400., 5600., 5600.0005, 5701., 5701.0005,
5771., 5771.0005, 5871., 5971., 5971.0005, 6061., 6151., 6151.0005,
6221., 6291., 6291.0005, 6346.6, 6346.6005, 6356., 6356.0005, 6368.,
6368.0005, 6371., 6371.0005, 8000.],
[13.08848, 13.07977, 13.05364, 13.01009, 12.94912, 12.87073, 12.77493,
12.76360, 12.16634, 12.06924, 11.94682, 11.80900, 11.65478, 11.48311,
11.29298, 11.08335, 10.85321, 10.60152, 10.32726, 10.02940, 9.90349,
5.56645, 5.50642, 5.49145, 5.49145, 5.40681, 5.30724, 5.20713,
5.10590, 5.00299, 4.89783, 4.78983, 4.67844, 4.56307, 4.44317,
4.44317, 4.38071, 3.99214, 3.97584, 3.97584, 3.84980, 3.72378,
3.54325, 3.48951, 3.43578, 3.35950, 3.36710, 3.37471, 3.37471,
3.38076, 2.90000, 2.90000, 2.60000, 2.60000, 1.02000, 1.02000,
0., 0.],
(0.4957, 0.4656, 0.4656),
6371.,
3480.,
1121.5)}
def __init__(self, name, y=None, rE=None, rOCore=None, rICore=None):
self.name = name
if name in self.models:
# see doc for models above
model = self.models[name]
self.y = model[2]
self.rE = model[3]
self.rOCore = model[4]
self.rICore = model[5]
# x (radius r) and y (matter density rho) values of the Earth density profile
profX = model[0]
profY = model[1]
elif path.isfile(name):
profFile = file(name, 'r')
profLines = profFile.readlines()
profFile.close()
# slightly unsafe, so don't accept nuCraft Earth profiles from untrusted sources ;)
self.y = eval(profLines[1])
self.rE = eval(profLines[2])
self.rOCore = eval(profLines[3])
self.rICore = eval(profLines[4])
profX = [float(prof.split()[0]) for prof in profLines[6:-1]]
profY = [float(prof.split()[1]) for prof in profLines[6:-1]]
else:
raise NotImplementedError("The Earth model name '%s' can not be found!" % name)
assert(len(profX) == len(profY))
# self.profInt = interpolate.interp1d(profX, profY) # same interface, but much slower
self.profInt = interpolate.InterpolatedUnivariateSpline(profX, profY, k=1)
# self.mod indicates whether this is the profile indicated by its name, or if it has
# been modified
self.mod = False
if not y == None:
if not self.y == y:
self.mod = True
self.y = y
assert(len(self.y) == 3 and logical_and(0 <= array(self.y), array(self.y) <= 1).all())
if not rE == None:
if not self.rE == rE:
self.mod = True
self.rE == rE
if not rOCore == None:
if not self.rOCore == rOCore:
self.mod = True
self.rOCore == rOCore
if not rICore == None:
if not self.rICore == rICore:
self.mod = True
self.rICore == rICore
assert(0. <= self.rICore <= self.rOCore <= self.rE)
# compute induced mass potentials for three flavors of neutrinos; see SetDim() below
self.SetDim(3)
def SetDim(self, dim):
"""
- Small function to set the number of neutrino flavors self._dim this instance of EarthModel
- should take into account; the dimension of A depends on this quantity. The function also
- computes the dim-dimensional A arrays.
+ Set the number of neutrino flavors (self._dim) this instance should take into account.
+
+ The dimension of A depends on this quantity. The function also computes the
+ dim-dimensional A arrays.
Calling this manually should not be needed because it is called by nuCraft automatically.
"""
assert(isinstance(dim, int) and dim >= 3)
self._dim = dim
# A and AxCore are constants of the squared mass potentials induced by matter
# effects, which still need to be multiplied by the neutrino energy, mass
# density, and by -1 if the particle is an anti-neutrino;
# the two versions take into account that the elementary composition of the
# inner and outer Earth cores is very different from that in the mantle, which
# is important because the relevant quantities are electron and neutron number
# densities instead of the mass densities given by the PREM profile
#
# A = 2*sqrt(2) * G_F * rho / m_N * (Y_e,0,0,1-Y_e,...) * E_nu
# = 2*sqrt(2)*1.16637e-5 / 0.939 * (Y_e,0,0,1-Y_e,...)
# * (1/1.783e-27)*(1.973e-15)**3 * rho/(kg/dm^3) * E/GeV * 1e18 * eV^2
self.A = array([15.256e-5*self.y[0], 0., 0.]+[7.6525e-5*(1-self.y[0])]*(dim-3))
self.AOCore = array([15.256e-5*self.y[1], 0., 0.]+[7.6525e-5*(1-self.y[1])]*(dim-3))
self.AICore = array([15.256e-5*self.y[2], 0., 0.]+[7.6525e-5*(1-self.y[2])]*(dim-3))
def __str__(self):
if self.mod:
return "EarthModel('%s', modified)" % self.name
else:
return "EarthModel('%s')" % self.name
def __repr__(self):
# assumes that EarthModel was imported into the main namespace (i.e., from NuCraft...),
# like numpy repr do; do nc = eval("NuCraft."+repr(nc)) otherwise
return "EarthModel('%s', y=(%.17f, %.17f, %.17f), rE=%.17f, rOCore=%.17f, rICore=%.17f)" \
% (self.name, self.y[0], self.y[1], self.y[2], self.rE, self.rOCore, self.rICore)
class NuCraft:
"""
- nuCraft is a Python class for calculation of neutrino oscillation probabilities
+ Main class that calculates atmospheric neutrino oscillation probabilities
by directly solving the Schroedinger equation. It includes matter effects (using the
PREM Earth density profile, modified for neutron abundance in the core), and supports
an arbitrary number of neutrino "flavors". The first flavor is interpreted as electron
neutrinos, the next two as non-electron non-sterile neutrinos, and all other flavors
as sterile neutrinos; to modify this, edit self.A and self.ACore.
An instance of the class is created by calling NuCraft(deltaMi1List, angleList),
where deltaMi1List is a list like described in ConstructMassMatrix, and angleList like
described in ConstructMixingMatrix.
The recommended method to calculate the oscillation probabilites is CalcWeights;
it solves the Schroedinger equation in the interaction basis (i.e., where vacuum
oscillations are flat), and it supports proper handling of the Earth's atmosphere,
which is enabled by default (atmMode 3; see docstring).
- The other method solves the Schroedinger equation in the flavor basis:
- CalcWeightsLegacy
- It is mostly kept for comparability, as it is much slower and less precise for
+ The other method CalcWeightsLegacy solves the Schroedinger equation in the flavor basis;
+ it is mostly kept for comparability, as it is much slower and less precise for
most problems. It can be suitable for high-energy sterile neutrino problems, where
vacuum oscillations are reeaaally slow. Its handling of the Earth's atmosphere is
very basic (atmMode == 0 of the other method).
The code should NEVER throw warnings like "excess work done"; if it does, adjust the
ODE solver parameters as specified in the warning, or contact the author.
+ Seldom warnings saying "the computed unitarity does not meet the specified precision"
+ are unproblematic if the inequalities thereafter are not off by much (e.g., 50%).
"""
def __init__(self, deltaMi1List, angleList, earthModel=EarthModel("prem")):
if len(deltaMi1List) == 2 and len(angleList) == 1 and angleList[0][1] == 2:
deltaMi1List = list(deltaMi1List)
deltaMi1List.append(0.)
angleList = list(angleList)
angleList.append((1,3,0.))
# mainly used for __repr__, do not modify in instances ("private")
self._deltaMi1List = deltaMi1List
self._angleList = angleList
self.M = self.ConstructMassMatrix(deltaMi1List)
self.U = self.ConstructMixingMatrix(angleList)
dim = len(self.M)
assert(self.M.shape == self.U.shape == eye(dim).shape)
if isinstance(earthModel, EarthModel):
earthModel.SetDim(dim)
self.earthModel = earthModel
else:
raise ValueError("The provided earth model '%s' is not of the EarthModel class." % earthModel)
# tuple entry 0 became redundant because of the transition from Geant-style particle IDs to
# PDG-style particle IDs; kept it for now because of performance reasons (prob. negligible)
self.mcTypeDict = {}
self.mcTypeDict[ 12] = ( 1, array([1.] +[0.j]*(dim-1)))
self.mcTypeDict[-12] = (-1, array([1.] +[0.j]*(dim-1)))
self.mcTypeDict[ 14] = ( 1, array([0.j,1.] +[0.j]*(dim-2)))
self.mcTypeDict[-14] = (-1, array([0.j,1.] +[0.j]*(dim-2)))
self.mcTypeDict[ 16] = ( 1, array([0.j,0.j,1.]+[0.j]*(dim-3)))
self.mcTypeDict[-16] = (-1, array([0.j,0.j,1.]+[0.j]*(dim-3)))
for i in range(1,dim-2):
self.mcTypeDict[ 80+i] = ( 1, array([0.j]*(2+i)+[1.]+[0.j]*(dim-3-i)))
self.mcTypeDict[-80-i] = (-1, array([0.j]*(2+i)+[1.]+[0.j]*(dim-3-i)))
# depth of the (center of the) detector below the surface of the Earth sphere, in km
self.detectorDepth = 2.
# extension of the atmosphere above the surface of the Earth sphere, in km; ignored in default
# atmosphere handling mode (uses Gaisser/Stanev model instead, see InteractionAlt method below)
self.atmHeight = 20.
def __repr__(self):
# assumes that nuCraft was imported into the main namespace (i.e., from NuCraft...),
# like numpy repr do; do nc = eval("NuCraft."+repr(nc)) otherwise;
# availability of the class EarthModel is essential
if ( (self.M == self.ConstructMassMatrix(self._deltaMi1List)).all()
and (self.U == self.ConstructMixingMatrix(self._angleList)).all()):
return "NuCraft(%s, %s, earthModel=%s)" % (self._deltaMi1List, self._angleList, repr(self.earthModel))
else:
return "NuCraft(%s, %s, mod)" % (self._deltaMi1List, self._angleList)
def __str__(self):
return "nuCraft(n=%d)" % len(self.M)
def ConstructMassMatrix(self, parList):
"""
- constructs and returns a mass matrix out of the input list; first parameter
- is the mass of mass state 1, the following parameters are the correctly
- ordered squared mass differences of the other states to state 1, i.e.,
+ Construct and return a squared-mass matrix out of the input list;
+ the first parameter is the mass of mass state 1, the following parameters are
+ the correctly ordered squared mass differences of the other states to state 1,
+ all given in units of eV or eV^2, respectively, i.e.,
parList[i] = m_i^2 - m_1^2 for i > 0; e.g.,
parList = [1., 7.50e-5, 7.50e-5+2.32e-3]
"""
# ensure that all masses are positive
assert(-parList[0]**2 <= min(parList[1:]))
return diag([parList[0]**2] + [parList[0]**2 + m for m in parList[1:]])
def ConstructMixingMatrix(self, parList):
"""
- constructs and returns a mixing matrix out of the input list; the list should
- consist of tuples of the format (i,j,theta_ij), i<j, with theta in degrees,
- and the mixing matrix will be constructed in reverse order, e.g.:
+ Construct and return a mixing matrix out of the input list;
+ the list should consist of tuples of the format (i,j,theta_ij), i<j, with theta
+ in degrees, and the mixing matrix will be constructed in reverse order, e.g.:
parList = [(1,2,33.89),(1,3,9.12),(2,3,45.00)]
=> U = R_23 . R_13 . R_12
- for CP-violating factors, use tuples like (i,j,theta_ij,delta_ij),
- with delta_ij in degrees
+ For CP-violating factors, use tuples like (i,j,theta_ij,delta_ij),
+ with delta_ij in degrees.
"""
dim = max([par[1] for par in parList])
def RotMat(dim, i, j, ang, cp): # actually not rotation matrices, but Gell-Mann-generated matrices
if not i < j <= dim:
raise Exception("Missconstructed rotation matrix: "+repr(dim)+", "+repr(i)+", "+repr(j)+", "+repr(ang))
if cp == 0:
R = eye(dim)
R[i,j] = sin(ang)
R[j,i] = -sin(ang)
else:
R = eye(dim, dtype='complex128')
R[i,j] = sin(ang) * exp(-1j*cp)
R[j,i] = -sin(ang) * exp(1j*cp)
R[i,i] = R[j,j] = cos(ang)
return R
degToRad = pi/180.
U = eye(dim)
for par in parList:
if len(par)>3:
U = dot(RotMat(dim,par[0]-1,par[1]-1,par[2]*degToRad,par[3]*degToRad), U)
else:
U = dot(RotMat(dim,par[0]-1,par[1]-1,par[2]*degToRad,0), U)
return U
def InteractionAlt(self, mcType, mcEn, mcZen, mode):
"""
- helper method; depending on mode, returns a list of one or more tuples of
- weights and altitudes in which the neutrinos should start to be propagated:
+ Return a list of weight-altitude tuples for atmospheric propagation.
+
+ Helper method; depending on mode, returns a list of one or more tuples of
+ weights and altitudes in which the neutrinos should start to be propagated.
+ The weights have to add up to one.
mode 0:
- returns self.rAtm = self.rE + 20. with weight 1., which means that interaction
- is expected to happen fixed 20 km above ground level
+ returns self.earthModel.rE + self.atmHeight with weight 1., which means that the
+ interaction is expected to happen at a fixed hight (default 20 km) above ground level
mode 1:
samples a single altitude from a parametrization to the atmospheric interaction
model presented in "Path length distributions of atmospheric neutrinos",
Gaisser and Stanev, PhysRevD.57.1977
mode 2 and mode 3:
returns eight equally probable altitudes from the whole range of values allowed
by the parametrization also used in mode 1
"""
if mode == 0:
return [(1., self.earthModel.rE + self.atmHeight)]
# get the cosine of the zenith angle, and compensate for not having a parametrization
# for the effective zenith angle (i.e., the zenith angle relative to the Earth's curvature
# averaged through the atmosphere) at cos(zen) < 0.05
cosZen = (fabs(scos(mcZen))**3 + 0.000125)**0.333333333
if mcType in (12, -12): # electron neutrinos
mu = 1.285e-9*(cosZen-4.677)**14. + 2.581
sigma = 0.6048*cosZen**0.7667 - 0.5308*cosZen + 0.1823
else:
mu = 1.546e-9*(cosZen-4.618)**14. + 2.553
sigma = 1.729*cosZen**0.8938 - 1.634*cosZen + 0.1844
logn = lognorm(sigma, scale=2*exp(mu), loc=-12)
if mode == 1:
z = logn.rvs()*cosZen
while z < 0:
z = logn.rvs()*cosZen
return [(1., z+self.earthModel.rE)]
elif mode in [2, 3]:
cdf0 = logn.cdf(0)
qList = cdf0 + array([0.0625,0.1875,0.3125,0.4375,0.5625,0.6875,0.8125,0.9375])*(1.-cdf0)
return list(zip(ones(8)/8., logn.ppf(qList)*cosZen + self.earthModel.rE))
else:
raise NotImplementedError("Unrecognized mode for neutrino interaction height estimation!")
def CalcWeights(self, inList, vacuum=False, atmMode=3, numPrec=5e-4):
"""
- calculates and returns the oscillation probabilities for
- NuE: mcType = 12
- NuEBar: mcType = -12
- NuMu: mcType = 14
- NuMuBar: mcType = -14
- NuTau: mcType = 16
- NuTauBar: mcType = -16
- NuSterile_i: mcType = 80+i i = 1,2,3...
- NuSterileBar_i: mcType = -80-i
+ Calculate neutrino oscillation probabilities in the interaction basis.
inList may be provided in three formats:
1 - as a tuple of three lists: mcType, neutrino energy and zenith angle,
i.e., ([type1, type2, ...],[energy1, energy2, ...],[zenith1, zenith2, ...])
2 - as a list of tuples (mcType, neutrinoEnergy, zenith angle),
i.e., [(type1, energy1, zenith1),(type2, energy2, zenith2),...]
- 3 - a list of particles of a customizable format
+ 3 - a list of particles of a customizable format:
from collections import namedtuple
SimPart = namedtuple("SimPart", (..., "zenMC", "eMC", "mcType", "oscProb", ...))
part = SimPart(..., zenith, energy, mctype, -1., ...)
+ mcType uses PDG conventions:
+ NuE: mcType = 12
+ NuEBar: mcType = -12
+ NuMu: mcType = 14
+ NuMuBar: mcType = -14
+ NuTau: mcType = 16
+ NuTauBar: mcType = -16
+ NuSterile_i: mcType = 80+i i = 1,2,3...
+ NuSterileBar_i: mcType = -80-i
+
+ energy is expected in units of GeV, zenith angles in radian
+
return format:
in input format cases 1 and 2, the return format is a list of [P_E, P_Mu, P_Tau, ...];
in case 3, the list of particles is returned with updated oscProb fields
atmMode controls the handling of the atmospheric interaction altitude:
- 0 assumes a fixed height of 20 km (old default)
+ 0 assumes a fixed height of self.atmHeight (20 km by default)
1 draws a single altitude from a parametrization described in the InteractionAlt method
2 calculates the average for eight altitudes distributed over the whole range according
- to the same parametrization as in mode 1; this is slow and only meant for debugging
+ to the same parametrization as in mode 1; this is SLOW and only meant for debugging
3 uses the same altitudes as mode 2, but only propagates the lowest-altitude neutrino
and adds the remaining lengths as vacuum oscillations afterwards; this is pretty fast
- and the new default mode
+ and used by default
numPrec governs the numerical precision with which the Schroedinger equation is solved;
the unitarity condition (i.e., the fact that the sum of the resulting probabilities
should be 1.) is used as a simple cross-check, a warning is issued if the deviation from
1. is larger than numPrec.
"""
# for the input format checks, assume that inList is homogeneous
assert(type(vacuum) is bool)
assert(atmMode in range(4))
partMode = False
if not len(inList): # empty input
return []
elif type(inList) is tuple and not type(inList[0]) is tuple: # input case 1
assert(len(inList[0]) == len(inList[1]) == len(inList[2]))
inList = zip(*inList)
elif not type(inList) is tuple and type(inList[0]) is tuple: # input case 2
assert(len(inList[0]) == 3)
elif type(inList[0]).__base__ is tuple: # input case 3
try:
if not set(("zenMC", "eMC", "mcType", "oscProb")).issubset(inList[0]._fields):
raise Exception("NamedTuple in inList does not have the required fields!")
except AttributeError:
raise NotImplementedError("Format of the input inList is not supported!")
partMode = True
else:
raise NotImplementedError("Wrong input format!")
VAC = dot(self.U, dot(self.M, conj(self.U).T))
(svdVAC0n, svdVAC1, svdVAC2n) = linalg.svd(VAC)
svdVAC0n = array(svdVAC0n, order='C')
svdVAC2n = array(svdVAC2n, order='C') # svd returns fortran-ordered arrays
svdVAC0a = conj(svdVAC0n)
svdVAC2a = conj(svdVAC2n)
# caching of some quantities that will be used often later on
rE = self.earthModel.rE
rOCore = self.earthModel.rOCore
rICore = self.earthModel.rICore
rDet = self.earthModel.rE - self.detectorDepth
rAtm = self.earthModel.rE + self.atmHeight
profInt = self.earthModel.profInt # matter density profile interpolation
global lCache
global aMSW
global modVAC
lCache = pi
- # function to calculate the oscillation probabilities of an individual nu
def calcProb(inTuple):
+ """
+ Calculate the oscillation probabilities of an individual nu.
+ """
# python 3 does not support tuple parameter unpacking anymore
mcType, mcEn, mcZen = inTuple
- # update the MSW-induced mass-squared matrix aMSW with the current density,
- # and update the state-transition matrix modVAC to the current time/position
def dscM(l, en, zen):
+ """
+ Update the MSW-induced mass-squared matrix aMSW with the current density,
+ and update the state-transition matrix modVAC to the current time/position.
+ """
global lCache, aMSW, modVAC
# if l did not change, the update is unnecessary
if l == lCache:
return
# modVAC is the time-dependent state-transition matrix that brings a
# state vector to the interaction basis, i.e., to the basis where
# the vacuum oscillations are flat
if isAnti:
modVAC = dot(svdVAC0a * exp(-2.533866j/en*svdVAC1*(L-l)), svdVAC2a)
else:
modVAC = dot(svdVAC0n * exp(-2.533866j/en*svdVAC1*(L-l)), svdVAC2n)
# <==> modVAC = dot(dot(svdVAC0, diag(exp(-2.533866j/mcEn*svdVAC1*(L-l)))), svdVAC2)
# distance from the center of Earth
r = ssqrt( l*l + rDet*rDet - 2*l*rDet*scos(pi - zen) )
if r <= rICore:
aMSW = aMSWWoRhoICore * profInt(r)
elif r <= rOCore:
aMSW = aMSWWoRhoOCore * profInt(r)
else:
aMSW = aMSWWoRho * profInt(r)
lCache = l
def f(l, psi, en, zen):
dscM(l, en, zen)
return -2.533866j/en * dot(modVAC, aMSW*dot(psi, conj(modVAC)))
# <==> return -2.533866j/en * dot(modVAC, dot(diag(aMSW), dot(conj(modVAC).T, psi)))
def jac(l, psi, en, zen):
dscM(l, en, zen)
return -2.533866j/en * dot(modVAC*aMSW, conj(modVAC).T)
# <==> return -2.533866j/en * dot(dot(modVAC, diag(aMSW)), conj(modVAC).T)
try:
mcType = self.mcTypeDict[mcType]
except KeyError:
raise KeyError("The mcType %d is not known to nuCraft!" % mcType)
isAnti = mcType[0] == -1
# inefficient performance-wise, but nicer code, and vacuum is fast enough anyway
if vacuum:
aMSWWoRho = zeros_like(self.earthModel.A)
aMSWWoRhoOCore = aMSWWoRho
aMSWWoRhoICore = aMSWWoRho
else:
aMSWWoRho = mcType[0] * self.earthModel.A * mcEn
aMSWWoRhoOCore = mcType[0] * self.earthModel.AOCore * mcEn
aMSWWoRhoICore = mcType[0] * self.earthModel.AICore * mcEn
# depending on the mode, get a list of interaction altitude tuples, see method doc string;
# the first number of the tuples is the weight, the second the distance of the interaction
# point to the center of the Earth; the weights have to add up to 1.
if atmMode == 3:
# first get the tuples and propagate only the lowest-altitude neutrino:
rAtmTuples = self.InteractionAlt(mcType, mcEn, mcZen, 2)
rAtm = rAtmTuples[0][1]
L = ssqrt( rAtm*rAtm + rDet*rDet - 2*rAtm*rDet*scos( mcZen - arcsin(sin(pi-mcZen)/rAtm*rDet) ) )
dscM(L, mcEn, mcZen)
solver = integrate.ode(f, jac).set_integrator('zvode', method='adams', order=4, with_jacobian=True,
nsteps=120000, min_step=0.0002, max_step=500.,
atol=numPrec*2e-3, rtol=numPrec*2e-3)
solver.set_initial_value(dot(modVAC, mcType[1]), L).set_f_params(mcEn, mcZen).set_jac_params(mcEn, mcZen)
solver.integrate(0.)
if isAnti:
endVAC = dot(svdVAC0a * exp(-2.533866j/mcEn*svdVAC1*L), svdVAC2a)
else:
endVAC = dot(svdVAC0n * exp(-2.533866j/mcEn*svdVAC1*L), svdVAC2n)
# <==> endVAC = dot(dot(svdVAC0, diag(exp(-2.533866j/mcEn*svdVAC1*L))), svdVAC2)
results = [rAtmTuples[0][0] * square(absolute( dot(conj(endVAC).T, solver.y) ))]
# now for all the other neutrinos, add the missing lengths as vacuum oscillations
# at the end of the track; keep in mind that atmosphere is always handled as vacuum
for rAtmWeight, rAtm in rAtmTuples[1:]:
L = ssqrt( rAtm*rAtm + rDet*rDet - 2*rAtm*rDet*scos( mcZen - arcsin(sin(pi-mcZen)/rAtm*rDet) ) )
if isAnti:
endVAC = dot(svdVAC0a * exp(-2.533866j/mcEn*svdVAC1*L), svdVAC2a)
else:
endVAC = dot(svdVAC0n * exp(-2.533866j/mcEn*svdVAC1*L), svdVAC2n)
results.append( rAtmWeight * square(absolute( dot(conj(endVAC).T, solver.y) )) )
else:
# in this case, just stupidly propagate every neutrino in the list...
results = []
for rAtmWeight, rAtm in self.InteractionAlt(mcType, mcEn, mcZen, atmMode):
L = ssqrt( rAtm*rAtm + rDet*rDet - 2*rAtm*rDet*scos( mcZen - arcsin(sin(pi-mcZen)/rAtm*rDet) ) )
dscM(L, mcEn, mcZen)
solver = integrate.ode(f, jac).set_integrator('zvode', method='adams', order=4, with_jacobian=True,
nsteps=120000, min_step=0.0002, max_step=500.,
atol=numPrec*2e-3, rtol=numPrec*2e-3)
solver.set_initial_value(dot(modVAC, mcType[1]), L).set_f_params(mcEn, mcZen).set_jac_params(mcEn, mcZen)
solver.integrate(0.)
if isAnti:
endVAC = dot(svdVAC0a * exp(-2.533866j/mcEn*svdVAC1*L), svdVAC2a)
else:
endVAC = dot(svdVAC0n * exp(-2.533866j/mcEn*svdVAC1*L), svdVAC2n)
# <==> endVAC = dot(dot(svdVAC0, diag(exp(-2.533866j/mcEn*svdVAC1*L))), svdVAC2)
results.append( rAtmWeight * square(absolute( dot(conj(endVAC).T, solver.y) )) )
prob = sum(results, 0)
if abs(1-sum(prob)) > numPrec:
warnings.warn("The computed unitarity does not meet the specified precision: %.2e > %.2e" % (abs(1-sum(prob)), numPrec))
return prob
if partMode:
# return the list of input particles with updated oscProb fields
return [p._replace(oscProb=calcProb((p.mcType, p.eMC, p.zenMC))) for p in inList]
else:
# return a list with one oscillation probability array per input tuple
return [calcProb(t) for t in inList]
def CalcWeightsLegacy(self, inList, vacuum=False):
"""
- legacy mode; solves the Schroedinger equation in the flavor basis; better
- use CalcWeights, unless you know what you are doing
+ Calculate neutrino oscillation probabilities in the flavor basis.
- does not properly take into account the atmosphere (only atmMode 0 is available),
- and does not support the numPrec keyword argument
-
- calculates and returns the oscillation probabilities for
- NuE: mcType = 12
- NuEBar: mcType = -12
- NuMu: mcType = 14
- NuMuBar: mcType = -14
- NuTau: mcType = 16
- NuTauBar: mcType = -16
- NuSterile_i: mcType = 80+i i = 1,2,3...
- NuSterileBar_i: mcType = -80-i
+ Legacy mode, solves the Schroedinger equation in the flavor basis; better
+ use CalcWeights, unless you know what you are doing.
inList may be provided in three formats:
1 - as a tuple of three lists: mcType, neutrino energy and zenith angle,
i.e., ([type1, type2, ...],[energy1, energy2, ...],[zenith1, zenith2, ...])
2 - as a list of tuples (mcType, neutrinoEnergy, zenith angle),
i.e., [(type1, energy1, zenith1),(type2, energy2, zenith2),...]
3 - a list of particles of a customizable format
from collections import namedtuple
SimPart = namedtuple("SimPart", (..., "zenMC", "eMC", "mcType", "oscProb", ...))
part = SimPart(..., zenith, energy, mctype, -1., ...)
+ mcType uses PDG conventions:
+ NuE: mcType = 12
+ NuEBar: mcType = -12
+ NuMu: mcType = 14
+ NuMuBar: mcType = -14
+ NuTau: mcType = 16
+ NuTauBar: mcType = -16
+ NuSterile_i: mcType = 80+i i = 1,2,3...
+ NuSterileBar_i: mcType = -80-i
+
+ energy is expected in units of GeV, zenith angles in radian
+
return format:
in input format cases 1 and 2, the return format is a list of [P_E, P_Mu, P_Tau, ...];
in case 3, the list of particles is returned with updated oscProb fields
+
+ Does not properly take into account the atmosphere (only atmMode 0 is available),
+ and does not support the numPrec keyword argument.
"""
warnings.warn("Legacy mode methods are often slower and less precise, use at your own risk")
# for the input format checks, assume that inList is homogeneous
assert(type(vacuum) is bool)
partMode = False
if not len(inList): # empty input
return []
elif type(inList) is tuple and not type(inList[0]) is tuple: # input case 1
assert(len(inList[0]) == len(inList[1]) == len(inList[2]))
inList = zip(*inList)
elif not type(inList) is tuple and type(inList[0]) is tuple: # input case 2
assert(len(inList[0]) == 3)
elif type(inList[0]).__base__ is tuple: # input case 3
try:
if not set(("zenMC", "eMC", "mcType", "oscProb")).issubset(inList[0]._fields):
raise Exception("NamedTuple in inList does not have the required fields!")
except AttributeError:
raise NotImplementedError("Format of the input inList is not supported!")
partMode = True
else:
raise NotImplementedError("Wrong input format!")
VACn = dot(self.U, dot(self.M, conj(self.U).T))
VACa = conj(VACn)
# radius of the Earth, radius to the center of the detector,
# and radius including the atmosphere, all in units of km
rE = self.earthModel.rE
rDet = self.earthModel.rE - self.detectorDepth
rAtm = self.earthModel.rE + self.atmHeight
global lCache
global aMSW
lCache = pi
- # function to calculate the oscillation probabilities of an individual nu
def calcProb(inTuple):
+ """
+ Calculate the oscillation probabilities of an individual nu.
+ """
# python 3 does not support tuple parameter unpacking anymore
mcType, mcEn, mcZen = inTuple
- # update the MSW-induced mass-squared matrix aMSW with the current density
def dscM(l, zen):
+ """
+ Update the MSW-induced mass-squared matrix aMSW with the current density.
+ """
global lCache, aMSW
# if l did not change, the update is unnecessary
if l == lCache:
return
# distance from the center of Earth
r = ssqrt( l*l + rDet*rDet - 2*l*rDet*scos(pi - zen) )
if r <= self.earthModel.rICore:
aMSW = aMSWWoRhoICore * self.earthModel.profInt(r)
elif r <= self.earthModel.rOCore:
aMSW = aMSWWoRhoOCore * self.earthModel.profInt(r)
else:
aMSW = aMSWWoRho * self.earthModel.profInt(r)
lCache = l
def f(l, psi, en, zen):
dscM(l, zen)
if isAnti:
return -2.533866j/en * dot((VACa + aMSW), psi)
else:
return -2.533866j/en * dot((VACn + aMSW), psi)
-
def jac(l, psi, en, zen):
dscM(l, zen)
if isAnti:
return -2.533866j/en * (VACa + aMSW)
else:
return -2.533866j/en * (VACn + aMSW)
try:
mcType = self.mcTypeDict[mcType]
except KeyError:
raise KeyError("The mcType %d is not known to nuCraft!" % mcType)
isAnti = mcType[0] == -1
# inefficient performance-wise, but nicer code, and vacuum is fast enough anyway
if vacuum:
aMSWWoRho = diag(zeros_like(self.earthModel.A))
aMSWWoRhoOCore = aMSWWoRho
aMSWWoRhoICore = aMSWWoRho
else:
aMSWWoRho = diag(mcType[0] * self.earthModel.A * mcEn)
aMSWWoRhoOCore = diag(mcType[0] * self.earthModel.AOCore * mcEn)
aMSWWoRhoICore = diag(mcType[0] * self.earthModel.AICore * mcEn)
L = ssqrt( rAtm*rAtm + rDet*rDet - 2*rAtm*rDet*scos( mcZen - arcsin(sin(pi-mcZen)/rAtm*rDet) ) )
dscM(L, mcZen)
solver = integrate.ode(f, jac).set_integrator('zvode', method='adams', order=4, with_jacobian=True,
nsteps=1200000, min_step=0.0002, max_step=500.,
atol=.1e-6, rtol=.1e-6)
solver.set_initial_value(mcType[1], L).set_f_params(mcEn, mcZen).set_jac_params(mcEn, mcZen)
solver.integrate(0.)
return square(absolute(solver.y))
if partMode:
# return the list of input particles with updated oscProb fields
return [p._replace(oscProb=calcProb((p.mcType, p.eMC, p.zenMC))) for p in inList]
else:
# return a list with one oscillation probability array per input tuple
return [calcProb(t) for t in inList]
# print( " \n Completed without fatal errors! \n " )
Index: trunk/__init__.py
===================================================================
--- trunk/__init__.py (revision 16)
+++ trunk/__init__.py (revision 17)
@@ -1,6 +1,11 @@
#! /usr/bin/env python
-
+"""
+The nuCraft package for calculation of atmospheric neutrino oscillation probabilities
+consists of the NuCraft module with the main class NuCraft and the auxiliary class
+EarthModel, which are both imported to the package namespace.
+For more information, see the module's and classes' docstrings.
+"""
__all__ = ["NuCraft","EarthModel"]
from NuCraft import EarthModel
from NuCraft import NuCraft