diff --git a/benchmarks/zdc_lambda/analysis/lambda_plots.py b/benchmarks/zdc_lambda/analysis/lambda_plots.py
index f81de0b7c321c8b00f0fc86aae1e4c4affb80c4a..c81159510a70e992c7976758c28597d7516586a3 100644
--- a/benchmarks/zdc_lambda/analysis/lambda_plots.py
+++ b/benchmarks/zdc_lambda/analysis/lambda_plots.py
@@ -60,216 +60,218 @@ for p in momenta:
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))
+ import sys
+ sys.exit(0)
+
+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))
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))
- 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))
- plt.sca(axs[0])
- plt.title(f"$E_{{\\Lambda}}=100-275$ GeV")
- x=[]
- y=[]
- for p in momenta:
- x+=list(ak.flatten(arrays_sim[p]['MCParticles.vertex.z'][:,3]+0*zvtx_recon[p])/1000)
- y+=list(ak.flatten(zvtx_recon[p])/1000)
- plt.scatter(x,y)
- #print(x,y)
- plt.xlabel("$z^{\\rm truth}_{\\rm vtx}$ [m]")
- plt.ylabel("$z^{\\rm recon}_{\\rm vtx}$ [m]")
- plt.xlim(0,40)
- plt.ylim(0,40)
-
- plt.sca(axs[1])
- plt.title(f"$E_{{\\Lambda}}=100-275$ GeV")
- y,x,_=plt.hist(y-np.array(x), bins=50, range=(-10,10))
+ 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))
+plt.sca(axs[0])
+plt.title(f"$E_{{\\Lambda}}=100-275$ GeV")
+x=[]
+y=[]
+for p in momenta:
+ x+=list(ak.flatten(arrays_sim[p]['MCParticles.vertex.z'][:,3]+0*zvtx_recon[p])/1000)
+ y+=list(ak.flatten(zvtx_recon[p])/1000)
+plt.scatter(x,y)
+#print(x,y)
+plt.xlabel("$z^{\\rm truth}_{\\rm vtx}$ [m]")
+plt.ylabel("$z^{\\rm recon}_{\\rm vtx}$ [m]")
+plt.xlim(0,40)
+plt.ylim(0,40)
+
+plt.sca(axs[1])
+plt.title(f"$E_{{\\Lambda}}=100-275$ GeV")
+y,x,_=plt.hist(y-np.array(x), bins=50, range=(-10,10))
+bc=(x[1:]+x[:-1])/2
+
+from scipy.optimize import curve_fit
+slc=abs(bc)<5
+fnc=gauss
+p0=[100, 0, 1]
+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(-5, 5)
+plt.plot(x, gauss(x, *coeff), color='tab:orange')
+print(coeff[2], np.sqrt(var_matrix[2][2]))
+plt.xlabel("$z^{*\\rm recon}_{\\rm vtx}-z^{*\\rm truth}_{\\rm vtx}$ [m]")
+plt.ylabel("events")
+
+plt.sca(axs[2])
+sigmas=[]
+dsigmas=[]
+xvals=[]
+for p in momenta:
+
+
+ a=ak.flatten((zvtx_recon[p]-arrays_sim[p]['MCParticles.vertex.z'][:,3])/1000)
+ y,x=np.histogram(a, bins=100, range=(-10,10))
bc=(x[1:]+x[:-1])/2
from scipy.optimize import curve_fit
slc=abs(bc)<5
fnc=gauss
- p0=[100, 0, 1]
- 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(-5, 5)
- plt.plot(x, gauss(x, *coeff), color='tab:orange')
- print(coeff[2], np.sqrt(var_matrix[2][2]))
- plt.xlabel("$z^{*\\rm recon}_{\\rm vtx}-z^{*\\rm truth}_{\\rm vtx}$ [m]")
- plt.ylabel("events")
-
- plt.sca(axs[2])
- sigmas=[]
- dsigmas=[]
- xvals=[]
- for p in momenta:
-
-
- a=ak.flatten((zvtx_recon[p]-arrays_sim[p]['MCParticles.vertex.z'][:,3])/1000)
- y,x=np.histogram(a, bins=100, range=(-10,10))
- bc=(x[1:]+x[:-1])/2
-
- from scipy.optimize import curve_fit
- slc=abs(bc)<5
- fnc=gauss
- p0=(100, 0, 1)
- #print(bc[slc],y[slc])
- 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, 2)
-
- plt.errorbar(xvals, sigmas, dsigmas, ls='', marker='o', color='k')
- x=np.linspace(100, 275, 100)
-
- avg=np.sum(sigmas/np.array(dsigmas)**2)/np.sum(1/np.array(dsigmas)**2)
- plt.axhline(avg, color='tab:orange')
- plt.text(150, 1.25,f"$\\sigma\\approx${avg:.1f} m")
-
- plt.xlabel("$E_{\\Lambda}$ [GeV]")
- plt.ylabel("$\\sigma[z_{\\rm vtx}]$ [m]")
- plt.tight_layout()
- plt.savefig(outdir+"zvtx_recon.pdf")
- #plt.show()
-
- p=100
- fig,axs=plt.subplots(1,2, figsize=(16, 8))
- plt.sca(axs[0])
- lambda_mass=1.115683
- vals=[]
- for p in momenta:
- vals+=list(ak.flatten(mass_recon[p]))
-
- y,x,_= plt.hist(vals, bins=100, range=(1.0, 1.25))
+ p0=(100, 0, 1)
+ #print(bc[slc],y[slc])
+ 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, 2)
+
+plt.errorbar(xvals, sigmas, dsigmas, ls='', marker='o', color='k')
+x=np.linspace(100, 275, 100)
+
+avg=np.sum(sigmas/np.array(dsigmas)**2)/np.sum(1/np.array(dsigmas)**2)
+plt.axhline(avg, color='tab:orange')
+plt.text(150, 1.25,f"$\\sigma\\approx${avg:.1f} m")
+
+plt.xlabel("$E_{\\Lambda}$ [GeV]")
+plt.ylabel("$\\sigma[z_{\\rm vtx}]$ [m]")
+plt.tight_layout()
+plt.savefig(outdir+"zvtx_recon.pdf")
+#plt.show()
+
+p=100
+fig,axs=plt.subplots(1,2, figsize=(16, 8))
+plt.sca(axs[0])
+lambda_mass=1.115683
+vals=[]
+for p in momenta:
+ vals+=list(ak.flatten(mass_recon[p]))
+
+y,x,_= plt.hist(vals, bins=100, range=(1.0, 1.25))
+bc=(x[1:]+x[:-1])/2
+plt.axvline(lambda_mass, ls='--', color='tab:green', lw=3)
+plt.text(lambda_mass+.01, np.max(y)*1.05, "PDG mass", color='tab:green')
+plt.xlabel("$m_{\\Lambda}^{\\rm recon}$ [GeV]")
+plt.ylim(0, np.max(y)*1.2)
+plt.xlim(1.0, 1.25)
+
+from scipy.optimize import curve_fit
+slc=abs(bc-lambda_mass)<0.05
+fnc=gauss
+p0=[100, lambda_mass, 0.04]
+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(0.8, 1.3, 200)
+plt.plot(x, gauss(x, *coeff), color='tab:orange')
+print(coeff[2], np.sqrt(var_matrix[2][2]))
+plt.xlabel("$m^{\\rm recon}_{\\Lambda}$ [GeV]")
+plt.ylabel("events")
+plt.title(f"$E_{{\\Lambda}}=100-275$ GeV")
+
+plt.sca(axs[1])
+xvals=[]
+sigmas=[]
+dsigmas=[]
+for p in momenta:
+ y,x= np.histogram(ak.flatten(mass_recon[p]), bins=100, range=(0.6,1.4))
bc=(x[1:]+x[:-1])/2
- plt.axvline(lambda_mass, ls='--', color='tab:green', lw=3)
- plt.text(lambda_mass+.01, np.max(y)*1.05, "PDG mass", color='tab:green')
- plt.xlabel("$m_{\\Lambda}^{\\rm recon}$ [GeV]")
- plt.ylim(0, np.max(y)*1.2)
- plt.xlim(1.0, 1.25)
from scipy.optimize import curve_fit
slc=abs(bc-lambda_mass)<0.05
fnc=gauss
- p0=[100, lambda_mass, 0.04]
- 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(0.8, 1.3, 200)
- plt.plot(x, gauss(x, *coeff), color='tab:orange')
- print(coeff[2], np.sqrt(var_matrix[2][2]))
- plt.xlabel("$m^{\\rm recon}_{\\Lambda}$ [GeV]")
- plt.ylabel("events")
- plt.title(f"$E_{{\\Lambda}}=100-275$ GeV")
-
- plt.sca(axs[1])
- xvals=[]
- sigmas=[]
- dsigmas=[]
- for p in momenta:
- y,x= np.histogram(ak.flatten(mass_recon[p]), bins=100, range=(0.6,1.4))
- bc=(x[1:]+x[:-1])/2
-
- from scipy.optimize import curve_fit
- slc=abs(bc-lambda_mass)<0.05
- fnc=gauss
- p0=[100, lambda_mass, 0.05]
- try:
- coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,
- sigma=list(np.sqrt(y[slc])+(y[slc]==0)), maxfev=10000)
- x=np.linspace(0.8, 1.3, 200)
- sigmas.append(coeff[2])
- dsigmas.append(np.sqrt(var_matrix[2][2]))
- xvals.append(p)
- except:
- print("fit failed")
-
- plt.errorbar(xvals, sigmas, dsigmas, ls='', marker='o', color='k')
- avg=np.sum(sigmas/np.array(dsigmas)**2)/np.sum(1/np.array(dsigmas)**2)
- plt.axhline(avg, color='tab:orange')
- plt.text(150, 0.01,f"$\\sigma\\approx${avg*1000:.0f} MeV")
- plt.xlabel("$E_{\\Lambda}$ [GeV]")
- plt.ylabel("$\\sigma[m_{\\Lambda}]$ [GeV]")
- plt.ylim(0, 0.02)
- plt.tight_layout()
- plt.savefig(outdir+"lambda_mass_rec.pdf")
+ p0=[100, lambda_mass, 0.05]
+ try:
+ coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,
+ sigma=list(np.sqrt(y[slc])+(y[slc]==0)), maxfev=10000)
+ x=np.linspace(0.8, 1.3, 200)
+ sigmas.append(coeff[2])
+ dsigmas.append(np.sqrt(var_matrix[2][2]))
+ xvals.append(p)
+ except:
+ print("fit failed")
+
+plt.errorbar(xvals, sigmas, dsigmas, ls='', marker='o', color='k')
+avg=np.sum(sigmas/np.array(dsigmas)**2)/np.sum(1/np.array(dsigmas)**2)
+plt.axhline(avg, color='tab:orange')
+plt.text(150, 0.01,f"$\\sigma\\approx${avg*1000:.0f} MeV")
+plt.xlabel("$E_{\\Lambda}$ [GeV]")
+plt.ylabel("$\\sigma[m_{\\Lambda}]$ [GeV]")
+plt.ylim(0, 0.02)
+plt.tight_layout()
+plt.savefig(outdir+"lambda_mass_rec.pdf")
#now for the CM stuff: