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 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 " )