diff --git a/benchmarks/zdc_lambda/analysis/lambda_plots.py b/benchmarks/zdc_lambda/analysis/lambda_plots.py index a083cc86a8dd819e3b63802c49b60c82342965db..f81de0b7c321c8b00f0fc86aae1e4c4affb80c4a 100644 --- a/benchmarks/zdc_lambda/analysis/lambda_plots.py +++ b/benchmarks/zdc_lambda/analysis/lambda_plots.py @@ -138,138 +138,138 @@ else: 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: - + #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) - a=ak.flatten((zvtx_recon[p]-arrays_sim[p]['MCParticles.vertex.z'][:,3])/1000) - y,x=np.histogram(a, bins=100, range=(-10,10)) + 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) - #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)) + 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)) 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.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.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") #now for the CM stuff: diff --git a/benchmarks/zdc_photon/analysis/zdc_photon_plots.py b/benchmarks/zdc_photon/analysis/zdc_photon_plots.py index bd6abf3c9f585c344cdbe6d7813ef72ea2d39ee8..a89f79589e71229b12ea7dcdb879f45dd41deedc 100644 --- a/benchmarks/zdc_photon/analysis/zdc_photon_plots.py +++ b/benchmarks/zdc_photon/analysis/zdc_photon_plots.py @@ -28,7 +28,7 @@ for p in momenta: for index in range(5) }) -if "ReconstructedFarForwardZDCNeutrals.PDG" in arrays_sim[p][momenta[0]].fields: +if "ReconstructedFarForwardZDCNeutrals.PDG" not in arrays_sim[p][momenta[0]].fields: print("ReconstructedFarForwardZDCNeutrals collection is not available (needs EICrecon 1.23)") import sys sys.exit(0)