From 2b5779ae02d41e941b822b0e9b4336fce833ffbe Mon Sep 17 00:00:00 2001 From: Barak Schmookler Date: Mon, 15 Jul 2024 19:46:46 -0400 Subject: [PATCH] Added benchmark --- Snakefile | 2 +- benchmarks/zdc_lyso/README.md | 6 + benchmarks/zdc_lyso/Snakefile | 88 +++++++++++++ benchmarks/zdc_lyso/analysis/analysis.py | 154 +++++++++++++++++++++++ benchmarks/zdc_lyso/gen_particles.cxx | 126 +++++++++++++++++++ 5 files changed, 375 insertions(+), 1 deletion(-) create mode 100644 benchmarks/zdc_lyso/README.md create mode 100644 benchmarks/zdc_lyso/Snakefile create mode 100644 benchmarks/zdc_lyso/analysis/analysis.py create mode 100644 benchmarks/zdc_lyso/gen_particles.cxx diff --git a/Snakefile b/Snakefile index e7a54350..41895ad0 100644 --- a/Snakefile +++ b/Snakefile @@ -3,7 +3,7 @@ include: "benchmarks/barrel_ecal/Snakefile" include: "benchmarks/ecal_gaps/Snakefile" include: "benchmarks/material_scan/Snakefile" include: "benchmarks/tracking_performances/Snakefile" - +include: "benchmarks/ rule fetch_epic: output: diff --git a/benchmarks/zdc_lyso/README.md b/benchmarks/zdc_lyso/README.md new file mode 100644 index 00000000..18dbcc5b --- /dev/null +++ b/benchmarks/zdc_lyso/README.md @@ -0,0 +1,6 @@ +Detector Benchmark for LYSO ZDC +=============================== + +## Overview +This benchmark generates events with single low-energy (5 MeV - 1 GeV) photons. The photons are generated with angles of 0-5.24 mRad with respect to the proton/ion beam direction. The benchmark creates performance plots of the LYSO ZDC acceptance and energy reconstruction resolution. + diff --git a/benchmarks/zdc_lyso/Snakefile b/benchmarks/zdc_lyso/Snakefile new file mode 100644 index 00000000..c34a3223 --- /dev/null +++ b/benchmarks/zdc_lyso/Snakefile @@ -0,0 +1,88 @@ +rule all: + 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."], + theta_min=["0"], + theta_max=["0.3"], + phi_min=["0"], + phi_max=["360"]) + + +rule ecal_lyso_sim_hepmc: + input: + gen_par = "gen_particles.cxx", + output: + hepmcfile="data/{particle}_{beam_energy}GeV_theta_{theta_min}deg_thru_{theta_max}deg.hepmc", + log: + "log/{particle}_{beam_energy}GeV_theta_{theta_min}deg_thru_{theta_max}deg.log", + wildcard_constraints: + + params: + physics_list="FTFP_BERT", + num_events=1000, + shell: + """ +root -l -b -q '{input.gen_par}({params.num_events}, "{output.hepmcfile}", "{wildcards.particle}", {wildcards.theta_min}, {wildcards.theta_max}, 0, 360, {wildcards.beam_energy})' +""" + + +rule ecal_lyso_sim: + input: + "data/{particle}_{beam_energy}GeV_theta_{theta_min}deg_thru_{theta_max}deg.hepmc", + output: + "data/{DETECTOR}_sim_{particle}_{beam_energy}GeV_theta_{theta_min}deg_thru_{theta_max}deg.edm4hep.root", + log: + "log/{DETECTOR}_sim_{particle}_{beam_energy}GeV_theta_{theta_min}deg_thru_{theta_max}deg.log", + wildcard_constraints: + + params: + physics_list="FTFP_BERT", + num_events=1000, + + shell: + """ +npsim \ + --compactFile $DETECTOR_PATH/{wildcards.DETECTOR}.xml \ + --numberOfEvents {params.num_events} \ + --physicsList {params.physics_list} \ + --inputFiles {input} \ + --outputFile {output} || exit +""" + + +rule ecal_lyso_reco: + input: + "data/{DETECTOR}_sim_{particle}_{beam_energy}GeV_theta_{theta_min}deg_thru_{theta_max}deg.edm4hep.root", + output: + "data/{DETECTOR}_reco_{particle}_{beam_energy}GeV_theta_{theta_min}deg_thru_{theta_max}deg.edm4hep.root", + log: + "log/{DETECTOR}_reco_{particle}_{beam_energy}GeV_theta_{theta_min}deg_thru_{theta_max}deg.log", + wildcard_constraints: + + params: + DETECTOR="epic_craterlake", + shell: + """ +eicrecon -Ppodio:output_include_collections=HcalFarForwardZDCRawHits,HcalFarForwardZDCRecHits,HcalFarForwardZDCClusters,EcalFarForwardZDCRawHits,EcalFarForwardZDCRecHits,EcalFarForwardZDCClusters,MCParticles {input} +mv podio_output.root {output} +""" + + +rule zdc_analysis: + input: + "data/epic_craterlake_reco_gamma_0.005GeV_theta_0deg_thru_0.3deg.edm4hep.root", + "data/epic_craterlake_reco_gamma_0.01GeV_theta_0deg_thru_0.3deg.edm4hep.root", + "data/epic_craterlake_reco_gamma_0.05GeV_theta_0deg_thru_0.3deg.edm4hep.root", + "data/epic_craterlake_reco_gamma_0.1GeV_theta_0deg_thru_0.3deg.edm4hep.root", + "data/epic_craterlake_reco_gamma_0.5GeV_theta_0deg_thru_0.3deg.edm4hep.root", + "data/epic_craterlake_reco_gamma_1GeV_theta_0deg_thru_0.3deg.edm4hep.root", + + script="analysis.py", + output: + "results/plots.pdf", + + shell: + """ +python {input.script} +""" diff --git a/benchmarks/zdc_lyso/analysis/analysis.py b/benchmarks/zdc_lyso/analysis/analysis.py new file mode 100644 index 00000000..80b489cc --- /dev/null +++ b/benchmarks/zdc_lyso/analysis/analysis.py @@ -0,0 +1,154 @@ +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 + +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) ) + +Energy = [0.005, 0.01, 0.05, 0.1, 0.5, 1] + +df = pd.DataFrame({}) +for eng in Energy: + tree = uproot.open(f'data/epic_craterlake_reco_gamma_{eng}GeV_theta_0deg_thru_0.3deg.edm4hep.root')['events'] + ecal_reco_energy = tree['EcalFarForwardZDCClusters/EcalFarForwardZDCClusters.energy'].array() + #hcal_reco_energy = tree['HcalFarForwardZDCClusters/HcalFarForwardZDCClusters.energy'].array() + + tree = uproot.open(f'data/epic_craterlake_sim_gamma_{eng}GeV_theta_0deg_thru_0.3deg.edm4hep.root')['events'] + ecal_sim_energy = tree['EcalFarForwardZDCHits/EcalFarForwardZDCHits.energy'].array() + hcal_sim_energy = tree['HcalFarForwardZDCHits/HcalFarForwardZDCHits.energy'].array() + + eng = int(eng*1000) + + ecal_reco_energy = pd.DataFrame({f'ecal_reco_energy_{eng}': np.array(ecal_reco_energy.tolist(), dtype=object)}) + #hcal_reco_energy = pd.DataFrame({f'hcal_reco_energy_{eng}': np.array(hcal_reco_energy.tolist(), dtype=object)}) + ecal_sim_energy = pd.DataFrame({f'ecal_sim_energy_{eng}': np.array(ecal_sim_energy.tolist(), dtype=object)}) + hcal_sim_energy = pd.DataFrame({f'hcal_sim_energy_{eng}': np.array(hcal_sim_energy.tolist(), dtype=object)}) + + df = pd.concat([df,ecal_reco_energy,ecal_sim_energy,hcal_sim_energy],axis=1) + + +mu = [] +sigma = [] +resolution = [] +fig1, ax = plt.subplots(3,2,figsize=(20,10)) +plt.tight_layout() +for i in range(6): + plt.sca(ax[i%3,i//3]) + eng = int(Energy[i]*1000) + plt.title(f'Gamma Energy: {eng} MeV') + temp = [item[0] for item in df[f'ecal_reco_energy_{eng}'] if item] + hist, x = np.histogram(temp,bins=15) + x = x[1:]/2 + x[:-1]/2 + plt.scatter(x,hist,marker='o') + coeff, covar = curve_fit(gaussian,x,hist,p0=(200,eng,eng/2)) + plt.plot(np.linspace(coeff[1]-5*coeff[2],coeff[1]+5*coeff[2],50),gaussian(np.linspace(coeff[1]-5*coeff[2],coeff[1]+5*coeff[2],50),*coeff) + ,label = f'$\mu$ = {coeff[1]:.3f} $\pm$ {covar[1][1]:.3f}\n$\sigma$ = {np.abs(coeff[2]):.3f} $\pm$ {covar[2][2]:.3f}\nResolution = {np.abs(coeff[2])*100/coeff[1]:.2f}%') + plt.xlabel('Energy (MeV)') + plt.legend() + mu.append(coeff[1]) + sigma.append(coeff[2]) + resolution.append(np.abs(coeff[2])*100/coeff[1]) +plt.savefig('results/Energy_reconstruction_cluster.pdf') + +plt.show() + +fig2, (ax1,ax2) = plt.subplots(2,1,figsize=(15,10),sharex=True) +plt.tight_layout() +# Plot data on primary axis +ax1.scatter(Energy, mu, c='b') +ax1.plot([0.0045,1],[0.0045,1],c='black',label='x=y') +ax1.set_ylabel('Reconstructed Energy (GeV)') +ax1.legend() +ax1.set_title('ECal Craterlake Energy Reconstruction') + +ax2.plot(Energy, resolution, c='r') +ax2.scatter(Energy, resolution, c='r') +ax2.set_ylabel('Resolution (%)') +ax2.set_xlabel('Gamma Energy (GeV)') +ax2.set_xscale('log') +plt.savefig('results/Energy_resolution.pdf') + +plt.show() + + +htower = [] +herr = [] +hmean = [] +hhits = [] +emean = [] +ehits = [] +etower = [] +eerr = [] +fig3, ax = plt.subplots(2,3,figsize=(20,10)) +for i in range(6): + plt.sca(ax[i//3,i%3]) + + eng = int(Energy[i]*1000) + + energy = np.array([sum(item) for item in df[f'hcal_sim_energy_{eng}'] if item])#df.eval(f'hcal_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="HCal") + hmean.append(sum(energy)*1000/len(energy)) + hhits.append(len(energy[energy!=0])) + energy = np.array([sum(item) for item in df[f'ecal_sim_energy_{eng}'] if item])#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)*1000/len(energy)) + ehits.append(len(energy[energy!=0])) + plt.legend() + plt.xscale('log') + plt.xlim(1e0,1e3) +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,hmean,label='HCal Active Layers',fmt='-o') +plt.errorbar(np.array(Energy)*1000,np.array(hmean)*47.619,label='HCal',fmt='-o') +plt.errorbar(np.array(Energy)*1000,np.array(hmean)*47.619+emean,label='HCal+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,np.array(hmean)/np.array(emean),label='HCal Active Layers/ECal',fmt='-o') +plt.errorbar(np.array(Energy)*1000,np.array(hmean)/np.array(emean)*47.619,label='HCal/ECal',fmt='-o') +plt.legend() +plt.ylabel('Leakage Ratio') +plt.xlabel('Truth Energy (MeV)') +plt.savefig('results/Energy_ratio_and_Leakage.pdf') +plt.tight_layout() +plt.show() + +fig5 = plt.figure() +plt.scatter(Energy,np.array(hhits)/1000*100,label='HCal Hits') +plt.scatter(Energy,np.array(ehits)/1000*100,label='ECal Hits') +plt.plot(Energy,np.array(hhits)/1000*100) +plt.plot(Energy,np.array(ehits)/1000*100) +plt.plot(Energy,np.array(hhits)/np.array(ehits)*100,label='HCal / ECal') +plt.scatter(Energy,np.array(hhits)/np.array(ehits)*100) + +plt.legend() +plt.xlabel('Simulation Truth Gamma Energy (GeV)') +plt.ylabel('Simulation Hits at ZDC (%)') +plt.savefig('results/Hits.pdf') +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("results/plots.pdf") as pdf: + pdf.savefig(fig1) + pdf.savefig(fig2) + pdf.savefig(fig3) + pdf.savefig(fig4) + pdf.savefig(fig5) diff --git a/benchmarks/zdc_lyso/gen_particles.cxx b/benchmarks/zdc_lyso/gen_particles.cxx new file mode 100644 index 00000000..b84b1170 --- /dev/null +++ b/benchmarks/zdc_lyso/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; +}