Page MenuHomeHEPForge
Files F7878961

No OneTemporary
Actions

View Options

Index: trunk/RELEASE_NOTES
===================================================================
--- trunk/RELEASE_NOTES (revision 21)
+++ trunk/RELEASE_NOTES (revision 22)
@@ -1,127 +1,135 @@
Release notes for nuCraft
+2020-05-24, rev22, Marius Wallraff
+-------------------------------------------------------------------------------
+
+Minor changes:
+- Fixed the example script to be compatible with various versions of matplotlib.
+- Updated README file.
+
+
2018-06-26, rev21, Marius Wallraff
-------------------------------------------------------------------------------
Minor changes:
- Fixed the incompatibility of EarthModel's construction from file with Python
3; to achieve this, the built-in function "file" was replaced by "open".
2017-01-09, rev20, Marius Wallraff
-------------------------------------------------------------------------------
Minor changes:
- Bug fix in EarthModel's init method; setting custom Earth radii by providing
init arguments did not work as intended because of "==" instead of "=".
Thanks to Justin Lanfranchi for finding and reporting this.
- Minor improvements in the documentation.
- Removal of matplotlib "usetex" in the standalone example for improved compat-
ibility; had been removed already in CPC version.
2015-03-09, rev19, Marius Wallraff
-------------------------------------------------------------------------------
Minor changes:
- Made detectorDepth and atmHeight more easily configurable by adding keyword
arguments with default values to the NuCraft class.
- Backported keyword argument numPrec to CalcWeightsLegacy method.
- Replaced the built-in eval function with ast.literal_eval for reading in
EarthModel files to eliminate security risks.
- Added text to the assertions and improved README, docstrings and comments;
cleaned up code.
- Improved the example script.
- Slightly adjusted ODE solver settings:
- Order was changed from 4 to 5, where 5 is the highest order that features
changes compared to lower orders; seldom requires less steps than before,
therefore a miniscule performance increase
- Removed min_steps and max_steps and increased nsteps; those parameters were/
are chosen not to be reached, anyway
2014-08-01, rev18, Marius Wallraff
-------------------------------------------------------------------------------
Minor changes:
- Added a small README file in compliance with the Computer Physics
Communications Programm Library guidelines.
- Fixed a small error in an inline documentation string.
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/example-standAlone.py
===================================================================
--- trunk/example-standAlone.py (revision 21)
+++ trunk/example-standAlone.py (revision 22)
@@ -1,171 +1,177 @@
#!/usr/bin/env python
from __future__ import print_function, division
from numpy import *
max = amax
min = amin
import os, sys
from time import time
from sys import stdout
from NuCraft import *
# Example script that reproduces the left side of figure 1 from the
# Akhmedov/Rassaque/Smirnov paper arxiv 1205.7071;
# for the core region, results differ because NuCraft takes into account
# the deviation of Y_e from 0.5 due to heavy elements in the Earth's core.
# number of energy bins
eBins = 800
# zenith angles for the four plots
zList = arccos([-1., -0.8, -0.6, -0.4])
# energy range in GeV
eList = linspace(1, 20, eBins)
# parameters from arxiv 1205.7071
theta23 = arcsin(sqrt(0.420))/pi*180.
theta13 = arcsin(sqrt(0.025))/pi*180.
theta12 = arcsin(sqrt(0.312))/pi*180.
DM21 = 7.60e-5
DM31 = 2.35e-3 + DM21
# Akhemdov is implicitely assuming an electron-to-neutron ratio of 0.5;
# he is also using the approximation DM31 = DM32;
# if you want to reproduce his numbers exactly, switch the lines below, and turn
# atmosphereMode to 0 (no handling of the atmosphere because of )
# AkhmedovOsci = NuCraft((1., DM21, DM31-DM21), [(1,2,theta12),(1,3,theta13,0),(2,3,theta23)],
# earthModel=EarthModel("prem", y=(0.5,0.5,0.5))
# detectorDepth=0., atmHeight=0.)
AkhmedovOsci = NuCraft((1., DM21, DM31), [(1,2,theta12),(1,3,theta13,0),(2,3,theta23)])
# To compute weights with a non-zero CP-violating phase, replace ^ this zero
# by the corresponding angle (in degrees); this will add the phase to the theta13 mixing matrix,
# as it is done in the standard parametrization; alternatively, you can also add CP-violating
# phases to the other matrices, but in the 3-flavor case more than one phase are redundant.
# atmosphereMode = 0 # use fixed atmsopheric depth (set to 0 km if line 42 is not commented out)
atmosphereMode = 3 # default: efficiently calculate eight path lenghts per neutrino and take the average
# This parameter governs the precision with which nuCraft computes the weights; it is the upper
# limit for the deviation of the sum of the resulting probabilities from unitarity.
# You can verify this by checking the output plot example-standAlone2.png.
numPrec = 5e-4
# 12, -12: NuE, NuEBar
# 14, -14: NuMu, NuMuBar
# 16, -16: NuTau, NuTauBar
pType = 14
print("Calculating...")
# two methods of using the nuCraft instance:
# saving the current time to measure the time needed for the execution of the following code
t = time()
### using particles
"""
from collections import namedtuple
DatPart = namedtuple("DatPart", ("zen", "eTrunc", "eMuex", "eMill", "eTrumi"))
SimPart = namedtuple("SimPart", DatPart._fields + ("atmWeight", "zenMC", "eMC", "eMuon", "mcType", "oscProb"))
prob = array([zeros([3,eBins]), zeros([3,eBins]), zeros([3,eBins]), zeros([3,eBins])])
for index, zenith in enumerate(zList):
for eIndex, energy in enumerate(eList):
p = SimPart(0.,0.,0.,0.,0.,0.,
zenith,energy,0.,pType, -1.)
# one call to nuCraft per particle; could also pass all particles in a list
pList = AkhmedovOsci.CalcWeights([p], numPrec=numPrec, atmMode=atmosphereMode)
prob[index][:, eIndex] = pList[0].oscProb
"""
### using lists (arrays, actually)
zListLong, eListLong = meshgrid(zList, eList)
zListLong = zListLong.flatten()
eListLong = eListLong.flatten()
tListLong = ones_like(eListLong)*pType
# actual call to nuCraft for weight calculations:
prob = AkhmedovOsci.CalcWeights((tListLong, eListLong, zListLong), numPrec=numPrec, atmMode=atmosphereMode)
prob = rollaxis( array(prob).reshape(len(eList), len(zList),-1), 0,3)
# rollaxis is only needed to get the same shape as prob from above,
# i.e., four elements for the different zenith angles, of which each is an
# array of 3 x eBins (three oscillation probabilities for every energy bin)
print("Calculating the probabilities took %f seconds." % (time()-t))
print("Plotting...")
-from matplotlib import rc
-rc('axes', grid=True, titlesize=14, labelsize=14, color_cycle=['b','r','k']) # only available in recent versions
-rc('xtick', labelsize=12)
-rc('ytick', labelsize=12)
-rc('lines', linewidth=2)
+try:
+ import matplotlib
+except ImportError:
+ raise ImportError( "Sorry, matplotlib does not seem to be installed. All calculations were finished, "
+ +"though, so nuCraft seems to work well.")
+if float(matplotlib.__version__[0:3]) < 1.5:
+ matplotlib.use('Agg') # to enable headless output, i.e., when no display is connected
+ matplotlib.rc('axes', grid=True, color_cycle=['b','r','k'])
+ # matplotlib.rc('text', usetex=True) # uncomment this line if plot labels are garbled
+else:
+ matplotlib.rc('axes', grid=True, prop_cycle=matplotlib.rcsetup.cycler(color=['b','r','k']))
from pylab import *
# plot the probabilities
fig = figure(figsize = (6,10))
ax1 = fig.add_subplot(411)
ax1.plot(eList, prob[0].T)
ax1.set_title('NuMu to NuE (blue), NuMu (red), and NuTau (black)')
ax1.set_ylabel(r'zenith angle $%d^\circ$' % (zList[0]/pi*180))
ax1.set_xlim([1,20])
ax2 = fig.add_subplot(412)
ax2.plot(eList, prob[1].T)
ax2.set_ylabel(r'zenith angle $%d^\circ$' % (zList[1]/pi*180))
ax2.set_xlim([1,20])
ax3 = fig.add_subplot(413)
ax3.plot(eList, prob[2].T)
ax3.set_ylabel(r'zenith angle $%d^\circ$' % (zList[2]/pi*180))
ax3.set_xlim([1,20])
ax4 = fig.add_subplot(414)
ax4.plot(eList, prob[3].T)
ax4.set_ylabel(r'zenith angle $%d^\circ$' % (zList[3]/pi*180))
ax4.set_xlim([1,20])
ax4.set_xlabel("neutrino energy / GeV")
savefig("example-standAlone1.png", dpi=160)
# plot the verification of unitarity (sum(P) = 1)
figa = figure(figsize = (6,10))
ax1a = figa.add_subplot(411)
ax1a.plot(eList, sum(prob[0],0))
ax1a.set_title('unitarity test')
ax1a.set_ylabel(r'zenith angle $%d^\circ$' % (zList[0]/pi*180))
ax1a.set_xlim([1,20])
ax2a = figa.add_subplot(412)
ax2a.plot(eList, sum(prob[1],0))
ax2a.set_ylabel(r'zenith angle $%d^\circ$' % (zList[1]/pi*180))
ax2a.set_xlim([1,20])
ax3a = figa.add_subplot(413)
ax3a.plot(eList, sum(prob[2],0))
ax3a.set_ylabel(r'zenith angle $%d^\circ$' % (zList[2]/pi*180))
ax3a.set_xlim([1,20])
ax4a = figa.add_subplot(414)
ax4a.plot(eList, sum(prob[3],0))
ax4a.set_ylabel(r'zenith angle $%d^\circ$' % (zList[3]/pi*180))
ax4a.set_xlim([1,20])
ax4a.set_xlabel("neutrino energy / GeV")
savefig("example-standAlone2.png", dpi=160)
print( " \n Completed without fatal errors! \n " )
Index: trunk/README
===================================================================
--- trunk/README (revision 21)
+++ trunk/README (revision 22)
@@ -1,123 +1,126 @@
================================================================================
nuCraft: Calculation of oscillation probabilities of atmospheric neutrinos
================================================================================
Welcome to nuCraft!
NuCraft is a Python project for the calculation of neutrino oscillation
probabilities, with an emphasis on atmospheric neutrinos. For more detailed
information and the most recent version, check out the official homepage:
http://nucraft.hepforge.org
As this is purely Python, installation is not necessary; you can just import
the project and get started.
Documentation can be found in the docstrings, with some additional explanations
in inline comments.
-A paper describing details of the implementation has been submitted to Computer
-Physics Communications for publication, with a preprint available at
+A paper describing details of the implementation has been published in Computer
+Physics Communications 197 (2015) 185-189, DOI: 10.1016/j.cpc.2015.07.010, with
+a preprint available at
http://arxiv.org/abs/1409.1387.
If any questions remain, please contact
mwallraff#physik.rwth-aachen.de
External dependencies:
------------------------
NuCraft heavily relies on NumPy and SciPy. It has been tested with several
versions of NumPy, beginning at 1.5.1, and SciPy, beginning at 0.8.0. It is
possible that it will work well with earlier versions of both packages.
-NuCraft is compatible with both Python 2 (2.6+) and Python 3 (tested with 3.3).
+NuCraft is compatible with both Python 2 (2.6+) and Python 3 (tested with 3.3
+and 3.7.7).
It will not work with versions earlier than 2.6 without modifications, but it
might work with earlier versions of Python 3.
-The example script example-standAlone.py also relies on matplotlib for plotting,
-and it uses TeX inside matplotlib so that dvipng needs to be installed (optional
-requirement of matplotlib).
+The example script example-standAlone.py also relies on matplotlib for plotting.
+Very early versions of matplotlib (<1.5) might produce garbled axis labeling and
+require the UseTeX flag in matplotlib so that dvipng needs to be installed
+(optional requirement of matplotlib).
Program Structure:
--------------------
The structure of nuCraft is fairly easy. The project consists of two classes,
NuCraft and EarthModel. For usage, an instance of NuCraft has to be created with
a set of oscillation parameters, which optionally accepts an instance of
EarthModel at creation for custom Earth density profiles (see next section), a
custom value of the depth of the detector below the ground, and a custom value
for the production height of the neutrino in the atmosphere that is only used if
advanced handling of the atmosphere is disabled later on.
For actual calculations, nuCraft has the two methods CalcWeights and
CalcWeightsLegacy. For the detailed difference between them, have a look at the
paper; in the mean time, do not use the legacy version, it is almost always
worse.
CalcWeights is the only method that should be relevant for normal users. It
accepts the actual particle quantities relevant for oscillations and returns the
oscillation probabilities. For supported input formats and parameters, please
have a look at the docstrings, or see the example script.
Using non-standard Earth models:
----------------------------------
The default Earth model used by nuCraft is the Preliminary Reference Earth Model
PREM (http://dx.doi.org/10.1016/0031-9201(81)90046-7), which is usually loaded
out of the dictionary 'models' in the class EarthModel. If you want to use other
density models, you can either add them to that dictionary, or you can load them
form a text file. The example file EarthModelPrem.txt is shipped with nuCraft as
a template.
Data from files is read in through the ast.literal_eval method. It works
similar to the built-in eval function, but is designed to be inherently safe.
Nonetheless, a check has been kept in place to raise an exception if there are
any underscores in the input file, as all known non-trivial exploits of the
standard eval function require double underscores. If you encounter this
exception, chances are good that there is merily an underscore in the comments;
remove it, and everything should work.
If you want to use nuCraft for reactor neutrino experiments or other not
spherically symmetric mass density distributions, you have to do some simple
modifications to the code:
First, you have create a suitable Earth model for your needs. In case of reactor
neutrinos, you might want to use a matter profile that does not depend on the
distance to the center of Earth, but on the travelled distance. In that case,
substitute profX with an array containing distances, profY with corresponding
matter densities, and set EarthModel.y to (0.4957, 0.4957, 0.4957) (assuming
that the neutrino travelled only through mantle and crust).
Now, EarthModel.profInt is not a function of radius anymore, so go to lines 582
to 586 and substitute profInt(r) by profInt(l).
Similarly, you can modify EarthModel such that EarthModel.profInt depends on
multiple parameters, e.g., radius and zenith angle, and modify the same lines to
call it correctly.
For reactor neutrinos, don't forget to turn the atmospheric handling to 0 (fixed
altitude) and set NuCraft.detectorDepth and NuCraft.atmHeight to 0.
Alternatively, you can modify NuCraft.InteractionAlt to model the neutrino
production region inside the reactor to smear out very-short-baseline effects.
Feel free to write a mail if you need assistance.
List of files:
----------------
- README
this file
- RELEASE_NOTES
comprehensive changelog
- NuCraft.py
actual nuCraft source code, containing the classes NuCraft and EarthModel
- EarthModelPrem.txt
example/template for files to be read in by the EarthModel class
- example-standAlone.py
example script demonstrating the usage of nuCraft by reproducing a plot from
the Akhmedov/Rassaque/Smirnov paper arxiv:1205.7071; can also be used as test
- __init__.py
package initiation file, required to facilitate proper importing

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 7:05 PM (1 d, 12 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805756
Default Alt Text
(17 KB)

Event Timeline