diff options
| -rw-r--r-- | .gitignore | 2 | ||||
| -rw-r--r-- | utils/.gitignore | 2 | ||||
| -rwxr-xr-x | utils/apply-atmospheric-oscillations | 322 | ||||
| -rwxr-xr-x | utils/calculate-atmospheric-oscillations | 116 | ||||
| -rwxr-xr-x | utils/plot-atmospheric-fluxes | 285 | ||||
| -rwxr-xr-x | utils/plot-atmospheric-oscillations | 114 |
6 files changed, 841 insertions, 0 deletions
@@ -1,5 +1,7 @@ *.swp *.swo +*.swm +*.swn *.o *.zdab *.txt diff --git a/utils/.gitignore b/utils/.gitignore new file mode 100644 index 0000000..e307f16 --- /dev/null +++ b/utils/.gitignore @@ -0,0 +1,2 @@ +*.eps +*.pdf diff --git a/utils/apply-atmospheric-oscillations b/utils/apply-atmospheric-oscillations new file mode 100755 index 0000000..3d9f32d --- /dev/null +++ b/utils/apply-atmospheric-oscillations @@ -0,0 +1,322 @@ +#!/usr/bin/env python +""" +Script to combine the low energy Battistoni atmospheric flux with the high +energy flux from Giles Barr and apply oscillations. + +Note: There are a couple of potential issues with this: + +1. The Battistoni fluxes give a total flux integrated over the zenith angle +whereas we need to account for the zenith dependence. Therefore, we just use +the zenith dependence of the lowest energy in the Barr flux files. + +2. I only have access to the Battistoni fluxes for electron neutrinos. +Therefore I will assume that the muon neutrino flux is twice as big which seems +to be roughly correct based on the paper "The Atmospheric Neutrino Flux Below +100 MeV: The FLUKA Results". + +To run it: + + $ ./apply-atmospheric-oscillations --flux-dir ~/atm-production/fluxes/irc01_plus_lowe \ + --osc-dir ~/atm-production/nucraft_osc +""" +from __future__ import print_function, division +import numpy as np + +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 + import matplotlib.pyplot as plt + from os.path import split, splitext, join + from scipy.interpolate import RectBivariateSpline + from scipy.integrate import dblquad, nquad + import sys + + parser = argparse.ArgumentParser("combine low energy and high energy atmospheric fluxes and apply oscillations") + parser.add_argument("--flux-dir",required=True,help="atmospheric production directory") + parser.add_argument("--osc-dir",required=True,help="directory with atmospheric oscillations") + args = parser.parse_args() + + # Energy and cos(zenith) bins for the fluxes from Barr + barr_cos_theta_bins = np.linspace(-1,1,21) + barr_energy_bins = np.logspace(-1,1,41) + + # Final energy bins + # We go here from 10 MeV -> 10 GeV with 20 points per decade like the Barr + # fluxes + energy_bins = np.logspace(-2,1,61) + cos_theta_bins = barr_cos_theta_bins.copy() + + # Read in Battistoni fluxes + nue_battistoni = np.genfromtxt(join(args.flux_dir,"lowe/nue.dat")) + nbe_battistoni = np.genfromtxt(join(args.flux_dir,"lowe/nbe.dat")) + + # Convert Battistoni fluxes from MeV -> GeV + nue_battistoni[:,0] /= 1000.0 + nue_battistoni[:,1] *= 1000.0 + nbe_battistoni[:,0] /= 1000.0 + nbe_battistoni[:,1] *= 1000.0 + + # Convert Battistoni fluxes from cm^-2 -> m^-2 + nue_battistoni[:,1] *= 100.0**2 + nbe_battistoni[:,1] *= 100.0**2 + + # Assume muon neutrino flux below 100 MeV is twice that of electron + # neutrino flux + num_battistoni = nue_battistoni.copy() + num_battistoni[:,1] *= 2 + nbm_battistoni = nbe_battistoni.copy() + nbm_battistoni[:,1] *= 2 + + # Read in Barr fluxes + nue_barr = np.genfromtxt(join(args.flux_dir,"irc01/fmax20_i0403z.sno_nue")) + nbe_barr = np.genfromtxt(join(args.flux_dir,"irc01/fmax20_i0403z.sno_nbe")) + num_barr = np.genfromtxt(join(args.flux_dir,"irc01/fmax20_i0403z.sno_num")) + nbm_barr = np.genfromtxt(join(args.flux_dir,"irc01/fmax20_i0403z.sno_nbm")) + + # Read in oscillation probabilities + nue_osc_prob = np.genfromtxt(join(args.osc_dir,"nue_osc_prob.txt")) + nbe_osc_prob = np.genfromtxt(join(args.osc_dir,"nbe_osc_prob.txt")) + num_osc_prob = np.genfromtxt(join(args.osc_dir,"num_osc_prob.txt")) + nbm_osc_prob = np.genfromtxt(join(args.osc_dir,"nbm_osc_prob.txt")) + + # Reshape oscillation probability arrays. They should now look like: + # + # [[[0.01, -1, P(nue), P(num), P(nut)], + # [0.01, -0.8, P(nue), P(num), P(nut)], + # ... + # [0.01, 1.0, P(nue), P(num), P(nut)]], + # [[0.02, -1, P(nue), P(num), P(nut)], + # [0.02, -0.8, P(nue), P(num), P(nut)], + # ... + # [0.02, 1.0, P(nue), P(num), P(nut)]], + # ... + # [[10.0, -1, P(nue), P(num), P(nut)], + # [10.0, -0.8, P(nue), P(num), P(nut)], + # ... + # [10.0, 1.0, P(nue), P(num), P(nut)]]] + shape0 = len(np.unique(nue_osc_prob[:,0])) + nue_osc_prob = nue_osc_prob.reshape((shape0,-1,5)) + nbe_osc_prob = nbe_osc_prob.reshape((shape0,-1,5)) + num_osc_prob = num_osc_prob.reshape((shape0,-1,5)) + nbm_osc_prob = nbm_osc_prob.reshape((shape0,-1,5)) + + # convert dn/dlnE -> dn/dE + nue_barr[:,2] /= nue_barr[:,0] + nbe_barr[:,2] /= nbe_barr[:,0] + num_barr[:,2] /= num_barr[:,0] + nbm_barr[:,2] /= nbm_barr[:,0] + + # Reshape Barr flux arrays. They should now look like: + # + # [[[0.1, -1, flux,...], + # [0.2, -1, flux,...], + # ... + # [10.0, -1, flux,...]], + # + # [[0.1, -0.9, flux,...], + # [0.2, -0.9, flux,...], + # ... + # [10.0, -0.9, flux,...]], + # ... + # ... + # [[0.1, 1.0, flux,...], + # [0.2, 1.0, flux,...], + # ... + # [10.0, 1.0, flux,...]]] + shape1 = len(np.unique(nue_barr[:,0])) + nue_barr = nue_barr.reshape((-1,shape1,5)) + nbe_barr = nbe_barr.reshape((-1,shape1,5)) + num_barr = num_barr.reshape((-1,shape1,5)) + nbm_barr = nbm_barr.reshape((-1,shape1,5)) + + nue_zenith_dist_x = nue_barr[:,0,1].copy() + nue_zenith_dist_y = nue_barr[:,0,2].copy() + nue_zenith_dist_y /= np.trapz(nue_zenith_dist_y,x=nue_zenith_dist_x) + nbe_zenith_dist_x = nbe_barr[:,0,1].copy() + nbe_zenith_dist_y = nbe_barr[:,0,2].copy() + nbe_zenith_dist_y /= np.trapz(nbe_zenith_dist_y,x=nbe_zenith_dist_x) + num_zenith_dist_x = num_barr[:,0,1].copy() + num_zenith_dist_y = num_barr[:,0,2].copy() + num_zenith_dist_y /= np.trapz(num_zenith_dist_y,x=num_zenith_dist_x) + nbm_zenith_dist_x = nbm_barr[:,0,1].copy() + nbm_zenith_dist_y = nbm_barr[:,0,2].copy() + nbm_zenith_dist_y /= np.trapz(nbm_zenith_dist_y,x=nbm_zenith_dist_x) + + nue_battistoniz = np.empty((nue_barr.shape[0],nue_battistoni.shape[0],4)) + nue_battistoniz[:,:,0] = nue_battistoni[:,0] + nue_battistoniz[:,:,1] = nue_zenith_dist_x[:,np.newaxis] + nue_battistoniz[:,:,2] = nue_battistoni[:,1]*nue_zenith_dist_y[:,np.newaxis]/(2*np.pi) + nue_battistoniz[:,:,3] = nue_battistoniz[:,:,2]*nue_battistoniz[:,:,0] + + nbe_battistoniz = np.empty((nbe_barr.shape[0],nbe_battistoni.shape[0],4)) + nbe_battistoniz[:,:,0] = nbe_battistoni[:,0] + nbe_battistoniz[:,:,1] = nbe_zenith_dist_x[:,np.newaxis] + nbe_battistoniz[:,:,2] = nbe_battistoni[:,1]*nbe_zenith_dist_y[:,np.newaxis]/(2*np.pi) + nbe_battistoniz[:,:,3] = nbe_battistoniz[:,:,2]*nbe_battistoniz[:,:,0] + + num_battistoniz = np.empty((num_barr.shape[0],num_battistoni.shape[0],4)) + num_battistoniz[:,:,0] = num_battistoni[:,0] + num_battistoniz[:,:,1] = num_zenith_dist_x[:,np.newaxis] + num_battistoniz[:,:,2] = num_battistoni[:,1]*num_zenith_dist_y[:,np.newaxis]/(2*np.pi) + num_battistoniz[:,:,3] = num_battistoniz[:,:,2]*num_battistoniz[:,:,0] + + nbm_battistoniz = np.empty((nbm_barr.shape[0],nbm_battistoni.shape[0],4)) + nbm_battistoniz[:,:,0] = nbm_battistoni[:,0] + nbm_battistoniz[:,:,1] = nbm_zenith_dist_x[:,np.newaxis] + nbm_battistoniz[:,:,2] = nbm_battistoni[:,1]*nbm_zenith_dist_y[:,np.newaxis]/(2*np.pi) + nbm_battistoniz[:,:,3] = nbm_battistoniz[:,:,2]*nbm_battistoniz[:,:,0] + + nue_total = np.empty((nue_barr.shape[0],nue_barr.shape[1]+nue_battistoniz.shape[1],4)) + nue_total[:,:nue_battistoniz.shape[1],:] = nue_battistoniz + nue_total[:,nue_battistoniz.shape[1]:,:] = nue_barr[:,:,:4] + + nbe_total = np.empty((nbe_barr.shape[0],nbe_barr.shape[1]+nbe_battistoniz.shape[1],4)) + nbe_total[:,:nbe_battistoniz.shape[1],:] = nbe_battistoniz + nbe_total[:,nbe_battistoniz.shape[1]:,:] = nbe_barr[:,:,:4] + + num_total = np.empty((num_barr.shape[0],num_barr.shape[1]+num_battistoniz.shape[1],4)) + num_total[:,:num_battistoniz.shape[1],:] = num_battistoniz + num_total[:,num_battistoniz.shape[1]:,:] = num_barr[:,:,:4] + + nbm_total = np.empty((nbm_barr.shape[0],nbm_barr.shape[1]+nbm_battistoniz.shape[1],4)) + nbm_total[:,:nbm_battistoniz.shape[1],:] = nbm_battistoniz + nbm_total[:,nbm_battistoniz.shape[1]:,:] = nbm_barr[:,:,:4] + + flux_arrays = {12: nue_total, -12: nbe_total, + 14: num_total, -14: nbm_total} + + # Save unoscillated fluxes + for p in [12,-12,14,-14]: + np.savetxt("fmax20_i0403z_plus_lowe.sno_%s" % pdg_code_to_string(p),flux_arrays[p].reshape((flux_arrays[p].shape[0]*flux_arrays[p].shape[1],-1)), + fmt='%10.4f %10.3f %20.8g %20.8g', + header="Energy (GeV), Cos(zenith), Flux (m^-2 s^-1 Gev^-1 steradian^-1), Flux (dn/dlnE)") + + osc_prob = {12: nue_osc_prob, -12: nbe_osc_prob, + 14: num_osc_prob, -14: nbm_osc_prob} + + def get_flux_interp(nu_type): + flux_array = flux_arrays[nu_type] + flux = flux_array[:,:,2] + e = flux_array[0,:,0] + cos_theta = flux_array[:,0,1] + + # Use first order splines here since we need to enforce the fact that + # the flux is not negative + return RectBivariateSpline(e,cos_theta,flux.T,kx=1,ky=1) + + def get_osc_interp(nu_type1,nu_type2): + prob = osc_prob[nu_type1] + e = prob[:,0,0] + cos_theta = prob[0,:,1] + + if np.sign(nu_type1) != np.sign(nu_type2): + raise ValueError("asking for oscillation probability from %s -> %s!" % tuple(map(pdg_code_to_string,(nu_type1,nu_type2)))) + + # Use second order splines here since oscillation probabilities are + # pretty smooth + if abs(nu_type2) == 12: + return RectBivariateSpline(e,cos_theta,prob[:,:,2],kx=2,ky=2) + elif abs(nu_type2) == 14: + return RectBivariateSpline(e,cos_theta,prob[:,:,3],kx=2,ky=2) + elif abs(nu_type2) == 16: + return RectBivariateSpline(e,cos_theta,prob[:,:,4],kx=2,ky=2) + + def get_oscillated_flux(nu_type1, nu_type2, elo, emid, ehi, zlo, zmid, zhi): + """ + Returns the average flux of neutrinos coming from `nu_type1` neutrinos + oscillating to `nu_type2` neutrinos for neutrinos in the bin between + energies elo and ehi and with a cosine of the zenith angle between zlo + and zhi. + + Note: We also pass the midpoint of the bin for both energy and + cos(theta) as `emid` and `zmid`. The reason for this is that at least + for the Barr fluxes, the midpoint of the bin was chosen as the midpoint + of the log of the left and right edge instead of the actual middle of + the bin. Therefore to be consistent when interpolating these fluxes we + should use the same value. + + The returned flux has units of 1/(m^2 sec steradian GeV). + """ + f_flux1 = get_flux_interp(nu_type1) + f_osc1 = get_osc_interp(nu_type1,nu_type2) + + def f_flux(e,z): + # Not sure exactly how scipy deals with boundary issues in + # RectBivariateSpline so we just make sure everything is within the + # range of the values given in the Barr and Battistoni fluxes + if e < 0.0106: + e = 0.0106 + if z < -0.95: + z = -0.95 + if z > 0.95: + z = 0.95 + return f_flux1(e,z) + + def f_osc(e,z): + # FIXME: Should remove this + if e < 0.02: + e = 0.02 + return f_osc1(e,z) + + # Here we have two options for calculating the oscillated flux: + # + # 1. We interpolate the value of the flux and then multiply by an + # average oscillation probability. This is more correct in the sense + # that the original flux values are actually a histogram and so it + # makes sense to just sample the flux at the middle of the bin. + # + # 2. We integrate the product of the flux and the oscillation + # probability over the whole bin. This is technically more correct + # since it will weight the oscillation probability by the flux, but it + # may have issues near the edge of the first and last bin since we have + # to extrapolate. + # + # For now, we do 1 since it is much faster. To do the second you can + # just comment the next line. + return f_flux(emid,zmid)*f_osc1.integral(elo,ehi,zlo, zhi)/(ehi-elo)/(zhi-zlo) + return nquad(lambda e, z: f_flux(e,z)*f_osc(e,z), [(elo,ehi), (zlo, zhi)], opts={'limit':1000,'epsrel':1e-3})[0]/(ehi-elo)/(zhi-zlo) + + data = {nu_type: np.zeros((cos_theta_bins.shape[0]-1,energy_bins.shape[0]-1,4)) for nu_type in [12,-12,14,-14,16,-16]} + for nu_type1 in [12,-12,14,-14]: + if nu_type1 > 0: + nu_type2s = [12,14,16] + else: + nu_type2s = [-12,-14,-16] + for nu_type2 in nu_type2s: + print("Calculating neutrino flux for %s -> %s" % tuple(map(pdg_code_to_string,(nu_type1,nu_type2)))) + for i, (elo, ehi) in enumerate(zip(energy_bins[:-1],energy_bins[1:])): + for j, (zlo, zhi) in enumerate(zip(cos_theta_bins[:-1],cos_theta_bins[1:])): + print("\r%i/%i" % (i*(len(cos_theta_bins)-1)+j+1,(len(cos_theta_bins)-1)*(len(energy_bins)-1)),end="") + sys.stdout.flush() + # We'll follow the Barr convention here and take the midpoint in log space + emid = np.exp((np.log(elo)+np.log(ehi))/2) + zmid = (zlo + zhi)/2 + data[nu_type2][j,i,0] = emid + data[nu_type2][j,i,1] = zmid + data[nu_type2][j,i,2] += get_oscillated_flux(nu_type1,nu_type2,elo,emid,ehi,zlo,zmid,zhi) + print() + + for nu_type2 in data: + for i, (elo, ehi) in enumerate(zip(energy_bins[:-1],energy_bins[1:])): + for j, (zlo, zhi) in enumerate(zip(cos_theta_bins[:-1],cos_theta_bins[1:])): + data[nu_type2][j,i,3] = data[nu_type2][j,i,2]*data[nu_type2][j,i,0] + + for p in [12,-12,14,-14,16,-16]: + np.savetxt("fmax20_i0403z_plus_lowe.sno_%s.osc" % pdg_code_to_string(p),data[p].reshape((data[p].shape[0]*data[p].shape[1],-1)), + fmt='%10.4f %10.3f %20.8g %20.8g', + header="Energy (GeV), Cos(zenith), Flux (m^-2 s^-1 Gev^-1 steradian^-1), Flux (dn/dlnE), Junk") 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)") diff --git a/utils/plot-atmospheric-fluxes b/utils/plot-atmospheric-fluxes new file mode 100755 index 0000000..fb4861d --- /dev/null +++ b/utils/plot-atmospheric-fluxes @@ -0,0 +1,285 @@ +#!/usr/bin/env python +""" +Script for plotting solar fluxes from +http://www-pnp.physics.ox.ac.uk/~barr/fluxfiles/0403i/index.html. +""" +from __future__ import print_function, division +import numpy as np +import os + +def splitext(path): + """ + Like os.path.splitext() except it returns the full extension if the + filename has multiple extensions, for example: + + splitext('foo.tar.gz') -> 'foo', '.tar.gz' + """ + full_root, full_ext = os.path.splitext(path) + while True: + root, ext = os.path.splitext(full_root) + if ext: + full_ext = ext + full_ext + full_root = root + else: + break + + return full_root, full_ext + +# Taken from https://raw.githubusercontent.com/mwaskom/seaborn/c73055b2a9d9830c6fbbace07127c370389d04dd/seaborn/utils.py +def despine(fig=None, ax=None, top=True, right=True, left=False, + bottom=False, offset=None, trim=False): + """Remove the top and right spines from plot(s). + + fig : matplotlib figure, optional + Figure to despine all axes of, default uses current figure. + ax : matplotlib axes, optional + Specific axes object to despine. + top, right, left, bottom : boolean, optional + If True, remove that spine. + offset : int or dict, optional + Absolute distance, in points, spines should be moved away + from the axes (negative values move spines inward). A single value + applies to all spines; a dict can be used to set offset values per + side. + trim : bool, optional + If True, limit spines to the smallest and largest major tick + on each non-despined axis. + + Returns + ------- + None + + """ + # Get references to the axes we want + if fig is None and ax is None: + axes = plt.gcf().axes + elif fig is not None: + axes = fig.axes + elif ax is not None: + axes = [ax] + + for ax_i in axes: + for side in ["top", "right", "left", "bottom"]: + # Toggle the spine objects + is_visible = not locals()[side] + ax_i.spines[side].set_visible(is_visible) + if offset is not None and is_visible: + try: + val = offset.get(side, 0) + except AttributeError: + val = offset + _set_spine_position(ax_i.spines[side], ('outward', val)) + + # Potentially move the ticks + if left and not right: + maj_on = any( + t.tick1line.get_visible() + for t in ax_i.yaxis.majorTicks + ) + min_on = any( + t.tick1line.get_visible() + for t in ax_i.yaxis.minorTicks + ) + ax_i.yaxis.set_ticks_position("right") + for t in ax_i.yaxis.majorTicks: + t.tick2line.set_visible(maj_on) + for t in ax_i.yaxis.minorTicks: + t.tick2line.set_visible(min_on) + + if bottom and not top: + maj_on = any( + t.tick1line.get_visible() + for t in ax_i.xaxis.majorTicks + ) + min_on = any( + t.tick1line.get_visible() + for t in ax_i.xaxis.minorTicks + ) + ax_i.xaxis.set_ticks_position("top") + for t in ax_i.xaxis.majorTicks: + t.tick2line.set_visible(maj_on) + for t in ax_i.xaxis.minorTicks: + t.tick2line.set_visible(min_on) + + if trim: + # clip off the parts of the spines that extend past major ticks + xticks = ax_i.get_xticks() + if xticks.size: + firsttick = np.compress(xticks >= min(ax_i.get_xlim()), + xticks)[0] + lasttick = np.compress(xticks <= max(ax_i.get_xlim()), + xticks)[-1] + ax_i.spines['bottom'].set_bounds(firsttick, lasttick) + ax_i.spines['top'].set_bounds(firsttick, lasttick) + newticks = xticks.compress(xticks <= lasttick) + newticks = newticks.compress(newticks >= firsttick) + ax_i.set_xticks(newticks) + + xticks_minor = [t for t in ax_i.get_xticks(minor=True) if firsttick < t < lasttick] + ax_i.set_xticks(xticks_minor,minor=True) + + yticks = ax_i.get_yticks() + if yticks.size: + firsttick = np.compress(yticks >= min(ax_i.get_ylim()), + yticks)[0] + lasttick = np.compress(yticks <= max(ax_i.get_ylim()), + yticks)[-1] + ax_i.spines['left'].set_bounds(firsttick, lasttick) + ax_i.spines['right'].set_bounds(firsttick, lasttick) + newticks = yticks.compress(yticks <= lasttick) + newticks = newticks.compress(newticks >= firsttick) + ax_i.set_yticks(newticks) + + yticks_minor = [t for t in ax_i.get_yticks(minor=True) if firsttick < t < lasttick] + ax_i.set_yticks(yticks_minor,minor=True) + +# Data is in the form: Fluxes in bins of neutrino energy (equally spaced bins +# in logE with 10 bins per decade with the low edge of the first bin at 100 +# MeV) and zenith angle (20 bins equally spaced in cos(zenith) with bin width +# 0.1), integrated over azimuth. Note logE means log to base e. Fluxes below 10 +# GeV are from the 3D calculation. + +# The files all have the same format. After the initial comment lines (starting +# with a # character), the files contain one line per bin. No smoothoing +# between bins has been done. The columns are: +# +# 1. Energy at centre of bin in GeV +# 2. Zenith files: Cos(zenith) at centre of bin +# Azimuth files: Azimuth at centre of bin (degrees) +# 3. Flux in dN/dlogE in /m**2/steradian/sec +# 4. MC Statistical error on the flux +# 5. Number of unweighted events entering the bin (not too useful) + +if __name__ == '__main__': + import argparse + import matplotlib + import glob + + parser = argparse.ArgumentParser("plot solar fluxes") + parser.add_argument("filenames", nargs='+', help="filenames of flux files") + parser.add_argument("--save", action="store_true", default=False, help="save plots") + args = parser.parse_args() + + if args.save: + # default \textwidth for a fullpage article in Latex is 16.50764 cm. + # You can figure this out by compiling the following TeX document: + # + # \documentclass{article} + # \usepackage{fullpage} + # \usepackage{layouts} + # \begin{document} + # textwidth in cm: \printinunitsof{cm}\prntlen{\textwidth} + # \end{document} + + width = 16.50764 + width /= 2.54 # cm -> inches + # According to this page: + # http://www-personal.umich.edu/~jpboyd/eng403_chap2_tuftegospel.pdf, + # Tufte suggests an aspect ratio of 1.5 - 1.6. + height = width/1.5 + FIGSIZE = (width,height) + + import matplotlib.pyplot as plt + + font = {'family':'serif', 'serif': ['computer modern roman']} + plt.rc('font',**font) + + plt.rc('text', usetex=True) + else: + # on retina screens, the default plots are way too small + # by using Qt5 and setting QT_AUTO_SCREEN_SCALE_FACTOR=1 + # Qt5 will scale everything using the dpi in ~/.Xresources + import matplotlib + matplotlib.use("Qt5Agg") + + import matplotlib.pyplot as plt + + # Default figure size. Currently set to my monitor width and height so that + # things are properly formatted + FIGSIZE = (13.78,7.48) + + font = {'family':'serif', 'serif': ['computer modern roman']} + plt.rc('font',**font) + + # Make the defalt font bigger + plt.rc('font', size=22) + + plt.rc('text', usetex=True) + + fig = plt.figure(figsize=FIGSIZE) + + colors = plt.rcParams["axes.prop_cycle"].by_key()["color"] + linestyles = ['-','--'] + + def key(filename): + head, tail = os.path.split(filename) + + k = 0 + if tail.startswith('fmax'): + k += 1 + if 'nue' in tail: + k += 10 + elif 'nbe' in tail: + k += 20 + elif 'num' in tail: + k += 30 + elif 'nbm' in tail: + k += 40 + elif 'nut' in tail: + k += 50 + elif 'nbt' in tail: + k += 60 + + return k + + for filename in sorted(args.filenames,key=key): + head, tail = os.path.split(filename) + + print(filename) + + data = np.genfromtxt(filename) + + shape1 = len(np.unique(data[:,0])) + + x = data[:,0].reshape((-1,shape1)) + y = data[:,1].reshape((-1,shape1)) + z = data[:,2].reshape((-1,shape1)) + + # Convert to MeV + x *= 1000.0 + z /= 1000.0 + + zbins = np.linspace(-1,1,21) + dz = zbins[1] - zbins[0] + + x = x[0] + # Integrate over cos(theta) and multiply by 2*pi to convert 3D flux to + # a total flux + y = np.sum(z*dz,axis=0)*2*np.pi + + if 'sno_nue' in tail: + plt.plot(x,y,color=colors[0],linestyle=linestyles[0],label=r'$\nu_e$') + elif 'sno_nbe' in tail: + plt.plot(x,y,color=colors[0],linestyle=linestyles[1],label=r'$\overline{\nu}_e$') + elif 'sno_num' in tail: + plt.plot(x,y,color=colors[1],linestyle=linestyles[0],label=r'$\nu_\mu$') + elif 'sno_nbm' in tail: + plt.plot(x,y,color=colors[1],linestyle=linestyles[1],label=r'$\overline{\nu}_\mu$') + elif 'sno_nut' in tail: + plt.plot(x,y,color=colors[2],linestyle=linestyles[0],label=r'$\nu_\tau$') + elif 'sno_nbt' in tail: + plt.plot(x,y,color=colors[2],linestyle=linestyles[1],label=r'$\overline{\nu}_\tau$') + + plt.gca().set_xscale("log") + plt.gca().set_yscale("log") + despine(fig,trim=True) + plt.xlabel("$E$ (MeV)") + plt.ylabel(r"$\mathrm{d}\Phi/\mathrm{d}E$ (1/$\mathrm{m}^2$/sec/MeV)") + plt.legend() + plt.tight_layout() + + if args.save: + plt.savefig("irc01_atmospheric_flux.pdf") + plt.savefig("irc01_atmospheric_flux.eps") + + plt.show() diff --git a/utils/plot-atmospheric-oscillations b/utils/plot-atmospheric-oscillations new file mode 100755 index 0000000..c24a776 --- /dev/null +++ b/utils/plot-atmospheric-oscillations @@ -0,0 +1,114 @@ +#!/usr/bin/env python +""" +Script to plot the probabilities for atmospheric neutrino oscillations. To run +it: + + $ ./plot-atmospheric-oscillations nue_osc_prob.txt +""" +from __future__ import print_function, division +import numpy as np + +if __name__ == '__main__': + import argparse + from os.path import split, splitext + + parser = argparse.ArgumentParser("script to plot atmospheric oscillations") + parser.add_argument("filenames", nargs='+', help="oscillation probability filenames") + parser.add_argument("--save", action="store_true", default=False, help="save plots") + args = parser.parse_args() + + if args.save: + # default \textwidth for a fullpage article in Latex is 16.50764 cm. + # You can figure this out by compiling the following TeX document: + # + # \documentclass{article} + # \usepackage{fullpage} + # \usepackage{layouts} + # \begin{document} + # textwidth in cm: \printinunitsof{cm}\prntlen{\textwidth} + # \end{document} + + width = 16.50764 + width /= 2.54 # cm -> inches + # According to this page: + # http://www-personal.umich.edu/~jpboyd/eng403_chap2_tuftegospel.pdf, + # Tufte suggests an aspect ratio of 1.5 - 1.6. + height = width/1.5 + FIGSIZE = (width,height) + + import matplotlib.pyplot as plt + + font = {'family':'serif', 'serif': ['computer modern roman']} + plt.rc('font',**font) + + plt.rc('text', usetex=True) + else: + # on retina screens, the default plots are way too small + # by using Qt5 and setting QT_AUTO_SCREEN_SCALE_FACTOR=1 + # Qt5 will scale everything using the dpi in ~/.Xresources + import matplotlib + matplotlib.use("Qt5Agg") + + import matplotlib.pyplot as plt + + # Default figure size. Currently set to my monitor width and height so that + # things are properly formatted + FIGSIZE = (13.78,7.48) + + # Make the defalt font bigger + plt.rc('font', size=22) + + for filename in args.filenames: + head, tail = split(filename) + root, ext = splitext(tail) + + e, z, pnue, pnum, pnut = np.genfromtxt(filename).T + + shape0 = len(np.unique(e)) + + ee = e.reshape((shape0,-1)) + zz = z.reshape((shape0,-1)) + pnue = pnue.reshape((shape0,-1)) + pnum = pnum.reshape((shape0,-1)) + pnut = pnut.reshape((shape0,-1)) + + levels = np.linspace(0,1,101) + + plt.figure(figsize=FIGSIZE) + plt.contourf(ee,zz,pnue,levels=levels) + plt.gca().set_xscale('log') + plt.xlabel("Energy (GeV)") + plt.ylabel("Cos(Zenith)") + plt.colorbar() + plt.tight_layout() + if args.save: + plt.savefig("%s_nue.pdf" % root) + plt.savefig("%s_nue.eps" % root) + else: + plt.title(r"Probability to oscillate to $\nu_e$") + plt.figure(figsize=FIGSIZE) + plt.contourf(ee,zz,pnum,levels=levels) + plt.gca().set_xscale('log') + plt.xlabel("Energy (GeV)") + plt.ylabel("Cos(Zenith)") + plt.colorbar() + plt.tight_layout() + if args.save: + plt.savefig("%s_num.pdf" % root) + plt.savefig("%s_num.eps" % root) + else: + plt.title(r"Probability to oscillate to $\nu_\mu$") + plt.figure(figsize=FIGSIZE) + plt.contourf(ee,zz,pnut,levels=levels) + plt.gca().set_xscale('log') + plt.xlabel("Energy (GeV)") + plt.ylabel("Cos(Zenith)") + plt.colorbar() + plt.tight_layout() + if args.save: + plt.savefig("%s_nut.pdf" % root) + plt.savefig("%s_nut.eps" % root) + else: + plt.title(r"Probability to oscillate to $\nu_\tau$") + + plt.show() |