From 697b44112b62c65174d0486cd92e74751d75288c Mon Sep 17 00:00:00 2001 From: Sebouh Paul Date: Fri, 23 Aug 2024 12:29:42 -0400 Subject: [PATCH 1/4] updated Snakefile --- Snakefile | 1 + 1 file changed, 1 insertion(+) diff --git a/Snakefile b/Snakefile index 1825ce1a..05997c6a 100644 --- a/Snakefile +++ b/Snakefile @@ -43,4 +43,5 @@ ddsim \ include: "benchmarks/diffractive_vm/Snakefile" include: "benchmarks/dis/Snakefile" include: "benchmarks/lambda/Snakefile" +include: "benchmarks/neutron/Snakefile" include: "benchmarks/demp/Snakefile" From b0ac6bcd6b92c4fda468fd4168692c29235a745a Mon Sep 17 00:00:00 2001 From: Sebouh Paul Date: Tue, 10 Sep 2024 21:33:40 -0400 Subject: [PATCH 2/4] added fit to the energy recon resolution for neutrons --- benchmarks/neutron/analysis/neutron_plots.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/benchmarks/neutron/analysis/neutron_plots.py b/benchmarks/neutron/analysis/neutron_plots.py index fc1d5b30..3bb8e2fb 100644 --- a/benchmarks/neutron/analysis/neutron_plots.py +++ b/benchmarks/neutron/analysis/neutron_plots.py @@ -203,6 +203,7 @@ def gauss(x, A,mu, sigma): print("making energy recon plot") fig, axs=plt.subplots(1,3, figsize=(24,8)) partitions=[3.2,3.4, 3.6, 3.8, 4.0] + for eta_min, eta_max in zip(partitions[:-1],partitions[1:]): pvals=[] resvals=[] @@ -213,7 +214,6 @@ def gauss(x, A,mu, sigma): best_res=1000 res_err=1000 - wrange=np.linspace(30, 70, 30)*0.0257 w=w_avg @@ -268,14 +268,21 @@ def gauss(x, A,mu, sigma): plt.ylabel("$\\mu[E]/E$") - - + if eta_min==3.4: + fnc=lambda E, a, b: np.hypot(a,b/np.sqrt(E)) + p0=[.1,.5] + coeff, var_matrix = curve_fit(fnc, pvals, resvals, p0=p0,sigma=reserrs) + xx=np.linspace(15, 85, 100) + axs[1].plot(xx, fnc(xx,*coeff), color='tab:purple',ls='--', + label=f'fit ${eta_min:.1f}<\\eta<{eta_max:.1f}$:\n'+\ + f'{coeff[0]*100:.1f}%$\\oplus\\frac{{{coeff[1]*100:.0f}\\%}}{{\\sqrt{{E}}}}$') + axs[2].set_xlabel("$p_n$ [GeV]") axs[2].axhline(1, ls='--', color='0.5', alpha=0.7) axs[0].set_ylim(0) axs[1].set_ylim(0, 0.35) axs[2].set_ylim(0) -axs[1].legend() -axs[2].legend() +axs[1].legend(fontsize=20) +axs[2].legend(fontsize=20) plt.tight_layout() plt.savefig(outdir+"neutron_energy_recon.pdf") From dfdcecad27784b32510728cb0202dc85ca27cd25 Mon Sep 17 00:00:00 2001 From: Sebouh Paul Date: Tue, 10 Sep 2024 21:58:34 -0400 Subject: [PATCH 3/4] added fit for theta res --- benchmarks/neutron/analysis/neutron_plots.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/benchmarks/neutron/analysis/neutron_plots.py b/benchmarks/neutron/analysis/neutron_plots.py index 3bb8e2fb..5d61d1ce 100644 --- a/benchmarks/neutron/analysis/neutron_plots.py +++ b/benchmarks/neutron/analysis/neutron_plots.py @@ -110,9 +110,17 @@ 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}$") + if eta_min==3.4: + fnc=lambda E, a, b: np.hypot(a,b/np.sqrt(E)) + p0=[.002,.05] + coeff, var_matrix = curve_fit(fnc, xvals, sigmas, p0=p0,sigma=dsigmas) + xx=np.linspace(15, 85, 100) + axs[1].plot(xx, fnc(xx,*coeff), color='tab:purple',ls='--', + label=f'fit ${eta_min:.1f}<\\eta<{eta_max:.1f}$:\n'+\ + f'({coeff[0]:.2f}$\\oplus\\frac{{{coeff[1]:.1f}}}{{\\sqrt{{E}}}}$) mrad') plt.xlabel("$p_{n}$ [GeV]") plt.ylabel("$\\sigma[\\theta]$ [mrad]") -plt.ylim(0) +plt.ylim(0, 10) plt.legend() plt.tight_layout() plt.savefig(outdir+"neutron_theta_recon.pdf") From c9f36f8c04bba2ccf0a5f775abfebce2c52b164a Mon Sep 17 00:00:00 2001 From: Sebouh Paul Date: Wed, 11 Sep 2024 10:03:04 -0400 Subject: [PATCH 4/4] refactoring of the correction factors; also removed constant term from the energy res fit --- benchmarks/neutron/analysis/neutron_plots.py | 28 ++++++++++---------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/benchmarks/neutron/analysis/neutron_plots.py b/benchmarks/neutron/analysis/neutron_plots.py index 5d61d1ce..5bbef717 100644 --- a/benchmarks/neutron/analysis/neutron_plots.py +++ b/benchmarks/neutron/analysis/neutron_plots.py @@ -120,7 +120,7 @@ def gauss(x, A,mu, sigma): f'({coeff[0]:.2f}$\\oplus\\frac{{{coeff[1]:.1f}}}{{\\sqrt{{E}}}}$) mrad') plt.xlabel("$p_{n}$ [GeV]") plt.ylabel("$\\sigma[\\theta]$ [mrad]") -plt.ylim(0, 10) +plt.ylim(0, 5) plt.legend() plt.tight_layout() plt.savefig(outdir+"neutron_theta_recon.pdf") @@ -140,7 +140,7 @@ def gauss(x, A,mu, sigma): best_res=1000 res_err=1000 best_s=1000 - wrange=np.linspace(30, 70, 41)*0.0257 + wrange=np.linspace(0.8, 1.2, 41) coeff_best=None wbest=0 @@ -149,7 +149,7 @@ def gauss(x, A,mu, sigma): 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']>eta_min)&(a['eta_truth']0)&(a['eta_truth']>3.4)&(a['eta_truth']<3.6)] + r=(e+h*wbest)[(h>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.xlabel("$E_{uncorr}=w\\times E_{Hcal}+E_{Ecal}$ [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) @@ -200,8 +200,8 @@ def gauss(x, A,mu, sigma): 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.xlabel("$E_{uncorr}=w\\times E_{Hcal}+E_{Ecal}$ [GeV]") +plt.title("$E_{n,recon}=s\\times(w\\times E_{Hcal}+E_{Ecal})$") plt.ylabel('parameter values') plt.legend() plt.ylim(0) @@ -229,8 +229,8 @@ def gauss(x, A,mu, sigma): h=np.sum(a[f'HcalEndcapPInsertClusters.energy'], axis=-1) e=np.sum(a[f'EcalEndcapPInsertClusters.energy'], axis=-1) #phi=a['phi_truth'] - uncorr=(e/w+h) - s=-0.0064*uncorr+1.80 + uncorr=(e+h*w) + s=-0.0047*uncorr+1.64 r=uncorr*s #reconstructed energy with correction r=r[(h>0)&(a['eta_truth']>eta_min)&(a['eta_truth']np.pi/2)] y,x=np.histogram(r,bins=50) @@ -277,13 +277,13 @@ def gauss(x, A,mu, sigma): plt.ylabel("$\\mu[E]/E$") if eta_min==3.4: - fnc=lambda E, a, b: np.hypot(a,b/np.sqrt(E)) - p0=[.1,.5] - coeff, var_matrix = curve_fit(fnc, pvals, resvals, p0=p0,sigma=reserrs) + fnc=lambda E, b: b/np.sqrt(E) + p0=[.5] + coeff, var_matrix = curve_fit(fnc, pvals, resvals, p0=p0,sigma=np.array(reserrs)) xx=np.linspace(15, 85, 100) axs[1].plot(xx, fnc(xx,*coeff), color='tab:purple',ls='--', - label=f'fit ${eta_min:.1f}<\\eta<{eta_max:.1f}$:\n'+\ - f'{coeff[0]*100:.1f}%$\\oplus\\frac{{{coeff[1]*100:.0f}\\%}}{{\\sqrt{{E}}}}$') + label=f'fit ${eta_min:.1f}<\\eta<{eta_max:.1f}$: '+\ + f'$\\frac{{{coeff[0]*100:.0f}\\%}}{{\\sqrt{{E}}}}$') axs[2].set_xlabel("$p_n$ [GeV]") axs[2].axhline(1, ls='--', color='0.5', alpha=0.7)