diff options
Diffstat (limited to 'utils/plot-muons')
| -rwxr-xr-x | utils/plot-muons | 29 |
1 files changed, 27 insertions, 2 deletions
diff --git a/utils/plot-muons b/utils/plot-muons index 69dda87..e80d0af 100755 --- a/utils/plot-muons +++ b/utils/plot-muons @@ -180,6 +180,9 @@ if __name__ == '__main__': stopping_muons = stopping_muons[stopping_muons.udotr < -0.5] stopping_muons_mc = stopping_muons_mc[stopping_muons_mc.udotr < -0.5] + 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) @@ -190,6 +193,27 @@ if __name__ == '__main__': 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") + plt.hist(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") + plt.hist(np.cos(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:") @@ -252,6 +276,7 @@ if __name__ == '__main__': 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 (\%)") @@ -268,8 +293,8 @@ if __name__ == '__main__': 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") hist = np.histogram(michel_low_nhit.ke.values,bins=bins)[0] hist_mc = np.histogram(michel_mc.ke.values,bins=bins)[0] - if len(michel_mc) > 0: - norm = len(michel)/len(michel_mc) + if hist_mc.sum() > 0: + norm = hist.sum()/hist_mc.sum() else: norm = 1.0 p = get_multinomial_prob(hist,hist_mc,norm) |