add scripts to apply neutrino oscillations to the atmospheric fluxes

This commit adds four scripts:

1. calculate-atmospheric-oscillations

This script uses an independent python package called nucraft to calculate the
neutrino oscillation probabilities for atmospheric neutrinos. These
probabilities are calculated as a function of energy and cosine of the zenith
angle and stored in text files.

2. apply-atmospheric-oscillations

This script combines the high energy 2D atmospheric neutrino flux from Barr and
the low energy 1D flux from Battistoni and then applies neutrino oscillations
to them. The results are then stored in new flux files that can be used with a
modified version of GENIE to simulate the oscillated atmospheric neutrino flux
from 10 MeV to 10 GeV.

3. plot-atmospheric-fluxes

This is a simple script to plot the atmospheric flux files produced by
apply-atmospheric-oscillations.

4. plot-atmospheric-oscillations

This is a simple script to plot the 2D neutrino oscillation probabilities.

1 files changed, 116 insertions, 0 deletions

diff --git a/utils/calculate-atmospheric-oscillations b/utils/calculate-atmospheric-oscillations
new file mode 100755
index 0000000..356fc1f
--- /dev/null
+++ b/utils/calculate-atmospheric-oscillations
@@ -0,0 +1,116 @@
+#!/usr/bin/env python
+"""
+Script to calculate probabilities for atmospheric neutrino oscillations. To run
+it:
+
+    $ ./calculate-atmospheric-oscillations --atmosphere-mode 3 --num-prec 5e-4
+
+which will then produce the following files:
+
+nue_osc_prob.txt: probability for nues to oscillate
+nbe_osc_prob.txt: probability for anti nues to oscillate
+num_osc_prob.txt: probability for numus to oscillate
+nbm_osc_prob.txt: probability for anti numus to oscillate
+
+The format of each text file is the same. It has the following columns:
+
+    1. Energy of the neutrino in GeV.
+    2. Cosine of the zenith angle.
+    3. Probability to oscillate into nue (or antinue)
+    4. Probability to oscillate into numu (or antinumu)
+    5. Probability to oscillate into nutau (or antinutau)
+"""
+from __future__ import print_function, division
+import numpy as np
+from time import time
+from NuCraft import *
+
+# depth of the SNO detector (km)
+# currently just converted 6800 feet -> km
+SNO_DEPTH = 2.072
+
+# number of energy bins
+eBins = 1000
+# number of zenith angle bins
+zBins = 100
+# zenith angle bins
+zList = np.arccos(np.linspace(-1,1,zBins))
+# energy range in GeV
+eList = np.logspace(np.log10(0.01), np.log10(10), eBins)
+
+# parameters from PDG 2019
+THETA23 = np.arcsin(np.sqrt(0.512))/np.pi*180.
+THETA13 = np.arcsin(np.sqrt(0.0218))/np.pi*180.
+THETA12 = np.arcsin(np.sqrt(0.307))/np.pi*180.
+DM21   = 7.53e-5
+DM31   = 2.444e-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)],detectorDepth=SNO_DEPTH)
+
+def pdg_code_to_string(pdg):
+    if pdg == 12:
+        return "nue"
+    elif pdg == -12:
+        return "nbe"
+    elif pdg == 14:
+        return "num"
+    elif pdg == -14:
+        return "nbm"
+    elif pdg == 16:
+        return "nut"
+    elif pdg == -16:
+        return "nbt"
+    else:
+        raise ValueError("Unknown pdg code %i" % pdg)
+
+if __name__ == '__main__':
+    import argparse
+
+    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
+
+    parser = argparse.ArgumentParser("script to calculate atmospheric neutrino oscillations")
+    parser.add_argument("--atmosphere-mode", type=int, default=atmosphereMode, help="controls the handling of the atmospheric interaction altitude")
+    parser.add_argument("--num-prec", type=float, default=numPrec, help="parameter which governs the precision with which nuCraft computs the weights")
+    args = parser.parse_args()
+
+
+    # 12, -12:  NuE, NuEBar
+    # 14, -14:  NuMu, NuMuBar
+    # 16, -16:  NuTau, NuTauBar
+    for pType in [12,-12,14,-14]:
+        print("Calculating oscillation probabilities for %s..." % pdg_code_to_string(pType))
+
+        # saving the current time to measure the time needed for the execution of the following code
+        t = time()
+
+        ### using lists (arrays, actually)
+        zz, ee = np.meshgrid(zList, eList)
+        zListLong = zz.flatten()
+        eListLong = ee.flatten()
+        tListLong = np.repeat(pType,len(eListLong))
+
+        # actual call to nuCraft for weight calculations:
+        prob = AkhmedovOsci.CalcWeights((tListLong, eListLong, zListLong), numPrec=args.num_prec, atmMode=args.atmosphere_mode)
+
+        prob = np.array(prob)
+
+        save = np.empty((prob.shape[0],5))
+        save[:,0] = eListLong
+        save[:,1] = np.cos(zListLong)
+        save[:,2:] = prob
+
+        print("Calculating the probabilities took %.0f seconds." % (time()-t))
+
+        np.savetxt("%s_osc_prob.txt" % pdg_code_to_string(pType),save,header="Energy (GeV), Cos(zenith), P(nue), P(num), P(nut)")