From 697b44112b62c65174d0486cd92e74751d75288c Mon Sep 17 00:00:00 2001
From: Sebouh Paul <sebouh.paul@gmail.com>
Date: Fri, 23 Aug 2024 12:29:42 -0400
Subject: [PATCH 1/3] 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 <sebouh.paul@gmail.com>
Date: Tue, 10 Sep 2024 21:33:40 -0400
Subject: [PATCH 2/3] 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 <sebouh.paul@gmail.com>
Date: Tue, 10 Sep 2024 21:58:34 -0400
Subject: [PATCH 3/3] 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")