Merge pull request #13 from JostMigenda/python3

XS scripts: Python 3 updates and improvements

xscns/scripts/README.md

@@ -0,0 +1,29 @@
+# Cross section scripts
+A collection of various Python scripts to handle cross section files.
+
+## Installation
+All scripts are tested under Python 3.x.
+
+Some scripts require the following additional libraries:
+* [matplotlib](https://matplotlib.org/stable/users/installing.html)
+* [NumPy](https://numpy.org/install/)
+* [PyROOT](https://root.cern/manual/python/)
+* [SciPy](https://scipy.org/install.html)
+
+## Usage information
+
+### compare_xs.py
+Plot multiple cross section files.
+Run `python compare_xs.py` for usage information.
+
+### compute_nu_e_cross_sections.py
+Generate file with cross section for neutrino-electron scattering, based on [DOI:10.1088/0954-3899/29/11/013](https://doi.org/10.1088/0954-3899/29/11/013).
+
+### populate_CEvNS_xs.py
+Generate cross-section files for CEvNS (`xs_coh_helm_{element}_{flavor}`) using the Helm form factor.
+
+### universal_xs.py
+Produce a cross section file for an arbitrary A, Z, and Q. This gives a rough approximation and may only be accurate to within a factor of a few.
+Run `python universal_xs.py` for usage information.
+
+The folder `parameterization/` contains scripts used to generate this universal XS fit.

xscns/generic_xs/compare_xs.py → xscns/scripts/compare_xs.py

@@ -1,10 +1,9 @@
 # Quick, plot a couple cross sections against each other!
-# Will work with python 2.6, but possibly not 3.x
+# Will work with Python 2.6 and ROOT 6.12 or Python 3.x and ROOT 6.22.
 
 import ROOT
 import array
 import sys
-import math
 import code
 
 ROOT.gROOT.SetStyle("Plain")
@@ -18,13 +17,12 @@
 
 flavornum = 1
 filelist = []
-for ifile in range(len(sys.argv)-1):
-   print sys.argv[ifile+1]
-   if str(sys.argv[ifile+1]).isdigit():
-      flavornum = int(sys.argv[ifile+1])
-      continue
+for arg in sys.argv[1:]:
+   print(arg)
+   if arg.isdigit():
+      flavornum = int(arg)
    else:
-      filelist.append(sys.argv[ifile+1])
+      filelist.append(arg)
 
 # quick code to skip yellow, the worst color in ROOT
 def SkipYellow(number):
@@ -36,7 +34,7 @@ def SkipYellow(number):
 # first, a script to read in a cross section file and create a histogram
 def ReadXS(filename, flavornum = 1, xsname = "unnamed"):
    # 1 = nue, 2 = numu, 3 = nutau, 4 = nuebar, etc.
-   print flavornum
+   print(flavornum)
    xs_data = []
    infile = open(filename, "r")
    for line in infile.readlines():
@@ -55,7 +53,6 @@ def ReadXS(filename, flavornum = 1, xsname = "unnamed"):
    for ibin in range (numbins+1):
       binlist.append(minbinlowedge + ibin * logstep)
    binarray = array.array("d", [pow(10, x) for x in binlist])
-   maxbinhighedge = (xs_data[-1][0] + logstep/2.0)
    flavornamelist = ["nue", "numu", "nutau", "nuebar", "numubar", "nutaubar"]
    hist = ROOT.TH1D("xs_" + xsname + "_" + flavornamelist[flavornum-1],
                      xsname + "_" + flavornamelist[flavornum-1] + ";E (GeV);10^{-38} cm^{2}",

xscns/generic_xs/parameterization/optimize_xs_param.py → xscns/scripts/parameterization/optimize_xs_param.py

@@ -1,6 +1,6 @@
 # Try to find a function to optimize a parameterization for neutrino cross sections
 # This optimizes cross section magnitudes for Michel spectrum neutrinos.
-# code initially written and tested with python 2.7.10
+# code initially written and tested with python 2.7.10, also works with 3.x
 # by jba 7/2018
 
 import scipy.optimize as optimize
@@ -42,7 +42,7 @@ def calc_xs(nuc_param_dict, fit_param_list):
    p4 = fit_param_list[4]
    p5 = fit_param_list[5]
 
-   return (p0 * pow(a, p1) + p2 * pow(a-z, p3)) * ((1 + p4 * pow(q, p5)))
+   return (p0 * pow(a, p1) + p2 * pow(a-z, p3)) * (1 + p4 * pow(q, p5))
 
 
 # Function to compute a chi2 for how well we fit a single cross section
@@ -55,7 +55,6 @@ def chi2(xs_datum, fit_param_list):
 # Sum the chi2 over all our data points
 def totalchi2(fit_param_list):
    xs_data = GetData("athar_xs.dat")
-   ndat = len(xs_data)
    sum_chi2 = 0
    for xs_datum in xs_data:
       sum_chi2 += chi2(xs_datum, fit_param_list)

xscns/generic_xs/parameterization/parameterize_shape.py → xscns/scripts/parameterization/parameterize_shape.py

@@ -1,7 +1,7 @@
 # Let's find a parameterization which does a good job of matching "known" cross sections
 # This code works on cross section shape
 # code initially written and tested with python 2.7.10 and ROOT 6.12
-# It may not work for python 3.x due to the ROOT import
+# It also works for Python 3.x and ROOT 6.22.
 # by jba 8/9/2018
 
 import ROOT

xscns/generic_xs/populate_CEvNS_xs.py → xscns/scripts/populate_CEvNS_xs.py

@@ -1,5 +1,5 @@
 ##############################################################################################################
-#	This script, hosted in /xscns/generic_xs/ is intended to populate the xs_coh_helm_{element}_{flavor}
+#	This script, hosted in /xscns/scripts/ is intended to populate the xs_coh_helm_{element}_{flavor}
 #   files which give the GLoBES-required cross-sections for CEvNS interactions. The total cross-section 
 #   of interaction (as a function of incident neutrino energy) is calculated for each isotope, and the
 #   weighted average is calculated using the molar fraction as the weighting factor. The weighted

xscns/generic_xs/universal_xs.py → xscns/scripts/universal_xs.py

@@ -18,11 +18,11 @@
 emin = 0.0002 # GeV
 emax = 0.1 # GeV
 nbins = 1001 # this one cannot change!
-if len(sys.argv) not in [4, 5, 7, 8]:
+if len(sys.argv) not in [4, 5, 6, 8]:
    print("Usage: python ./universal_xs.py Z A Q [+/- for nue/nuebar] [<filename>] [emin emax]")
    print("Q must be in MeV, the optional emin and emax must be in GeV")
    print("These are very approximate cross sections, use with bucket of salt")
-if len(sys.argv) in [4,5,7,8]:
+if len(sys.argv) in [4,5,6,8]:
    Z = int(sys.argv[1])
    A = int(sys.argv[2])
    Q = float(sys.argv[3])