From 0ef75c01fdbc60c10b4cf88545f45558b3b9ba3f Mon Sep 17 00:00:00 2001
From: Dmitry Kalinkin <dmitry.kalinkin@gmail.com>
Date: Sun, 2 Mar 2025 20:40:23 -0500
Subject: [PATCH] zdc_{lambda,photon}: disable parts that don't work for
EICrecon 1.22 (#140)
---
.../zdc_lambda/analysis/lambda_plots.py | 145 +++++++++---------
.../zdc_photon/analysis/zdc_photon_plots.py | 6 +
2 files changed, 80 insertions(+), 71 deletions(-)
diff --git a/benchmarks/zdc_lambda/analysis/lambda_plots.py b/benchmarks/zdc_lambda/analysis/lambda_plots.py
index 2ac8d33..a083cc8 100644
--- a/benchmarks/zdc_lambda/analysis/lambda_plots.py
+++ b/benchmarks/zdc_lambda/analysis/lambda_plots.py
@@ -58,82 +58,85 @@ for p in momenta:
pt_truth[p]=np.hypot(px*np.cos(tilt)-pz*np.sin(tilt), py)
theta_truth[p]=np.arctan2(pt_truth[p],pz*np.cos(tilt)+px*np.sin(tilt))
-theta_recon={}
-E_recon={}
-zvtx_recon={}
-mass_recon={}
-print(arrays_sim[p].fields)
-for p in momenta:
-
- px,py,pz,m=(arrays_sim[p][f"ReconstructedFarForwardZDCLambdas.{a}"] for a in "momentum.x momentum.y momentum.z mass".split())
- theta_recon[p]=np.arctan2(np.hypot(px*np.cos(tilt)-pz*np.sin(tilt), py),pz*np.cos(tilt)+px*np.sin(tilt))
- E_recon[p]=np.sqrt(px**2+py**2+pz**2+m**2)
- zvtx_recon[p]=arrays_sim[p][f"ReconstructedFarForwardZDCLambdas.referencePoint.z"]*np.cos(tilt)+arrays_sim[p][f"ReconstructedFarForwardZDCLambdas.referencePoint.x"]*np.sin(tilt)
- mass_recon[p]=m
-
-#theta plots
-fig,axs=plt.subplots(1,3, figsize=(24, 8))
-plt.sca(axs[0])
-plt.title(f"$E_{{\\Lambda}}=100-275$ GeV")
-x=[]
-y=[]
-import awkward as ak
-for p in momenta:
- x+=list(ak.flatten(theta_truth[p]+0*theta_recon[p])*1000)
- y+=list(ak.flatten(theta_recon[p]*1000))
-plt.scatter(x,y)
-plt.xlabel("$\\theta^{*\\rm truth}_{\\Lambda}$ [mrad]")
-plt.ylabel("$\\theta^{*\\rm recon}_{\\Lambda}$ [mrad]")
-plt.xlim(0,3.2)
-plt.ylim(0,3.2)
-
-plt.sca(axs[1])
-plt.title(f"$E_{{\\Lambda}}=100-275$ GeV")
-y,x,_=plt.hist(y-np.array(x), bins=50, range=(-1,1))
-bc=(x[1:]+x[:-1])/2
-
-from scipy.optimize import curve_fit
-slc=abs(bc)<0.3
-fnc=gauss
-p0=[100, 0, 0.05]
-coeff, var_matrix = curve_fit(fnc, bc[slc], y[slc], p0=p0,
- sigma=np.sqrt(y[slc])+(y[slc]==0), maxfev=10000)
-x=np.linspace(-1, 1)
-plt.plot(x, gauss(x, *coeff), color='tab:orange')
-plt.xlabel("$\\theta^{*\\rm recon}_{\\Lambda}-\\theta^{*\\rm truth}_{\\Lambda}$ [mrad]")
-plt.ylabel("events")
-
-plt.sca(axs[2])
-sigmas=[]
-dsigmas=[]
-xvals=[]
-for p in momenta:
- y,x=np.histogram(ak.flatten(theta_recon[p]-theta_truth[p])*1000, bins=100, range=(-1,1))
+if "ReconstructedFarForwardZDCLambdas.momentum.x" not in arrays_sim[momenta[0]].fields:
+ print("ReconstructedFarForwardZDCLambdas collection is not available (needs EICrecon 1.23)")
+else:
+ theta_recon={}
+ E_recon={}
+ zvtx_recon={}
+ mass_recon={}
+ print(arrays_sim[p].fields)
+ for p in momenta:
+
+ px,py,pz,m=(arrays_sim[p][f"ReconstructedFarForwardZDCLambdas.{a}"] for a in "momentum.x momentum.y momentum.z mass".split())
+ theta_recon[p]=np.arctan2(np.hypot(px*np.cos(tilt)-pz*np.sin(tilt), py),pz*np.cos(tilt)+px*np.sin(tilt))
+ E_recon[p]=np.sqrt(px**2+py**2+pz**2+m**2)
+ zvtx_recon[p]=arrays_sim[p][f"ReconstructedFarForwardZDCLambdas.referencePoint.z"]*np.cos(tilt)+arrays_sim[p][f"ReconstructedFarForwardZDCLambdas.referencePoint.x"]*np.sin(tilt)
+ mass_recon[p]=m
+
+ #theta plots
+ fig,axs=plt.subplots(1,3, figsize=(24, 8))
+ plt.sca(axs[0])
+ plt.title(f"$E_{{\\Lambda}}=100-275$ GeV")
+ x=[]
+ y=[]
+ import awkward as ak
+ for p in momenta:
+ x+=list(ak.flatten(theta_truth[p]+0*theta_recon[p])*1000)
+ y+=list(ak.flatten(theta_recon[p]*1000))
+ plt.scatter(x,y)
+ plt.xlabel("$\\theta^{*\\rm truth}_{\\Lambda}$ [mrad]")
+ plt.ylabel("$\\theta^{*\\rm recon}_{\\Lambda}$ [mrad]")
+ plt.xlim(0,3.2)
+ plt.ylim(0,3.2)
+
+ plt.sca(axs[1])
+ plt.title(f"$E_{{\\Lambda}}=100-275$ GeV")
+ y,x,_=plt.hist(y-np.array(x), bins=50, range=(-1,1))
bc=(x[1:]+x[:-1])/2
from scipy.optimize import curve_fit
slc=abs(bc)<0.3
fnc=gauss
- p0=(100, 0, 0.06)
- sigma=np.sqrt(y[slc])+(y[slc]==0)
- try:
- coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000)
- sigmas.append(coeff[2])
- dsigmas.append(np.sqrt(var_matrix[2][2]))
- xvals.append(p)
- except:
- print("fit failed")
-plt.ylim(0, 0.3)
-
-plt.errorbar(xvals, sigmas, dsigmas, ls='', marker='o', color='k')
-x=np.linspace(100, 275, 100)
-plt.plot(x, 3/np.sqrt(x), color='tab:orange')
-plt.text(170, .23, "YR requirement:\n 3 mrad/$\\sqrt{E}$")
-plt.xlabel("$E_{\\Lambda}$ [GeV]")
-plt.ylabel("$\\sigma[\\theta^*_{\\Lambda}]$ [mrad]")
-plt.tight_layout()
-plt.savefig(outdir+"thetastar_recon.pdf")
-#plt.show()
+ p0=[100, 0, 0.05]
+ coeff, var_matrix = curve_fit(fnc, bc[slc], y[slc], p0=p0,
+ sigma=np.sqrt(y[slc])+(y[slc]==0), maxfev=10000)
+ x=np.linspace(-1, 1)
+ plt.plot(x, gauss(x, *coeff), color='tab:orange')
+ plt.xlabel("$\\theta^{*\\rm recon}_{\\Lambda}-\\theta^{*\\rm truth}_{\\Lambda}$ [mrad]")
+ plt.ylabel("events")
+
+ plt.sca(axs[2])
+ sigmas=[]
+ dsigmas=[]
+ xvals=[]
+ for p in momenta:
+ y,x=np.histogram(ak.flatten(theta_recon[p]-theta_truth[p])*1000, bins=100, range=(-1,1))
+ bc=(x[1:]+x[:-1])/2
+
+ from scipy.optimize import curve_fit
+ slc=abs(bc)<0.3
+ fnc=gauss
+ p0=(100, 0, 0.06)
+ sigma=np.sqrt(y[slc])+(y[slc]==0)
+ try:
+ coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000)
+ sigmas.append(coeff[2])
+ dsigmas.append(np.sqrt(var_matrix[2][2]))
+ xvals.append(p)
+ except:
+ print("fit failed")
+ plt.ylim(0, 0.3)
+
+ plt.errorbar(xvals, sigmas, dsigmas, ls='', marker='o', color='k')
+ x=np.linspace(100, 275, 100)
+ plt.plot(x, 3/np.sqrt(x), color='tab:orange')
+ plt.text(170, .23, "YR requirement:\n 3 mrad/$\\sqrt{E}$")
+ plt.xlabel("$E_{\\Lambda}$ [GeV]")
+ plt.ylabel("$\\sigma[\\theta^*_{\\Lambda}]$ [mrad]")
+ plt.tight_layout()
+ plt.savefig(outdir+"thetastar_recon.pdf")
+ #plt.show()
#vtx z
fig,axs=plt.subplots(1,3, figsize=(24, 8))
diff --git a/benchmarks/zdc_photon/analysis/zdc_photon_plots.py b/benchmarks/zdc_photon/analysis/zdc_photon_plots.py
index 009e21c..bd6abf3 100644
--- a/benchmarks/zdc_photon/analysis/zdc_photon_plots.py
+++ b/benchmarks/zdc_photon/analysis/zdc_photon_plots.py
@@ -27,6 +27,12 @@ for p in momenta:
f'sim_output/zdc_photon/{config}_rec_zdc_photon_{p}GeV_{index}.edm4eic.root': 'events'
for index in range(5)
})
+
+if "ReconstructedFarForwardZDCNeutrals.PDG" in arrays_sim[p][momenta[0]].fields:
+ print("ReconstructedFarForwardZDCNeutrals collection is not available (needs EICrecon 1.23)")
+ import sys
+ sys.exit(0)
+
import awkward as ak
fig,axs=plt.subplots(1,3, figsize=(24, 8))
pvals=[]
--
GitLab