update chi2 to correct for energy bias

This commit updates the chi2 script to correct for the energy bias of
the reconstruction relative to Monte Carlo.

1 files changed, 14 insertions, 2 deletions

diff --git a/utils/plot-fit-results b/utils/plot-fit-results
index c90c13a..20c3c68 100755
--- a/utils/plot-fit-results
+++ b/utils/plot-fit-results
@@ -28,6 +28,7 @@ if __name__ == '__main__':
 
     parser = argparse.ArgumentParser("plot fit results")
     parser.add_argument("filenames", nargs='+', help="input files")
+    parser.add_argument("-o", "--output", default=None, help="output file")
     parser.add_argument("--save", action="store_true", default=False, help="save plots")
     args = parser.parse_args()
 
@@ -110,11 +111,15 @@ if __name__ == '__main__':
     # 100 bins between 50 MeV and 1 GeV
     bins = np.arange(50,1000,100)
 
+    T = (bins[1:] + bins[:-1])/2
+
     markers = itertools.cycle(('o', 'v')) 
 
     fig3, ax3 = plt.subplots(3,1,num=3,sharex=True)
     fig4, ax4 = plt.subplots(3,1,num=4,sharex=True)
 
+    output = pd.DataFrame({'T':T})
+
     for id in [IDP_E_MINUS, IDP_MU_MINUS]:
         events = data_true[data_true['mcgn_id'] == id]
 
@@ -128,8 +133,6 @@ if __name__ == '__main__':
 
         label = 'Muon' if id == IDP_MU_MINUS else 'Electron'
 
-        T = (bins[1:] + bins[:-1])/2
-
         marker = markers.next()
 
         plt.figure(1)
@@ -142,6 +145,11 @@ if __name__ == '__main__':
         ax3[1].errorbar(T,dy['median'],yerr=dy['median_err'],fmt=marker,label=label)
         ax3[2].errorbar(T,dz['median'],yerr=dz['median_err'],fmt=marker,label=label)
 
+        if id == IDP_E_MINUS:
+            output['e_bias'] = (dT['median']/T).values
+        else:
+            output['u_bias'] = (dT['median']/T).values
+
         ax4[0].errorbar(T,dx['iqr_std'],yerr=dx['iqr_std_err'],fmt=marker,label=label)
         ax4[1].errorbar(T,dy['iqr_std'],yerr=dy['iqr_std_err'],fmt=marker,label=label)
         ax4[2].errorbar(T,dz['iqr_std'],yerr=dz['iqr_std_err'],fmt=marker,label=label)
@@ -152,6 +160,10 @@ if __name__ == '__main__':
         plt.figure(6)
         plt.scatter(events['Te'],events['ratio'],marker=marker,label=label)
 
+    if args.output:
+        with h5py.File(args.output,"w") as f:
+            f.create_dataset('energy_bias',data=output.to_records())
+
     fig = plt.figure(1)
     despine(fig,trim=True)
     plt.xlabel("Kinetic Energy (MeV)")