From be3b6e19510ddee930df1d9435368d66c18da551 Mon Sep 17 00:00:00 2001
From: tlatorre <tlatorre@uchicago.edu>
Date: Tue, 15 Dec 2020 14:09:26 -0600
Subject: add code to reweight the tau neutrino events

This commit updates the code to reweight the MC data from tau neutrinos
since I stupidly simulated the muon neutrino flux instead of the tau
neutrino flux.
---
 utils/print-event-rate | 44 ++++++++++++--------------------------------
 1 file changed, 12 insertions(+), 32 deletions(-)

(limited to 'utils/print-event-rate')

diff --git a/utils/print-event-rate b/utils/print-event-rate
index 691e806..6350ea8 100755
--- a/utils/print-event-rate
+++ b/utils/print-event-rate
@@ -28,27 +28,7 @@ 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
+from sddm.renormalize import *
 
 def get_mcgn(filename):
     ev = read_hdf(filename, "ev")
@@ -76,28 +56,28 @@ if __name__ == '__main__':
     for filename in args.filenames:
         head, tail = split(filename)
         try:
-            mcpl = read_mcpl(filename)
+            mcpl = read_mcpl_df(filename)
         except Exception as e:
             print_warning("failed to read mcpl file '%s'" % filename)
+            continue
         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])
+        for _, row in mcpl.iterrows():
+            nu = row['nu_id']
+            r = row['nu_r']
             if r < AV_RADIUS:
-                if nu in ev[:,1]:
+                if not row['cc']:
                     # nc event
-                    nc_av[nu] += 1
+                    nc_av[nu] += row['flux_weight']
                 else:
                     # cc event
-                    cc_av[nu] += 1
+                    cc_av[nu] += row['flux_weight']
             if r < PSUP_RADIUS:
-                if nu in ev[:,1]:
+                if not row['cc']:
                     # nc event
-                    nc_psup[nu] += 1
+                    nc_psup[nu] += row['flux_weight']
                 else:
                     # cc event
-                    cc_psup[nu] += 1
+                    cc_psup[nu] += row['flux_weight']
         for i in range(run_info.shape[0]):
             if run_info[i][0] == run:
                 livetime += run_info[i][3]
-- 
cgit