From 266b8a5be4cf1da9f5df4ffd40fb00123149b959 Mon Sep 17 00:00:00 2001 From: Sebouh Paul Date: Mon, 26 Aug 2024 15:20:43 -0400 Subject: [PATCH 1/4] Update neutron_plots.py import uproot as ur --- benchmarks/neutron/analysis/neutron_plots.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/neutron/analysis/neutron_plots.py b/benchmarks/neutron/analysis/neutron_plots.py index a7b919b6..40316cc5 100644 --- a/benchmarks/neutron/analysis/neutron_plots.py +++ b/benchmarks/neutron/analysis/neutron_plots.py @@ -1,4 +1,4 @@ -import numpy as np, pandas as pd, matplotlib.pyplot as plt, matplotlib as mpl, awkward as ak, sys +import numpy as np, pandas as pd, matplotlib.pyplot as plt, matplotlib as mpl, awkward as ak, sys, uproot as ur import mplhep as hep hep.style.use("CMS") From 1fd87ea74c43a72e0819874960016ec3cb238508 Mon Sep 17 00:00:00 2001 From: Sebouh Paul Date: Mon, 26 Aug 2024 15:28:31 -0400 Subject: [PATCH 2/4] fixed another bug in the neutron_plots.py --- benchmarks/neutron/analysis/neutron_plots.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/neutron/analysis/neutron_plots.py b/benchmarks/neutron/analysis/neutron_plots.py index 40316cc5..41583b9c 100644 --- a/benchmarks/neutron/analysis/neutron_plots.py +++ b/benchmarks/neutron/analysis/neutron_plots.py @@ -57,7 +57,7 @@ def gauss(x, A,mu, sigma): w=E array['theta_recon']=np.sum(np.arccos(zp/r)*w, axis=-1)/np.sum(w, axis=-1) - array['eta_recon']=-np.log(np.tan(['theta_recon']/2)) + array['eta_recon']=-np.log(np.tan(array['theta_recon']/2)) array['E_Hcal']=np.sum(array['HcalEndcapPInsertClusters.energy'], axis=-1)#*20/12.5 From ad66e982057f0573c9f031bf701c2fb03e975cb4 Mon Sep 17 00:00:00 2001 From: Sebouh Paul Date: Tue, 27 Aug 2024 16:07:32 -0400 Subject: [PATCH 3/4] different fit for the neutron energy corrections; also add YR requirements --- benchmarks/neutron/analysis/neutron_plots.py | 38 +++++++++++++++----- 1 file changed, 29 insertions(+), 9 deletions(-) diff --git a/benchmarks/neutron/analysis/neutron_plots.py b/benchmarks/neutron/analysis/neutron_plots.py index 41583b9c..6eedf615 100644 --- a/benchmarks/neutron/analysis/neutron_plots.py +++ b/benchmarks/neutron/analysis/neutron_plots.py @@ -62,11 +62,11 @@ def gauss(x, A,mu, sigma): array['E_Hcal']=np.sum(array['HcalEndcapPInsertClusters.energy'], axis=-1)#*20/12.5 array['E_Ecal']=np.sum(array['EcalEndcapPInsertClusters.energy'], axis=-1)#*27/20 - - array['E_recon']=array['E_Hcal']/(0.70-.30/np.sqrt(array['E_Hcal']))\ - +array['E_Ecal']/(0.82-0.35/np.sqrt(array['E_Ecal'])) + coeffs=[1.0540361, 1.12617385, 2.87336987, 1.9086172 ] + array['E_recon']=coeffs[0]*array['E_Hcal']+coeffs[1]*array['E_Ecal']+coeffs[2]*array['E_Hcal']/np.sqrt(array['E_Hcal']+array['E_Ecal'])+coeffs[3]*array['E_Ecal']/np.sqrt(array['E_Hcal']+array['E_Ecal']) #plot theta residuals: +print("making theta recon plot") from scipy.optimize import curve_fit fig, axs=plt.subplots(1,2, figsize=(16,8)) @@ -115,6 +115,8 @@ def gauss(x, A,mu, sigma): pass plt.sca(axs[1]) plt.errorbar(xvals, sigmas, dsigmas, ls='', marker='o', label=f"${eta_min}<\\eta<{eta_max}$") +xx=np.linspace(15,85, 100) +plt.plot(xx, 3/np.sqrt(xx), ls='--', label='YR req.: (3 mrad)/$\\sqrt{E}$') plt.xlabel("$p_{n}$ [GeV]") plt.ylabel("$\\sigma[\\theta]$ [mrad]") plt.ylim(0) @@ -122,7 +124,9 @@ def gauss(x, A,mu, sigma): plt.tight_layout() plt.savefig(outdir+"neutron_theta_recon.pdf") -fig, axs=plt.subplots(1,2, figsize=(16,8)) +#now make the energy plot +print("making energy recon plot") +fig, axs=plt.subplots(1,3, figsize=(24,8)) plt.sca(axs[0]) p=50 eta_min=3.4; eta_max=3.6 @@ -130,11 +134,12 @@ def gauss(x, A,mu, sigma): [(arrays_sim[p]['eta_truth']>eta_min)&(arrays_sim[p]['eta_truth']0)], bins=50, range=(-.5,.5), histtype='step') bc=(x[1:]+x[:-1])/2 -slc=abs(bc)<0.4 +slc=abs(bc)<.2 # try: -p0=(100, 0, 0.3) +p0=(100, 0, 0.15) fnc=gauss sigma=np.sqrt(y[slc])+(y[slc]==0) + coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,sigma=list(sigma)) xx=np.linspace(-.5,.5,100) plt.plot(xx,fnc(xx,*coeff)) @@ -149,27 +154,42 @@ def gauss(x, A,mu, sigma): xvals=[] sigmas=[] dsigmas=[] + mus=[] + dmus=[] for p in 20,30,40, 50, 60, 70, 80: y,x=np.histogram(((arrays_sim[p]['E_recon']-arrays_sim[p]['E_truth'])/arrays_sim[p]['E_truth'])\ [(arrays_sim[p]['eta_truth']>eta_min)&(arrays_sim[p]['eta_truth'] Date: Tue, 3 Sep 2024 16:10:32 -0400 Subject: [PATCH 4/4] removed YR requirements from neutron plots; also used a different energy correction formula --- benchmarks/neutron/analysis/neutron_plots.py | 222 +++++++++++++------ 1 file changed, 154 insertions(+), 68 deletions(-) diff --git a/benchmarks/neutron/analysis/neutron_plots.py b/benchmarks/neutron/analysis/neutron_plots.py index 6eedf615..fc1d5b30 100644 --- a/benchmarks/neutron/analysis/neutron_plots.py +++ b/benchmarks/neutron/analysis/neutron_plots.py @@ -21,7 +21,6 @@ arrays_sim[p] = ur.open(f'sim_output/neutron/{config}_rec_neutron_{p}GeV.edm4hep.root:events')\ .arrays() -#get reconstructed theta, eta, and E def gauss(x, A,mu, sigma): return A * np.exp(-(x-mu)**2/(2*sigma**2)) @@ -41,7 +40,8 @@ def gauss(x, A,mu, sigma): array['eta_truth']=1/2*np.log((p+pzp)/(p-pzp)) array['theta_truth']=np.arccos(pzp/p) -#weighted avg of theta of cluster centers, by energy +# +# get reconstructed theta as avg of theta of cluster centers, weighted by energy for array in arrays_sim.values(): tilt=-0.025 x=array['HcalEndcapPInsertClusters.position.x'] @@ -59,11 +59,6 @@ def gauss(x, A,mu, sigma): array['theta_recon']=np.sum(np.arccos(zp/r)*w, axis=-1)/np.sum(w, axis=-1) array['eta_recon']=-np.log(np.tan(array['theta_recon']/2)) - - array['E_Hcal']=np.sum(array['HcalEndcapPInsertClusters.energy'], axis=-1)#*20/12.5 - array['E_Ecal']=np.sum(array['EcalEndcapPInsertClusters.energy'], axis=-1)#*27/20 - coeffs=[1.0540361, 1.12617385, 2.87336987, 1.9086172 ] - array['E_recon']=coeffs[0]*array['E_Hcal']+coeffs[1]*array['E_Ecal']+coeffs[2]*array['E_Hcal']/np.sqrt(array['E_Hcal']+array['E_Ecal'])+coeffs[3]*array['E_Ecal']/np.sqrt(array['E_Hcal']+array['E_Ecal']) #plot theta residuals: print("making theta recon plot") @@ -91,7 +86,7 @@ def gauss(x, A,mu, sigma): plt.ylabel("events") plt.title(f"$p={p}$ GeV, ${eta_min}<\\eta<{eta_max}$") -r=[3.0,3.2,3.4,3.6,3.8] +r=[3.2,3.4,3.6,3.8,4.0] for eta_min, eta_max in zip(r[:-1],r[1:]): xvals=[] sigmas=[] @@ -115,8 +110,6 @@ def gauss(x, A,mu, sigma): pass plt.sca(axs[1]) plt.errorbar(xvals, sigmas, dsigmas, ls='', marker='o', label=f"${eta_min}<\\eta<{eta_max}$") -xx=np.linspace(15,85, 100) -plt.plot(xx, 3/np.sqrt(xx), ls='--', label='YR req.: (3 mrad)/$\\sqrt{E}$') plt.xlabel("$p_{n}$ [GeV]") plt.ylabel("$\\sigma[\\theta]$ [mrad]") plt.ylim(0) @@ -124,72 +117,165 @@ def gauss(x, A,mu, sigma): plt.tight_layout() plt.savefig(outdir+"neutron_theta_recon.pdf") +#now determine the energy recon parameters +pvals=[] +resvals=[] +reserrs=[] +wvals=[] +svals=[] +Euncorr=[] + +print("determining the energy recon correction parameters") +fig,axs=plt.subplots(1,2, figsize=(20,10)) +eta_min=3.4;eta_max=3.6 +for p in 20, 30,40,50,60,70, 80: + best_res=1000 + res_err=1000 + best_s=1000 + wrange=np.linspace(30, 70, 41)*0.0257 + coeff_best=None + + wbest=0 + a=arrays_sim[p] + h=np.sum(a[f'HcalEndcapPInsertClusters.energy'], axis=-1) + e=np.sum(a[f'EcalEndcapPInsertClusters.energy'], axis=-1) + for w in wrange: + + r=(e/w+h)[(h>0)&(a['eta_truth']>eta_min)&(a['eta_truth']0)&(a['eta_truth']>3.4)&(a['eta_truth']<3.6)] + plt.sca(axs[0]) + y, x, _= plt.hist(r, histtype='step', bins=50) + xx=np.linspace(20, 55, 100) + plt.plot(xx,fnc(xx, *coeff_best), ls='-') + plt.xlabel("$E_{uncorr}=E_{Hcal}+E_{Ecal}/w$ [GeV]") + plt.title(f"p=50 GeV, ${eta_min}<\\eta<{eta_max}$, w={wbest:.2f}") + plt.axvline(np.sqrt(50**2+.9406**2), color='g', ls=':') + plt.text(40, max(y)*0.9, "generated\nenergy", color='g', fontsize=20) + + Euncorr.append(Euncorr_best) + resvals.append(best_res) + reserrs.append(res_err) + pvals.append(p) + svals.append(best_s) + wvals.append(wbest) + +pvals=np.array(pvals) +svals=np.array(svals) +Euncorr=np.array(Euncorr) +plt.sca(axs[1]) +plt.plot(Euncorr,wvals, label="w") +w_avg=np.mean(wvals) +plt.axhline(w_avg, label=f'w avg: {w_avg:.2f}', ls=':') +plt.plot(Euncorr,svals, label="s") +m=(np.sum(svals*Euncorr)*len(Euncorr)-np.sum(Euncorr)*np.sum(svals))/(np.sum(Euncorr**2)*len(Euncorr)-np.sum(Euncorr)**2) +b=np.mean(svals)-np.mean(Euncorr)*m +plt.plot(Euncorr,Euncorr*m+b, label=f"s fit: ${m:.4f}E_{{uncorr}}+{b:.2f}$", ls=':') + +plt.xlabel("$E_{uncorr}=E_{Hcal}+E_{Ecal}/w$ [GeV]") +plt.title("$E_{n,recon}=s\\times(E_{Hcal}+E_{Ecal}/w)$") +plt.ylabel('parameter values') +plt.legend() +plt.ylim(0) +plt.savefig(outdir+"neutron_energy_params.pdf") + #now make the energy plot print("making energy recon plot") fig, axs=plt.subplots(1,3, figsize=(24,8)) -plt.sca(axs[0]) -p=50 -eta_min=3.4; eta_max=3.6 -y,x,_=plt.hist(((arrays_sim[p]['E_recon']-arrays_sim[p]['E_truth'])/arrays_sim[p]['E_truth'])\ - [(arrays_sim[p]['eta_truth']>eta_min)&(arrays_sim[p]['eta_truth']0)], bins=50, - range=(-.5,.5), histtype='step') -bc=(x[1:]+x[:-1])/2 -slc=abs(bc)<.2 -# try: -p0=(100, 0, 0.15) -fnc=gauss -sigma=np.sqrt(y[slc])+(y[slc]==0) +partitions=[3.2,3.4, 3.6, 3.8, 4.0] +for eta_min, eta_max in zip(partitions[:-1],partitions[1:]): + pvals=[] + resvals=[] + reserrs=[] + scalevals=[] + scaleerrs=[] + for p in 20, 30,40,50,60,70, 80: + best_res=1000 + res_err=1000 -coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,sigma=list(sigma)) -xx=np.linspace(-.5,.5,100) -plt.plot(xx,fnc(xx,*coeff)) -# except: -# pass -plt.xlabel("$(E_{rec}-E_{truth})/E_{truth}$") -plt.ylabel("events") -plt.title(f"$p={p}$ GeV, ${eta_min}<\\eta<{eta_max}$") -r=[3.0,3.2,3.4,3.6,3.8] -for eta_min, eta_max in zip(r[:-1],r[1:]): - xvals=[] - sigmas=[] - dsigmas=[] - mus=[] - dmus=[] - for p in 20,30,40, 50, 60, 70, 80: - y,x=np.histogram(((arrays_sim[p]['E_recon']-arrays_sim[p]['E_truth'])/arrays_sim[p]['E_truth'])\ - [(arrays_sim[p]['eta_truth']>eta_min)&(arrays_sim[p]['eta_truth']0)&(a['eta_truth']>eta_min)&(a['eta_truth']np.pi/2)] + y,x=np.histogram(r,bins=50) + bcs=(x[1:]+x[:-1])/2 fnc=gauss - p0=(100, 0, 0.15) - #print(bc[slc],y[slc]) - sigma=np.sqrt(y[slc])+(y[slc]==0) + slc=abs(bcs-np.mean(r)*1.25)<2*np.std(r) + sigma=np.sqrt(y[slc])+0.5*(y[slc]==0) + p0=(100, np.mean(r), np.std(r)) try: - coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,sigma=list(sigma)) - sigmas.append(np.abs(coeff[2])) - dsigmas.append(np.sqrt(var_matrix[2][2])) - mus.append(coeff[1]) - dmus.append(np.sqrt(var_matrix[2][2])) - xvals.append(p) - except: - pass + coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma)) + res=np.abs(coeff[2]/coeff[1]) + + if res