add a script to calculate the cc, nc, and total event rate

This commit adds a script called print-event-rate which calculates the
event rate per year from the GENIE MCPL files. The livetime comes from
the autosno run_info.log file. The output of the script is a latex table
that I included in my thesis.

1 files changed, 116 insertions, 0 deletions

diff --git a/utils/print-event-rate b/utils/print-event-rate
new file mode 100755
index 0000000..691e806
--- /dev/null
+++ b/utils/print-event-rate
@@ -0,0 +1,116 @@
+#!/usr/bin/env python
+# Copyright (c) 2019, Anthony Latorre <tlatorre at uchicago>
+#
+# This program is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the Free
+# Software Foundation, either version 3 of the License, or (at your option)
+# any later version.
+#
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
+# more details.
+#
+# You should have received a copy of the GNU General Public License along with
+# this program. If not, see <https://www.gnu.org/licenses/>.
+"""
+Script to print the event rate per year for different neutrino interactions. It
+works by reading the MCPL files produced by GENIE to determine the number of
+neutrino interactions and dividing by the livetime in the autosno run_info.log
+file. To run it, just run:
+
+    $ ./print-event-rate /path/to/mcpl/files/*.dat.gz --run-info /path/to/autosno/log/run_info.log
+
+which will print out a latex style table of the event rates.
+"""
+from __future__ import print_function, division
+import numpy as np
+import pandas as pd
+from sddm import read_hdf, AV_RADIUS, PSUP_RADIUS
+import gzip
+
+def read_mcpl(filename):
+    data = np.genfromtxt(filename,skip_header=1)
+    rv = {}
+    n = 1
+    i = 1
+    with gzip.open(filename,"r") as f:
+        lines = f.readlines()
+        nevents = int(lines[0].split()[0])
+        while True:
+            nparticles = int(data[i-1,0])
+            rv[n] = np.genfromtxt(lines[i:i+nparticles])
+            n += 1
+            i += nparticles
+
+            if i - 1 >= data.shape[0]:
+                break
+
+    assert nevents == len(rv)
+
+    return rv
+
+def get_mcgn(filename):
+    ev = read_hdf(filename, "ev")
+    mcgn = read_hdf(filename, "mcgn").assign(run=ev.run[0])
+    mcgn = mcgn.groupby(['run','evn']).first().reset_index()
+    return mcgn
+
+if __name__ == '__main__':
+    import argparse
+    import sys
+    from os.path import join, split
+    from sddm import *
+
+    parser = argparse.ArgumentParser("recreate MCPL files")
+    parser.add_argument("filenames", nargs='+', help="mcpl filenames")
+    parser.add_argument("--run-info", required=True, help="run_info.log autosno file")
+    args = parser.parse_args()
+
+    cc_av = {nu: 0 for nu in range(30,36)}
+    cc_psup = {nu: 0 for nu in range(30,36)}
+    nc_av = {nu: 0 for nu in range(30,36)}
+    nc_psup = {nu: 0 for nu in range(30,36)}
+    livetime = 0.0
+    run_info = np.genfromtxt(args.run_info,usecols=range(4),dtype=(np.int,np.int,np.double,np.double))
+    for filename in args.filenames:
+        head, tail = split(filename)
+        try:
+            mcpl = read_mcpl(filename)
+        except Exception as e:
+            print_warning("failed to read mcpl file '%s'" % filename)
+        run = int(tail.split('_')[3])
+        for ev in mcpl.values():
+            ev = np.atleast_2d(ev)
+            nu = ev[0,14]
+            r = np.linalg.norm(ev[0,16:19])
+            if r < AV_RADIUS:
+                if nu in ev[:,1]:
+                    # nc event
+                    nc_av[nu] += 1
+                else:
+                    # cc event
+                    cc_av[nu] += 1
+            if r < PSUP_RADIUS:
+                if nu in ev[:,1]:
+                    # nc event
+                    nc_psup[nu] += 1
+                else:
+                    # cc event
+                    cc_psup[nu] += 1
+        for i in range(run_info.shape[0]):
+            if run_info[i][0] == run:
+                livetime += run_info[i][3]
+
+    for nu in range(30,36):
+        cc_av[nu] *= 3600*24*365/livetime/100
+        cc_psup[nu] *= 3600*24*365/livetime/100
+        nc_av[nu] *= 3600*24*365/livetime/100
+        nc_psup[nu] *= 3600*24*365/livetime/100
+
+    for nu in range(30,36):
+        line = [cc_av[nu],nc_av[nu],
+                cc_av[nu]+nc_av[nu],
+                cc_psup[nu],nc_psup[nu],
+                cc_psup[nu]+nc_psup[nu]]
+        print(' & '.join(["%.1f" % rate for rate in line]))