From 660df3252ee7f05fd472df091433e8fb17288a8a Mon Sep 17 00:00:00 2001 From: steinber Date: Wed, 21 Aug 2024 12:34:17 -0400 Subject: [PATCH 01/30] first commit of LFHCAL benchmarks, with analysis.py still unmodified from zdc_lyso example. Not working just yet --- benchmarks/lfhcal/README.md | 11 + benchmarks/lfhcal/Snakefile | 87 ++++++++ benchmarks/lfhcal/analysis/analysis.py | 281 +++++++++++++++++++++++++ benchmarks/lfhcal/config.yml | 17 ++ benchmarks/lfhcal/gen_particles.cxx | 126 +++++++++++ 5 files changed, 522 insertions(+) create mode 100644 benchmarks/lfhcal/README.md create mode 100644 benchmarks/lfhcal/Snakefile create mode 100644 benchmarks/lfhcal/analysis/analysis.py create mode 100644 benchmarks/lfhcal/config.yml create mode 100644 benchmarks/lfhcal/gen_particles.cxx diff --git a/benchmarks/lfhcal/README.md b/benchmarks/lfhcal/README.md new file mode 100644 index 00000000..5f56298b --- /dev/null +++ b/benchmarks/lfhcal/README.md @@ -0,0 +1,11 @@ +Detector Benchmark for LFHCAL +=============================== + +## Overview +This benchmark generates events with a range of particle species at fixed energies, all aimed at, or forward of, the LFHCAL, and produces basic plots of clustering performance, both in isolation and combined with the FEMCAL. + + +## Contacts +[Peter Steinberg](steinberg@bnl.gov) + + diff --git a/benchmarks/lfhcal/Snakefile b/benchmarks/lfhcal/Snakefile new file mode 100644 index 00000000..d0204ce3 --- /dev/null +++ b/benchmarks/lfhcal/Snakefile @@ -0,0 +1,87 @@ +import os + + +rule lfhcal_sim_hepmc: + input: + script = "benchmarks/lfhcal/gen_particles.cxx", + output: + hepmcfile="data/{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc", + log: + "data/{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc.log", + params: + num_events=1000, + shell: + """ +root -l -b -q '{input.script}({params.num_events}, "{output.hepmcfile}", "{wildcards.PARTICLE}", {wildcards.THETA_MIN}, {wildcards.THETA_MAX}, 0, 360, {wildcards.BEAM_ENERGY})' +""" + + +rule lfhcal_sim: + input: + hepmcfile="data/{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc", + warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root", + output: + "sim_output/lfhcal/{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.edm4hep.root", + log: + "sim_output/lfhcal/{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.edm4hep.root.log", + params: + num_events=1000, + shell: + """ +npsim \ + --runType batch \ + -v WARNING \ + --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \ + --numberOfEvents {params.num_events} \ + --inputFiles {input.hepmcfile} \ + --outputFile {output} +""" + + +rule lfhcal_reco: + input: + "sim_output/lfhcal/{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.edm4hep.root", + output: + "sim_output/lfhcal/{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.eicrecon.tree.edm4eic.root", + log: + "sim_output/lfhcal/{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.eicrecon.tree.edm4eic.root.log", + shell: + """ +eicrecon -Ppodio:output_collections=MCParticles,ReconstructedParticles,LFHCALTruthClusters,LFHCALClusters,LFHCALHits,EcalEndcapPTruthClusters,EcalEndcapPClusters,EcalEnd\ +capPHits {input} +mv podio_output.root {output} +""" + + +rule lfhcal_analysis: + input: + expand("sim_output/lfhcal/{{DETECTOR_CONFIG}}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.eicrecon.tree.edm4eic.root", + PARTICLE=["neutron","gamma"], + BEAM_ENERGY=["2", "5", "10", "20", "40"], + THETA_MIN=["0"], + THETA_MAX=["50"]), + script="benchmarks/lfhcal/analysis/analysis.py", + output: + "results/{DETECTOR_CONFIG}/lfhcal/plots.pdf", + shell: + """ +python {input.script} +""" + + +# Examples of invocation +rule lfhcal_create_all_hepmc: + input: + expand("data/{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc", + PARTICLE=["gamma"], + BEAM_ENERGY=["0.005", "0.01", "0.05", "0.1", "0.5", "1.0"], + THETA_MIN=["0"], + THETA_MAX=["0.3"]) + + +rule lfhcal_run_all_locally: + input: + "results/" + os.environ["DETECTOR_CONFIG"] + "/lfhcal/plots.pdf", + message: + "See output in {input[0]}" + diff --git a/benchmarks/lfhcal/analysis/analysis.py b/benchmarks/lfhcal/analysis/analysis.py new file mode 100644 index 00000000..1ef1f828 --- /dev/null +++ b/benchmarks/lfhcal/analysis/analysis.py @@ -0,0 +1,281 @@ +import numpy as np +import matplotlib.pyplot as plt +import mplhep as hep +import uproot +import pandas as pd +from scipy.optimize import curve_fit +from matplotlib.backends.backend_pdf import PdfPages +import os +import awkward as ak + +plt.figure() +hep.set_style(hep.style.CMS) +hep.set_style("CMS") + +def gaussian(x, amp, mean, sigma): + return amp * np.exp( -(x - mean)**2 / (2*sigma**2) ) + +def rotateY(xdata, zdata, angle): + s = np.sin(angle) + c = np.cos(angle) + rotatedz = c*zdata - s*xdata + rotatedx = s*zdata + c*xdata + return rotatedx, rotatedz + +Energy = [0.005, 0.01, 0.05, 0.1, 0.5, 1.0] + + +df = pd.DataFrame({}) +for eng in Energy: + tree = uproot.open(f'sim_output/lfhcal/{os.environ["DETECTOR_CONFIG"]}_gamma_{eng}GeV_theta_0deg_thru_0.3deg.eicrecon.tree.edm4eic.root')['events'] + ecal_reco_energy = list(map(sum, tree['EcalFarForwardZDCClusters/EcalFarForwardZDCClusters.energy'].array())) + hcal_reco_energy = list(map(sum, tree['HcalFarForwardZDCClusters/HcalFarForwardZDCClusters.energy'].array())) + ecal_rec_energy = list(map(sum, tree['EcalFarForwardZDCRecHits/EcalFarForwardZDCRecHits.energy'].array())) + hcal_rec_energy = list(map(sum, tree['HcalFarForwardZDCRecHits/HcalFarForwardZDCRecHits.energy'].array())) + ecal_reco_clusters = [len(row) if len(row)>=1 else 0 for row in tree['EcalFarForwardZDCClusters/EcalFarForwardZDCClusters.nhits'].array()] + ecal_reco_nhits = [row[0] if len(row)>=1 else 0 for row in tree['EcalFarForwardZDCClusters/EcalFarForwardZDCClusters.nhits'].array()] + + tree = uproot.open(f'sim_output/lfhcal/{os.environ["DETECTOR_CONFIG"]}_gamma_{eng}GeV_theta_0deg_thru_0.3deg.edm4hep.root')['events'] + ecal_sim_energy = list(map(sum, tree['EcalFarForwardZDCHits/EcalFarForwardZDCHits.energy'].array())) + hcal_sim_energy = list(map(sum, tree['HcalFarForwardZDCHits/HcalFarForwardZDCHits.energy'].array())) + + par_x = tree['MCParticles/MCParticles.momentum.x'].array()[:,2] + par_y = tree['MCParticles/MCParticles.momentum.y'].array()[:,2] + par_z = tree['MCParticles/MCParticles.momentum.z'].array()[:,2] + + eng = int(eng*1000) + + ecal_reco_energy = pd.DataFrame({f'ecal_reco_energy_{eng}': np.array(ecal_reco_energy, dtype=object)}) + hcal_reco_energy = pd.DataFrame({f'hcal_reco_energy_{eng}': np.array(hcal_reco_energy, dtype=object)}) + ecal_rec_energy = pd.DataFrame({f'ecal_rec_energy_{eng}': np.array(ecal_rec_energy, dtype=object)}) + hcal_rec_energy = pd.DataFrame({f'hcal_rec_energy_{eng}': np.array(hcal_rec_energy, dtype=object)}) + ecal_sim_energy = pd.DataFrame({f'ecal_sim_energy_{eng}': np.array(ecal_sim_energy, dtype=object)}) + hcal_sim_energy = pd.DataFrame({f'hcal_sim_energy_{eng}': np.array(hcal_sim_energy, dtype=object)}) + ecal_reco_nhits = pd.DataFrame({f'ecal_reco_nhits_{eng}': np.array(ecal_reco_nhits, dtype=object)}) + ecal_reco_clusters = pd.DataFrame({f'ecal_reco_clusters_{eng}': np.array(ecal_reco_clusters, dtype=object)}) + par_x = pd.DataFrame({f'par_x_{eng}': np.array(par_x.tolist(), dtype=object)}) + par_y = pd.DataFrame({f'par_y_{eng}': np.array(par_y.tolist(), dtype=object)}) + par_z = pd.DataFrame({f'par_z_{eng}': np.array(par_z.tolist(), dtype=object)}) + + + 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) + + +mu = [] +sigma = [] +fig1, ax = plt.subplots(3,2,figsize=(20,10)) +fig1.suptitle('ZDC ECal Cluster Energy Reconstruction') + +plt.tight_layout() +for i in range(6): + x = df[f'par_x_{eng}'].astype(float).to_numpy() + y = df[f'par_y_{eng}'].astype(float).to_numpy() + z = df[f'par_z_{eng}'].astype(float).to_numpy() + x, z = rotateY(x,z, 0.025) + theta = np.arccos(z/np.sqrt((x**2+y**2+z**2)))*1000 + condition = theta <= 3.5 + + plt.sca(ax[i%3,i//3]) + eng = int(Energy[i]*1000) + plt.title(f'Gamma Energy: {eng} MeV') + temp = np.array(df[f'ecal_reco_energy_{eng}'].astype(float).to_numpy()[condition])*1000 + hist, x = np.histogram(temp,bins=np.linspace(min(temp),max(temp)+np.std(abs(temp)),2*int(np.sqrt(len(temp))))) + x = x[1:]/2 + x[:-1]/2 + plt.errorbar(x,hist,yerr=np.sqrt(hist),fmt='-o',label='Cluster') + 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]))) + #plt.plot(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),gaussian(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),*coeff)) + mu.append(coeff[1]) + sigma.append(coeff[2]) + + temp = np.array(df[f'ecal_rec_energy_{eng}'].astype(float).to_numpy()[condition])*1000 + hist, x = np.histogram(temp,bins=np.linspace(min(temp),max(temp)+np.std(abs(temp)),2*int(np.sqrt(len(temp))))) + x = x[1:]/2 + x[:-1]/2 + plt.errorbar(x,hist,yerr=np.sqrt(hist),fmt='-o',label='Digitization') + 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]))) + #plt.plot(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),gaussian(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),*coeff)) + mu.append(coeff[1]) + sigma.append(coeff[2]) + + temp = np.array(df[f'ecal_sim_energy_{eng}'].astype(float).to_numpy()[condition])*1000 + hist, x = np.histogram(temp,bins=np.linspace(min(temp),max(temp)+np.std(abs(temp)),2*int(np.sqrt(len(temp))))) + x = x[1:]/2 + x[:-1]/2 + plt.errorbar(x,hist,yerr=np.sqrt(hist),fmt='-o',label='Simulation') + 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]))) + #plt.plot(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),gaussian(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),*coeff)) + mu.append(coeff[1]) + sigma.append(coeff[2]) + + plt.xlabel('Energy (MeV)') + plt.legend() + +#plt.savefig('results/Energy_reconstruction_cluster.pdf') + +mu = np.array(mu) +sigma = np.array(sigma) + +plt.show() + +fig2, (ax1,ax2) = plt.subplots(2,1,figsize=(15,10),sharex=True) + +plt.tight_layout() +# Plot data on primary axis +ax1.scatter(np.array(Energy)*1000, mu[::3], label='cluster') +ax1.scatter(np.array(Energy)*1000, mu[1::3], label='digitization') +ax1.scatter(np.array(Energy)*1000, mu[2::3], label='simulation') + +ax1.plot([4.5,1000],[4.5,1000],c='black',label='x=y') +ax1.set_ylabel('Reconstructed Energy (MeV)') +ax1.set_yscale('log') +ax1.legend() +ax1.set_title('ECal Craterlake Cluster Energy Reconstruction') + +ax2.errorbar(np.array(Energy)*1000, abs(sigma[::3]/mu[::3])*100, fmt='-o', label='cluster') +ax2.errorbar(np.array(Energy)*1000, abs(sigma[1::3]/mu[1::3])*100, fmt='-o', label='digitization') +ax2.errorbar(np.array(Energy)*1000, abs(sigma[2::3]/mu[2::3])*100, fmt='-o', label='simulation') + +ax2.set_ylabel('Resolution (%)') +ax2.set_xlabel('Gamma Energy (MeV)') +ax2.set_xscale('log') +ax2.legend() + +#plt.savefig('results/Energy_resolution.pdf') + +plt.show() + + +htower = [] +herr = [] +hmean = [] +hhits = [] +hhits_cut = [] +emean = [] +ehits = [] +etower = [] +eerr = [] +ehits_cut = [] + +fig3, ax = plt.subplots(2,3,figsize=(20,10)) +fig3.suptitle('ZDC Simulation Energy Reconstruction') +for i in range(6): + plt.sca(ax[i//3,i%3]) + eng = int(Energy[i]*1000) + + x = df[f'par_x_{eng}'].astype(float).to_numpy() + y = df[f'par_y_{eng}'].astype(float).to_numpy() + z = df[f'par_z_{eng}'].astype(float).to_numpy() + x, z = rotateY(x,z, 0.025) + theta = np.arccos(z/np.sqrt((x**2+y**2+z**2)))*1000 + condition = theta <= 3.5 + + plt.title(f'Gamma Energy: {eng} MeV') + energy1 = df[f'hcal_sim_energy_{eng}'].astype(float).to_numpy()#df.eval(f'hcal_sim_energy_{eng}').apply(lambda row: sum(row)) + hist, x = np.histogram(energy1*1000,bins=np.logspace(0,3,200)) + x = x[1:]/2 + x[:-1]/2 + plt.plot(x,hist,marker='o',label="HCal") + hhits.append(len(energy1[energy1!=0])) + condition1 = energy1!=0 + hhits_cut.append(len(energy1[condition & condition1])/len(condition[condition==True])) + energy = df[f'ecal_sim_energy_{eng}'].astype(float).to_numpy()#df.eval(f'ecal_sim_energy_{eng}').apply(lambda row: sum(row)) + hist, x = np.histogram(energy*1000,bins=np.logspace(0,3,200)) + x = x[1:]/2 + x[:-1]/2 + plt.plot(x,hist,marker='o',label="ECal") + emean.append(sum(energy[energy!=0])*1000/len(energy[energy!=0])) + hmean.append(sum(energy1[energy!=0])*1000/len(energy[energy!=0])) + condition1 = energy!=0 + ehits_cut.append(len(energy[condition & condition1])/len(condition[condition==True])) + ehits.append(len(energy[energy!=0])) + plt.legend() + plt.xscale('log') + plt.xlim(1e0,1e3) + + + + + +plt.xlabel('Energy (MeV)') + +#plt.savefig('results/Energy_deposition.pdf') +plt.show() + +fig4, ax = plt.subplots(2,1,sharex=True,gridspec_kw={'height_ratios': [2,1]}) +plt.sca(ax[0]) +plt.errorbar(np.array(Energy)*1000,np.array(hmean)*47.619+np.array(emean),label='HCal/sf+ECal',fmt='-o') +plt.errorbar(np.array(Energy)*1000,emean,label='ECal',fmt='-o') +plt.legend() +plt.yscale('log') +plt.xscale('log') +plt.ylabel('Simulation Energy (MeV)') +plt.sca(ax[1]) +plt.errorbar(np.array(Energy)*1000,(1 - np.array(emean)/(np.array(hmean)*47.619+np.array(emean)))*100,label='Total/ECal',fmt='-o') +plt.legend() +plt.ylabel('Fraction of energy\n deposited in Hcal (%)') +plt.xlabel('Truth Energy (MeV)') +#plt.savefig('results/Energy_ratio_and_Leakage.pdf') +plt.tight_layout() +plt.show() + +fig5 = plt.figure() +plt.errorbar(np.array(Energy)*1000,np.array(hhits)/1000*100,label='HCal Hits',fmt='-o') +plt.errorbar(np.array(Energy)*1000,np.array(ehits)/1000*100,label='ECal Hits',fmt='-o') +#plt.errorbar(np.array(Energy)*1000,np.array(hhits)/np.array(ehits)*100,label='HCal / ECal',fmt='-o',c='b') + +plt.errorbar(np.array(Energy)*1000,np.array(hhits_cut)*100,label='HCal Hits with 3.5 mRad cut',fmt='-^') +plt.errorbar(np.array(Energy)*1000,np.array(ehits_cut)*100,label='ECal Hits with 3.5 mRad cut',fmt='-^') +#plt.errorbar(np.array(Energy)*1000,np.array(hhits_cut)/np.array(ehits_cut)*100,label='HCal / ECal with 3.5 mRad cut',fmt='-^',c='b') +### 3mrad cuts + +plt.legend() +plt.xlabel('Simulation Truth Gamma Energy (MeV)') +plt.ylabel('Fraction of Events with non-zero energy (%)') +#plt.savefig('results/Hits.pdf') +plt.xscale('log') +plt.show() + +fig6, ax = plt.subplots(2,3,figsize=(20,10)) +fig6.suptitle('ZDC Clustering') +fig6.tight_layout(pad=1.8) +for i in range(6): + plt.sca(ax[i//3,i%3]) + eng = int(Energy[i]*1000) + + x = df[f'par_x_{eng}'].astype(float).to_numpy() + y = df[f'par_y_{eng}'].astype(float).to_numpy() + z = df[f'par_z_{eng}'].astype(float).to_numpy() + x, z = rotateY(x,z, 0.025) + theta = np.arccos(z/np.sqrt((x**2+y**2+z**2)))*1000 + condition = theta <= 3.5 + + plt.hist(df[f'ecal_reco_clusters_{eng}'][condition],bins=np.linspace(0,5,6)) + plt.xlabel('Number of Clusters') + plt.title(f'Gamma Energy: {eng} MeV') +plt.show() + +fig7, ax = plt.subplots(2,3,figsize=(20,10)) +fig7.suptitle('ZDC Towering in Clusters') +fig7.tight_layout(pad=1.8) +for i in range(6): + plt.sca(ax[i//3,i%3]) + eng = int(Energy[i]*1000) + + x = df[f'par_x_{eng}'].astype(float).to_numpy() + y = df[f'par_y_{eng}'].astype(float).to_numpy() + z = df[f'par_z_{eng}'].astype(float).to_numpy() + x, z = rotateY(x,z, 0.025) + theta = np.arccos(z/np.sqrt((x**2+y**2+z**2)))*1000 + condition = theta <= 3.5 + + 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)) + plt.xlabel('Number of tower in Clusters') + plt.title(f'Gamma Energy: {eng} MeV') +plt.show() + + +#pdfs = ['results/Energy_reconstruction_cluster.pdf','results/Energy_resolution.pdf','results/Energy_deposition.pdf','results/Energy_ratio_and_Leakage.pdf','results/Hits.pdf'] +with PdfPages(f'results/{os.environ["DETECTOR_CONFIG"]}/lfhcal/plots.pdf') as pdf: + pdf.savefig(fig1) + pdf.savefig(fig2) + pdf.savefig(fig3) + pdf.savefig(fig4) + pdf.savefig(fig5) + pdf.savefig(fig6) + pdf.savefig(fig7) diff --git a/benchmarks/lfhcal/config.yml b/benchmarks/lfhcal/config.yml new file mode 100644 index 00000000..5f589a26 --- /dev/null +++ b/benchmarks/lfhcal/config.yml @@ -0,0 +1,17 @@ +sim:lfhcal: + extends: .det_benchmark + stage: simulate + script: + - snakemake --cores 1 run_all_locally + retry: + max: 2 + when: + - runner_system_failure + +collect_results:lfhcal: + extends: .det_benchmark + stage: collect + needs: + - "sim:lfhcal" + script: + - ls -lrht diff --git a/benchmarks/lfhcal/gen_particles.cxx b/benchmarks/lfhcal/gen_particles.cxx new file mode 100644 index 00000000..b84b1170 --- /dev/null +++ b/benchmarks/lfhcal/gen_particles.cxx @@ -0,0 +1,126 @@ +#include "HepMC3/GenEvent.h" +#include "HepMC3/ReaderAscii.h" +#include "HepMC3/WriterAscii.h" +#include "HepMC3/Print.h" + +#include "TRandom3.h" +#include "TVector3.h" + +#include +#include +#include +#include +#include +#include +#include + +using namespace HepMC3; + +// Generate single electron as input to the Insert simulation. +// -- +// We generate events with a constant polar angle with respect to +// the proton direction and then rotate so that the events are given +// in normal lab coordinate system +// -- +void gen_particles( + int n_events = 10, + const char* out_fname = "gen_particles.hepmc", + TString particle_name = "e-", + double th_min = 3., // Minimum polar angle, in degrees + double th_max = 3., // Maximum polar angle, in degrees + double phi_min = 0., // Minimum azimuthal angle, in degrees + double phi_max = 360., // Maximum azimuthal angle, in degrees + double p = 10., // Momentum in GeV/c + int dist = 0 //Momentum distribution: 0=fixed, 1=uniform, 2=Gaussian + ) +{ + WriterAscii hepmc_output(out_fname); + int events_parsed = 0; + GenEvent evt(Units::GEV, Units::MM); + + // Random number generator + TRandom3 *r1 = new TRandom3(0); //Use time as random seed + + // Getting generated particle information + TDatabasePDG *pdg = new TDatabasePDG(); + TParticlePDG *particle = pdg->GetParticle(particle_name); + const double mass = particle->Mass(); + const int pdgID = particle->PdgCode(); + + for (events_parsed = 0; events_parsed < n_events; events_parsed++) { + + //Set the event number + evt.set_event_number(events_parsed); + + // FourVector(px,py,pz,e,pdgid,status) + // type 4 is beam + // pdgid 11 - electron + // pdgid 111 - pi0 + // pdgid 2212 - proton + GenParticlePtr p1 = + std::make_shared(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4); + GenParticlePtr p2 = std::make_shared( + FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4); + + // Define momentum with respect to proton direction + double phi = r1->Uniform(phi_min*TMath::DegToRad(),phi_max*TMath::DegToRad()); + double th = r1->Uniform(th_min*TMath::DegToRad(),th_max*TMath::DegToRad()); + + //Total momentum distribution + double pevent = -1; + if(dist==0){ //fixed + pevent = p; + } + else if(dist==1){ //Uniform: +-50% variation + pevent = p*(1. + r1->Uniform(-0.5,0.5) ); + } + else if(dist==2){ //Gaussian: Sigma = 0.1*mean + while(pevent<0) //Avoid negative values + pevent = r1->Gaus(p,0.1*p); + } + + double px = pevent * std::cos(phi) * std::sin(th); + double py = pevent * std::sin(phi) * std::sin(th); + double pz = pevent * std::cos(th); + TVector3 pvec(px,py,pz); + + //Rotate to lab coordinate system + double cross_angle = -25./1000.; //in Rad + TVector3 pbeam_dir(sin(cross_angle),0,cos(cross_angle)); //proton beam direction + pvec.RotateY(-pbeam_dir.Theta()); // Theta is returned positive, beam in negative X + // type 1 is final state + // pdgid 11 - electron 0.510 MeV/c^2 + GenParticlePtr p3 = std::make_shared( + FourVector( + pvec.X(), pvec.Y(), pvec.Z(), + sqrt(pevent*pevent + (mass * mass))), + pdgID, 1); + + //If wanted, set non-zero vertex + double vx = 0.; + double vy = 0.; + double vz = 0.; + double vt = 0.; + + GenVertexPtr v1 = std::make_shared(); + evt.shift_position_by(FourVector(vx, vy, vz, vt)); + v1->add_particle_in(p1); + v1->add_particle_in(p2); + + v1->add_particle_out(p3); + evt.add_vertex(v1); + + if (events_parsed == 0) { + std::cout << "First event: " << std::endl; + Print::listing(evt); + } + + hepmc_output.write_event(evt); + if (events_parsed % 100 == 0) { + std::cout << "Event: " << events_parsed << std::endl; + } + evt.clear(); + } + hepmc_output.close(); + std::cout << "Events parsed and written: " << events_parsed << std::endl; +} From 217f6a93da81c7e019625df5d59f082cf6b4feec Mon Sep 17 00:00:00 2001 From: steinber Date: Wed, 21 Aug 2024 12:38:44 -0400 Subject: [PATCH 02/30] added LFHCAL benchmarks --- Snakefile | 1 + 1 file changed, 1 insertion(+) diff --git a/Snakefile b/Snakefile index 26cfc5d0..464a981a 100644 --- a/Snakefile +++ b/Snakefile @@ -4,6 +4,7 @@ include: "benchmarks/ecal_gaps/Snakefile" include: "benchmarks/material_scan/Snakefile" include: "benchmarks/tracking_performances/Snakefile" include: "benchmarks/zdc_lyso/Snakefile" +include: "benchmarks/lfhcal/Snakefile" rule fetch_epic: output: From 87f5a6e4cb2d71006e3b441ae60a2b2ff42bd9b8 Mon Sep 17 00:00:00 2001 From: steinber Date: Wed, 21 Aug 2024 13:45:01 -0400 Subject: [PATCH 03/30] tried to resolve rules, inputs and outputs relative to zdc_lyso --- benchmarks/lfhcal/Snakefile | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/benchmarks/lfhcal/Snakefile b/benchmarks/lfhcal/Snakefile index d0204ce3..0e176c48 100644 --- a/benchmarks/lfhcal/Snakefile +++ b/benchmarks/lfhcal/Snakefile @@ -1,13 +1,12 @@ import os - rule lfhcal_sim_hepmc: input: script = "benchmarks/lfhcal/gen_particles.cxx", output: - hepmcfile="data/{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc", + hepmcfile="data/lfhcal_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc", log: - "data/{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc.log", + "data/lfhcal_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc.log", params: num_events=1000, shell: @@ -18,12 +17,12 @@ root -l -b -q '{input.script}({params.num_events}, "{output.hepmcfile}", "{wildc rule lfhcal_sim: input: - hepmcfile="data/{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc", + hepmcfile="data/lfhcal_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc", warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root", output: - "sim_output/lfhcal/{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.edm4hep.root", + "sim_output/lfhcal/lfhcal_{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.edm4hep.root", log: - "sim_output/lfhcal/{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.edm4hep.root.log", + "sim_output/lfhcal/lfhcal_{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.edm4hep.root.log", params: num_events=1000, shell: @@ -40,11 +39,11 @@ npsim \ rule lfhcal_reco: input: - "sim_output/lfhcal/{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.edm4hep.root", + "sim_output/lfhcal/lfhcal_{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.edm4hep.root", output: - "sim_output/lfhcal/{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.eicrecon.tree.edm4eic.root", + "sim_output/lfhcal/lfhcal_{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.eicrecon.tree.edm4eic.root", log: - "sim_output/lfhcal/{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.eicrecon.tree.edm4eic.root.log", + "sim_output/lfhcal/lfhcal_{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.eicrecon.tree.edm4eic.root.log", shell: """ eicrecon -Ppodio:output_collections=MCParticles,ReconstructedParticles,LFHCALTruthClusters,LFHCALClusters,LFHCALHits,EcalEndcapPTruthClusters,EcalEndcapPClusters,EcalEnd\ @@ -55,14 +54,14 @@ mv podio_output.root {output} rule lfhcal_analysis: input: - expand("sim_output/lfhcal/{{DETECTOR_CONFIG}}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.eicrecon.tree.edm4eic.root", + expand("sim_output/lfhcal/lfhcal_{{DETECTOR_CONFIG}}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.eicrecon.tree.edm4eic.root", PARTICLE=["neutron","gamma"], BEAM_ENERGY=["2", "5", "10", "20", "40"], THETA_MIN=["0"], THETA_MAX=["50"]), script="benchmarks/lfhcal/analysis/analysis.py", output: - "results/{DETECTOR_CONFIG}/lfhcal/plots.pdf", + "results/{DETECTOR_CONFIG}/lfhcal/lfhcal_plots.pdf", shell: """ python {input.script} @@ -72,7 +71,7 @@ python {input.script} # Examples of invocation rule lfhcal_create_all_hepmc: input: - expand("data/{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc", + expand("data/lfhcal_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc", PARTICLE=["gamma"], BEAM_ENERGY=["0.005", "0.01", "0.05", "0.1", "0.5", "1.0"], THETA_MIN=["0"], @@ -81,7 +80,7 @@ rule lfhcal_create_all_hepmc: rule lfhcal_run_all_locally: input: - "results/" + os.environ["DETECTOR_CONFIG"] + "/lfhcal/plots.pdf", + "results/" + os.environ["DETECTOR_CONFIG"] + "/lfhcal/lfhcal_plots.pdf", message: "See output in {input[0]}" From 9434d7a10e21b87a289752d6890534fca629cedb Mon Sep 17 00:00:00 2001 From: steinber Date: Wed, 21 Aug 2024 13:55:40 -0400 Subject: [PATCH 04/30] resolved zdc_lyso ambiguitiy by changing output directory name, as changing filename seemed to just confuse the zdc_lyso snakemake... --- benchmarks/lfhcal/Snakefile | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/benchmarks/lfhcal/Snakefile b/benchmarks/lfhcal/Snakefile index 0e176c48..75bc9015 100644 --- a/benchmarks/lfhcal/Snakefile +++ b/benchmarks/lfhcal/Snakefile @@ -4,9 +4,9 @@ rule lfhcal_sim_hepmc: input: script = "benchmarks/lfhcal/gen_particles.cxx", output: - hepmcfile="data/lfhcal_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc", + hepmcfile="data_lfhcal/{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc", log: - "data/lfhcal_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc.log", + "data_lfhcal/lfhcal_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc.log", params: num_events=1000, shell: @@ -17,7 +17,7 @@ root -l -b -q '{input.script}({params.num_events}, "{output.hepmcfile}", "{wildc rule lfhcal_sim: input: - hepmcfile="data/lfhcal_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc", + hepmcfile="data_lfhcal/{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc", warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root", output: "sim_output/lfhcal/lfhcal_{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.edm4hep.root", @@ -72,7 +72,7 @@ python {input.script} rule lfhcal_create_all_hepmc: input: expand("data/lfhcal_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc", - PARTICLE=["gamma"], + PARTICLE=["neutron","gamma"], BEAM_ENERGY=["0.005", "0.01", "0.05", "0.1", "0.5", "1.0"], THETA_MIN=["0"], THETA_MAX=["0.3"]) From 7bc36351e18699618291f647c7c7aba11e5713f0 Mon Sep 17 00:00:00 2001 From: Peter Steinberg Date: Tue, 24 Sep 2024 11:30:05 -0400 Subject: [PATCH 05/30] Update benchmarks/lfhcal/config.yml Using the proper "namespaced" call. Co-authored-by: Dmitry Kalinkin --- benchmarks/lfhcal/config.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/lfhcal/config.yml b/benchmarks/lfhcal/config.yml index 5f589a26..226d2a4a 100644 --- a/benchmarks/lfhcal/config.yml +++ b/benchmarks/lfhcal/config.yml @@ -2,7 +2,7 @@ sim:lfhcal: extends: .det_benchmark stage: simulate script: - - snakemake --cores 1 run_all_locally + - snakemake --cores 1 lfhcal_run_all_locally retry: max: 2 when: From b5470523d5bcb2be79a174315cb27149800db0f6 Mon Sep 17 00:00:00 2001 From: steinber Date: Wed, 25 Sep 2024 13:39:05 -0400 Subject: [PATCH 06/30] adapting the tracking_performances approach, using the standard steering --- benchmarks/lfhcal/LFHCAL_Performance.C | 113 ++++++++++++++++++ benchmarks/lfhcal/Snakefile | 151 ++++++++++++++++--------- 2 files changed, 212 insertions(+), 52 deletions(-) create mode 100644 benchmarks/lfhcal/LFHCAL_Performance.C diff --git a/benchmarks/lfhcal/LFHCAL_Performance.C b/benchmarks/lfhcal/LFHCAL_Performance.C new file mode 100644 index 00000000..563506ae --- /dev/null +++ b/benchmarks/lfhcal/LFHCAL_Performance.C @@ -0,0 +1,113 @@ +// Code to extract the Tracking Performances +// Shyam Kumar; INFN Bari, Italy +// shyam.kumar@ba.infn.it; shyam.kumar@cern.ch + +#include "TGraphErrors.h" +#include "TF1.h" +#include "TRandom.h" +#include "TCanvas.h" +#include "TLegend.h" +#include "TMath.h" +#define mpi 0.139 // 1.864 GeV/c^2 + +void Tracking_Performances(TString filename="tracking_output",TString particle="pi-", double mom=0.1, Double_t pTcut = 0.15, TString name = "") +{ + + // style of the plot + gStyle->SetPalette(1); + gStyle->SetOptTitle(1); + gStyle->SetTitleOffset(.85,"X");gStyle->SetTitleOffset(.85,"Y"); + gStyle->SetTitleSize(.05,"X");gStyle->SetTitleSize(.05,"Y"); + gStyle->SetLabelSize(.04,"X");gStyle->SetLabelSize(.04,"Y"); + gStyle->SetHistLineWidth(2); + gStyle->SetOptFit(1); + gStyle->SetOptStat(1); + + TString dir = ""; + TString dist_dir_mom = ""; TString dist_dir_dca = ""; + if (name=="") {dist_dir_mom = "mom_resol_truth"; dist_dir_dca = "dca_resol_truth"; dir = "truthseed";} + else {dist_dir_mom = "mom_resol_realseed"; dist_dir_dca = "dca_resol_realseed"; dir = "realseed";} + + bool debug=true; + // Tree with reconstructed tracks + const int nbins_eta = 5; + int theta_val[nbins_eta+1] ={3,50,45,135,130,177}; + int nfiles = 100; + double eta[nbins_eta+1]={-3.5,-2.5,-1.0,1.0,2.5,3.5}; + double pt[nbins_eta+1]={0.5,1.0,2.0,5.0,10.0,20.1}; + TH1D *histp[nbins_eta]; + + TH3D *h_d0xy_3d= new TH3D("h_d0xy_3d","Transverse Pointing Resolution",500,-0.1,0.1,70,-3.5,3.5,201,0.,20.1); + TH3D *h_d0z_3d= new TH3D("h_d0z_3d","Longitudinal Pointing Resolution",500,-0.1,0.1,70,-3.5,3.5,201,0.,20.1); + + for (int i=0; iSetTitle(Form("%1.1f < #eta < %1.1f && p = %1.1f ",eta[i],eta[i+1],mom)); + histp[i]->SetName(Form("hist_mom_%1.1f_%1.1f_pmax_%1.1f",mom,eta[i],eta[i+1])); + } + + TFile* file = TFile::Open(filename.Data()); + if (!file) {printf("file not found !!!"); return;} + TTreeReader myReader("events", file); // name of tree and file + if (debug) cout<<"Filename: "<GetName()<<"\t NEvents: "< charge(myReader, "MCParticles.charge"); + TTreeReaderArray vx_mc(myReader, "MCParticles.vertex.x"); + TTreeReaderArray vy_mc(myReader, "MCParticles.vertex.y"); + TTreeReaderArray vz_mc(myReader, "MCParticles.vertex.z"); + TTreeReaderArray px_mc(myReader, "MCParticles.momentum.x"); + TTreeReaderArray py_mc(myReader, "MCParticles.momentum.y"); + TTreeReaderArray pz_mc(myReader, "MCParticles.momentum.z"); + TTreeReaderArray status(myReader, "MCParticles.generatorStatus"); + TTreeReaderArray pdg(myReader, "MCParticles.PDG"); + TTreeReaderArray match_flag(myReader, Form("CentralCKF%sTrackParameters.type",name.Data())); + TTreeReaderArray d0xy(myReader, Form("CentralCKF%sTrackParameters.loc.a",name.Data())); + TTreeReaderArray d0z(myReader, Form("CentralCKF%sTrackParameters.loc.b",name.Data())); + TTreeReaderArray theta(myReader, Form("CentralCKF%sTrackParameters.theta",name.Data())); + TTreeReaderArray phi(myReader, Form("CentralCKF%sTrackParameters.phi",name.Data())); + TTreeReaderArray qoverp(myReader, Form("CentralCKF%sTrackParameters.qOverP",name.Data())); + + int count =0; + int matchId = 1; // Always matched track assigned the index 0 + while (myReader.Next()) + { + if (match_flag.GetSize()==0) continue; // Events with no reco tracks skip them + for (int i = 0; i < matchId; ++i){ + + for (int j = 0; j < pdg.GetSize(); ++j){ + + if (status[j] !=1 && pdg.GetSize()!=1) continue; + Double_t ptmc = sqrt(px_mc[j]*px_mc[j]+py_mc[j]*py_mc[j]); + + if (fabs(ptmc) < pTcut) continue; + + Double_t pmc = (1./charge[j])*sqrt(px_mc[j]*px_mc[j]+py_mc[j]*py_mc[j]+pz_mc[j]*pz_mc[j]); // 1./(q/p); similar to prec + Double_t prec = 1./qoverp[j]; + + Double_t pzrec = prec*TMath::Cos(theta[j]); Double_t pt_rec = sqrt(prec*prec-pzrec*pzrec); + Double_t pzmc = pz_mc[j]; + + Double_t etamc = -1.0*TMath::Log(TMath::Tan((TMath::ACos(pzmc/fabs(pmc)))/2)); + Double_t p_resol = (prec-pmc)/pmc; + + for (int ibin=0; ibineta[ibin] && etamcFill(p_resol); + } + h_d0xy_3d->Fill(d0xy[j]*0.1, etamc, ptmc); // cm + h_d0z_3d->Fill(d0z[j]*0.1, etamc, ptmc); // cm + } // Generated Tracks + } // Reco Tracks + + }// event loop ends + + TFile *fout_mom = new TFile(Form("%s/mom/lfhcal_mom_%1.1f_%s_%s.root",particle.Data(),mom,dist_dir_mom.Data(),particle.Data()),"recreate"); + fout_mom->cd(); + for (int ibin=0; ibinWrite(); + fout_mom->Close(); + +} + + + + diff --git a/benchmarks/lfhcal/Snakefile b/benchmarks/lfhcal/Snakefile index 75bc9015..eaa84f2e 100644 --- a/benchmarks/lfhcal/Snakefile +++ b/benchmarks/lfhcal/Snakefile @@ -1,86 +1,133 @@ -import os - -rule lfhcal_sim_hepmc: - input: - script = "benchmarks/lfhcal/gen_particles.cxx", - output: - hepmcfile="data_lfhcal/{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc", - log: - "data_lfhcal/lfhcal_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc.log", - params: - num_events=1000, - shell: - """ -root -l -b -q '{input.script}({params.num_events}, "{output.hepmcfile}", "{wildcards.PARTICLE}", {wildcards.THETA_MIN}, {wildcards.THETA_MAX}, 0, 360, {wildcards.BEAM_ENERGY})' -""" - - rule lfhcal_sim: input: - hepmcfile="data_lfhcal/{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc", + steering_file="EPIC/EVGEN/SINGLE/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.steer", warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root", output: - "sim_output/lfhcal/lfhcal_{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.edm4hep.root", + "sim_output/lfhcal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.edm4hep.root", log: - "sim_output/lfhcal/lfhcal_{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.edm4hep.root.log", + "sim_output/lfhcal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.edm4hep.root.log", + wildcard_constraints: + PARTICLE="pi-", + ENERGY="[0-9]+[kMG]eV", + PHASE_SPACE="3to50deg", + INDEX="\d{4}", params: - num_events=1000, + N_EVENTS=10000 shell: """ -npsim \ +ddsim \ --runType batch \ + --enableGun \ + --steeringFile "{input.steering_file}" \ + --random.seed 1{wildcards.INDEX} \ + --filter.tracker edep0 \ -v WARNING \ + --numberOfEvents {params.N_EVENTS} \ --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \ - --numberOfEvents {params.num_events} \ - --inputFiles {input.hepmcfile} \ --outputFile {output} """ -rule lfhcal_reco: +rule lfhcal_recon: input: - "sim_output/lfhcal/lfhcal_{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.edm4hep.root", + "sim_output/lfhcal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.edm4hep.root", output: - "sim_output/lfhcal/lfhcal_{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.eicrecon.tree.edm4eic.root", + "sim_output/lfhcal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.eicrecon.tree.edm4eic.root", log: - "sim_output/lfhcal/lfhcal_{DETECTOR_CONFIG}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.eicrecon.tree.edm4eic.root.log", + "sim_output/lfhcal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.eicrecon.tree.edm4eic.root.log", + wildcard_constraints: + INDEX="\d{4}", + shell: """ +env DETECTOR_CONFIG={wildcards.DETECTOR_CONFIG} \ + eicrecon {input} -Ppodio:output_file={output} \ + -Ppodio:output_include_collections=MCParticles,ReconstructedParticles,LFHCALTruthClusters,LFHCALClusters,LFHCALHits,EcalEndcapPTruthClusters,EcalEndcapPClusters,EcalEndcapPHits +""" + + +rule lfhcal_at_momentum: + input: + script="benchmarks/lfhcal/LFHCAL_performance.C", + # TODO pass as a file list? + sim=lambda wildcards: + expand( + "sim_output/lfhcal/{DETECTOR_CONFIG}/{{PARTICLE}}/{ENERGY}/{PHASE_SPACE}/{{PARTICLE}}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root", + DETECTOR_CONFIG="epic_craterlake", + ENERGY=f"{float(wildcards.MOMENTUM):.0f}GeV" if float(wildcards.MOMENTUM) >= 1 else f"{float(wildcards.MOMENTUM) * 1000:.0f}MeV", + PHASE_SPACE=["3to50deg"], + INDEX=range(1), + ) + if wildcards.CAMPAIGN == "local" else + expand( + "EPIC/RECO/{{CAMPAIGN}}/epic_craterlake/SINGLE/{{PARTICLE}}/{ENERGY}/{PHASE_SPACE}/{{PARTICLE}}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root", + DETECTOR_CONFIG="epic_craterlake", + ENERGY=f"{float(wildcards.MOMENTUM):.0f}GeV" if float(wildcards.MOMENTUM) >= 1 else f"{float(wildcards.MOMENTUM) * 1000:.0f}MeV", + PHASE_SPACE=["3to50deg"], + INDEX=range(1), + ), + output: + "{CAMPAIGN}/pi-/mom/lfhcal_mom_{MOMENTUM}_mom_resol_{PARTICLE}.root", + combined_root="{CAMPAIGN}/pi-/mom/lfhcal_mom_{MOMENTUM}_mom_resol_{PARTICLE}.root", shell: """ -eicrecon -Ppodio:output_collections=MCParticles,ReconstructedParticles,LFHCALTruthClusters,LFHCALClusters,LFHCALHits,EcalEndcapPTruthClusters,EcalEndcapPClusters,EcalEnd\ -capPHits {input} -mv podio_output.root {output} +hadd {output} {input.sim} +cd {wildcards.CAMPAIGN} +root -l -b -q ../{input.script}'("../{output.combined_root}", "{wildcards.PARTICLE}", {wildcards.MOMENTUM}, 0.15)' """ -rule lfhcal_analysis: + +rule lfhcal_debug_montage: input: - expand("sim_output/lfhcal/lfhcal_{{DETECTOR_CONFIG}}_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.eicrecon.tree.edm4eic.root", - PARTICLE=["neutron","gamma"], - BEAM_ENERGY=["2", "5", "10", "20", "40"], - THETA_MIN=["0"], - THETA_MAX=["50"]), - script="benchmarks/lfhcal/analysis/analysis.py", + expand( + [ + "{{CAMPAIGN}}/Debug_Plots/pi-/mom/lfhcal_mom_resol_mom{MOMENTUM:.1f}_{{ETA_BIN}}.png", + ], + MOMENTUM=[0.5, 1.0, 2.0, 5.0, 10.0, 20.0], + ), output: - "results/{DETECTOR_CONFIG}/lfhcal/lfhcal_plots.pdf", + "{CAMPAIGN}/Debug_Plots/pi-/mom/lfhcal_mom_resol_debug_{ETA_BIN}.png", shell: """ -python {input.script} +montage -mode concatenate {input} {output} || true +ls {output} """ -# Examples of invocation -rule lfhcal_create_all_hepmc: +LFHCAL_ETA_BINS = [1.2,1.5,2.0,2.5,3.0,3.5] + +rule lfhcal: input: - expand("data/lfhcal_{PARTICLE}_{BEAM_ENERGY}GeV_theta_{THETA_MIN}deg_thru_{THETA_MAX}deg.hepmc", - PARTICLE=["neutron","gamma"], - BEAM_ENERGY=["0.005", "0.01", "0.05", "0.1", "0.5", "1.0"], - THETA_MIN=["0"], - THETA_MAX=["0.3"]) + expand( + [ + "{{CAMPAIGN}}/Final_Results/pi-/mom/mom_resol_{ETA_BIN}.png", + "{{CAMPAIGN}}/Final_Results/pi-/mom/mom_resol_{ETA_BIN}.root", + "{{CAMPAIGN}}/Debug_Plots/pi-/mom/mom_resol_debug_{ETA_BIN}.png", + ], + ETA_BIN=[f"{eta_min:.1f}_eta_{eta_max:.1f}" for eta_min, eta_max in zip(LFHCAL_ETA_BINS[:-1], LFHCAL_ETA_BINS[1:])], + ), + output: + directory("results/lfhcal/{CAMPAIGN}/") + shell: + """ +mkdir {output} +cp {input} {output} +""" -rule lfhcal_run_all_locally: - input: - "results/" + os.environ["DETECTOR_CONFIG"] + "/lfhcal/lfhcal_plots.pdf", - message: - "See output in {input[0]}" +rule lfhcal_local: + input: + "results/lfhcal/local", + +rule lfhcal_campaigns: + input: + expand( + "results/lfhcal/{CAMPAIGN}", + CAMPAIGN=[ + "23.10.0", + "23.11.0", + "23.12.0", + "24.03.1", + "24.04.0", + ], + ) From 8b6eac91e5a682469f2aca11f2f53eb947a1fbb3 Mon Sep 17 00:00:00 2001 From: steinber Date: Wed, 25 Sep 2024 13:41:58 -0400 Subject: [PATCH 07/30] moved to tracking_performances approach --- benchmarks/lfhcal/config.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/lfhcal/config.yml b/benchmarks/lfhcal/config.yml index 5f589a26..3da3cf52 100644 --- a/benchmarks/lfhcal/config.yml +++ b/benchmarks/lfhcal/config.yml @@ -2,7 +2,7 @@ sim:lfhcal: extends: .det_benchmark stage: simulate script: - - snakemake --cores 1 run_all_locally + - snakemake --cores 1 lfhcal retry: max: 2 when: From 0ac2c08c278d4b2b5065b60ca5a44cc1bb4d8682 Mon Sep 17 00:00:00 2001 From: steinber Date: Sat, 5 Oct 2024 09:53:33 -0400 Subject: [PATCH 08/30] Working version with pi- only --- benchmarks/lfhcal/Snakefile | 67 +++++++++++++++++++++++-------------- 1 file changed, 41 insertions(+), 26 deletions(-) diff --git a/benchmarks/lfhcal/Snakefile b/benchmarks/lfhcal/Snakefile index eaa84f2e..840ec52c 100644 --- a/benchmarks/lfhcal/Snakefile +++ b/benchmarks/lfhcal/Snakefile @@ -7,12 +7,12 @@ rule lfhcal_sim: log: "sim_output/lfhcal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.edm4hep.root.log", wildcard_constraints: - PARTICLE="pi-", ENERGY="[0-9]+[kMG]eV", + PARTICLE="pi-", PHASE_SPACE="3to50deg", INDEX="\d{4}", params: - N_EVENTS=10000 + N_EVENTS=100 shell: """ ddsim \ @@ -21,7 +21,7 @@ ddsim \ --steeringFile "{input.steering_file}" \ --random.seed 1{wildcards.INDEX} \ --filter.tracker edep0 \ - -v WARNING \ + -v INFO \ --numberOfEvents {params.N_EVENTS} \ --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \ --outputFile {output} @@ -43,68 +43,83 @@ env DETECTOR_CONFIG={wildcards.DETECTOR_CONFIG} \ -Ppodio:output_include_collections=MCParticles,ReconstructedParticles,LFHCALTruthClusters,LFHCALClusters,LFHCALHits,EcalEndcapPTruthClusters,EcalEndcapPClusters,EcalEndcapPHits """ - rule lfhcal_at_momentum: input: - script="benchmarks/lfhcal/LFHCAL_performance.C", + script="benchmarks/lfhcal/LFHCAL_Performance.C", # TODO pass as a file list? sim=lambda wildcards: expand( - "sim_output/lfhcal/{DETECTOR_CONFIG}/{{PARTICLE}}/{ENERGY}/{PHASE_SPACE}/{{PARTICLE}}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root", + "sim_output/lfhcal/{DETECTOR_CONFIG}/{{PARTICLE}}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root", DETECTOR_CONFIG="epic_craterlake", ENERGY=f"{float(wildcards.MOMENTUM):.0f}GeV" if float(wildcards.MOMENTUM) >= 1 else f"{float(wildcards.MOMENTUM) * 1000:.0f}MeV", PHASE_SPACE=["3to50deg"], + PARTICLE=["pi-"], INDEX=range(1), ) if wildcards.CAMPAIGN == "local" else - expand( - "EPIC/RECO/{{CAMPAIGN}}/epic_craterlake/SINGLE/{{PARTICLE}}/{ENERGY}/{PHASE_SPACE}/{{PARTICLE}}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root", + ancient(expand( + "EPIC/RECO/{{CAMPAIGN}}/epic_craterlake/SINGLE/{{PARTICLE}}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root", DETECTOR_CONFIG="epic_craterlake", ENERGY=f"{float(wildcards.MOMENTUM):.0f}GeV" if float(wildcards.MOMENTUM) >= 1 else f"{float(wildcards.MOMENTUM) * 1000:.0f}MeV", PHASE_SPACE=["3to50deg"], + PARTICLE=["pi-"], INDEX=range(1), - ), + )), output: - "{CAMPAIGN}/pi-/mom/lfhcal_mom_{MOMENTUM}_mom_resol_{PARTICLE}.root", - combined_root="{CAMPAIGN}/pi-/mom/lfhcal_mom_{MOMENTUM}_mom_resol_{PARTICLE}.root", + "{CAMPAIGN}/{PARTICLE}/mom/lfhcal_mom_{MOMENTUM}_mom_resol_{PARTICLE}.root", + combined_root=temp("{CAMPAIGN}/lfhcal_sim_{MOMENTUM}_{PARTICLE}.root"), shell: """ -hadd {output} {input.sim} +hadd {output.combined_root} {input.sim} cd {wildcards.CAMPAIGN} root -l -b -q ../{input.script}'("../{output.combined_root}", "{wildcards.PARTICLE}", {wildcards.MOMENTUM}, 0.15)' """ - - -rule lfhcal_debug_montage: +rule lfhcal_summary_at_eta: input: expand( [ - "{{CAMPAIGN}}/Debug_Plots/pi-/mom/lfhcal_mom_resol_mom{MOMENTUM:.1f}_{{ETA_BIN}}.png", + "{{CAMPAIGN}}/{PARTICLE}/mom/lfhcal_mom_{MOMENTUM:.1f}_mom_resol_{PARTICLE}.root", ], MOMENTUM=[0.5, 1.0, 2.0, 5.0, 10.0, 20.0], + PARTICLE=["pi-"] ), + script="benchmarks/lfhcal/doCompare_widebins_mom.C", output: - "{CAMPAIGN}/Debug_Plots/pi-/mom/lfhcal_mom_resol_debug_{ETA_BIN}.png", + "{CAMPAIGN}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{ETA_MIN}_eta_{ETA_MAX}.png", + "{CAMPAIGN}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{ETA_MIN}_eta_{ETA_MAX}.root", shell: """ -montage -mode concatenate {input} {output} || true -ls {output} -""" +if [[ "{wildcards.CAMPAIGN}" == "local" ]]; then + set +e + EPIC_VERSION="${{DETECTOR_VERSION:-}}" + EICRECON_VERSION="$(eicrecon -v | sed -n -e 's/.*\(v[0-9\.]\+\).*/\\1/p')" + # Legacy detection + : ${{EPIC_VERSION:="$(echo $DETECTOR_PATH | sed -n -e 's/.*epic-\([^-/]\+\).*/\\1/p')"}} + set -e + echo "ePIC version: $EPIC_VERSION" + echo "EICrecon version: $EICRECON_VERSION" + EXTRA_LEGEND="ePIC $EPIC_VERSION / EICrecon $EICRECON_VERSION" +else + EXTRA_LEGEND="ePIC Simulation {wildcards.CAMPAIGN}" +fi +cd {wildcards.CAMPAIGN} +root -l -b -q ../{input.script}'("pi-", {wildcards.ETA_MIN}, {wildcards.ETA_MAX}, 1., true, "'"$EXTRA_LEGEND"'")' +""" -LFHCAL_ETA_BINS = [1.2,1.5,2.0,2.5,3.0,3.5] +LFHCAL_ETA_BINS = [1.2,1.5,2,2.5,3,3.5] rule lfhcal: input: expand( [ - "{{CAMPAIGN}}/Final_Results/pi-/mom/mom_resol_{ETA_BIN}.png", - "{{CAMPAIGN}}/Final_Results/pi-/mom/mom_resol_{ETA_BIN}.root", - "{{CAMPAIGN}}/Debug_Plots/pi-/mom/mom_resol_debug_{ETA_BIN}.png", + "{{CAMPAIGN}}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{ETA_BIN}.png", + "{{CAMPAIGN}}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{ETA_BIN}.root", ], - ETA_BIN=[f"{eta_min:.1f}_eta_{eta_max:.1f}" for eta_min, eta_max in zip(LFHCAL_ETA_BINS[:-1], LFHCAL_ETA_BINS[1:])], - ), + ETA_BIN=[f"{eta_min:.1f}_eta_{eta_max:.1f}" for eta_min, eta_max in zip(LFHCAL_ETA_BINS[:-1], LFHCAL_ETA_BINS[1:])], + PARTICLE="pi-", + ) output: directory("results/lfhcal/{CAMPAIGN}/") shell: From 9c9802ddde150efab5d0f7ea031a8b7145bb9004 Mon Sep 17 00:00:00 2001 From: steinber Date: Sat, 5 Oct 2024 09:53:46 -0400 Subject: [PATCH 09/30] new per-sample analysis --- benchmarks/lfhcal/LFHCAL_Performance.C | 125 ++++++++++++++++--------- 1 file changed, 80 insertions(+), 45 deletions(-) diff --git a/benchmarks/lfhcal/LFHCAL_Performance.C b/benchmarks/lfhcal/LFHCAL_Performance.C index 563506ae..fc4dcbbb 100644 --- a/benchmarks/lfhcal/LFHCAL_Performance.C +++ b/benchmarks/lfhcal/LFHCAL_Performance.C @@ -8,9 +8,11 @@ #include "TCanvas.h" #include "TLegend.h" #include "TMath.h" +#include "TVector3.h" + #define mpi 0.139 // 1.864 GeV/c^2 -void Tracking_Performances(TString filename="tracking_output",TString particle="pi-", double mom=0.1, Double_t pTcut = 0.15, TString name = "") +void LFHCAL_Performance(TString filename="tracking_output",TString particle="pi-", double mom=0.1, Double_t pTcut = 0.0, TString name = "") { // style of the plot @@ -24,21 +26,16 @@ void Tracking_Performances(TString filename="tracking_output",TString particle=" gStyle->SetOptStat(1); TString dir = ""; - TString dist_dir_mom = ""; TString dist_dir_dca = ""; - if (name=="") {dist_dir_mom = "mom_resol_truth"; dist_dir_dca = "dca_resol_truth"; dir = "truthseed";} - else {dist_dir_mom = "mom_resol_realseed"; dist_dir_dca = "dca_resol_realseed"; dir = "realseed";} - + TString dist_dir_mom = "mom_resol"; + bool debug=true; // Tree with reconstructed tracks const int nbins_eta = 5; - int theta_val[nbins_eta+1] ={3,50,45,135,130,177}; int nfiles = 100; - double eta[nbins_eta+1]={-3.5,-2.5,-1.0,1.0,2.5,3.5}; + double eta[nbins_eta+1]={1.2,1.5,2,2.5,3,3.5}; double pt[nbins_eta+1]={0.5,1.0,2.0,5.0,10.0,20.1}; TH1D *histp[nbins_eta]; - TH3D *h_d0xy_3d= new TH3D("h_d0xy_3d","Transverse Pointing Resolution",500,-0.1,0.1,70,-3.5,3.5,201,0.,20.1); - TH3D *h_d0z_3d= new TH3D("h_d0z_3d","Longitudinal Pointing Resolution",500,-0.1,0.1,70,-3.5,3.5,201,0.,20.1); for (int i=0; i vx_mc(myReader, "MCParticles.vertex.x"); TTreeReaderArray vy_mc(myReader, "MCParticles.vertex.y"); - TTreeReaderArray vz_mc(myReader, "MCParticles.vertex.z"); + TTreeReaderArray vz_mc(myReader, "MCParticles.vertex.z"); TTreeReaderArray px_mc(myReader, "MCParticles.momentum.x"); TTreeReaderArray py_mc(myReader, "MCParticles.momentum.y"); - TTreeReaderArray pz_mc(myReader, "MCParticles.momentum.z"); + TTreeReaderArray pz_mc(myReader, "MCParticles.momentum.z"); TTreeReaderArray status(myReader, "MCParticles.generatorStatus"); - TTreeReaderArray pdg(myReader, "MCParticles.PDG"); - TTreeReaderArray match_flag(myReader, Form("CentralCKF%sTrackParameters.type",name.Data())); - TTreeReaderArray d0xy(myReader, Form("CentralCKF%sTrackParameters.loc.a",name.Data())); - TTreeReaderArray d0z(myReader, Form("CentralCKF%sTrackParameters.loc.b",name.Data())); - TTreeReaderArray theta(myReader, Form("CentralCKF%sTrackParameters.theta",name.Data())); - TTreeReaderArray phi(myReader, Form("CentralCKF%sTrackParameters.phi",name.Data())); - TTreeReaderArray qoverp(myReader, Form("CentralCKF%sTrackParameters.qOverP",name.Data())); + TTreeReaderArray pdg(myReader, "MCParticles.PDG"); + + TTreeReaderArray pe_lc(myReader, "LFHCALClusters.energy"); + TTreeReaderArray px_lc(myReader, "LFHCALClusters.position.x"); + TTreeReaderArray py_lc(myReader, "LFHCALClusters.position.y"); + TTreeReaderArray pz_lc(myReader, "LFHCALClusters.position.z"); + TTreeReaderArray pe_ec(myReader, "EcalEndcapPClusters.energy"); + TTreeReaderArray px_ec(myReader, "EcalEndcapPClusters.position.x"); + TTreeReaderArray py_ec(myReader, "EcalEndcapPClusters.position.y"); + TTreeReaderArray pz_ec(myReader, "EcalEndcapPClusters.position.z"); int count =0; int matchId = 1; // Always matched track assigned the index 0 while (myReader.Next()) - { - if (match_flag.GetSize()==0) continue; // Events with no reco tracks skip them - for (int i = 0; i < matchId; ++i){ - - for (int j = 0; j < pdg.GetSize(); ++j){ - - if (status[j] !=1 && pdg.GetSize()!=1) continue; - Double_t ptmc = sqrt(px_mc[j]*px_mc[j]+py_mc[j]*py_mc[j]); - - if (fabs(ptmc) < pTcut) continue; + { + for (int j = 0; j < pdg.GetSize(); ++j) + { + if (status[j] !=1 && pdg.GetSize()!=1) continue; + Double_t pzmc = pz_mc[j]; + Double_t ptmc = sqrt(px_mc[j]*px_mc[j]+py_mc[j]*py_mc[j]); + Double_t pmc = sqrt(px_mc[j]*px_mc[j]+py_mc[j]*py_mc[j]+pz_mc[j]*pz_mc[j]); // 1./(q/p); similar to prec + Double_t etamc = -1.0*TMath::Log(TMath::Tan((TMath::ACos(pzmc/fabs(pmc)))/2)); + Double_t phimc = TMath::ATan2(py_mc[j],px_mc[j]); + if (fabs(ptmc) < pTcut) continue; - Double_t pmc = (1./charge[j])*sqrt(px_mc[j]*px_mc[j]+py_mc[j]*py_mc[j]+pz_mc[j]*pz_mc[j]); // 1./(q/p); similar to prec - Double_t prec = 1./qoverp[j]; + Double_t l_px_tot; + Double_t l_py_tot; + Double_t l_pz_tot; + Double_t l_e_tot; - Double_t pzrec = prec*TMath::Cos(theta[j]); Double_t pt_rec = sqrt(prec*prec-pzrec*pzrec); - Double_t pzmc = pz_mc[j]; - - Double_t etamc = -1.0*TMath::Log(TMath::Tan((TMath::ACos(pzmc/fabs(pmc)))/2)); - Double_t p_resol = (prec-pmc)/pmc; - - for (int ibin=0; ibineta[ibin] && etamcFill(p_resol); - } - h_d0xy_3d->Fill(d0xy[j]*0.1, etamc, ptmc); // cm - h_d0z_3d->Fill(d0z[j]*0.1, etamc, ptmc); // cm - } // Generated Tracks - } // Reco Tracks - - }// event loop ends + for (int jl = 0;jleta[ibin] && etamcFill(p_resol); + } + } // Generated Tracks + + }// event loop ends TFile *fout_mom = new TFile(Form("%s/mom/lfhcal_mom_%1.1f_%s_%s.root",particle.Data(),mom,dist_dir_mom.Data(),particle.Data()),"recreate"); fout_mom->cd(); From 14dfb3fd1266f144358623d0b0b9b3a70c92dbad Mon Sep 17 00:00:00 2001 From: steinber Date: Sat, 5 Oct 2024 09:53:57 -0400 Subject: [PATCH 10/30] new integration step --- benchmarks/lfhcal/doCompare_widebins_mom.C | 173 +++++++++++++++++++++ 1 file changed, 173 insertions(+) create mode 100644 benchmarks/lfhcal/doCompare_widebins_mom.C diff --git a/benchmarks/lfhcal/doCompare_widebins_mom.C b/benchmarks/lfhcal/doCompare_widebins_mom.C new file mode 100644 index 00000000..41285e03 --- /dev/null +++ b/benchmarks/lfhcal/doCompare_widebins_mom.C @@ -0,0 +1,173 @@ +// Code to compare the tracking performances: Truth seeding vs real seeding +// Shyam Kumar; shyam.kumar@ba.infn.it; shyam055119@gmail.com + +#include "TGraphErrors.h" +#include "TF1.h" +#include "TRandom.h" +#include "TCanvas.h" +#include "TLegend.h" +#include "TMath.h" +#define mpi 0.139 // 1.864 GeV/c^2 + +//void draw_req_Mom(double etamin, double etamax, double xmin=0., double xmax=0.); +void doCompare_widebins_mom(TString particle = "pi-",double etamin=-1.0, double etamax=1.0, double range =0.3, Bool_t drawreq=1, TString extra_legend = "") // name = p, pt for getting p or pt dependence fitted results +{ + + //=== style of the plot========= + gStyle->SetPalette(1); + gStyle->SetOptTitle(1); + gStyle->SetTitleOffset(1.0,"XY"); + gStyle->SetTitleSize(.04,"XY"); + gStyle->SetLabelSize(.04,"XY"); + gStyle->SetHistLineWidth(2); + gStyle->SetOptFit(1); + gStyle->SetOptStat(1); + + const Int_t nfiles = 6; + double mom[nfiles] ={0.5,1.0,2.0,5.0,10.0,20.0}; + std::vector momV, momresolV, err_momresolV; + momV.clear(); momresolV.clear(); err_momresolV.clear(); + TString symbolname = ""; + if (particle == "pi-") symbolname = "#pi^{-}"; + else symbolname = particle; + ofstream outfile; + outfile.open ("Mom_resol.txt",ios_base::app); + + TF1 *f1=new TF1("f1","FitMomentumResolution",0.,30.0,2); + f1->SetParLimits(0,0.,0.1); + f1->SetParLimits(1,0.,5.0); + + TCanvas *c_mom = new TCanvas("cmom","cmom",1400,1000); + c_mom->SetMargin(0.10, 0.05 ,0.1,0.05); + c_mom->SetGridy(); + + //Reading the root file + TFile *fmom[nfiles]; + TGraphErrors *gr_mom; + TMultiGraph *mgMom; + TLegend *lmom; + mgMom = new TMultiGraph("mgMom",";p (GeV/c); #sigmap/p %"); + + lmom = new TLegend(0.65,0.80,0.90,0.93); + lmom->SetTextSize(0.03); + lmom->SetBorderSize(0); + lmom->SetHeader(extra_legend.Data(), "C"); + lmom->AddEntry((TObject*)0, Form("%s, %1.1f < #eta < %1.1f", symbolname.Data(), etamin, etamax), "C"); + + TF1 *func = new TF1("func","gaus",-0.5,0.5); + + for (int i =0; iSetMargin(0.10, 0.05 ,0.1,0.07); + + //pi-/mom/lfhcal_mom_20.0_mom_resol_pi-.root + fmom[i] = TFile::Open(Form("./pi-/mom/lfhcal_mom_%1.1f_mom_resol_%s.root",mom[i],particle.Data())); + + TH1D *hist = (TH1D*) fmom[i]->Get(Form("hist_mom_%1.1f_%1.1f_pmax_%1.1f",mom[i],etamin,etamax)); + hist->Rebin(2); + hist->SetName(Form("hist_mom_%1.1f_%1.1f_eta_%1.1f_%s",mom[i],etamin,etamax,particle.Data())); + hist->SetTitle(Form("Momentum = %1.1f && %1.1f<#eta<%1.1f;#Delta p/p; Entries(a.u.)",particle.Data(),mom[i],etamin,etamax)); + + double mu = hist->GetMean(); + double sigma = hist->GetStdDev(); + hist->GetXaxis()->SetRangeUser(-1.0*range,1.0*range); + func->SetRange(mu-2.0*sigma,mu+2.0*sigma); // fit with in 2 sigma range + hist->Fit(func,"NR+"); + mu = func->GetParameter(1); + sigma = func->GetParameter(2); + func->SetRange(mu-2.0*sigma,mu+2.0*sigma); + hist->Fit(func,"R+"); + float par2 = func->GetParameter(2)*100; + float par2_err = func->GetParError(2)*100; + momV.push_back(mom[i]); + momresolV.push_back(par2); + err_momresolV.push_back(par2_err); + + cp->cd(); + hist->Draw(); + //cp->SaveAs(Form("Debug_Plots/%s/mom/mom_resol_mom%1.1f_%1.1f_eta_%1.1f.png",particle.Data(),mom[i],etamin,etamax)); + } // all files + + const int size = momV.size(); + double p[size], err_p[size], sigma_p[size], err_sigma_p[size]; + + for (int i=0; iSetName("grseed"); + gr1->SetMarkerStyle(25); + gr1->SetMarkerColor(kBlue); + gr1->SetMarkerSize(2.0); + gr1->SetTitle(";p (GeV/c);#sigmap/p"); + gr1->GetXaxis()->CenterTitle(); + gr1->GetYaxis()->CenterTitle(); + + + mgMom->Add(gr1); + c_mom->cd(); + mgMom->GetXaxis()->SetRangeUser(0.40,20.2); + mgMom->GetYaxis()->SetRangeUser(0.0,1.50*TMath::MaxElement(gr1->GetN(),gr1->GetY())); // 50% more of the maximum value on yaxis + mgMom->Draw("AP"); + lmom->AddEntry(gr1,"Nominal"); + lmom->Draw("same"); + //draw_req_Mom(etamin,etamax,0.,mgMom->GetXaxis()->GetXmax()); + c_mom->SaveAs(Form("Final_Results/%s/mom/lfhcal_mom_resol_%1.1f_eta_%1.1f.png",particle.Data(),etamin,etamax)); + + // Write the numbers in output file for comparisons + outfile << extra_legend << endl; + outfile<<"Etamin"<GetN(); ++i){ + double x,y; + gr1->GetPoint(i,x,y); + outfile<cd(); + mgMom->SetName(Form("mom_resol_%1.1f_eta_%1.1f",etamin,etamax)); + mgMom->Write(); + fout->Close(); +} + +//===Fit Momentum Resolution +float FitMomentumResolution(Double_t *x, Double_t *par) +{ + float func = sqrt(par[0]*par[0]*x[0]*x[0]+par[1]*par[1]); + return func; +} + +//From Yellow report from section 11.2.2 + +// 1.2,1.5,2,2.5,3,3.5 + +/* +void draw_req_Mom(double etamin, double etamax, double xmin=0., double xmax=0.) +{ + + TF1 *dd4hep_p; + if (etamin >= -3.5 && etamax <= -2.5) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.1*x)^2+2.0^2)",xmin,xmax); + else if (etamin >= -2.5 && etamax <= -1.0) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.05*x)^2+1.0^2)",xmin,xmax); + else if (etamin >= -1.0 && etamax <= 1.0) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.05*x)^2+0.5^2)",xmin,xmax); + else if (etamin >= 1.0 && etamax <= 2.5) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.05*x)^2+1.0^2)",xmin,xmax); + else if (etamin >= 2.5 && etamax <= 3.5) dd4hep_p = new TF1("dd4hep_p", "TMath::Sqrt((0.1*x)^2+2.0^2)",xmin,xmax); + else return; + dd4hep_p->SetLineStyle(7); + dd4hep_p->SetLineColor(kMagenta); + dd4hep_p->SetLineWidth(3.0); + dd4hep_p->Draw("same"); + + TLegend *l= new TLegend(0.70,0.75,0.90,0.80); + l->SetTextSize(0.03); + l->SetBorderSize(0); + l->AddEntry(dd4hep_p,"PWGReq","l"); + l->Draw("same"); + } +*/ From 79df73ac8680aaf886fa8085095aea175f15a3b4 Mon Sep 17 00:00:00 2001 From: steinber Date: Sat, 5 Oct 2024 09:55:59 -0400 Subject: [PATCH 11/30] clean up of unused files --- benchmarks/lfhcal/analysis/analysis.py | 281 ------------------------- benchmarks/lfhcal/gen_particles.cxx | 126 ----------- 2 files changed, 407 deletions(-) delete mode 100644 benchmarks/lfhcal/analysis/analysis.py delete mode 100644 benchmarks/lfhcal/gen_particles.cxx diff --git a/benchmarks/lfhcal/analysis/analysis.py b/benchmarks/lfhcal/analysis/analysis.py deleted file mode 100644 index 1ef1f828..00000000 --- a/benchmarks/lfhcal/analysis/analysis.py +++ /dev/null @@ -1,281 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -import mplhep as hep -import uproot -import pandas as pd -from scipy.optimize import curve_fit -from matplotlib.backends.backend_pdf import PdfPages -import os -import awkward as ak - -plt.figure() -hep.set_style(hep.style.CMS) -hep.set_style("CMS") - -def gaussian(x, amp, mean, sigma): - return amp * np.exp( -(x - mean)**2 / (2*sigma**2) ) - -def rotateY(xdata, zdata, angle): - s = np.sin(angle) - c = np.cos(angle) - rotatedz = c*zdata - s*xdata - rotatedx = s*zdata + c*xdata - return rotatedx, rotatedz - -Energy = [0.005, 0.01, 0.05, 0.1, 0.5, 1.0] - - -df = pd.DataFrame({}) -for eng in Energy: - tree = uproot.open(f'sim_output/lfhcal/{os.environ["DETECTOR_CONFIG"]}_gamma_{eng}GeV_theta_0deg_thru_0.3deg.eicrecon.tree.edm4eic.root')['events'] - ecal_reco_energy = list(map(sum, tree['EcalFarForwardZDCClusters/EcalFarForwardZDCClusters.energy'].array())) - hcal_reco_energy = list(map(sum, tree['HcalFarForwardZDCClusters/HcalFarForwardZDCClusters.energy'].array())) - ecal_rec_energy = list(map(sum, tree['EcalFarForwardZDCRecHits/EcalFarForwardZDCRecHits.energy'].array())) - hcal_rec_energy = list(map(sum, tree['HcalFarForwardZDCRecHits/HcalFarForwardZDCRecHits.energy'].array())) - ecal_reco_clusters = [len(row) if len(row)>=1 else 0 for row in tree['EcalFarForwardZDCClusters/EcalFarForwardZDCClusters.nhits'].array()] - ecal_reco_nhits = [row[0] if len(row)>=1 else 0 for row in tree['EcalFarForwardZDCClusters/EcalFarForwardZDCClusters.nhits'].array()] - - tree = uproot.open(f'sim_output/lfhcal/{os.environ["DETECTOR_CONFIG"]}_gamma_{eng}GeV_theta_0deg_thru_0.3deg.edm4hep.root')['events'] - ecal_sim_energy = list(map(sum, tree['EcalFarForwardZDCHits/EcalFarForwardZDCHits.energy'].array())) - hcal_sim_energy = list(map(sum, tree['HcalFarForwardZDCHits/HcalFarForwardZDCHits.energy'].array())) - - par_x = tree['MCParticles/MCParticles.momentum.x'].array()[:,2] - par_y = tree['MCParticles/MCParticles.momentum.y'].array()[:,2] - par_z = tree['MCParticles/MCParticles.momentum.z'].array()[:,2] - - eng = int(eng*1000) - - ecal_reco_energy = pd.DataFrame({f'ecal_reco_energy_{eng}': np.array(ecal_reco_energy, dtype=object)}) - hcal_reco_energy = pd.DataFrame({f'hcal_reco_energy_{eng}': np.array(hcal_reco_energy, dtype=object)}) - ecal_rec_energy = pd.DataFrame({f'ecal_rec_energy_{eng}': np.array(ecal_rec_energy, dtype=object)}) - hcal_rec_energy = pd.DataFrame({f'hcal_rec_energy_{eng}': np.array(hcal_rec_energy, dtype=object)}) - ecal_sim_energy = pd.DataFrame({f'ecal_sim_energy_{eng}': np.array(ecal_sim_energy, dtype=object)}) - hcal_sim_energy = pd.DataFrame({f'hcal_sim_energy_{eng}': np.array(hcal_sim_energy, dtype=object)}) - ecal_reco_nhits = pd.DataFrame({f'ecal_reco_nhits_{eng}': np.array(ecal_reco_nhits, dtype=object)}) - ecal_reco_clusters = pd.DataFrame({f'ecal_reco_clusters_{eng}': np.array(ecal_reco_clusters, dtype=object)}) - par_x = pd.DataFrame({f'par_x_{eng}': np.array(par_x.tolist(), dtype=object)}) - par_y = pd.DataFrame({f'par_y_{eng}': np.array(par_y.tolist(), dtype=object)}) - par_z = pd.DataFrame({f'par_z_{eng}': np.array(par_z.tolist(), dtype=object)}) - - - 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) - - -mu = [] -sigma = [] -fig1, ax = plt.subplots(3,2,figsize=(20,10)) -fig1.suptitle('ZDC ECal Cluster Energy Reconstruction') - -plt.tight_layout() -for i in range(6): - x = df[f'par_x_{eng}'].astype(float).to_numpy() - y = df[f'par_y_{eng}'].astype(float).to_numpy() - z = df[f'par_z_{eng}'].astype(float).to_numpy() - x, z = rotateY(x,z, 0.025) - theta = np.arccos(z/np.sqrt((x**2+y**2+z**2)))*1000 - condition = theta <= 3.5 - - plt.sca(ax[i%3,i//3]) - eng = int(Energy[i]*1000) - plt.title(f'Gamma Energy: {eng} MeV') - temp = np.array(df[f'ecal_reco_energy_{eng}'].astype(float).to_numpy()[condition])*1000 - hist, x = np.histogram(temp,bins=np.linspace(min(temp),max(temp)+np.std(abs(temp)),2*int(np.sqrt(len(temp))))) - x = x[1:]/2 + x[:-1]/2 - plt.errorbar(x,hist,yerr=np.sqrt(hist),fmt='-o',label='Cluster') - 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]))) - #plt.plot(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),gaussian(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),*coeff)) - mu.append(coeff[1]) - sigma.append(coeff[2]) - - temp = np.array(df[f'ecal_rec_energy_{eng}'].astype(float).to_numpy()[condition])*1000 - hist, x = np.histogram(temp,bins=np.linspace(min(temp),max(temp)+np.std(abs(temp)),2*int(np.sqrt(len(temp))))) - x = x[1:]/2 + x[:-1]/2 - plt.errorbar(x,hist,yerr=np.sqrt(hist),fmt='-o',label='Digitization') - 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]))) - #plt.plot(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),gaussian(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),*coeff)) - mu.append(coeff[1]) - sigma.append(coeff[2]) - - temp = np.array(df[f'ecal_sim_energy_{eng}'].astype(float).to_numpy()[condition])*1000 - hist, x = np.histogram(temp,bins=np.linspace(min(temp),max(temp)+np.std(abs(temp)),2*int(np.sqrt(len(temp))))) - x = x[1:]/2 + x[:-1]/2 - plt.errorbar(x,hist,yerr=np.sqrt(hist),fmt='-o',label='Simulation') - 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]))) - #plt.plot(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),gaussian(np.linspace(coeff[1]-3*coeff[2],coeff[1]+3*coeff[2],50),*coeff)) - mu.append(coeff[1]) - sigma.append(coeff[2]) - - plt.xlabel('Energy (MeV)') - plt.legend() - -#plt.savefig('results/Energy_reconstruction_cluster.pdf') - -mu = np.array(mu) -sigma = np.array(sigma) - -plt.show() - -fig2, (ax1,ax2) = plt.subplots(2,1,figsize=(15,10),sharex=True) - -plt.tight_layout() -# Plot data on primary axis -ax1.scatter(np.array(Energy)*1000, mu[::3], label='cluster') -ax1.scatter(np.array(Energy)*1000, mu[1::3], label='digitization') -ax1.scatter(np.array(Energy)*1000, mu[2::3], label='simulation') - -ax1.plot([4.5,1000],[4.5,1000],c='black',label='x=y') -ax1.set_ylabel('Reconstructed Energy (MeV)') -ax1.set_yscale('log') -ax1.legend() -ax1.set_title('ECal Craterlake Cluster Energy Reconstruction') - -ax2.errorbar(np.array(Energy)*1000, abs(sigma[::3]/mu[::3])*100, fmt='-o', label='cluster') -ax2.errorbar(np.array(Energy)*1000, abs(sigma[1::3]/mu[1::3])*100, fmt='-o', label='digitization') -ax2.errorbar(np.array(Energy)*1000, abs(sigma[2::3]/mu[2::3])*100, fmt='-o', label='simulation') - -ax2.set_ylabel('Resolution (%)') -ax2.set_xlabel('Gamma Energy (MeV)') -ax2.set_xscale('log') -ax2.legend() - -#plt.savefig('results/Energy_resolution.pdf') - -plt.show() - - -htower = [] -herr = [] -hmean = [] -hhits = [] -hhits_cut = [] -emean = [] -ehits = [] -etower = [] -eerr = [] -ehits_cut = [] - -fig3, ax = plt.subplots(2,3,figsize=(20,10)) -fig3.suptitle('ZDC Simulation Energy Reconstruction') -for i in range(6): - plt.sca(ax[i//3,i%3]) - eng = int(Energy[i]*1000) - - x = df[f'par_x_{eng}'].astype(float).to_numpy() - y = df[f'par_y_{eng}'].astype(float).to_numpy() - z = df[f'par_z_{eng}'].astype(float).to_numpy() - x, z = rotateY(x,z, 0.025) - theta = np.arccos(z/np.sqrt((x**2+y**2+z**2)))*1000 - condition = theta <= 3.5 - - plt.title(f'Gamma Energy: {eng} MeV') - energy1 = df[f'hcal_sim_energy_{eng}'].astype(float).to_numpy()#df.eval(f'hcal_sim_energy_{eng}').apply(lambda row: sum(row)) - hist, x = np.histogram(energy1*1000,bins=np.logspace(0,3,200)) - x = x[1:]/2 + x[:-1]/2 - plt.plot(x,hist,marker='o',label="HCal") - hhits.append(len(energy1[energy1!=0])) - condition1 = energy1!=0 - hhits_cut.append(len(energy1[condition & condition1])/len(condition[condition==True])) - energy = df[f'ecal_sim_energy_{eng}'].astype(float).to_numpy()#df.eval(f'ecal_sim_energy_{eng}').apply(lambda row: sum(row)) - hist, x = np.histogram(energy*1000,bins=np.logspace(0,3,200)) - x = x[1:]/2 + x[:-1]/2 - plt.plot(x,hist,marker='o',label="ECal") - emean.append(sum(energy[energy!=0])*1000/len(energy[energy!=0])) - hmean.append(sum(energy1[energy!=0])*1000/len(energy[energy!=0])) - condition1 = energy!=0 - ehits_cut.append(len(energy[condition & condition1])/len(condition[condition==True])) - ehits.append(len(energy[energy!=0])) - plt.legend() - plt.xscale('log') - plt.xlim(1e0,1e3) - - - - - -plt.xlabel('Energy (MeV)') - -#plt.savefig('results/Energy_deposition.pdf') -plt.show() - -fig4, ax = plt.subplots(2,1,sharex=True,gridspec_kw={'height_ratios': [2,1]}) -plt.sca(ax[0]) -plt.errorbar(np.array(Energy)*1000,np.array(hmean)*47.619+np.array(emean),label='HCal/sf+ECal',fmt='-o') -plt.errorbar(np.array(Energy)*1000,emean,label='ECal',fmt='-o') -plt.legend() -plt.yscale('log') -plt.xscale('log') -plt.ylabel('Simulation Energy (MeV)') -plt.sca(ax[1]) -plt.errorbar(np.array(Energy)*1000,(1 - np.array(emean)/(np.array(hmean)*47.619+np.array(emean)))*100,label='Total/ECal',fmt='-o') -plt.legend() -plt.ylabel('Fraction of energy\n deposited in Hcal (%)') -plt.xlabel('Truth Energy (MeV)') -#plt.savefig('results/Energy_ratio_and_Leakage.pdf') -plt.tight_layout() -plt.show() - -fig5 = plt.figure() -plt.errorbar(np.array(Energy)*1000,np.array(hhits)/1000*100,label='HCal Hits',fmt='-o') -plt.errorbar(np.array(Energy)*1000,np.array(ehits)/1000*100,label='ECal Hits',fmt='-o') -#plt.errorbar(np.array(Energy)*1000,np.array(hhits)/np.array(ehits)*100,label='HCal / ECal',fmt='-o',c='b') - -plt.errorbar(np.array(Energy)*1000,np.array(hhits_cut)*100,label='HCal Hits with 3.5 mRad cut',fmt='-^') -plt.errorbar(np.array(Energy)*1000,np.array(ehits_cut)*100,label='ECal Hits with 3.5 mRad cut',fmt='-^') -#plt.errorbar(np.array(Energy)*1000,np.array(hhits_cut)/np.array(ehits_cut)*100,label='HCal / ECal with 3.5 mRad cut',fmt='-^',c='b') -### 3mrad cuts - -plt.legend() -plt.xlabel('Simulation Truth Gamma Energy (MeV)') -plt.ylabel('Fraction of Events with non-zero energy (%)') -#plt.savefig('results/Hits.pdf') -plt.xscale('log') -plt.show() - -fig6, ax = plt.subplots(2,3,figsize=(20,10)) -fig6.suptitle('ZDC Clustering') -fig6.tight_layout(pad=1.8) -for i in range(6): - plt.sca(ax[i//3,i%3]) - eng = int(Energy[i]*1000) - - x = df[f'par_x_{eng}'].astype(float).to_numpy() - y = df[f'par_y_{eng}'].astype(float).to_numpy() - z = df[f'par_z_{eng}'].astype(float).to_numpy() - x, z = rotateY(x,z, 0.025) - theta = np.arccos(z/np.sqrt((x**2+y**2+z**2)))*1000 - condition = theta <= 3.5 - - plt.hist(df[f'ecal_reco_clusters_{eng}'][condition],bins=np.linspace(0,5,6)) - plt.xlabel('Number of Clusters') - plt.title(f'Gamma Energy: {eng} MeV') -plt.show() - -fig7, ax = plt.subplots(2,3,figsize=(20,10)) -fig7.suptitle('ZDC Towering in Clusters') -fig7.tight_layout(pad=1.8) -for i in range(6): - plt.sca(ax[i//3,i%3]) - eng = int(Energy[i]*1000) - - x = df[f'par_x_{eng}'].astype(float).to_numpy() - y = df[f'par_y_{eng}'].astype(float).to_numpy() - z = df[f'par_z_{eng}'].astype(float).to_numpy() - x, z = rotateY(x,z, 0.025) - theta = np.arccos(z/np.sqrt((x**2+y**2+z**2)))*1000 - condition = theta <= 3.5 - - 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)) - plt.xlabel('Number of tower in Clusters') - plt.title(f'Gamma Energy: {eng} MeV') -plt.show() - - -#pdfs = ['results/Energy_reconstruction_cluster.pdf','results/Energy_resolution.pdf','results/Energy_deposition.pdf','results/Energy_ratio_and_Leakage.pdf','results/Hits.pdf'] -with PdfPages(f'results/{os.environ["DETECTOR_CONFIG"]}/lfhcal/plots.pdf') as pdf: - pdf.savefig(fig1) - pdf.savefig(fig2) - pdf.savefig(fig3) - pdf.savefig(fig4) - pdf.savefig(fig5) - pdf.savefig(fig6) - pdf.savefig(fig7) diff --git a/benchmarks/lfhcal/gen_particles.cxx b/benchmarks/lfhcal/gen_particles.cxx deleted file mode 100644 index b84b1170..00000000 --- a/benchmarks/lfhcal/gen_particles.cxx +++ /dev/null @@ -1,126 +0,0 @@ -#include "HepMC3/GenEvent.h" -#include "HepMC3/ReaderAscii.h" -#include "HepMC3/WriterAscii.h" -#include "HepMC3/Print.h" - -#include "TRandom3.h" -#include "TVector3.h" - -#include -#include -#include -#include -#include -#include -#include - -using namespace HepMC3; - -// Generate single electron as input to the Insert simulation. -// -- -// We generate events with a constant polar angle with respect to -// the proton direction and then rotate so that the events are given -// in normal lab coordinate system -// -- -void gen_particles( - int n_events = 10, - const char* out_fname = "gen_particles.hepmc", - TString particle_name = "e-", - double th_min = 3., // Minimum polar angle, in degrees - double th_max = 3., // Maximum polar angle, in degrees - double phi_min = 0., // Minimum azimuthal angle, in degrees - double phi_max = 360., // Maximum azimuthal angle, in degrees - double p = 10., // Momentum in GeV/c - int dist = 0 //Momentum distribution: 0=fixed, 1=uniform, 2=Gaussian - ) -{ - WriterAscii hepmc_output(out_fname); - int events_parsed = 0; - GenEvent evt(Units::GEV, Units::MM); - - // Random number generator - TRandom3 *r1 = new TRandom3(0); //Use time as random seed - - // Getting generated particle information - TDatabasePDG *pdg = new TDatabasePDG(); - TParticlePDG *particle = pdg->GetParticle(particle_name); - const double mass = particle->Mass(); - const int pdgID = particle->PdgCode(); - - for (events_parsed = 0; events_parsed < n_events; events_parsed++) { - - //Set the event number - evt.set_event_number(events_parsed); - - // FourVector(px,py,pz,e,pdgid,status) - // type 4 is beam - // pdgid 11 - electron - // pdgid 111 - pi0 - // pdgid 2212 - proton - GenParticlePtr p1 = - std::make_shared(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4); - GenParticlePtr p2 = std::make_shared( - FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4); - - // Define momentum with respect to proton direction - double phi = r1->Uniform(phi_min*TMath::DegToRad(),phi_max*TMath::DegToRad()); - double th = r1->Uniform(th_min*TMath::DegToRad(),th_max*TMath::DegToRad()); - - //Total momentum distribution - double pevent = -1; - if(dist==0){ //fixed - pevent = p; - } - else if(dist==1){ //Uniform: +-50% variation - pevent = p*(1. + r1->Uniform(-0.5,0.5) ); - } - else if(dist==2){ //Gaussian: Sigma = 0.1*mean - while(pevent<0) //Avoid negative values - pevent = r1->Gaus(p,0.1*p); - } - - double px = pevent * std::cos(phi) * std::sin(th); - double py = pevent * std::sin(phi) * std::sin(th); - double pz = pevent * std::cos(th); - TVector3 pvec(px,py,pz); - - //Rotate to lab coordinate system - double cross_angle = -25./1000.; //in Rad - TVector3 pbeam_dir(sin(cross_angle),0,cos(cross_angle)); //proton beam direction - pvec.RotateY(-pbeam_dir.Theta()); // Theta is returned positive, beam in negative X - // type 1 is final state - // pdgid 11 - electron 0.510 MeV/c^2 - GenParticlePtr p3 = std::make_shared( - FourVector( - pvec.X(), pvec.Y(), pvec.Z(), - sqrt(pevent*pevent + (mass * mass))), - pdgID, 1); - - //If wanted, set non-zero vertex - double vx = 0.; - double vy = 0.; - double vz = 0.; - double vt = 0.; - - GenVertexPtr v1 = std::make_shared(); - evt.shift_position_by(FourVector(vx, vy, vz, vt)); - v1->add_particle_in(p1); - v1->add_particle_in(p2); - - v1->add_particle_out(p3); - evt.add_vertex(v1); - - if (events_parsed == 0) { - std::cout << "First event: " << std::endl; - Print::listing(evt); - } - - hepmc_output.write_event(evt); - if (events_parsed % 100 == 0) { - std::cout << "Event: " << events_parsed << std::endl; - } - evt.clear(); - } - hepmc_output.close(); - std::cout << "Events parsed and written: " << events_parsed << std::endl; -} From 4e20fefb876005550f7c640083e648eb6a7fc2f3 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 5 Oct 2024 21:27:38 -0400 Subject: [PATCH 12/30] .gitlab-ci.yml: add lfhcal --- .gitlab-ci.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index a56c79d2..f25c74e9 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -132,6 +132,7 @@ include: - local: 'benchmarks/tracking_performances_dis/config.yml' - local: 'benchmarks/barrel_ecal/config.yml' - local: 'benchmarks/barrel_hcal/config.yml' + - local: 'benchmarks/lfhcal/config.yml' - local: 'benchmarks/zdc/config.yml' - local: 'benchmarks/zdc_lyso/config.yml' - local: 'benchmarks/zdc_photon/config.yml' From ffdfae931e8a8cd218bb651f310a21bf5fa6fa1f Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 5 Oct 2024 21:47:07 -0400 Subject: [PATCH 13/30] DONTMERGE disable other benchmarks --- .gitlab-ci.yml | 25 +------------------------ 1 file changed, 1 insertion(+), 24 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index f25c74e9..5e61024a 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -123,31 +123,8 @@ get_data: when: - runner_system_failure -include: - - local: 'benchmarks/backgrounds/config.yml' - - local: 'benchmarks/backwards_ecal/config.yml' - - local: 'benchmarks/ecal_gaps/config.yml' - - local: 'benchmarks/tracking_detectors/config.yml' - - local: 'benchmarks/tracking_performances/config.yml' - - local: 'benchmarks/tracking_performances_dis/config.yml' - - local: 'benchmarks/barrel_ecal/config.yml' - - local: 'benchmarks/barrel_hcal/config.yml' +include: - local: 'benchmarks/lfhcal/config.yml' - - local: 'benchmarks/zdc/config.yml' - - local: 'benchmarks/zdc_lyso/config.yml' - - local: 'benchmarks/zdc_photon/config.yml' - - local: 'benchmarks/zdc_pi0/config.yml' - - local: 'benchmarks/material_scan/config.yml' - - local: 'benchmarks/pid/config.yml' - - local: 'benchmarks/timing/config.yml' - - local: 'benchmarks/b0_tracker/config.yml' - - local: 'benchmarks/insert_muon/config.yml' - - local: 'benchmarks/zdc_sigma/config.yml' - - local: 'benchmarks/zdc_lambda/config.yml' - - local: 'benchmarks/insert_neutron/config.yml' - - local: 'benchmarks/femc_electron/config.yml' - - local: 'benchmarks/femc_photon/config.yml' - - local: 'benchmarks/femc_pi0/config.yml' deploy_results: allow_failure: true stage: deploy From 9e24471e569d574b53f6d85be436c4f60b698221 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 5 Oct 2024 21:55:22 -0400 Subject: [PATCH 14/30] Snakemake: copy more stuff from tracking becnhmark, run lfhcal_local rule --- benchmarks/lfhcal/config.yml | 27 +++++++++++++++++++++------ 1 file changed, 21 insertions(+), 6 deletions(-) diff --git a/benchmarks/lfhcal/config.yml b/benchmarks/lfhcal/config.yml index 3da3cf52..3cac2538 100644 --- a/benchmarks/lfhcal/config.yml +++ b/benchmarks/lfhcal/config.yml @@ -1,17 +1,32 @@ sim:lfhcal: extends: .det_benchmark stage: simulate + parallel: + matrix: + - PARTICLE: ["pi-"] + MOMENTUM: ["500MeV", "1GeV", "2GeV", "5GeV", "10GeV", "20GeV"] script: - - snakemake --cores 1 lfhcal - retry: - max: 2 - when: - - runner_system_failure + - | + snakemake --cache --cores 5 \ + sim_output/lfhcal/epic_craterlake/${PARTICLE}/${MOMENTUM}/3to50deg/${PARTICLE}_${MOMENTUM}_3to50deg.0000.eicrecon.tree.edm4eic.root \ + sim_output/lfhcal/epic_craterlake/${PARTICLE}/${MOMENTUM}/45to135deg/${PARTICLE}_${MOMENTUM}_45to135deg.0000.eicrecon.tree.edm4eic.root \ + sim_output/lfhcal/epic_craterlake/${PARTICLE}/${MOMENTUM}/130to177deg/${PARTICLE}_${MOMENTUM}_130to177deg.0000.eicrecon.tree.edm4eic.root + +bench:lfhcal: + extends: .det_benchmark + stage: benchmarks + needs: + - ["sim:lfhcal"] + script: + - snakemake --cores 3 lfhcal_local collect_results:lfhcal: extends: .det_benchmark stage: collect needs: - - "sim:lfhcal" + - "bench:lfhcal" script: - ls -lrht + - mv results{,_save}/ # move results directory out of the way to preserve it + - snakemake --cores 1 --delete-all-output lfhcal_local + - mv results{_save,}/ From 64ae88b19b9b70c0cf75c8a8eeadc63f4640a933 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 5 Oct 2024 21:55:53 -0400 Subject: [PATCH 15/30] .gitlab-ci.yml: collect_results:lfhcal --- .gitlab-ci.yml | 22 +--------------------- 1 file changed, 1 insertion(+), 21 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 5e61024a..c4120c46 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -129,27 +129,7 @@ deploy_results: allow_failure: true stage: deploy needs: - - "collect_results:backgrounds" - - "collect_results:backwards_ecal" - - "collect_results:barrel_ecal" - - "collect_results:barrel_hcal" - - "collect_results:ecal_gaps" - - "collect_results:material_scan" - - "collect_results:pid" - - "collect_results:tracking_performance" - - "collect_results:tracking_performance_campaigns" - - "collect_results:zdc_sigma" - - "collect_results:zdc_lambda" - - "collect_results:insert_neutron" - - "collect_results:tracking_performances_dis" - - "collect_results:zdc" - - "collect_results:zdc_lyso" - - "collect_results:insert_muon" - - "collect_results:zdc_photon" - - "collect_results:zdc_pi0" - - "collect_results:femc_electron" - - "collect_results:femc_photon" - - "collect_results:femc_pi0" + - "collect_results:lfhcal" script: - echo "deploy results!" - find results -print | sort | tee summary.txt From 51c3a4d8c5b50a6f248f4496c1101d87a73b47ce Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 5 Oct 2024 22:04:23 -0400 Subject: [PATCH 16/30] only simulate 3to50deg --- benchmarks/lfhcal/config.yml | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/benchmarks/lfhcal/config.yml b/benchmarks/lfhcal/config.yml index 3cac2538..21159308 100644 --- a/benchmarks/lfhcal/config.yml +++ b/benchmarks/lfhcal/config.yml @@ -8,9 +8,7 @@ sim:lfhcal: script: - | snakemake --cache --cores 5 \ - sim_output/lfhcal/epic_craterlake/${PARTICLE}/${MOMENTUM}/3to50deg/${PARTICLE}_${MOMENTUM}_3to50deg.0000.eicrecon.tree.edm4eic.root \ - sim_output/lfhcal/epic_craterlake/${PARTICLE}/${MOMENTUM}/45to135deg/${PARTICLE}_${MOMENTUM}_45to135deg.0000.eicrecon.tree.edm4eic.root \ - sim_output/lfhcal/epic_craterlake/${PARTICLE}/${MOMENTUM}/130to177deg/${PARTICLE}_${MOMENTUM}_130to177deg.0000.eicrecon.tree.edm4eic.root + sim_output/lfhcal/epic_craterlake/${PARTICLE}/${MOMENTUM}/3to50deg/${PARTICLE}_${MOMENTUM}_3to50deg.0000.eicrecon.tree.edm4eic.root bench:lfhcal: extends: .det_benchmark From 64cba694e6b3e6f6f2725e256bf7cf1f8618eb46 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 5 Oct 2024 22:23:01 -0400 Subject: [PATCH 17/30] Snakefile: add neutron --- benchmarks/lfhcal/Snakefile | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/benchmarks/lfhcal/Snakefile b/benchmarks/lfhcal/Snakefile index 840ec52c..d83a2d51 100644 --- a/benchmarks/lfhcal/Snakefile +++ b/benchmarks/lfhcal/Snakefile @@ -8,7 +8,7 @@ rule lfhcal_sim: "sim_output/lfhcal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.edm4hep.root.log", wildcard_constraints: ENERGY="[0-9]+[kMG]eV", - PARTICLE="pi-", + PARTICLE="(neutron|pi-)", PHASE_SPACE="3to50deg", INDEX="\d{4}", params: @@ -53,7 +53,7 @@ rule lfhcal_at_momentum: DETECTOR_CONFIG="epic_craterlake", ENERGY=f"{float(wildcards.MOMENTUM):.0f}GeV" if float(wildcards.MOMENTUM) >= 1 else f"{float(wildcards.MOMENTUM) * 1000:.0f}MeV", PHASE_SPACE=["3to50deg"], - PARTICLE=["pi-"], + PARTICLE=["neutron", "pi-"], INDEX=range(1), ) if wildcards.CAMPAIGN == "local" else @@ -62,7 +62,7 @@ rule lfhcal_at_momentum: DETECTOR_CONFIG="epic_craterlake", ENERGY=f"{float(wildcards.MOMENTUM):.0f}GeV" if float(wildcards.MOMENTUM) >= 1 else f"{float(wildcards.MOMENTUM) * 1000:.0f}MeV", PHASE_SPACE=["3to50deg"], - PARTICLE=["pi-"], + PARTICLE=["neutron", "pi-"], INDEX=range(1), )), output: @@ -82,7 +82,7 @@ rule lfhcal_summary_at_eta: "{{CAMPAIGN}}/{PARTICLE}/mom/lfhcal_mom_{MOMENTUM:.1f}_mom_resol_{PARTICLE}.root", ], MOMENTUM=[0.5, 1.0, 2.0, 5.0, 10.0, 20.0], - PARTICLE=["pi-"] + PARTICLE=["neutron", "pi-"] ), script="benchmarks/lfhcal/doCompare_widebins_mom.C", output: @@ -105,7 +105,7 @@ else EXTRA_LEGEND="ePIC Simulation {wildcards.CAMPAIGN}" fi cd {wildcards.CAMPAIGN} -root -l -b -q ../{input.script}'("pi-", {wildcards.ETA_MIN}, {wildcards.ETA_MAX}, 1., true, "'"$EXTRA_LEGEND"'")' +root -l -b -q ../{input.script}'("{wildcards.PARTICLE}", {wildcards.ETA_MIN}, {wildcards.ETA_MAX}, 1., true, "'"$EXTRA_LEGEND"'")' """ LFHCAL_ETA_BINS = [1.2,1.5,2,2.5,3,3.5] @@ -118,7 +118,7 @@ rule lfhcal: "{{CAMPAIGN}}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{ETA_BIN}.root", ], ETA_BIN=[f"{eta_min:.1f}_eta_{eta_max:.1f}" for eta_min, eta_max in zip(LFHCAL_ETA_BINS[:-1], LFHCAL_ETA_BINS[1:])], - PARTICLE="pi-", + PARTICLE=["neutron", "pi-"], ) output: directory("results/lfhcal/{CAMPAIGN}/") From 83e069a5cb4108ffc4232950e7ee0c4c0820e922 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 5 Oct 2024 22:23:33 -0400 Subject: [PATCH 18/30] config.yml: add neutron --- benchmarks/lfhcal/config.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/lfhcal/config.yml b/benchmarks/lfhcal/config.yml index 21159308..3bff69f4 100644 --- a/benchmarks/lfhcal/config.yml +++ b/benchmarks/lfhcal/config.yml @@ -3,7 +3,7 @@ sim:lfhcal: stage: simulate parallel: matrix: - - PARTICLE: ["pi-"] + - PARTICLE: ["neutron", "pi-"] MOMENTUM: ["500MeV", "1GeV", "2GeV", "5GeV", "10GeV", "20GeV"] script: - | From 0d48166e7e5c0ac6d672ff292d0145d5f8b791a8 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 5 Oct 2024 22:44:54 -0400 Subject: [PATCH 19/30] Snakefile: fix sim input --- benchmarks/lfhcal/Snakefile | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/benchmarks/lfhcal/Snakefile b/benchmarks/lfhcal/Snakefile index d83a2d51..4b366627 100644 --- a/benchmarks/lfhcal/Snakefile +++ b/benchmarks/lfhcal/Snakefile @@ -49,20 +49,18 @@ rule lfhcal_at_momentum: # TODO pass as a file list? sim=lambda wildcards: expand( - "sim_output/lfhcal/{DETECTOR_CONFIG}/{{PARTICLE}}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root", + "sim_output/lfhcal/{DETECTOR_CONFIG}/{{PARTICLE}}/{ENERGY}/{PHASE_SPACE}/{{PARTICLE}}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root", DETECTOR_CONFIG="epic_craterlake", ENERGY=f"{float(wildcards.MOMENTUM):.0f}GeV" if float(wildcards.MOMENTUM) >= 1 else f"{float(wildcards.MOMENTUM) * 1000:.0f}MeV", PHASE_SPACE=["3to50deg"], - PARTICLE=["neutron", "pi-"], INDEX=range(1), ) if wildcards.CAMPAIGN == "local" else ancient(expand( - "EPIC/RECO/{{CAMPAIGN}}/epic_craterlake/SINGLE/{{PARTICLE}}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root", + "EPIC/RECO/{{CAMPAIGN}}/epic_craterlake/SINGLE/{{PARTICLE}}/{ENERGY}/{PHASE_SPACE}/{{PARTICLE}}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root", DETECTOR_CONFIG="epic_craterlake", ENERGY=f"{float(wildcards.MOMENTUM):.0f}GeV" if float(wildcards.MOMENTUM) >= 1 else f"{float(wildcards.MOMENTUM) * 1000:.0f}MeV", PHASE_SPACE=["3to50deg"], - PARTICLE=["neutron", "pi-"], INDEX=range(1), )), output: From ee5c94b9b20794f7cb0c25621961880e73995789 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 5 Oct 2024 22:58:56 -0400 Subject: [PATCH 20/30] doCompare_widebins_mom.C: fix pi- --- benchmarks/lfhcal/doCompare_widebins_mom.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/lfhcal/doCompare_widebins_mom.C b/benchmarks/lfhcal/doCompare_widebins_mom.C index 41285e03..0e8bd0b9 100644 --- a/benchmarks/lfhcal/doCompare_widebins_mom.C +++ b/benchmarks/lfhcal/doCompare_widebins_mom.C @@ -62,7 +62,7 @@ void doCompare_widebins_mom(TString particle = "pi-",double etamin=-1.0, double cp->SetMargin(0.10, 0.05 ,0.1,0.07); //pi-/mom/lfhcal_mom_20.0_mom_resol_pi-.root - fmom[i] = TFile::Open(Form("./pi-/mom/lfhcal_mom_%1.1f_mom_resol_%s.root",mom[i],particle.Data())); + fmom[i] = TFile::Open(Form("./%s/mom/lfhcal_mom_%1.1f_mom_resol_%s.root",particle.Data(),mom[i],particle.Data())); TH1D *hist = (TH1D*) fmom[i]->Get(Form("hist_mom_%1.1f_%1.1f_pmax_%1.1f",mom[i],etamin,etamax)); hist->Rebin(2); From 0fde4a5cb20e9cd94cda71399667417301a34ea3 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Sat, 5 Oct 2024 23:23:33 -0400 Subject: [PATCH 21/30] rename output files to include particle name --- benchmarks/lfhcal/Snakefile | 11 +++++------ benchmarks/lfhcal/doCompare_widebins_mom.C | 4 ++-- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/benchmarks/lfhcal/Snakefile b/benchmarks/lfhcal/Snakefile index 4b366627..2992cfda 100644 --- a/benchmarks/lfhcal/Snakefile +++ b/benchmarks/lfhcal/Snakefile @@ -77,15 +77,14 @@ rule lfhcal_summary_at_eta: input: expand( [ - "{{CAMPAIGN}}/{PARTICLE}/mom/lfhcal_mom_{MOMENTUM:.1f}_mom_resol_{PARTICLE}.root", + "{{CAMPAIGN}}/{{PARTICLE}}/mom/lfhcal_mom_{MOMENTUM:.1f}_mom_resol_{{PARTICLE}}.root", ], MOMENTUM=[0.5, 1.0, 2.0, 5.0, 10.0, 20.0], - PARTICLE=["neutron", "pi-"] ), script="benchmarks/lfhcal/doCompare_widebins_mom.C", output: - "{CAMPAIGN}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{ETA_MIN}_eta_{ETA_MAX}.png", - "{CAMPAIGN}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{ETA_MIN}_eta_{ETA_MAX}.root", + "{CAMPAIGN}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{PARTICLE}_{ETA_MIN}_eta_{ETA_MAX}.png", + "{CAMPAIGN}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{PARTICLE}_{ETA_MIN}_eta_{ETA_MAX}.root", shell: """ if [[ "{wildcards.CAMPAIGN}" == "local" ]]; then @@ -112,8 +111,8 @@ rule lfhcal: input: expand( [ - "{{CAMPAIGN}}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{ETA_BIN}.png", - "{{CAMPAIGN}}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{ETA_BIN}.root", + "{{CAMPAIGN}}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{PARTICLE}_{ETA_BIN}.png", + "{{CAMPAIGN}}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{PARTICLE}_{ETA_BIN}.root", ], ETA_BIN=[f"{eta_min:.1f}_eta_{eta_max:.1f}" for eta_min, eta_max in zip(LFHCAL_ETA_BINS[:-1], LFHCAL_ETA_BINS[1:])], PARTICLE=["neutron", "pi-"], diff --git a/benchmarks/lfhcal/doCompare_widebins_mom.C b/benchmarks/lfhcal/doCompare_widebins_mom.C index 0e8bd0b9..5b4e7a47 100644 --- a/benchmarks/lfhcal/doCompare_widebins_mom.C +++ b/benchmarks/lfhcal/doCompare_widebins_mom.C @@ -100,7 +100,7 @@ void doCompare_widebins_mom(TString particle = "pi-",double etamin=-1.0, double } - TFile *fout = new TFile(Form("Final_Results/%s/mom/lfhcal_mom_resol_%1.1f_eta_%1.1f.root",particle.Data(),etamin,etamax),"recreate"); + TFile *fout = new TFile(Form("Final_Results/%s/mom/lfhcal_mom_resol_%s_%1.1f_eta_%1.1f.root",particle.Data(),particle.Data(),etamin,etamax),"recreate"); TGraphErrors *gr1 = new TGraphErrors(size,p,sigma_p,err_p,err_sigma_p); gr1->SetName("grseed"); gr1->SetMarkerStyle(25); @@ -119,7 +119,7 @@ void doCompare_widebins_mom(TString particle = "pi-",double etamin=-1.0, double lmom->AddEntry(gr1,"Nominal"); lmom->Draw("same"); //draw_req_Mom(etamin,etamax,0.,mgMom->GetXaxis()->GetXmax()); - c_mom->SaveAs(Form("Final_Results/%s/mom/lfhcal_mom_resol_%1.1f_eta_%1.1f.png",particle.Data(),etamin,etamax)); + c_mom->SaveAs(Form("Final_Results/%s/mom/lfhcal_mom_resol_%s_%1.1f_eta_%1.1f.png",particle.Data(),particle.Data(),etamin,etamax)); // Write the numbers in output file for comparisons outfile << extra_legend << endl; From 3d5f9c07856acbadf47b02793808a6062dc63d3a Mon Sep 17 00:00:00 2001 From: steinber Date: Sun, 6 Oct 2024 09:59:20 -0400 Subject: [PATCH 22/30] bugfixes to basic analysis script --- benchmarks/lfhcal/LFHCAL_Performance.C | 44 +++++++++++++++----------- 1 file changed, 26 insertions(+), 18 deletions(-) diff --git a/benchmarks/lfhcal/LFHCAL_Performance.C b/benchmarks/lfhcal/LFHCAL_Performance.C index fc4dcbbb..dcfb7595 100644 --- a/benchmarks/lfhcal/LFHCAL_Performance.C +++ b/benchmarks/lfhcal/LFHCAL_Performance.C @@ -38,7 +38,7 @@ void LFHCAL_Performance(TString filename="tracking_output",TString particle="pi- for (int i=0; iSetTitle(Form("%1.1f < #eta < %1.1f && p = %1.1f ",eta[i],eta[i+1],mom)); histp[i]->SetName(Form("hist_mom_%1.1f_%1.1f_pmax_%1.1f",mom,eta[i],eta[i+1])); } @@ -72,6 +72,7 @@ void LFHCAL_Performance(TString filename="tracking_output",TString particle="pi- int matchId = 1; // Always matched track assigned the index 0 while (myReader.Next()) { + std::cout << "events = " << count++ << std::endl; for (int j = 0; j < pdg.GetSize(); ++j) { if (status[j] !=1 && pdg.GetSize()!=1) continue; @@ -80,30 +81,35 @@ void LFHCAL_Performance(TString filename="tracking_output",TString particle="pi- Double_t pmc = sqrt(px_mc[j]*px_mc[j]+py_mc[j]*py_mc[j]+pz_mc[j]*pz_mc[j]); // 1./(q/p); similar to prec Double_t etamc = -1.0*TMath::Log(TMath::Tan((TMath::ACos(pzmc/fabs(pmc)))/2)); Double_t phimc = TMath::ATan2(py_mc[j],px_mc[j]); + std::cout << "neutron p=" << pmc << " pt=" << ptmc << std::endl; + if (fabs(ptmc) < pTcut) continue; - Double_t l_px_tot; - Double_t l_py_tot; - Double_t l_pz_tot; - Double_t l_e_tot; + float l_px_tot=0; + float l_py_tot=0; + float l_pz_tot=0; + float l_e_tot=0; + std::cout << "LFHCAL nclus=" << px_lc.GetSize() << " ECAL nclus=" << pe_ec.GetSize() << std::endl; for (int jl = 0;jl Date: Mon, 7 Oct 2024 13:26:50 -0400 Subject: [PATCH 24/30] Revert ".gitlab-ci.yml: collect_results:lfhcal" This reverts commit 64ae88b19b9b70c0cf75c8a8eeadc63f4640a933. --- .gitlab-ci.yml | 22 +++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index c4120c46..5e61024a 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -129,7 +129,27 @@ deploy_results: allow_failure: true stage: deploy needs: - - "collect_results:lfhcal" + - "collect_results:backgrounds" + - "collect_results:backwards_ecal" + - "collect_results:barrel_ecal" + - "collect_results:barrel_hcal" + - "collect_results:ecal_gaps" + - "collect_results:material_scan" + - "collect_results:pid" + - "collect_results:tracking_performance" + - "collect_results:tracking_performance_campaigns" + - "collect_results:zdc_sigma" + - "collect_results:zdc_lambda" + - "collect_results:insert_neutron" + - "collect_results:tracking_performances_dis" + - "collect_results:zdc" + - "collect_results:zdc_lyso" + - "collect_results:insert_muon" + - "collect_results:zdc_photon" + - "collect_results:zdc_pi0" + - "collect_results:femc_electron" + - "collect_results:femc_photon" + - "collect_results:femc_pi0" script: - echo "deploy results!" - find results -print | sort | tee summary.txt From 95d7939957dd7b00df2df25f70e2839acc3b9935 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Mon, 7 Oct 2024 13:26:59 -0400 Subject: [PATCH 25/30] Revert "DONTMERGE disable other benchmarks" This reverts commit ffdfae931e8a8cd218bb651f310a21bf5fa6fa1f. --- .gitlab-ci.yml | 25 ++++++++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 5e61024a..f25c74e9 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -123,8 +123,31 @@ get_data: when: - runner_system_failure -include: +include: + - local: 'benchmarks/backgrounds/config.yml' + - local: 'benchmarks/backwards_ecal/config.yml' + - local: 'benchmarks/ecal_gaps/config.yml' + - local: 'benchmarks/tracking_detectors/config.yml' + - local: 'benchmarks/tracking_performances/config.yml' + - local: 'benchmarks/tracking_performances_dis/config.yml' + - local: 'benchmarks/barrel_ecal/config.yml' + - local: 'benchmarks/barrel_hcal/config.yml' - local: 'benchmarks/lfhcal/config.yml' + - local: 'benchmarks/zdc/config.yml' + - local: 'benchmarks/zdc_lyso/config.yml' + - local: 'benchmarks/zdc_photon/config.yml' + - local: 'benchmarks/zdc_pi0/config.yml' + - local: 'benchmarks/material_scan/config.yml' + - local: 'benchmarks/pid/config.yml' + - local: 'benchmarks/timing/config.yml' + - local: 'benchmarks/b0_tracker/config.yml' + - local: 'benchmarks/insert_muon/config.yml' + - local: 'benchmarks/zdc_sigma/config.yml' + - local: 'benchmarks/zdc_lambda/config.yml' + - local: 'benchmarks/insert_neutron/config.yml' + - local: 'benchmarks/femc_electron/config.yml' + - local: 'benchmarks/femc_photon/config.yml' + - local: 'benchmarks/femc_pi0/config.yml' deploy_results: allow_failure: true stage: deploy From cb4897c25706f0c15f5ef2d4453580c828ff9297 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Mon, 7 Oct 2024 16:09:11 -0400 Subject: [PATCH 26/30] Snakefile: disable cache --- Snakefile | 1 - 1 file changed, 1 deletion(-) diff --git a/Snakefile b/Snakefile index 9c36bc2b..c1114a90 100644 --- a/Snakefile +++ b/Snakefile @@ -35,7 +35,6 @@ def get_remote_path(path): rule fetch_epic: output: filepath="EPIC/{PATH}" - cache: True retries: 3 shell: """ xrdcp --debug 2 root://dtn-eic.jlab.org//work/eic2/{output.filepath} {output.filepath} From e002b390cafd09cb6f741507eae29c8dda056173 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Mon, 7 Oct 2024 18:54:05 -0400 Subject: [PATCH 27/30] add bench:lfhcal_campaigns --- benchmarks/lfhcal/Snakefile | 4 ++-- benchmarks/lfhcal/config.yml | 18 ++++++++++++++++++ 2 files changed, 20 insertions(+), 2 deletions(-) diff --git a/benchmarks/lfhcal/Snakefile b/benchmarks/lfhcal/Snakefile index 3526ec39..b0177433 100644 --- a/benchmarks/lfhcal/Snakefile +++ b/benchmarks/lfhcal/Snakefile @@ -109,13 +109,13 @@ LFHCAL_ETA_BINS = [1.2,1.5,2,2.5,3,3.5] rule lfhcal: input: - expand( + lambda wildcards: expand( [ "{{CAMPAIGN}}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{PARTICLE}_{ETA_BIN}.png", "{{CAMPAIGN}}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{PARTICLE}_{ETA_BIN}.root", ], ETA_BIN=[f"{eta_min:.1f}_eta_{eta_max:.1f}" for eta_min, eta_max in zip(LFHCAL_ETA_BINS[:-1], LFHCAL_ETA_BINS[1:])], - PARTICLE=["neutron", "pi-"], + PARTICLE=["neutron", "pi-"] if wildcards.CAMPAIGN == "local" else ["pi-"], ) output: directory("results/lfhcal/{CAMPAIGN}/") diff --git a/benchmarks/lfhcal/config.yml b/benchmarks/lfhcal/config.yml index 3bff69f4..f7b9ad3a 100644 --- a/benchmarks/lfhcal/config.yml +++ b/benchmarks/lfhcal/config.yml @@ -28,3 +28,21 @@ collect_results:lfhcal: - mv results{,_save}/ # move results directory out of the way to preserve it - snakemake --cores 1 --delete-all-output lfhcal_local - mv results{_save,}/ + +bench:lfhcal_campaigns: + extends: .det_benchmark + stage: benchmarks + #when: manual + script: + - snakemake --cores 1 lfhcal_campaigns + +collect_results:lfhcal_campaigns: + extends: .det_benchmark + stage: collect + needs: + - "bench:lfhcal_campaigns" + script: + - ls -lrht + - mv results{,_save}/ # move results directory out of the way to preserve it + - snakemake --cores 1 --delete-all-output lfhcal_campaigns + - mv results{_save,}/ From 7381745cc26effdcc64264029796d98e60af38ab Mon Sep 17 00:00:00 2001 From: steinber Date: Mon, 7 Oct 2024 20:28:12 -0400 Subject: [PATCH 28/30] added gammas to particle list --- benchmarks/lfhcal/Snakefile | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/benchmarks/lfhcal/Snakefile b/benchmarks/lfhcal/Snakefile index b0177433..5afbf17e 100644 --- a/benchmarks/lfhcal/Snakefile +++ b/benchmarks/lfhcal/Snakefile @@ -8,7 +8,7 @@ rule lfhcal_sim: "sim_output/lfhcal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.edm4hep.root.log", wildcard_constraints: ENERGY="[0-9]+[kMG]eV", - PARTICLE="(neutron|pi-)", + PARTICLE="(neutron|pi-|gamma)", PHASE_SPACE="3to50deg", INDEX="\d{4}", params: @@ -115,7 +115,7 @@ rule lfhcal: "{{CAMPAIGN}}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{PARTICLE}_{ETA_BIN}.root", ], ETA_BIN=[f"{eta_min:.1f}_eta_{eta_max:.1f}" for eta_min, eta_max in zip(LFHCAL_ETA_BINS[:-1], LFHCAL_ETA_BINS[1:])], - PARTICLE=["neutron", "pi-"] if wildcards.CAMPAIGN == "local" else ["pi-"], + PARTICLE=["neutron", "pi-", "gamma"] if wildcards.CAMPAIGN == "local" else ["pi-"], ) output: directory("results/lfhcal/{CAMPAIGN}/") From 8bc65611613449deb6afd4aac7643f4f03e57ad9 Mon Sep 17 00:00:00 2001 From: steinber Date: Mon, 7 Oct 2024 20:28:36 -0400 Subject: [PATCH 29/30] replaced fit with straight RMS --- benchmarks/lfhcal/doCompare_widebins_mom.C | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/benchmarks/lfhcal/doCompare_widebins_mom.C b/benchmarks/lfhcal/doCompare_widebins_mom.C index 5b4e7a47..5258e5ef 100644 --- a/benchmarks/lfhcal/doCompare_widebins_mom.C +++ b/benchmarks/lfhcal/doCompare_widebins_mom.C @@ -70,7 +70,9 @@ void doCompare_widebins_mom(TString particle = "pi-",double etamin=-1.0, double hist->SetTitle(Form("Momentum = %1.1f && %1.1f<#eta<%1.1f;#Delta p/p; Entries(a.u.)",particle.Data(),mom[i],etamin,etamax)); double mu = hist->GetMean(); - double sigma = hist->GetStdDev(); + double sigma = hist->GetRMS(); + double sigma_err = hist->GetRMSError(); + /* hist->GetXaxis()->SetRangeUser(-1.0*range,1.0*range); func->SetRange(mu-2.0*sigma,mu+2.0*sigma); // fit with in 2 sigma range hist->Fit(func,"NR+"); @@ -80,9 +82,10 @@ void doCompare_widebins_mom(TString particle = "pi-",double etamin=-1.0, double hist->Fit(func,"R+"); float par2 = func->GetParameter(2)*100; float par2_err = func->GetParError(2)*100; + */ momV.push_back(mom[i]); - momresolV.push_back(par2); - err_momresolV.push_back(par2_err); + momresolV.push_back(sigma); + err_momresolV.push_back(sigma_err); cp->cd(); hist->Draw(); From 3e37b66f98a5fea8c831e01df01f3fa2257ad2c3 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Tue, 8 Oct 2024 10:36:17 -0400 Subject: [PATCH 30/30] Snakefile: avoid using deprecated eicrecon option --- benchmarks/lfhcal/Snakefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmarks/lfhcal/Snakefile b/benchmarks/lfhcal/Snakefile index 5afbf17e..459392c0 100644 --- a/benchmarks/lfhcal/Snakefile +++ b/benchmarks/lfhcal/Snakefile @@ -40,7 +40,7 @@ rule lfhcal_recon: shell: """ env DETECTOR_CONFIG={wildcards.DETECTOR_CONFIG} \ eicrecon {input} -Ppodio:output_file={output} \ - -Ppodio:output_include_collections=MCParticles,ReconstructedParticles,LFHCALTruthClusters,LFHCALClusters,LFHCALHits,EcalEndcapPTruthClusters,EcalEndcapPClusters,EcalEndcapPHits + -Ppodio:output_collections=MCParticles,ReconstructedParticles,LFHCALTruthClusters,LFHCALClusters,LFHCALHits,EcalEndcapPTruthClusters,EcalEndcapPClusters,EcalEndcapPHits """ rule lfhcal_at_momentum: