Skip to content
Snippets Groups Projects
Unverified Commit 7ba375c3 authored by Sebouh Paul's avatar Sebouh Paul Committed by GitHub
Browse files

ZDC benchmarks use reconstructed particles from EICrecon, rather than... 

ZDC benchmarks use reconstructed particles from EICrecon, rather than attempting to reconstruct them directly (#138)
parent 1ad5cea7
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
benchmarks/zdc_lambda/Snakefile
+ 1
− 1
View file @ 7ba375c3
...@@ -52,7 +52,7 @@ rule zdc_lambda_recon: ...@@ -52,7 +52,7 @@ rule zdc_lambda_recon:
shell: shell:
""" """
env DETECTOR_CONFIG={params.DETECTOR_CONFIG} \ env DETECTOR_CONFIG={params.DETECTOR_CONFIG} \
eicrecon {input.SIM_FILE} -Ppodio:output_file={output.REC_FILE} -Ppodio:output_collections=MCParticles,HcalFarForwardZDCClusters,HcalFarForwardZDCRecHits,HcalFarForwardZDCSubcellHits -Pjana:nevents={params.N_EVENTS} eicrecon {input.SIM_FILE} -Ppodio:output_file={output.REC_FILE} -Ppodio:output_collections=MCParticles,HcalFarForwardZDCClusters,HcalFarForwardZDCRecHits,HcalFarForwardZDCSubcellHits,ReconstructedFarForwardZDCLambdas,ReconstructedFarForwardZDCLambdaDecayProductsCM -Pjana:nevents={params.N_EVENTS}
""" """
rule zdc_lambda_analysis: rule zdc_lambda_analysis:
......
benchmarks/zdc_lambda/analysis/lambda_plots.py
+ 66
− 125
View file @ 7ba375c3
...@@ -39,7 +39,7 @@ for p in momenta: ...@@ -39,7 +39,7 @@ for p in momenta:
plt.hist(nclusters[p],bins=20, range=(0,20)) plt.hist(nclusters[p],bins=20, range=(0,20))
plt.xlabel("number of clusters") plt.xlabel("number of clusters")
plt.yscale('log') plt.yscale('log')
plt.title(f"$p_\Lambda={p}$ GeV") plt.title(f"$p_\\Lambda={p}$ GeV")
plt.ylim(1) plt.ylim(1)
plt.savefig(outdir+f"nclust_{p}GeV_recon.pdf") plt.savefig(outdir+f"nclust_{p}GeV_recon.pdf")
print("saved file ", outdir+f"nclust_{p}GeV_recon.pdf") print("saved file ", outdir+f"nclust_{p}GeV_recon.pdf")
...@@ -58,123 +58,29 @@ for p in momenta: ...@@ -58,123 +58,29 @@ for p in momenta:
pt_truth[p]=np.hypot(px*np.cos(tilt)-pz*np.sin(tilt), py) 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_truth[p]=np.arctan2(pt_truth[p],pz*np.cos(tilt)+px*np.sin(tilt))
theta_recon={}
#create an array with the same shape as the cluster-level arrays E_recon={}
is_neutron_cand={}
for p in momenta:
is_neutron_cand[p]=(0*arrays_sim[p][f"HcalFarForwardZDCClusters.energy"]).to_list()
#largest_eigenvalues
for i in range(len(arrays_sim[p])):
pars=arrays_sim[p]["_HcalFarForwardZDCClusters_shapeParameters"][i]
index_of_max=-1
max_val=0
eigs=[]
shape_params_per_cluster=7
for j in range(len(pars)//shape_params_per_cluster):
largest_eigenvalue=max(pars[shape_params_per_cluster*j+4:shape_params_per_cluster*j+7])
eigs.append(largest_eigenvalue)
if(largest_eigenvalue>max_val):
max_val=largest_eigenvalue
index_of_max=j
if index_of_max >=0:
is_neutron_cand[p][i][index_of_max]=1
eigs.sort()
is_neutron_cand[p]=ak.Array(is_neutron_cand[p])
#with the position of the vertex determined by assuming the mass of the pi0
#corrected pt* and theta* recon
pt_recon_corr={}
theta_recon_corr={}
mass_recon_corr={}
mass_pi0_recon_corr={}
pi0_converged={}
zvtx_recon={} zvtx_recon={}
mass_recon={}
maxZ=35800 print(arrays_sim[p].fields)
for p in momenta: for p in momenta:
xvtx=0
yvtx=0
zvtx=0
for iteration in range(20):
#compute the value of theta* using the clusters in the ZDC
xc=arrays_sim[p][f"HcalFarForwardZDCClusters.position.x"]
yc=arrays_sim[p][f"HcalFarForwardZDCClusters.position.y"]
zc=arrays_sim[p][f"HcalFarForwardZDCClusters.position.z"]
E=arrays_sim[p][f"HcalFarForwardZDCClusters.energy"]
#apply correction to the neutron candidates only
A,B,C=-0.0756, -1.91, 2.30
neutron_corr=(1-is_neutron_cand[p])+is_neutron_cand[p]/(1+A+B/np.sqrt(E)+C/E)
E=E*neutron_corr
E_recon=np.sum(E, axis=-1)
pabs=np.sqrt(E**2-is_neutron_cand[p]*.9406**2)
tilt=-0.025
xcp=xc*np.cos(tilt)-zc*np.sin(tilt)
ycp=yc
zcp=zc*np.cos(tilt)+xc*np.sin(tilt)
rcp=np.sqrt(xcp**2+ycp**2+zcp**2)
ux=(xcp-xvtx)
uy=(ycp-yvtx)
uz=(zcp-zvtx)
norm=np.sqrt(ux**2+uy**2+uz**2)
ux=ux/norm
uy=uy/norm
uz=uz/norm
px_recon,py_recon,pz_recon=np.sum(pabs*ux, axis=-1),np.sum(pabs*uy, axis=-1),np.sum(pabs*uz, axis=-1)
pt_recon_corr[p]=np.hypot(px_recon,py_recon) px,py,pz,m=(arrays_sim[p][f"ReconstructedFarForwardZDCLambdas.{a}"] for a in "momentum.x momentum.y momentum.z mass".split())
theta_recon_corr[p]=np.arctan2(pt_recon_corr[p], pz_recon) 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)
mass_recon_corr[p]=np.sqrt((E_recon)**2\ zvtx_recon[p]=arrays_sim[p][f"ReconstructedFarForwardZDCLambdas.referencePoint.z"]*np.cos(tilt)+arrays_sim[p][f"ReconstructedFarForwardZDCLambdas.referencePoint.x"]*np.sin(tilt)
-(px_recon)**2\ mass_recon[p]=m
-(py_recon)**2\
-(pz_recon)**2) #theta plots
mass_pi0_recon_corr[p]=np.sqrt(np.sum(pabs*(1-is_neutron_cand[p]), axis=-1)**2\
-np.sum(pabs*ux*(1-is_neutron_cand[p]), axis=-1)**2\
-np.sum(pabs*uy*(1-is_neutron_cand[p]), axis=-1)**2\
-np.sum(pabs*uz*(1-is_neutron_cand[p]), axis=-1)**2)
alpha=1
if iteration==0:
u=np.sqrt(px_recon**2+py_recon**2+pz_recon**2)
ux=px_recon/u
uy=py_recon/u
uz=pz_recon/u
zeta=1/2
zvtx=maxZ*np.power(zeta,alpha)
xvtx=ux/uz*zvtx
yvtx=uy/uz*zvtx
else :
u=np.sqrt(px_recon**2+py_recon**2+pz_recon**2)
ux=px_recon/u
uy=py_recon/u
uz=pz_recon/u
s=2*(mass_pi0_recon_corr[p]<0.135)-1
zeta=np.power(zvtx/maxZ, 1/alpha)
zeta=zeta+s*1/2**(1+iteration)
zvtx=maxZ*np.power(zeta,alpha)
xvtx=ux/uz*zvtx
yvtx=uy/uz*zvtx
#print(zvtx)
pi0_converged[p]=np.abs(mass_pi0_recon_corr[p]-0.135)<0.01
zvtx_recon[p]=zvtx
fig,axs=plt.subplots(1,3, figsize=(24, 8)) fig,axs=plt.subplots(1,3, figsize=(24, 8))
plt.sca(axs[0]) plt.sca(axs[0])
plt.title(f"$E_{{\\Lambda}}=100-275$ GeV") plt.title(f"$E_{{\\Lambda}}=100-275$ GeV")
x=[] x=[]
y=[] y=[]
import awkward as ak
for p in momenta: for p in momenta:
accept=(nclusters[p]==3) &(pi0_converged[p]) x+=list(ak.flatten(theta_truth[p]+0*theta_recon[p])*1000)
x+=list(theta_truth[p][accept]*1000) y+=list(ak.flatten(theta_recon[p]*1000))
y+=list(theta_recon_corr[p][accept]*1000)
plt.scatter(x,y) plt.scatter(x,y)
plt.xlabel("$\\theta^{*\\rm truth}_{\\Lambda}$ [mrad]") plt.xlabel("$\\theta^{*\\rm truth}_{\\Lambda}$ [mrad]")
plt.ylabel("$\\theta^{*\\rm recon}_{\\Lambda}$ [mrad]") plt.ylabel("$\\theta^{*\\rm recon}_{\\Lambda}$ [mrad]")
...@@ -202,16 +108,13 @@ sigmas=[] ...@@ -202,16 +108,13 @@ sigmas=[]
dsigmas=[] dsigmas=[]
xvals=[] xvals=[]
for p in momenta: for p in momenta:
y,x=np.histogram(ak.flatten(theta_recon[p]-theta_truth[p])*1000, bins=100, range=(-1,1))
accept=(nclusters[p]==3) &(pi0_converged[p])
y,x=np.histogram((theta_recon_corr[p]-theta_truth[p])[accept]*1000, bins=100, range=(-1,1))
bc=(x[1:]+x[:-1])/2 bc=(x[1:]+x[:-1])/2
from scipy.optimize import curve_fit from scipy.optimize import curve_fit
slc=abs(bc)<0.3 slc=abs(bc)<0.3
fnc=gauss fnc=gauss
p0=(100, 0, 0.06) p0=(100, 0, 0.06)
#print(bc[slc],y[slc])
sigma=np.sqrt(y[slc])+(y[slc]==0) sigma=np.sqrt(y[slc])+(y[slc]==0)
try: try:
coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000) coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000)
...@@ -232,16 +135,15 @@ plt.tight_layout() ...@@ -232,16 +135,15 @@ plt.tight_layout()
plt.savefig(outdir+"thetastar_recon.pdf") plt.savefig(outdir+"thetastar_recon.pdf")
#plt.show() #plt.show()
#vtx z
fig,axs=plt.subplots(1,3, figsize=(24, 8)) fig,axs=plt.subplots(1,3, figsize=(24, 8))
plt.sca(axs[0]) plt.sca(axs[0])
plt.title(f"$E_{{\\Lambda}}=100-275$ GeV") plt.title(f"$E_{{\\Lambda}}=100-275$ GeV")
x=[] x=[]
y=[] y=[]
for p in momenta: for p in momenta:
accept=(nclusters[p]==3) &(pi0_converged[p]) x+=list(ak.flatten(arrays_sim[p]['MCParticles.vertex.z'][:,3]+0*zvtx_recon[p])/1000)
x+=list(arrays_sim[p]['MCParticles.vertex.z'][:,3][accept]/1000) y+=list(ak.flatten(zvtx_recon[p])/1000)
y+=list(zvtx_recon[p][accept]/1000)
plt.scatter(x,y) plt.scatter(x,y)
#print(x,y) #print(x,y)
plt.xlabel("$z^{\\rm truth}_{\\rm vtx}$ [m]") plt.xlabel("$z^{\\rm truth}_{\\rm vtx}$ [m]")
...@@ -272,8 +174,8 @@ dsigmas=[] ...@@ -272,8 +174,8 @@ dsigmas=[]
xvals=[] xvals=[]
for p in momenta: for p in momenta:
accept=(nclusters[p]==3) &(pi0_converged[p])
a=list((zvtx_recon[p]-arrays_sim[p]['MCParticles.vertex.z'][:,3])[accept]/1000) a=ak.flatten((zvtx_recon[p]-arrays_sim[p]['MCParticles.vertex.z'][:,3])/1000)
y,x=np.histogram(a, bins=100, range=(-10,10)) y,x=np.histogram(a, bins=100, range=(-10,10))
bc=(x[1:]+x[:-1])/2 bc=(x[1:]+x[:-1])/2
...@@ -311,8 +213,7 @@ plt.sca(axs[0]) ...@@ -311,8 +213,7 @@ plt.sca(axs[0])
lambda_mass=1.115683 lambda_mass=1.115683
vals=[] vals=[]
for p in momenta: for p in momenta:
accept=(nclusters[p]==3) &(pi0_converged[p]) vals+=list(ak.flatten(mass_recon[p]))
vals+=list(mass_recon_corr[p][accept])
y,x,_= plt.hist(vals, bins=100, range=(1.0, 1.25)) y,x,_= plt.hist(vals, bins=100, range=(1.0, 1.25))
bc=(x[1:]+x[:-1])/2 bc=(x[1:]+x[:-1])/2
...@@ -323,7 +224,7 @@ plt.ylim(0, np.max(y)*1.2) ...@@ -323,7 +224,7 @@ plt.ylim(0, np.max(y)*1.2)
plt.xlim(1.0, 1.25) plt.xlim(1.0, 1.25)
from scipy.optimize import curve_fit from scipy.optimize import curve_fit
slc=abs(bc-lambda_mass)<0.07 slc=abs(bc-lambda_mass)<0.05
fnc=gauss fnc=gauss
p0=[100, lambda_mass, 0.04] p0=[100, lambda_mass, 0.04]
coeff, var_matrix = curve_fit(fnc, bc[slc], y[slc], p0=p0, coeff, var_matrix = curve_fit(fnc, bc[slc], y[slc], p0=p0,
...@@ -340,12 +241,11 @@ xvals=[] ...@@ -340,12 +241,11 @@ xvals=[]
sigmas=[] sigmas=[]
dsigmas=[] dsigmas=[]
for p in momenta: for p in momenta:
accept=(nclusters[p]==3) &(pi0_converged[p]) y,x= np.histogram(ak.flatten(mass_recon[p]), bins=100, range=(0.6,1.4))
y,x= np.histogram(mass_recon_corr[p][accept], bins=100, range=(0.6,1.4))
bc=(x[1:]+x[:-1])/2 bc=(x[1:]+x[:-1])/2
from scipy.optimize import curve_fit from scipy.optimize import curve_fit
slc=abs(bc-lambda_mass)<0.07 slc=abs(bc-lambda_mass)<0.05
fnc=gauss fnc=gauss
p0=[100, lambda_mass, 0.05] p0=[100, lambda_mass, 0.05]
try: try:
...@@ -367,3 +267,44 @@ plt.ylabel("$\\sigma[m_{\\Lambda}]$ [GeV]") ...@@ -367,3 +267,44 @@ plt.ylabel("$\\sigma[m_{\\Lambda}]$ [GeV]")
plt.ylim(0, 0.02) plt.ylim(0, 0.02)
plt.tight_layout() plt.tight_layout()
plt.savefig(outdir+"lambda_mass_rec.pdf") plt.savefig(outdir+"lambda_mass_rec.pdf")
#now for the CM stuff:
phi_residuals=[]
theta_residuals=[]
for p in momenta:
isNeutron=arrays_sim[p]['ReconstructedFarForwardZDCLambdaDecayProductsCM.PDG']==2112
pxcm=arrays_sim[p]['ReconstructedFarForwardZDCLambdaDecayProductsCM.momentum.x']
pycm=arrays_sim[p]['ReconstructedFarForwardZDCLambdaDecayProductsCM.momentum.y']
pzcm=arrays_sim[p]['ReconstructedFarForwardZDCLambdaDecayProductsCM.momentum.z']
import ROOT
px,py,pz,E=arrays_sim[p]['MCParticles.momentum.x'], arrays_sim[p]['MCParticles.momentum.y'], arrays_sim[p]['MCParticles.momentum.z'],np.sqrt(arrays_sim[p]['MCParticles.momentum.x']**2+arrays_sim[p]['MCParticles.momentum.y']**2+arrays_sim[p]['MCParticles.momentum.z']**2\
+arrays_sim[p]['MCParticles.mass']**2)
phicmtruth=list(np.repeat(-9999, len(arrays_sim[p])))
thetacmtruth=list(np.repeat(-9999, len(arrays_sim[p])))
for i in range(len(arrays_sim[p])):
l=ROOT.TLorentzVector(px[i,2], py[i,2], pz[i,2], E[i,2])
n=ROOT.TLorentzVector(px[i,3], py[i,3], pz[i,3], E[i,3])
ncm=n.Clone();
ncm.Boost(-l.BoostVector())
phicmtruth[i]=ncm.Phi()
thetacmtruth[i]=ncm.Theta()
arrays_sim[p]["phicmtruth"]=phicmtruth
arrays_sim[p]["thetacmtruth"]=thetacmtruth
phicmtruth=arrays_sim[p]["phicmtruth"]
thetacmtruth=arrays_sim[p]["thetacmtruth"]
phi_residuals=np.concatenate((phi_residuals, ak.flatten((np.arctan2(pycm,pxcm)[isNeutron]-phicmtruth)*np.sin(thetacmtruth))))
theta_residuals=np.concatenate((theta_residuals, ak.flatten(np.arctan2(np.hypot(pycm,pxcm),pzcm)[isNeutron]-thetacmtruth)))
plt.figure()
plt.hist(phi_residuals*1000, bins=100, range=(-300, 300))
plt.xlabel("$(\\phi^n_{cm,rec}-\\phi^n_{cm,truth})\\times\\sin\\theta^n_{cm,truth} [mrad]$")
plt.savefig(outdir+"neutron_phi_cm_res.pdf")
plt.figure()
plt.hist(1000*theta_residuals, bins=100, range=(-1000, 1000))
plt.xlabel("$\\theta^n_{cm,rec}-\\theta^n_{cm,truth}$ [mrad]")
plt.savefig(outdir+"neutron_theta_cm_res.pdf")
benchmarks/zdc_neutron/Snakefile
+ 1
− 1
View file @ 7ba375c3
...@@ -44,7 +44,7 @@ rule zdc_neutron_reco: ...@@ -44,7 +44,7 @@ rule zdc_neutron_reco:
set -m # monitor mode to prevent lingering processes set -m # monitor mode to prevent lingering processes
exec env DETECTOR_CONFIG={wildcards.DETECTOR_CONFIG} \ exec env DETECTOR_CONFIG={wildcards.DETECTOR_CONFIG} \
eicrecon {input} -Ppodio:output_file={output} \ eicrecon {input} -Ppodio:output_file={output} \
-Ppodio:output_collections=MCParticles,EcalFarForwardZDCRawHits,EcalFarForwardZDCRecHits,EcalFarForwardZDCClusters,HcalFarForwardZDCRawHits,HcalFarForwardZDCRecHits,HcalFarForwardZDCClusters,ReconstructedFarForwardZDCNeutrons -Ppodio:output_collections=MCParticles,EcalFarForwardZDCRawHits,EcalFarForwardZDCRecHits,EcalFarForwardZDCClusters,HcalFarForwardZDCRawHits,HcalFarForwardZDCRecHits,HcalFarForwardZDCClusters,ReconstructedFarForwardZDCNeutrals
""" """
......
benchmarks/zdc_neutron/analysis/fwd_neutrons_recon.C
+ 7
− 5
View file @ 7ba375c3
...@@ -112,11 +112,11 @@ void fwd_neutrons_recon(std::string inputfile, std::string outputfile){ ...@@ -112,11 +112,11 @@ void fwd_neutrons_recon(std::string inputfile, std::string outputfile){
TTreeReaderArray<float> hcal_cluster_energy(tr, "HcalFarForwardZDCClusters.energy"); TTreeReaderArray<float> hcal_cluster_energy(tr, "HcalFarForwardZDCClusters.energy");
//Reconstructed neutron quantity //Reconstructed neutron quantity
TTreeReaderArray<float> rec_neutron_energy(tr,"ReconstructedFarForwardZDCNeutrons.energy"); TTreeReaderArray<float> rec_neutron_energy(tr,"ReconstructedFarForwardZDCNeutrals.energy");
TTreeReaderArray<float> rec_neutron_px(tr,"ReconstructedFarForwardZDCNeutrons.momentum.x"); TTreeReaderArray<float> rec_neutron_px(tr,"ReconstructedFarForwardZDCNeutrals.momentum.x");
TTreeReaderArray<float> rec_neutron_py(tr,"ReconstructedFarForwardZDCNeutrons.momentum.y"); TTreeReaderArray<float> rec_neutron_py(tr,"ReconstructedFarForwardZDCNeutrals.momentum.y");
TTreeReaderArray<float> rec_neutron_pz(tr,"ReconstructedFarForwardZDCNeutrons.momentum.z"); TTreeReaderArray<float> rec_neutron_pz(tr,"ReconstructedFarForwardZDCNeutrals.momentum.z");
TTreeReaderArray<int> rec_neutron_PDG(tr,"ReconstructedFarForwardZDCNeutrals.PDG");
//Other variables //Other variables
int counter(0); int counter(0);
...@@ -195,6 +195,8 @@ void fwd_neutrons_recon(std::string inputfile, std::string outputfile){ ...@@ -195,6 +195,8 @@ void fwd_neutrons_recon(std::string inputfile, std::string outputfile){
//Reconstructed neutron(s) //Reconstructed neutron(s)
for(int irec=0;irec<rec_neutron_energy.GetSize();irec++){ for(int irec=0;irec<rec_neutron_energy.GetSize();irec++){
if (rec_neutron_PDG[irec] != 2112)
continue;
if(neut_true_rot.Theta()*1000. < 3.5 && hcal_clus_size==1){ if(neut_true_rot.Theta()*1000. < 3.5 && hcal_clus_size==1){
h1_neut_rec->Fill(rec_neutron_energy[irec]); h1_neut_rec->Fill(rec_neutron_energy[irec]);
......
benchmarks/zdc_photon/Snakefile
+ 1
− 1
View file @ 7ba375c3
...@@ -55,7 +55,7 @@ rule zdc_photon_recon: ...@@ -55,7 +55,7 @@ rule zdc_photon_recon:
shell: shell:
""" """
env DETECTOR_CONFIG={params.DETECTOR_CONFIG} \ env DETECTOR_CONFIG={params.DETECTOR_CONFIG} \
eicrecon {input.SIM_FILE} -Ppodio:output_file={output.REC_FILE} -Ppodio:output_collections=MCParticles,HcalFarForwardZDCRecHits,HcalFarForwardZDCClusters,HcalFarForwardZDCSubcellHits -Pjana:nevents={params.N_EVENTS} eicrecon {input.SIM_FILE} -Ppodio:output_file={output.REC_FILE} -Ppodio:output_collections=MCParticles,HcalFarForwardZDCRecHits,HcalFarForwardZDCClusters,HcalFarForwardZDCSubcellHits,ReconstructedFarForwardZDCNeutrals -Pjana:nevents={params.N_EVENTS}
""" """
rule zdc_photon_analysis: rule zdc_photon_analysis:
......
benchmarks/zdc_photon/analysis/zdc_photon_plots.py
+ 27
− 34
View file @ 7ba375c3
...@@ -27,7 +27,7 @@ for p in momenta: ...@@ -27,7 +27,7 @@ for p in momenta:
f'sim_output/zdc_photon/{config}_rec_zdc_photon_{p}GeV_{index}.edm4eic.root': 'events' f'sim_output/zdc_photon/{config}_rec_zdc_photon_{p}GeV_{index}.edm4eic.root': 'events'
for index in range(5) for index in range(5)
}) })
import awkward as ak
fig,axs=plt.subplots(1,3, figsize=(24, 8)) fig,axs=plt.subplots(1,3, figsize=(24, 8))
pvals=[] pvals=[]
resvals=[] resvals=[]
...@@ -35,27 +35,22 @@ dresvals=[] ...@@ -35,27 +35,22 @@ dresvals=[]
scalevals=[] scalevals=[]
dscalevals=[] dscalevals=[]
for p in momenta: for p in momenta:
selection=[len(arrays_sim[p]["HcalFarForwardZDCClusters.energy"][i])==1 for i in range(len(arrays_sim[p]))] selection=arrays_sim[p]["ReconstructedFarForwardZDCNeutrals.PDG"]==22
E=arrays_sim[p][selection]["HcalFarForwardZDCClusters.energy"] Etot=arrays_sim[p]["ReconstructedFarForwardZDCNeutrals.energy"][selection]
Etot=np.sum(E, axis=-1)
#print(len(Etot))
#print(p, res, mrecon)
if p==100: if p==100:
plt.sca(axs[0]) plt.sca(axs[0])
y, x, _=plt.hist(Etot, bins=100, range=(p*.75, p*1.25), histtype='step') y, x, _=plt.hist(ak.flatten(Etot), bins=100, range=(p*.75, p*1.25), histtype='step')
plt.ylabel("events") plt.ylabel("events")
plt.title(f"$p_{{\gamma}}$={p} GeV") plt.title(f"$p_{{\\gamma}}$={p} GeV")
plt.xlabel("$E^{\\gamma}_{recon}$ [GeV]") plt.xlabel("$E^{\\gamma}_{recon}$ [GeV]")
else: else:
y, x = np.histogram(Etot, bins=100, range=(p*.75, p*1.25)) y, x = np.histogram(ak.flatten(Etot), bins=100, range=(p*.75, p*1.25))
bc=(x[1:]+x[:-1])/2 bc=(x[1:]+x[:-1])/2
from scipy.optimize import curve_fit from scipy.optimize import curve_fit
slc=abs(bc-p)<10 slc=abs(bc-p)<10
fnc=gauss fnc=gauss
p0=[100, p, 10] p0=[100, p, 10]
#print(list(y), list(x))
coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0, 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) sigma=list(np.sqrt(y[slc])+(y[slc]==0)), maxfev=10000)
if p==100: if p==100:
...@@ -98,47 +93,45 @@ pvals=[] ...@@ -98,47 +93,45 @@ pvals=[]
resvals=[] resvals=[]
dresvals=[] dresvals=[]
for p in momenta: for p in momenta:
selection=[len(arrays_sim[p]["HcalFarForwardZDCClusters.energy"][i])==1 for i in range(len(arrays_sim[p]))] selection=arrays_sim[p]["ReconstructedFarForwardZDCNeutrals.PDG"]==22
x=arrays_sim[p][selection]["HcalFarForwardZDCClusters.position.x"][::,0] px_recon=arrays_sim[p]["ReconstructedFarForwardZDCNeutrals.momentum.x"][selection]
y=arrays_sim[p][selection]["HcalFarForwardZDCClusters.position.y"][::,0] py_recon=arrays_sim[p]["ReconstructedFarForwardZDCNeutrals.momentum.y"][selection]
z=arrays_sim[p][selection]["HcalFarForwardZDCClusters.position.z"][::,0] pz_recon=arrays_sim[p]["ReconstructedFarForwardZDCNeutrals.momentum.z"][selection]
theta_recon=np.arctan2(np.hypot(x*np.cos(-.025)-z*np.sin(-.025), y), z*np.cos(-.025)+x*np.sin(-.025)) theta_recon=np.arctan2(np.hypot(px_recon*np.cos(-.025)-pz_recon*np.sin(-.025), py_recon), pz_recon*np.cos(-.025)+px_recon*np.sin(-.025))
px=arrays_sim[p][selection]["MCParticles.momentum.x"][::,2] px=arrays_sim[p]["MCParticles.momentum.x"][::,2]
py=arrays_sim[p][selection]["MCParticles.momentum.y"][::,2] py=arrays_sim[p]["MCParticles.momentum.y"][::,2]
pz=arrays_sim[p][selection]["MCParticles.momentum.z"][::,2] pz=arrays_sim[p]["MCParticles.momentum.z"][::,2]
theta_truth=np.arctan2(np.hypot(px*np.cos(-.025)-pz*np.sin(-.025), py), pz*np.cos(-.025)+px*np.sin(-.025)) theta_truth=np.arctan2(np.hypot(px*np.cos(-.025)-pz*np.sin(-.025), py), pz*np.cos(-.025)+px*np.sin(-.025))
Etot=np.sum(E, axis=-1)
#print(p, res, mrecon)
if p==100: if p==100:
plt.sca(axs[0]) plt.sca(axs[0])
y, x, _=plt.hist(1000*(theta_recon-theta_truth), bins=100, range=(-0.5, 0.5), histtype='step') y, x, _=plt.hist(1000*ak.flatten(theta_recon-theta_truth), bins=100, range=(-0.5, 0.5), histtype='step')
plt.ylabel("events") plt.ylabel("events")
plt.title(f"$p_{{\gamma}}$={p} GeV") plt.title(f"$p_{{\\gamma}}$={p} GeV")
plt.xlabel("$\\theta^{\\gamma}_{recon}$ [mrad]") plt.xlabel("$\\theta^{\\gamma}_{recon}$ [mrad]")
else: else:
y, x = np.histogram(1000*(theta_recon-theta_truth), bins=100, range=(-0.5, 0.5)) y, x = np.histogram(1000*ak.flatten(theta_recon-theta_truth), bins=100, range=(-0.5, 0.5))
bc=(x[1:]+x[:-1])/2 bc=(x[1:]+x[:-1])/2
from scipy.optimize import curve_fit from scipy.optimize import curve_fit
slc=abs(bc)<0.2#1.5*np.std(1000*(theta_recon-theta_truth)) slc=abs(bc)<0.2#1.5*np.std(1000*(theta_recon-theta_truth))
fnc=gauss fnc=gauss
p0=[100, 0, 0.1] p0=[100, 0, 0.1]
#print(list(y), list(x)) try:
coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0, 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) sigma=list(np.sqrt(y[slc])+(y[slc]==0)), maxfev=10000)
if p==100: if p==100:
xx=np.linspace(-0.5,0.5, 100) xx=np.linspace(-0.5,0.5, 100)
plt.plot(xx, fnc(xx,*coeff)) plt.plot(xx, fnc(xx,*coeff))
pvals.append(p) pvals.append(p)
resvals.append(np.abs(coeff[2])) resvals.append(np.abs(coeff[2]))
dresvals.append(np.sqrt(var_matrix[2][2])) dresvals.append(np.sqrt(var_matrix[2][2]))
except:
pass
plt.sca(axs[1]) plt.sca(axs[1])
plt.errorbar(pvals, resvals, dresvals, ls='', marker='o') plt.errorbar(pvals, resvals, dresvals, ls='', marker='o')
#print(dresvals)
fnc=lambda E,a, b: np.hypot(a/np.sqrt(E), b) fnc=lambda E,a, b: np.hypot(a/np.sqrt(E), b)
#pvals, resvals, dresvals #pvals, resvals, dresvals
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment