only plot Michels from stopping muons

1 files changed, 129 insertions, 129 deletions

diff --git a/utils/plot-muons b/utils/plot-muons
index 128eaa2..fc8427c 100755
--- a/utils/plot-muons
+++ b/utils/plot-muons
@@ -115,29 +115,6 @@ if __name__ == '__main__':
                Line2D([0], [0], color='C1')]
     labels = ('Data','Monte Carlo')
 
-    # For the Michel energy plot, we only look at events where the
-    # corresponding muon had less than 2500 nhit.  The reason for only looking
-    # at Michel electrons from muons with less than 2500 nhit is because there
-    # is significant ringing and afterpulsing after a large muon which can
-    # cause the reconstruction to overestimate the energy.
-    michel_low_nhit = michel[michel.muon_nhit < 2500]
-
-    print("Particle ID probability for Michel electrons:")
-    print("Data")
-    print_particle_probs(michel_low_nhit)
-    print("Monte Carlo")
-    print_particle_probs(michel_mc)
-
-    fig = plt.figure()
-    plot_hist2_data_mc(michel_low_nhit,michel_mc)
-    despine(fig,trim=True)
-    fig.legend(handles,labels,loc='upper right')
-    if args.save:
-        plt.savefig("michel_electrons.pdf")
-        plt.savefig("michel_electrons.eps")
-    else:
-        plt.suptitle("Michel Electrons")
-
     fig = plt.figure()
     plot_hist2_data_mc(muons,muons_mc)
     despine(fig,trim=True)
@@ -185,118 +162,141 @@ if __name__ == '__main__':
     stopping_muons = stopping_muons[stopping_muons.cos_theta < -0.5]
     stopping_muons_mc = stopping_muons_mc[stopping_muons_mc.cos_theta < -0.5]
 
