diff options
Diffstat (limited to 'utils/plot-muons')
| -rwxr-xr-x | utils/plot-muons | 258 |
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: |