File indexing completed on 2024-09-27 07:02:39
0001 import numpy as np
0002 import matplotlib.pyplot as plt
0003 import mplhep as hep
0004 import uproot
0005 import pandas as pd
0006 from scipy.optimize import curve_fit
0007 from matplotlib.backends.backend_pdf import PdfPages
0008 import os
0009 import awkward as ak
0010
0011 plt.figure()
0012 hep.set_style(hep.style.CMS)
0013 hep.set_style("CMS")
0014
0015 def gaussian(x, amp, mean, sigma):
0016 return amp * np.exp( -(x - mean)**2 / (2*sigma**2) )
0017
0018 def rotateY(xdata, zdata, angle):
0019 s = np.sin(angle)
0020 c = np.cos(angle)
0021 rotatedz = c*zdata - s*xdata
0022 rotatedx = s*zdata + c*xdata
0023 return rotatedx, rotatedz
0024
0025 Energy = [0.005, 0.01, 0.05, 0.1, 0.5, 1.0]
0026
0027
0028 df = pd.DataFrame({})
0029 for eng in Energy:
0030 tree = uproot.open(f'sim_output/zdc_lyso/{os.environ["DETECTOR_CONFIG"]}_gamma_{eng}GeV_theta_0deg_thru_0.3deg.eicrecon.tree.edm4eic.root')['events']
0031 ecal_reco_energy = list(map(sum, tree['EcalFarForwardZDCClusters/EcalFarForwardZDCClusters.energy'].array()))
0032 hcal_reco_energy = list(map(sum, tree['HcalFarForwardZDCClusters/HcalFarForwardZDCClusters.energy'].array()))
0033 ecal_rec_energy = list(map(sum, tree['EcalFarForwardZDCRecHits/EcalFarForwardZDCRecHits.energy'].array()))
0034 hcal_rec_energy = list(map(sum, tree['HcalFarForwardZDCRecHits/HcalFarForwardZDCRecHits.energy'].array()))
0035 ecal_reco_clusters = [len(row) if len(row)>=1 else 0 for row in tree['EcalFarForwardZDCClusters/EcalFarForwardZDCClusters.nhits'].array()]
0036 ecal_reco_nhits = [row[0] if len(row)>=1 else 0 for row in tree['EcalFarForwardZDCClusters/EcalFarForwardZDCClusters.nhits'].array()]
0037
0038 tree = uproot.open(f'sim_output/zdc_lyso/{os.environ["DETECTOR_CONFIG"]}_gamma_{eng}GeV_theta_0deg_thru_0.3deg.edm4hep.root')['events']
0039 ecal_sim_energy = list(map(sum, tree['EcalFarForwardZDCHits/EcalFarForwardZDCHits.energy'].array()))
0040 hcal_sim_energy = list(map(sum, tree['HcalFarForwardZDCHits/HcalFarForwardZDCHits.energy'].array()))
0041
0042 par_x = tree['MCParticles/MCParticles.momentum.x'].array()[:,2]
0043 par_y = tree['MCParticles/MCParticles.momentum.y'].array()[:,2]
0044 par_z = tree['MCParticles/MCParticles.momentum.z'].array()[:,2]
0045
0046 eng = int(eng*1000)
0047
0048 ecal_reco_energy = pd.DataFrame({f'ecal_reco_energy_{eng}': np.array(ecal_reco_energy, dtype=object)})
0049 hcal_reco_energy = pd.DataFrame({f'hcal_reco_energy_{eng}': np.array(hcal_reco_energy, dtype=object)})
0050 ecal_rec_energy = pd.DataFrame({f'ecal_rec_energy_{eng}': np.array(ecal_rec_energy, dtype=object)})
0051 hcal_rec_energy = pd.DataFrame({f'hcal_rec_energy_{eng}': np.array(hcal_rec_energy, dtype=object)})
0052 ecal_sim_energy = pd.DataFrame({f'ecal_sim_energy_{eng}': np.array(ecal_sim_energy, dtype=object)})
0053 hcal_sim_energy = pd.DataFrame({f'hcal_sim_energy_{eng}': np.array(hcal_sim_energy, dtype=object)})
0054 ecal_reco_nhits = pd.DataFrame({f'ecal_reco_nhits_{eng}': np.array(ecal_reco_nhits, dtype=object)})
0055 ecal_reco_clusters = pd.DataFrame({f'ecal_reco_clusters_{eng}': np.array(ecal_reco_clusters, dtype=object)})
0056 par_x = pd.DataFrame({f'par_x_{eng}': np.array(par_x.tolist(), dtype=object)})
0057 par_y = pd.DataFrame({f'par_y_{eng}': np.array(par_y.tolist(), dtype=object)})
0058 par_z = pd.DataFrame({f'par_z_{eng}': np.array(par_z.tolist(), dtype=object)})
0059
0060
0061 df = pd.concat([df,ecal_reco_energy,ecal_rec_energy,ecal_sim_energy,hcal_reco_energy,hcal_rec_energy,hcal_sim_energy,ecal_reco_clusters,ecal_reco_nhits,par_x,par_y,par_z],axis=1)
0062
0063
0064 mu = []
0065 sigma = []
0066 fig1, ax = plt.subplots(3,2,figsize=(20,10))
0067 fig1.suptitle('ZDC ECal Cluster Energy Reconstruction')
0068
0069 plt.tight_layout()
0070 for i in range(6):
0071 x = df[f'par_x_{eng}'].astype(float).to_numpy()
0072 y = df[f'par_y_{eng}'].astype(float).to_numpy()
0073 z = df[f'par_z_{eng}'].astype(float).to_numpy()
0074 x, z = rotateY(x,z, 0.025)
0075 theta = np.arccos(z/np.sqrt((x**2+y**2+z**2)))*1000
0076 condition = theta <= 3.5
0077
0078 plt.sca(ax[i%3,i//3])
0079 eng = int(Energy[i]*1000)
0080 plt.title(f'Gamma Energy: {eng} MeV')
0081 temp = np.array(df[f'ecal_reco_energy_{eng}'].astype(float).to_numpy()[condition])*1000
0082 hist, x = np.histogram(temp,bins=np.linspace(min(temp),max(temp)+np.std(abs(temp)),2*int(np.sqrt(len(temp)))))
0083 x = x[1:]/2 + x[:-1]/2
0084 plt.errorbar(x,hist,yerr=np.sqrt(hist),fmt='-o',label='Cluster')
0085 coeff, covar = curve_fit(gaussian,x[1:],hist[1:],p0=(max(hist[x>=np.std(abs(temp))]),np.mean(temp[temp!=0]),np.std(temp[temp!=0])))
0086
0087 mu.append(coeff[1])
0088 sigma.append(coeff[2])
0089
0090 temp = np.array(df[f'ecal_rec_energy_{eng}'].astype(float).to_numpy()[condition])*1000
0091 hist, x = np.histogram(temp,bins=np.linspace(min(temp),max(temp)+np.std(abs(temp)),2*int(np.sqrt(len(temp)))))
0092 x = x[1:]/2 + x[:-1]/2
0093 plt.errorbar(x,hist,yerr=np.sqrt(hist),fmt='-o',label='Digitization')
0094 coeff, covar = curve_fit(gaussian,x[1:],hist[1:],p0=(max(hist[x>=np.std(abs(temp))]),np.mean(temp[temp!=0]),np.std(temp[temp!=0])))
0095
0096 mu.append(coeff[1])
0097 sigma.append(coeff[2])
0098
0099 temp = np.array(df[f'ecal_sim_energy_{eng}'].astype(float).to_numpy()[condition])*1000
0100 hist, x = np.histogram(temp,bins=np.linspace(min(temp),max(temp)+np.std(abs(temp)),2*int(np.sqrt(len(temp)))))
0101 x = x[1:]/2 + x[:-1]/2
0102 plt.errorbar(x,hist,yerr=np.sqrt(hist),fmt='-o',label='Simulation')
0103 coeff, covar = curve_fit(gaussian,x[1:],hist[1:],p0=(max(hist[x>=np.std(abs(temp))]),np.mean(temp[temp!=0]),np.std(temp[temp!=0])))
0104
0105 mu.append(coeff[1])
0106 sigma.append(coeff[2])
0107
0108 plt.xlabel('Energy (MeV)')
0109 plt.legend()
0110
0111
0112
0113 mu = np.array(mu)
0114 sigma = np.array(sigma)
0115
0116 plt.show()
0117
0118 fig2, (ax1,ax2) = plt.subplots(2,1,figsize=(15,10),sharex=True)
0119
0120 plt.tight_layout()
0121
0122 ax1.scatter(np.array(Energy)*1000, mu[::3], label='cluster')
0123 ax1.scatter(np.array(Energy)*1000, mu[1::3], label='digitization')
0124 ax1.scatter(np.array(Energy)*1000, mu[2::3], label='simulation')
0125
0126 ax1.plot([4.5,1000],[4.5,1000],c='black',label='x=y')
0127 ax1.set_ylabel('Reconstructed Energy (MeV)')
0128 ax1.set_yscale('log')
0129 ax1.legend()
0130 ax1.set_title('ECal Craterlake Cluster Energy Reconstruction')
0131
0132 ax2.errorbar(np.array(Energy)*1000, abs(sigma[::3]/mu[::3])*100, fmt='-o', label='cluster')
0133 ax2.errorbar(np.array(Energy)*1000, abs(sigma[1::3]/mu[1::3])*100, fmt='-o', label='digitization')
0134 ax2.errorbar(np.array(Energy)*1000, abs(sigma[2::3]/mu[2::3])*100, fmt='-o', label='simulation')
0135
0136 ax2.set_ylabel('Resolution (%)')
0137 ax2.set_xlabel('Gamma Energy (MeV)')
0138 ax2.set_xscale('log')
0139 ax2.legend()
0140
0141
0142
0143 plt.show()
0144
0145
0146 htower = []
0147 herr = []
0148 hmean = []
0149 hhits = []
0150 hhits_cut = []
0151 emean = []
0152 ehits = []
0153 etower = []
0154 eerr = []
0155 ehits_cut = []
0156
0157 fig3, ax = plt.subplots(2,3,figsize=(20,10))
0158 fig3.suptitle('ZDC Simulation Energy Reconstruction')
0159 for i in range(6):
0160 plt.sca(ax[i//3,i%3])
0161 eng = int(Energy[i]*1000)
0162
0163 x = df[f'par_x_{eng}'].astype(float).to_numpy()
0164 y = df[f'par_y_{eng}'].astype(float).to_numpy()
0165 z = df[f'par_z_{eng}'].astype(float).to_numpy()
0166 x, z = rotateY(x,z, 0.025)
0167 theta = np.arccos(z/np.sqrt((x**2+y**2+z**2)))*1000
0168 condition = theta <= 3.5
0169
0170 plt.title(f'Gamma Energy: {eng} MeV')
0171 energy1 = df[f'hcal_sim_energy_{eng}'].astype(float).to_numpy()
0172 hist, x = np.histogram(energy1*1000,bins=np.logspace(0,3,200))
0173 x = x[1:]/2 + x[:-1]/2
0174 plt.plot(x,hist,marker='o',label="HCal")
0175 hhits.append(len(energy1[energy1!=0]))
0176 condition1 = energy1!=0
0177 hhits_cut.append(len(energy1[condition & condition1])/len(condition[condition==True]))
0178 energy = df[f'ecal_sim_energy_{eng}'].astype(float).to_numpy()
0179 hist, x = np.histogram(energy*1000,bins=np.logspace(0,3,200))
0180 x = x[1:]/2 + x[:-1]/2
0181 plt.plot(x,hist,marker='o',label="ECal")
0182 emean.append(sum(energy[energy!=0])*1000/len(energy[energy!=0]))
0183 hmean.append(sum(energy1[energy!=0])*1000/len(energy[energy!=0]))
0184 condition1 = energy!=0
0185 ehits_cut.append(len(energy[condition & condition1])/len(condition[condition==True]))
0186 ehits.append(len(energy[energy!=0]))
0187 plt.legend()
0188 plt.xscale('log')
0189 plt.xlim(1e0,1e3)
0190
0191
0192
0193
0194
0195 plt.xlabel('Energy (MeV)')
0196
0197
0198 plt.show()
0199
0200 fig4, ax = plt.subplots(2,1,sharex=True,gridspec_kw={'height_ratios': [2,1]})
0201 plt.sca(ax[0])
0202 plt.errorbar(np.array(Energy)*1000,np.array(hmean)*47.619+np.array(emean),label='HCal/sf+ECal',fmt='-o')
0203 plt.errorbar(np.array(Energy)*1000,emean,label='ECal',fmt='-o')
0204 plt.legend()
0205 plt.yscale('log')
0206 plt.xscale('log')
0207 plt.ylabel('Simulation Energy (MeV)')
0208 plt.sca(ax[1])
0209 plt.errorbar(np.array(Energy)*1000,(1 - np.array(emean)/(np.array(hmean)*47.619+np.array(emean)))*100,label='Total/ECal',fmt='-o')
0210 plt.legend()
0211 plt.ylabel('Fraction of energy\n deposited in Hcal (%)')
0212 plt.xlabel('Truth Energy (MeV)')
0213
0214 plt.tight_layout()
0215 plt.show()
0216
0217 fig5 = plt.figure()
0218 plt.errorbar(np.array(Energy)*1000,np.array(hhits)/1000*100,label='HCal Hits',fmt='-o')
0219 plt.errorbar(np.array(Energy)*1000,np.array(ehits)/1000*100,label='ECal Hits',fmt='-o')
0220
0221
0222 plt.errorbar(np.array(Energy)*1000,np.array(hhits_cut)*100,label='HCal Hits with 3.5 mRad cut',fmt='-^')
0223 plt.errorbar(np.array(Energy)*1000,np.array(ehits_cut)*100,label='ECal Hits with 3.5 mRad cut',fmt='-^')
0224
0225
0226
0227 plt.legend()
0228 plt.xlabel('Simulation Truth Gamma Energy (MeV)')
0229 plt.ylabel('Fraction of Events with non-zero energy (%)')
0230
0231 plt.xscale('log')
0232 plt.show()
0233
0234 fig6, ax = plt.subplots(2,3,figsize=(20,10))
0235 fig6.suptitle('ZDC Clustering')
0236 fig6.tight_layout(pad=1.8)
0237 for i in range(6):
0238 plt.sca(ax[i//3,i%3])
0239 eng = int(Energy[i]*1000)
0240
0241 x = df[f'par_x_{eng}'].astype(float).to_numpy()
0242 y = df[f'par_y_{eng}'].astype(float).to_numpy()
0243 z = df[f'par_z_{eng}'].astype(float).to_numpy()
0244 x, z = rotateY(x,z, 0.025)
0245 theta = np.arccos(z/np.sqrt((x**2+y**2+z**2)))*1000
0246 condition = theta <= 3.5
0247
0248 plt.hist(df[f'ecal_reco_clusters_{eng}'][condition],bins=np.linspace(0,5,6))
0249 plt.xlabel('Number of Clusters')
0250 plt.title(f'Gamma Energy: {eng} MeV')
0251 plt.show()
0252
0253 fig7, ax = plt.subplots(2,3,figsize=(20,10))
0254 fig7.suptitle('ZDC Towering in Clusters')
0255 fig7.tight_layout(pad=1.8)
0256 for i in range(6):
0257 plt.sca(ax[i//3,i%3])
0258 eng = int(Energy[i]*1000)
0259
0260 x = df[f'par_x_{eng}'].astype(float).to_numpy()
0261 y = df[f'par_y_{eng}'].astype(float).to_numpy()
0262 z = df[f'par_z_{eng}'].astype(float).to_numpy()
0263 x, z = rotateY(x,z, 0.025)
0264 theta = np.arccos(z/np.sqrt((x**2+y**2+z**2)))*1000
0265 condition = theta <= 3.5
0266
0267 plt.hist(df[f'ecal_reco_nhits_{eng}'][condition],bins=np.linspace(0,max(df[f'ecal_reco_nhits_{eng}'][condition]),max(df[f'ecal_reco_nhits_{eng}'][condition])+1))
0268 plt.xlabel('Number of tower in Clusters')
0269 plt.title(f'Gamma Energy: {eng} MeV')
0270 plt.show()
0271
0272
0273
0274 with PdfPages(f'results/{os.environ["DETECTOR_CONFIG"]}/zdc_lyso/plots.pdf') as pdf:
0275 pdf.savefig(fig1)
0276 pdf.savefig(fig2)
0277 pdf.savefig(fig3)
0278 pdf.savefig(fig4)
0279 pdf.savefig(fig5)
0280 pdf.savefig(fig6)
0281 pdf.savefig(fig7)