-    if len(stopping_muons):
-        # project muon to PSUP
-        stopping_muons['dx'] = stopping_muons.apply(get_dx,axis=1)
-        stopping_muons_mc['dx'] = stopping_muons_mc.apply(get_dx,axis=1)
-        # energy based on distance travelled
-        stopping_muons['T_dx'] = dx_to_energy(stopping_muons.dx)
-        stopping_muons_mc['T_dx'] = dx_to_energy(stopping_muons_mc.dx)
-        stopping_muons['dT'] = stopping_muons['ke'] - stopping_muons['T_dx']
-        stopping_muons_mc['dT'] = stopping_muons_mc['ke'] - stopping_muons_mc['T_dx']
-
-        # Plot the energy and angular distribution for external muons
-        fig = plt.figure()
-        plt.subplot(2,1,1)
-        plt.hist(stopping_muons.ke.values, bins=np.logspace(3,7,100), histtype='step', color='C0', label="Data")
-        norm = len(stopping_muons.ke.values)/len(stopping_muons_mc.ke.values)
-        plt.hist(stopping_muons_mc.ke.values, weights=np.tile(norm,len(stopping_muons_mc.ke.values)), bins=np.logspace(3,7,100), histtype='step', color='C1', label="Monte Carlo")
-        plt.legend()
-        plt.xlabel("Energy (MeV)")
-        plt.gca().set_xscale("log")
-        plt.subplot(2,1,2)
-        plt.hist(np.cos(stopping_muons.theta.values), bins=np.linspace(-1,1,100), histtype='step', color='C0', label="Data")
-        norm = len(stopping_muons.theta.values)/len(stopping_muons_mc.theta.values)
-        plt.hist(np.cos(stopping_muons_mc.theta.values), weights=np.tile(norm,len(stopping_muons_mc.theta.values)), bins=np.linspace(-1,1,100), histtype='step', color='C1', label="Monte Carlo")
-        plt.legend()
-        despine(fig,trim=True)
-        plt.xlabel(r"$\cos(\theta)$")
-        plt.tight_layout()
-        if args.save:
-            plt.savefig("stopping_muon_energy_cos_theta.pdf")
-            plt.savefig("stopping_muon_energy_cos_theta.eps")
-        else:
-            plt.suptitle("Stopping Muons")
-
-        print(stopping_muons[['run','gtid','ke','T_dx','dT','gtid_michel','r_michel','ftp_r_michel','id','r']])
-
-        print("Particle ID probability for Stopping Muons:")
-        print("Data")
-        print_particle_probs(stopping_muons)
-        print("Monte Carlo")
-        print_particle_probs(stopping_muons_mc)
-
-        fig = plt.figure()
-        plot_hist2_data_mc(stopping_muons,stopping_muons_mc)
-        despine(fig,trim=True)
-        if len(muons):
-            plt.tight_layout()
-        fig.legend(handles,labels,loc='upper right')
-        if args.save:
-            plt.savefig("stopping_muons.pdf")
-            plt.savefig("stopping_muons.eps")
-        else:
-            plt.suptitle("Stopping Muons")
-
-        fig = plt.figure()
-        plt.hist((stopping_muons['ke']-stopping_muons['T_dx'])*100/stopping_muons['T_dx'], bins=np.linspace(-100,100,200), histtype='step', color='C0', label="Data")
-        plt.hist((stopping_muons_mc['ke']-stopping_muons_mc['T_dx'])*100/stopping_muons_mc['T_dx'], bins=np.linspace(-100,100,200), histtype='step', color='C1', label="Monte Carlo")
-        plt.legend()
-        despine(fig,trim=True)
-        plt.xlabel("Fractional energy difference (\%)")
-        plt.title("Fractional energy difference for Stopping Muons")
-        plt.tight_layout()
-        if args.save:
-            plt.savefig("stopping_muon_fractional_energy_difference.pdf")
-            plt.savefig("stopping_muon_fractional_energy_difference.eps")
-        else:
-            plt.title("Stopping Muon Fractional Energy Difference")
-
-        # 100 bins between 50 MeV and 10 GeV
-        bins = np.linspace(50,2000,10)
-
-        pd_bins = pd.cut(stopping_muons['T_dx'],bins)
-        pd_bins_mc = pd.cut(stopping_muons_mc['T_dx'],bins)
-
-        T = (bins[1:] + bins[:-1])/2
-
-        dT = stopping_muons.groupby(pd_bins)['dT'].agg(['mean','sem','std',std_err,median,median_err,iqr_std,iqr_std_err])
-        dT_mc = stopping_muons_mc.groupby(pd_bins_mc)['dT'].agg(['mean','sem','std',std_err,median,median_err,iqr_std,iqr_std_err])
-
-        fig = plt.figure()
-        plt.errorbar(T, dT['median']*100/T, yerr=dT['median_err']*100/T, color='C0', label="Data")
-        plt.errorbar(T, dT_mc['median']*100/T, yerr=dT_mc['median_err']*100/T, color='C1', label="Monte Carlo")
-        plt.legend()
-        despine(fig,trim=True)
-        plt.xlabel("Kinetic Energy (MeV)")
-        plt.ylabel(r"Energy bias (\%)")
-        plt.tight_layout()
-        if args.save:
-            plt.savefig("stopping_muon_energy_bias.pdf")
-            plt.savefig("stopping_muon_energy_bias.eps")
-        else:
-            plt.title("Stopping Muon Energy Bias")
-
-        fig = plt.figure()
-        plt.errorbar(T, dT['iqr_std']*100/T, yerr=dT['iqr_std_err']*100/T, color='C0', label="Data")
-        plt.errorbar(T, dT_mc['iqr_std']*100/T, yerr=dT_mc['iqr_std_err']*100/T, color='C1', label="Monte Carlo")
-        plt.legend()
-        despine(fig,trim=True)
-        plt.xlabel("Kinetic Energy (MeV)")
-        plt.ylabel(r"Energy resolution (\%)")
+    # project muon to PSUP
+    stopping_muons['dx'] = stopping_muons.apply(get_dx,axis=1)
+    stopping_muons_mc['dx'] = stopping_muons_mc.apply(get_dx,axis=1)
+    # energy based on distance travelled
+    stopping_muons['T_dx'] = dx_to_energy(stopping_muons.dx)
+    stopping_muons_mc['T_dx'] = dx_to_energy(stopping_muons_mc.dx)
+    stopping_muons['dT'] = stopping_muons['ke'] - stopping_muons['T_dx']
+    stopping_muons_mc['dT'] = stopping_muons_mc['ke'] - stopping_muons_mc['T_dx']
+
+    # Plot the energy and angular distribution for external muons
+    fig = plt.figure()
+    plt.subplot(2,1,1)
+    plt.hist(stopping_muons.ke.values, bins=np.logspace(3,7,100), histtype='step', color='C0', label="Data")
+    norm = len(stopping_muons.ke.values)/len(stopping_muons_mc.ke.values)
+    plt.hist(stopping_muons_mc.ke.values, weights=np.tile(norm,len(stopping_muons_mc.ke.values)), bins=np.logspace(3,7,100), histtype='step', color='C1', label="Monte Carlo")
+    plt.legend()
+    plt.xlabel("Energy (MeV)")
+    plt.gca().set_xscale("log")
+    plt.subplot(2,1,2)
+    plt.hist(np.cos(stopping_muons.theta.values), bins=np.linspace(-1,1,100), histtype='step', color='C0', label="Data")
+    norm = len(stopping_muons.theta.values)/len(stopping_muons_mc.theta.values)
+    plt.hist(np.cos(stopping_muons_mc.theta.values), weights=np.tile(norm,len(stopping_muons_mc.theta.values)), bins=np.linspace(-1,1,100), histtype='step', color='C1', label="Monte Carlo")
+    plt.legend()
+    despine(fig,trim=True)
+    plt.xlabel(r"$\cos(\theta)$")
+    plt.tight_layout()
+    if args.save:
+        plt.savefig("stopping_muon_energy_cos_theta.pdf")
+        plt.savefig("stopping_muon_energy_cos_theta.eps")
+    else:
+        plt.suptitle("Stopping Muons")
+
+    print(stopping_muons[['run','gtid','ke','T_dx','dT','gtid_michel','r_michel','ftp_r_michel','id','r']])
+
+    print("Particle ID probability for Stopping Muons:")
+    print("Data")
+    print_particle_probs(stopping_muons)
+    print("Monte Carlo")
+    print_particle_probs(stopping_muons_mc)
+
+    fig = plt.figure()
+    plot_hist2_data_mc(stopping_muons,stopping_muons_mc)
+    despine(fig,trim=True)
+    if len(muons):
         plt.tight_layout()
-        if args.save:
-            plt.savefig("stopping_muon_energy_resolution.pdf")
-            plt.savefig("stopping_muon_energy_resolution.eps")
-        else:
-            plt.title("Stopping Muon Energy Resolution")
+    fig.legend(handles,labels,loc='upper right')
+    if args.save:
+        plt.savefig("stopping_muons.pdf")
+        plt.savefig("stopping_muons.eps")
+    else:
+        plt.suptitle("Stopping Muons")
+
+    fig = plt.figure()
+    plt.hist((stopping_muons['ke']-stopping_muons['T_dx'])*100/stopping_muons['T_dx'], bins=np.linspace(-100,100,200), histtype='step', color='C0', label="Data")
+    plt.hist((stopping_muons_mc['ke']-stopping_muons_mc['T_dx'])*100/stopping_muons_mc['T_dx'], bins=np.linspace(-100,100,200), histtype='step', color='C1', label="Monte Carlo")
+    plt.legend()
+    despine(fig,trim=True)
+    plt.xlabel("Fractional energy difference (\%)")
+    plt.title("Fractional energy difference for Stopping Muons")
+    plt.tight_layout()
+    if args.save:
+        plt.savefig("stopping_muon_fractional_energy_difference.pdf")
+        plt.savefig("stopping_muon_fractional_energy_difference.eps")
+    else:
+        plt.title("Stopping Muon Fractional Energy Difference")
+
+    # 100 bins between 50 MeV and 10 GeV
+    bins = np.linspace(50,2000,10)
+
+    pd_bins = pd.cut(stopping_muons['T_dx'],bins)
+    pd_bins_mc = pd.cut(stopping_muons_mc['T_dx'],bins)
+
+    T = (bins[1:] + bins[:-1])/2
+
+    dT = stopping_muons.groupby(pd_bins)['dT'].agg(['mean','sem','std',std_err,median,median_err,iqr_std,iqr_std_err])
+    dT_mc = stopping_muons_mc.groupby(pd_bins_mc)['dT'].agg(['mean','sem','std',std_err,median,median_err,iqr_std,iqr_std_err])
+
+    fig = plt.figure()
+    plt.errorbar(T, dT['median']*100/T, yerr=dT['median_err']*100/T, color='C0', label="Data")
+    plt.errorbar(T, dT_mc['median']*100/T, yerr=dT_mc['median_err']*100/T, color='C1', label="Monte Carlo")
+    plt.legend()
+    despine(fig,trim=True)
+    plt.xlabel("Kinetic Energy (MeV)")
+    plt.ylabel(r"Energy bias (\%)")
+    plt.tight_layout()
+    if args.save:
+        plt.savefig("stopping_muon_energy_bias.pdf")
+        plt.savefig("stopping_muon_energy_bias.eps")
+    else:
+        plt.title("Stopping Muon Energy Bias")
+
+    fig = plt.figure()
+    plt.errorbar(T, dT['iqr_std']*100/T, yerr=dT['iqr_std_err']*100/T, color='C0', label="Data")
+    plt.errorbar(T, dT_mc['iqr_std']*100/T, yerr=dT_mc['iqr_std_err']*100/T, color='C1', label="Monte Carlo")
+    plt.legend()
+    despine(fig,trim=True)
+    plt.xlabel("Kinetic Energy (MeV)")
+    plt.ylabel(r"Energy resolution (\%)")
+    plt.tight_layout()
+    if args.save:
+        plt.savefig("stopping_muon_energy_resolution.pdf")
+        plt.savefig("stopping_muon_energy_resolution.eps")
+    else:
+        plt.title("Stopping Muon Energy Resolution")
+
+    # For the Michel energy plot, we only look at events where the
+    # corresponding muon had less than 2500 nhit.  The reason for only looking
+    # at Michel electrons from muons with less than 2500 nhit is because there
+    # is significant ringing and afterpulsing after a large muon which can
+    # cause the reconstruction to overestimate the energy.
+    michel_low_nhit = michel[michel.muon_gtid.isin(stopping_muons.gtid.values) & (michel.muon_nhit < 2500)]
+    michel_low_nhit_mc = michel_mc[michel_mc.muon_gtid.isin(stopping_muons_mc.gtid.values) & (michel_mc.muon_nhit < 2500)]
+
+    print("Particle ID probability for Michel electrons:")
+    print("Data")
+    print_particle_probs(michel_low_nhit)
+    print("Monte Carlo")
+    print_particle_probs(michel_low_nhit_mc)
+
+    fig = plt.figure()
+    plot_hist2_data_mc(michel_low_nhit,michel_low_nhit_mc)
+    despine(fig,trim=True)
+    fig.legend(handles,labels,loc='upper right')
+    if args.save:
+        plt.savefig("michel_electrons.pdf")
+        plt.savefig("michel_electrons.eps")
+    else:
+        plt.suptitle("Michel Electrons")
 
     fig = plt.figure()
     bins = np.linspace(0,100,41)
     plt.hist(michel_low_nhit.ke.values, bins=bins, histtype='step', color='C0', label="Data")
-    plt.hist(michel_mc.ke.values, weights=np.tile(len(michel_low_nhit)/len(michel_mc),len(michel_mc.ke.values)), bins=bins, histtype='step', color='C1', label="Monte Carlo")
+    plt.hist(michel_low_nhit_mc.ke.values, weights=np.tile(len(michel_low_nhit)/len(michel_low_nhit_mc),len(michel_low_nhit_mc.ke.values)), bins=bins, histtype='step', color='C1', label="Monte Carlo")
     hist = np.histogram(michel_low_nhit.ke.values,bins=bins)[0]
-    hist_mc = np.histogram(michel_mc.ke.values,bins=bins)[0]
+    hist_mc = np.histogram(michel_low_nhit_mc.ke.values,bins=bins)[0]
     if hist_mc.sum() > 0:
         norm = hist.sum()/hist_mc.sum()
     else: