From 7246e0d3b0018ec3968d307da2e6c4b264f50772 Mon Sep 17 00:00:00 2001
From: Sebouh Paul <sebouh.paul@gmail.com>
Date: Mon, 23 Sep 2024 13:21:21 -0400
Subject: [PATCH] move zdc pi0 and photon benchmarks from physics benchmarks
(#70)
Co-authored-by: Dmitry Kalinkin <dmitry.kalinkin@gmail.com>
---
.gitlab-ci.yml | 4 +
Snakefile | 4 +-
benchmarks/zdc_photon/Snakefile | 70 ++++++
.../zdc_photon/analysis/gen_particles.cxx | 127 ++++++++++
.../zdc_photon/analysis/zdc_photon_plots.py | 161 ++++++++++++
benchmarks/zdc_photon/config.yml | 36 +++
benchmarks/zdc_pi0/Snakefile | 56 +++++
benchmarks/zdc_pi0/analysis/gen_pi0_decay.cxx | 182 ++++++++++++++
benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py | 232 ++++++++++++++++++
benchmarks/zdc_pi0/config.yml | 34 +++
10 files changed, 905 insertions(+), 1 deletion(-)
create mode 100644 benchmarks/zdc_photon/Snakefile
create mode 100644 benchmarks/zdc_photon/analysis/gen_particles.cxx
create mode 100644 benchmarks/zdc_photon/analysis/zdc_photon_plots.py
create mode 100644 benchmarks/zdc_photon/config.yml
create mode 100644 benchmarks/zdc_pi0/Snakefile
create mode 100644 benchmarks/zdc_pi0/analysis/gen_pi0_decay.cxx
create mode 100644 benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py
create mode 100644 benchmarks/zdc_pi0/config.yml
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 2e836ec6..4b3c9c8e 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -141,6 +141,8 @@ include:
- local: 'benchmarks/barrel_hcal/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'
@@ -169,6 +171,8 @@ deploy_results:
- "collect_results:tracking_performances_dis"
- "collect_results:zdc"
- "collect_results:zdc_lyso"
+ - "collect_results:zdc_photon"
+ - "collect_results:zdc_pi0"
script:
- echo "deploy results!"
- find results -print | sort | tee summary.txt
diff --git a/Snakefile b/Snakefile
index 4a12021b..58357199 100644
--- a/Snakefile
+++ b/Snakefile
@@ -7,8 +7,10 @@ include: "benchmarks/ecal_gaps/Snakefile"
include: "benchmarks/material_scan/Snakefile"
include: "benchmarks/tracking_performances/Snakefile"
include: "benchmarks/tracking_performances_dis/Snakefile"
-include: "benchmarks/zdc_lyso/Snakefile"
include: "benchmarks/zdc_lambda/Snakefile"
+include: "benchmarks/zdc_lyso/Snakefile"
+include: "benchmarks/zdc_photon/Snakefile"
+include: "benchmarks/zdc_pi0/Snakefile"
include: "benchmarks/zdc_sigma/Snakefile"
include: "benchmarks/insert_neutron/Snakefile"
diff --git a/benchmarks/zdc_photon/Snakefile b/benchmarks/zdc_photon/Snakefile
new file mode 100644
index 00000000..cfb51242
--- /dev/null
+++ b/benchmarks/zdc_photon/Snakefile
@@ -0,0 +1,70 @@
+rule zdc_photon_generate:
+ input:
+ script="benchmarks/zdc_photon/analysis/gen_particles.cxx",
+ params:
+ th_max=0.23,
+ th_min=0
+ output:
+ GEN_FILE="sim_output/zdc_photon/zdc_photon_{P}GeV.hepmc"
+ shell:
+ """
+if [[ {wildcards.P} -gt 140 ]]; then
+ NEVENTS_GEN=300
+else
+ NEVENTS_GEN=1000
+fi
+mkdir -p sim_output/zdc_photon
+root -l -b -q '{input.script}('$NEVENTS_GEN',"{output.GEN_FILE}", "gamma", {params.th_min}, {params.th_max}, 0., 360., {wildcards.P})'
+"""
+
+rule zdc_photon_simulate:
+ input:
+ GEN_FILE="sim_output/zdc_photon/zdc_photon_{P}GeV.hepmc"
+ params:
+ PHYSICS_LIST="FTFP_BERT"
+ output:
+ SIM_FILE="sim_output/zdc_photon/{DETECTOR_CONFIG}_sim_zdc_photon_{P}GeV.edm4hep.root"
+ shell:
+ """
+# Running simulation
+if [[ {wildcards.P} -gt 140 ]]; then
+ NEVENTS_SIM=300
+else
+ NEVENTS_SIM=1000
+fi
+npsim \
+ --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \
+ --physicsList {params.PHYSICS_LIST} \
+ --numberOfEvents $NEVENTS_SIM \
+ --inputFiles {input.GEN_FILE} \
+ --outputFile {output.SIM_FILE}
+"""
+
+rule zdc_photon_recon:
+ input:
+ SIM_FILE="sim_output/zdc_photon/{DETECTOR_CONFIG}_sim_zdc_photon_{P}GeV.edm4hep.root"
+ output:
+ REC_FILE="sim_output/zdc_photon/{DETECTOR_CONFIG}_rec_zdc_photon_{P}GeV.edm4hep.root"
+ shell:
+ """
+if [[ {wildcards.P} -gt 140 ]]; then
+ NEVENTS_REC=300
+else
+ NEVENTS_REC=1000
+fi
+eicrecon {input.SIM_FILE} -Ppodio:output_file={output.REC_FILE} -Pdd4hep:xml_files=$DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml -Ppodio:output_collections=MCParticles,HcalFarForwardZDCRecHits,HcalFarForwardZDCClusters,HcalFarForwardZDCSubcellHits -Pjana:nevents=$NEVENTS_REC
+"""
+
+rule zdc_photon_analysis:
+ input:
+ expand("sim_output/zdc_photon/{DETECTOR_CONFIG}_rec_zdc_photon_{P}GeV.edm4hep.root",
+ P=[20, 30, 50, 70, 100, 150, 200, 275],
+ DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]),
+ script="benchmarks/zdc_photon/analysis/zdc_photon_plots.py",
+ output:
+ results_dir=directory("results/{DETECTOR_CONFIG}/zdc_photon"),
+ shell:
+ """
+mkdir -p {output.results_dir}
+python {input.script} {output.results_dir}
+"""
diff --git a/benchmarks/zdc_photon/analysis/gen_particles.cxx b/benchmarks/zdc_photon/analysis/gen_particles.cxx
new file mode 100644
index 00000000..0877a852
--- /dev/null
+++ b/benchmarks/zdc_photon/analysis/gen_particles.cxx
@@ -0,0 +1,127 @@
+#include "HepMC3/GenEvent.h"
+#include "HepMC3/ReaderAscii.h"
+#include "HepMC3/WriterAscii.h"
+#include "HepMC3/Print.h"
+
+#include "TRandom3.h"
+#include "TVector3.h"
+
+#include <iostream>
+#include <random>
+#include <cmath>
+#include <math.h>
+#include <TMath.h>
+#include <TDatabasePDG.h>
+#include <TParticlePDG.h>
+
+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 = 1000,
+ 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
+ int useCrossingAngle=1 // 0= no rotation, 1 = -25 mrad
+ )
+{
+ 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<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
+ GenParticlePtr p2 = std::make_shared<GenParticle>(
+ 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.*useCrossingAngle; //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<GenParticle>(
+ 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<GenVertex>();
+ 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;
+}
diff --git a/benchmarks/zdc_photon/analysis/zdc_photon_plots.py b/benchmarks/zdc_photon/analysis/zdc_photon_plots.py
new file mode 100644
index 00000000..b06cf434
--- /dev/null
+++ b/benchmarks/zdc_photon/analysis/zdc_photon_plots.py
@@ -0,0 +1,161 @@
+import numpy as np, pandas as pd, matplotlib.pyplot as plt, matplotlib as mpl, awkward as ak, sys
+import mplhep as hep
+hep.style.use("CMS")
+
+plt.rcParams['figure.facecolor']='white'
+plt.rcParams['savefig.facecolor']='white'
+plt.rcParams['savefig.bbox']='tight'
+
+
+outdir=sys.argv[1]+"/"
+config=outdir.split("/")[1]
+try:
+ import os
+ os.mkdir(outdir[:-1])
+except:
+ pass
+
+def gauss(x, A,mu, sigma):
+ return A * np.exp(-(x-mu)**2/(2*sigma**2))
+
+#load the files
+import uproot as ur
+arrays_sim={}
+momenta=20, 30, 50, 70, 100, 150, 200, 275
+for p in momenta:
+ filename=f'sim_output/zdc_photon/{config}_rec_zdc_photon_{p}GeV.edm4hep.root'
+ print("opening file", filename)
+ events = ur.open(filename+':events')
+ arrays_sim[p] = events.arrays()#[:-1] #remove last event, which for some reason is blank
+ import gc
+ gc.collect()
+ print("read", filename)
+
+fig,axs=plt.subplots(1,3, figsize=(24, 8))
+pvals=[]
+resvals=[]
+dresvals=[]
+scalevals=[]
+dscalevals=[]
+for p in momenta:
+ selection=[len(arrays_sim[p]["HcalFarForwardZDCClusters.energy"][i])==1 for i in range(len(arrays_sim[p]))]
+ E=arrays_sim[p][selection]["HcalFarForwardZDCClusters.energy"]
+
+ Etot=np.sum(E, axis=-1)
+ #print(len(Etot))
+ #print(p, res, mrecon)
+ if p==100:
+ plt.sca(axs[0])
+ y, x, _=plt.hist(Etot, bins=100, range=(p*.75, p*1.25), histtype='step')
+ plt.ylabel("events")
+ plt.title(f"$p_{{\gamma}}$={p} GeV")
+ plt.xlabel("$E^{\\gamma}_{recon}$ [GeV]")
+ else:
+ y, x = np.histogram(Etot, bins=100, range=(p*.75, p*1.25))
+
+ bc=(x[1:]+x[:-1])/2
+ from scipy.optimize import curve_fit
+ slc=abs(bc-p)<10
+ fnc=gauss
+ p0=[100, p, 10]
+ #print(list(y), list(x))
+ coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,
+ sigma=list(np.sqrt(y[slc])+(y[slc]==0)))
+ if p==100:
+ xx=np.linspace(p*0.75,p*1.25, 100)
+ plt.plot(xx, fnc(xx,*coeff))
+ pvals.append(p)
+ resvals.append(np.abs(coeff[2])/coeff[1])
+ dresvals.append(np.sqrt(var_matrix[2][2])/coeff[1])
+ scalevals.append(np.abs(coeff[1])/p)
+ dscalevals.append(np.sqrt(var_matrix[2][2])/p)
+
+plt.sca(axs[1])
+plt.errorbar(pvals, resvals, dresvals, ls='', marker='o')
+plt.ylim(0)
+plt.ylabel("$\\sigma[E_{\\gamma}]/\\mu[E_{\\gamma}]$")
+plt.xlabel("$p_{\\gamma}$ [GeV]")
+
+xx=np.linspace(15, 275, 100)
+
+fnc=lambda E,a: a/np.sqrt(E)
+#pvals, resvals, dresvals
+coeff, var_matrix = curve_fit(fnc, pvals, resvals, p0=(1,),
+ sigma=dresvals)
+
+xx=np.linspace(15, 275, 100)
+plt.plot(xx, fnc(xx, *coeff), label=f'fit: $\\frac{{{coeff[0]*100:.0f}\\%}}{{\\sqrt{{E}}}}$')
+plt.legend()
+plt.sca(axs[2])
+plt.errorbar(pvals, scalevals, dscalevals, ls='', marker='o')
+plt.ylim(0.8, 1.2)
+plt.ylabel("$\\mu[E_{\\gamma}]/E_{\\gamma}$")
+plt.xlabel("$p_{\\gamma}$ [GeV]")
+plt.axhline(1, ls='--', alpha=0.7, color='0.5')
+plt.tight_layout()
+plt.savefig(outdir+"photon_energy_res.pdf")
+
+#theta res
+fig,axs=plt.subplots(1,2, figsize=(16, 8))
+pvals=[]
+resvals=[]
+dresvals=[]
+for p in momenta:
+ selection=[len(arrays_sim[p]["HcalFarForwardZDCClusters.energy"][i])==1 for i in range(len(arrays_sim[p]))]
+ x=arrays_sim[p][selection]["HcalFarForwardZDCClusters.position.x"][::,0]
+ y=arrays_sim[p][selection]["HcalFarForwardZDCClusters.position.y"][::,0]
+ z=arrays_sim[p][selection]["HcalFarForwardZDCClusters.position.z"][::,0]
+
+ theta_recon=np.arctan2(np.hypot(x*np.cos(-.025)-z*np.sin(-.025), y), z*np.cos(-.025)+x*np.sin(-.025))
+
+ px=arrays_sim[p][selection]["MCParticles.momentum.x"][::,2]
+ py=arrays_sim[p][selection]["MCParticles.momentum.y"][::,2]
+ pz=arrays_sim[p][selection]["MCParticles.momentum.z"][::,2]
+
+ theta_truth=np.arctan2(np.hypot(px*np.cos(-.025)-pz*np.sin(-.025), py), pz*np.cos(-.025)+px*np.sin(-.025))
+
+ Etot=np.sum(E, axis=-1)
+ #print(p, res, mrecon)
+ if p==100:
+ plt.sca(axs[0])
+ y, x, _=plt.hist(1000*(theta_recon-theta_truth), bins=100, range=(-0.5, 0.5), histtype='step')
+ plt.ylabel("events")
+ plt.title(f"$p_{{\gamma}}$={p} GeV")
+ plt.xlabel("$\\theta^{\\gamma}_{recon}$ [mrad]")
+ else:
+ y, x = np.histogram(1000*(theta_recon-theta_truth), bins=100, range=(-0.5, 0.5))
+
+ bc=(x[1:]+x[:-1])/2
+ from scipy.optimize import curve_fit
+ slc=abs(bc)<0.2#1.5*np.std(1000*(theta_recon-theta_truth))
+ fnc=gauss
+ p0=[100, 0, 0.1]
+ #print(list(y), list(x))
+ coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,
+ sigma=list(np.sqrt(y[slc])+(y[slc]==0)))
+ if p==100:
+ xx=np.linspace(-0.5,0.5, 100)
+ plt.plot(xx, fnc(xx,*coeff))
+ pvals.append(p)
+ resvals.append(np.abs(coeff[2]))
+ dresvals.append(np.sqrt(var_matrix[2][2]))
+plt.sca(axs[1])
+plt.errorbar(pvals, resvals, dresvals, ls='', marker='o')
+#print(dresvals)
+
+fnc=lambda E,a, b: np.hypot(a/np.sqrt(E), b)
+#pvals, resvals, dresvals
+coeff, var_matrix = curve_fit(fnc, pvals, resvals, p0=(1,.1),
+ sigma=dresvals)
+
+xx=np.linspace(15, 275, 100)
+
+plt.plot(xx, fnc(xx, *coeff), label=f'fit: $\\frac{{{coeff[0]:.2f}}}{{\\sqrt{{E}}}}\\oplus {coeff[1]:.3f}$ mrad')
+
+plt.ylabel("$\\sigma[\\theta_{\\gamma}]$ [mrad]")
+plt.xlabel("$p_{\\gamma}$ [GeV]")
+
+plt.ylim(0, 0.1)
+plt.legend()
+plt.tight_layout()
+plt.savefig(outdir+"photon_theta_res.pdf")
diff --git a/benchmarks/zdc_photon/config.yml b/benchmarks/zdc_photon/config.yml
new file mode 100644
index 00000000..c8e8bdb5
--- /dev/null
+++ b/benchmarks/zdc_photon/config.yml
@@ -0,0 +1,36 @@
+sim:zdc_photon:
+ stage: simulate
+ extends: .det_benchmark
+ needs: ["common:detector"]
+ parallel:
+ matrix:
+ - P: 20
+ - P: 30
+ - P: 50
+ - P: 70
+ - P: 100
+ - P: 150
+ - P: 200
+ - P: 275
+ timeout: 1 hours
+ script:
+ - snakemake --cores 1 sim_output/zdc_photon/epic_zdc_sipm_on_tile_only_rec_zdc_photon_${P}GeV.edm4hep.root
+ retry:
+ max: 2
+ when:
+ - runner_system_failure
+
+bench:zdc_photon:
+ stage: benchmarks
+ extends: .det_benchmark
+ needs: ["sim:zdc_photon"]
+ script:
+ - mkdir -p results/epic_zdc_sipm_on_tile_only
+ - python benchmarks/zdc_photon/analysis/zdc_photon_plots.py results/epic_zdc_sipm_on_tile_only/zdc_photon
+
+collect_results:zdc_photon:
+ stage: collect
+ extends: .det_benchmark
+ needs: ["bench:zdc_photon"]
+ script:
+ - ls -al
diff --git a/benchmarks/zdc_pi0/Snakefile b/benchmarks/zdc_pi0/Snakefile
new file mode 100644
index 00000000..57fd2c90
--- /dev/null
+++ b/benchmarks/zdc_pi0/Snakefile
@@ -0,0 +1,56 @@
+rule zdc_pi0_generate:
+ input:
+ script="benchmarks/zdc_pi0/analysis/gen_pi0_decay.cxx",
+ params:
+ NEVENTS_GEN=1000,
+ output:
+ GEN_FILE="sim_output/zdc_pi0/zdc_pi0_{P}GeV.hepmc"
+ shell:
+ """
+mkdir -p sim_output/zdc_pi0
+root -l -b -q '{input.script}({params.NEVENTS_GEN},0,"{output.GEN_FILE}",{wildcards.P},{wildcards.P})'
+"""
+
+rule zdc_pi0_simulate:
+ input:
+ GEN_FILE="sim_output/zdc_pi0/zdc_pi0_{P}GeV.hepmc"
+ params:
+ PHYSICS_LIST="FTFP_BERT"
+ output:
+ SIM_FILE="sim_output/zdc_pi0/{DETECTOR_CONFIG}_sim_zdc_pi0_{P}GeV.edm4hep.root"
+ shell:
+ """
+NEVENTS_SIM=1000
+# Running simulation
+npsim \
+ --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \
+ --numberOfEvents $NEVENTS_SIM \
+ --physicsList {params.PHYSICS_LIST} \
+ --inputFiles {input.GEN_FILE} \
+ --outputFile {output.SIM_FILE}
+"""
+
+rule zdc_pi0_recon:
+ input:
+ SIM_FILE="sim_output/zdc_pi0/{DETECTOR_CONFIG}_sim_zdc_pi0_{P}GeV.edm4hep.root"
+ output:
+ REC_FILE="sim_output/zdc_pi0/{DETECTOR_CONFIG}_rec_zdc_pi0_{P}GeV.edm4hep.root"
+ shell:
+ """
+NEVENTS_REC=1000
+eicrecon {input.SIM_FILE} -Ppodio:output_file={output.REC_FILE} -Pdd4hep:xml_files=$DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml -Ppodio:output_collections=MCParticles,HcalFarForwardZDCRecHits,HcalFarForwardZDCClusters,HcalFarForwardZDCSubcellHits -Pjana:nevents=$NEVENTS_REC
+"""
+
+rule zdc_pi0_analysis:
+ input:
+ expand("sim_output/zdc_pi0/{DETECTOR_CONFIG}_rec_zdc_pi0_{P}GeV.edm4hep.root",
+ P=[60, 80, 100, 130, 160],
+ DETECTOR_CONFIG=["{DETECTOR_CONFIG}"]),
+ script="benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py",
+ output:
+ results_dir=directory("results/{DETECTOR_CONFIG}/zdc_pi0"),
+ shell:
+ """
+mkdir -p {output.results_dir}
+python {input.script} {output.results_dir}
+"""
diff --git a/benchmarks/zdc_pi0/analysis/gen_pi0_decay.cxx b/benchmarks/zdc_pi0/analysis/gen_pi0_decay.cxx
new file mode 100644
index 00000000..13b3d807
--- /dev/null
+++ b/benchmarks/zdc_pi0/analysis/gen_pi0_decay.cxx
@@ -0,0 +1,182 @@
+#include "HepMC3/GenEvent.h"
+#include "HepMC3/ReaderAscii.h"
+#include "HepMC3/WriterAscii.h"
+#include "HepMC3/Print.h"
+
+#include "TRandom3.h"
+#include "TVector3.h"
+
+#include <TDatabasePDG.h>
+#include <TParticlePDG.h>
+
+#include <iostream>
+#include <random>
+#include <TMath.h>
+
+using namespace HepMC3;
+
+std::tuple<double, int, double> GetParticleInfo(TDatabasePDG* pdg, TString particle_name)
+{
+ TParticlePDG *particle = pdg->GetParticle(particle_name);
+ const double mass = particle->Mass();
+ const int pdgID = particle->PdgCode();
+ const double lifetime = particle->Lifetime();
+ return std::make_tuple(mass, pdgID, lifetime);
+}
+// Calculates the decay length of a particle. Samples from an exponential decay.
+double GetDecayLength(TRandom3* r1, double lifetime, double mass, double momentum_magnitude)
+{
+ double c_speed = TMath::C() * 1000.; // speed of light im mm/sec
+ double average_decay_length = (momentum_magnitude/mass) * lifetime * c_speed;
+ return r1->Exp(average_decay_length);
+}
+
+// Generate single pi0 mesons and decay them to 2 photons
+void gen_pi0_decay(int n_events = 100000, UInt_t seed = 0, char* out_fname = "pi0_decay.hepmc",
+ double p_min = 60., // in GeV/c
+ double p_max = 275.) // in GeV/c
+{
+
+ const double theta_min = 0.0; // in mRad
+ const double theta_max = 4; // in mRad
+ //const double p_min = 100.; // in GeV/c
+ //const double p_max = 275.; // in GeV/c
+
+ WriterAscii hepmc_output(out_fname);
+ int events_parsed = 0;
+ GenEvent evt(Units::GEV, Units::MM);
+
+ // Random number generator
+ TRandom3 *r1 = new TRandom3(seed); //Default = 0, which uses clock to set seed
+ cout<<"Random number seed is "<<r1->GetSeed()<<"!"<<endl;
+
+ // Getting generated particle information
+ TDatabasePDG *pdg = new TDatabasePDG();
+
+ auto pi0_info = GetParticleInfo(pdg, "pi0");
+ double pi0_mass = std::get<0>(pi0_info);
+ int pi0_pdgID = std::get<1>(pi0_info);
+ double pi0_lifetime = std::get<2>(pi0_info);
+
+ auto photon_info = GetParticleInfo(pdg, "gamma");
+ double photon_mass = std::get<0>(photon_info);
+ int photon_pdgID = std::get<1>(photon_info);
+
+ 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 2212 - proton
+ GenParticlePtr p1 =
+ std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
+ GenParticlePtr p2 = std::make_shared<GenParticle>(
+ FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
+
+ // Define momentum with respect to EIC proton beam direction
+ Double_t pi0_p = r1->Uniform(p_min, p_max);
+ Double_t pi0_phi = r1->Uniform(0.0, 2.0 * M_PI);
+ Double_t pi0_th = r1->Uniform(theta_min/1000., theta_max/1000.); // Divide by 1000 for radians
+ Double_t pi0_px = pi0_p * TMath::Cos(pi0_phi) * TMath::Sin(pi0_th);
+ Double_t pi0_py = pi0_p * TMath::Sin(pi0_phi) * TMath::Sin(pi0_th);
+ Double_t pi0_pz = pi0_p * TMath::Cos(pi0_th);
+ Double_t pi0_E = TMath::Sqrt(pi0_p*pi0_p + pi0_mass*pi0_mass);
+
+ // Rotate to lab coordinate system
+ TVector3 pi0_pvec(pi0_px, pi0_py, pi0_pz);
+ double cross_angle = -25./1000.; // in Rad
+ TVector3 pbeam_dir(TMath::Sin(cross_angle), 0, TMath::Cos(cross_angle)); //proton beam direction
+ pi0_pvec.RotateY(cross_angle); // Theta is returned positive, beam in negative X
+
+ // type 2 is state that will decay
+ GenParticlePtr p_pi0 = std::make_shared<GenParticle>(
+ FourVector(pi0_pvec.X(), pi0_pvec.Y(), pi0_pvec.Z(), pi0_E), pi0_pdgID, 2 );
+
+ // Generating pi0 particle, will be generated at origin
+ // Must have input electron + proton for vertex
+ GenVertexPtr pi0_initial_vertex = std::make_shared<GenVertex>();
+ pi0_initial_vertex->add_particle_in(p1);
+ pi0_initial_vertex->add_particle_in(p2);
+ pi0_initial_vertex->add_particle_out(p_pi0);
+ evt.add_vertex(pi0_initial_vertex);
+
+ // Generate gamma1 + gamma1 in pi0 rest frame
+ TLorentzVector gamma1_rest, gamma2_rest;
+
+ // Generating uniformly along a sphere
+ double cost_gamma1_rest = r1->Uniform(-1,1);
+ double th_gamma1_rest = TMath::ACos(cost_gamma1_rest);
+ double sint_gamma1_rest = TMath::Sin(th_gamma1_rest);
+
+ double phi_gamma1_rest = r1->Uniform(-1.*TMath::Pi(),1.*TMath::Pi());
+ double cosp_gamma1_rest = TMath::Cos(phi_gamma1_rest);
+ double sinp_gamma1_rest = TMath::Sin(phi_gamma1_rest);
+
+ // Calculate energy of each particle in the pi0 rest frame
+ // This is half the pi0 mass
+ double E_gamma1_rest = pi0_mass/2;
+ double E_gamma2_rest = pi0_mass/2;
+
+ // Both particles will have the same momentum, and are massless
+ double momentum_rest = pi0_mass/2;
+
+ gamma1_rest.SetE(E_gamma1_rest);
+ gamma1_rest.SetPx( momentum_rest * sint_gamma1_rest * cosp_gamma1_rest );
+ gamma1_rest.SetPy( momentum_rest * sint_gamma1_rest * sinp_gamma1_rest );
+ gamma1_rest.SetPz( momentum_rest * cost_gamma1_rest );
+
+ gamma2_rest.SetE(E_gamma2_rest);
+ gamma2_rest.SetPx( -gamma1_rest.Px() );
+ gamma2_rest.SetPy( -gamma1_rest.Py() );
+ gamma2_rest.SetPz( -gamma1_rest.Pz() );
+
+ // Boost both gammas to lab frame
+ TLorentzVector pi0_lab(pi0_pvec.X(), pi0_pvec.Y(), pi0_pvec.Z(), pi0_E);
+ TVector3 pi0_boost = pi0_lab.BoostVector();
+ TLorentzVector gamma1_lab, gamma2_lab;
+ gamma1_lab = gamma1_rest;
+ gamma1_lab.Boost(pi0_boost);
+ gamma2_lab = gamma2_rest;
+ gamma2_lab.Boost(pi0_boost);
+
+ // Calculating position for pi0 decay
+ TVector3 pi0_unit = pi0_lab.Vect().Unit();
+ double pi0_decay_length = GetDecayLength(r1, pi0_lifetime, pi0_mass, pi0_lab.P());
+ TVector3 pi0_decay_position = pi0_unit * pi0_decay_length;
+ double pi0_decay_time = pi0_decay_length / pi0_lab.Beta() ; // Decay time in lab frame in length units (mm)
+
+ // Generating vertex for pi0 decay
+ GenParticlePtr p_gamma1 = std::make_shared<GenParticle>(
+ FourVector(gamma1_lab.Px(), gamma1_lab.Py(), gamma1_lab.Pz(), gamma1_lab.E()), photon_pdgID, 1 );
+
+ GenParticlePtr p_gamma2 = std::make_shared<GenParticle>(
+ FourVector(gamma2_lab.Px(), gamma2_lab.Py(), gamma2_lab.Pz(), gamma2_lab.E()), photon_pdgID, 1 );
+
+ GenVertexPtr v_pi0_decay = std::make_shared<GenVertex>(FourVector(pi0_decay_position.X(), pi0_decay_position.Y(), pi0_decay_position.Z(), pi0_decay_time));
+ v_pi0_decay->add_particle_in(p_pi0);
+ v_pi0_decay->add_particle_out(p_gamma1);
+ v_pi0_decay->add_particle_out(p_gamma2);
+
+ evt.add_vertex(v_pi0_decay);
+
+ if (events_parsed == 0) {
+ std::cout << "First event: " << std::endl;
+ Print::listing(evt);
+ }
+ double zdc_z=35800;
+ TVector3 extrap_gamma1=pi0_decay_position+gamma1_lab.Vect()*((zdc_z-pbeam_dir.Dot(pi0_decay_position))/(pbeam_dir.Dot(gamma1_lab.Vect())));
+ TVector3 extrap_gamma2=pi0_decay_position+gamma2_lab.Vect()*((zdc_z-pbeam_dir.Dot(pi0_decay_position))/(pbeam_dir.Dot(gamma2_lab.Vect())));
+ if (extrap_gamma1.Angle(pbeam_dir)<0.004 && extrap_gamma2.Angle(pbeam_dir)<0.004 && pi0_decay_position.Dot(pbeam_dir)<zdc_z)
+ hepmc_output.write_event(evt);
+ if (events_parsed % 1000 == 0) {
+ std::cout << "Event: " << events_parsed << std::endl;
+ }
+ evt.clear();
+ }
+ hepmc_output.close();
+
+ std::cout << "Events parsed and written: " << events_parsed << std::endl;
+}
diff --git a/benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py b/benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py
new file mode 100644
index 00000000..92adac55
--- /dev/null
+++ b/benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py
@@ -0,0 +1,232 @@
+import numpy as np, pandas as pd, matplotlib.pyplot as plt, matplotlib as mpl, awkward as ak, sys
+import mplhep as hep
+hep.style.use("CMS")
+
+plt.rcParams['figure.facecolor']='white'
+plt.rcParams['savefig.facecolor']='white'
+plt.rcParams['savefig.bbox']='tight'
+
+outdir=sys.argv[1]+"/"
+config=outdir.split("/")[1]
+try:
+ import os
+ os.mkdir(outdir[:-1])
+except:
+ pass
+
+def gauss(x, A,mu, sigma):
+ return A * np.exp(-(x-mu)**2/(2*sigma**2))
+
+import uproot as ur
+arrays_sim={}
+momenta=60, 80, 100, 130, 160,
+for p in momenta:
+ filename=f'sim_output/zdc_pi0/{config}_rec_zdc_pi0_{p}GeV.edm4hep.root'
+ print("opening file", filename)
+ events = ur.open(filename+':events')
+ arrays_sim[p] = events.arrays()
+ import gc
+ gc.collect()
+ print("read", filename)
+
+#energy res plot
+fig,axs=plt.subplots(1,3, figsize=(24, 8))
+pvals=[]
+resvals=[]
+dresvals=[]
+scalevals=[]
+dscalevals=[]
+for p in momenta:
+ selection=[len(arrays_sim[p]["HcalFarForwardZDCClusters.energy"][i])==2 for i in range(len(arrays_sim[p]))]
+ E=arrays_sim[p][selection]["HcalFarForwardZDCClusters.energy"]
+
+ Etot=np.sum(E, axis=-1)
+ if len(Etot)<25:
+ continue
+ #print(p, res, mrecon)
+ if p==100:
+ plt.sca(axs[0])
+ y, x, _=plt.hist(Etot, bins=100, range=(p*.5, p*1.5), histtype='step')
+ plt.ylabel("events")
+ plt.title(f"$p_{{\pi^0}}$={p} GeV")
+ plt.xlabel("$E^{\\pi^{0}}_{recon}$ [GeV]")
+ else:
+ y, x = np.histogram(Etot, bins=100, range=(p*.5, p*1.5))
+
+ bc=(x[1:]+x[:-1])/2
+ from scipy.optimize import curve_fit
+ slc=abs(bc-p)<10
+ fnc=gauss
+ p0=[100, p, 10]
+ #print(list(y), list(x))
+ coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,
+ sigma=list(np.sqrt(y[slc])+(y[slc]==0)))
+ if p==100:
+ xx=np.linspace(p*0.5,p*1.5, 100)
+ plt.plot(xx, fnc(xx,*coeff))
+ pvals.append(p)
+ resvals.append(np.abs(coeff[2])/coeff[1])
+ dresvals.append(np.sqrt(var_matrix[2][2])/coeff[1])
+ scalevals.append(np.abs(coeff[1])/p)
+ dscalevals.append(np.sqrt(var_matrix[2][2])/p)
+
+plt.sca(axs[1])
+plt.errorbar(pvals, resvals, dresvals, ls='', marker='o')
+
+plt.ylabel("$\\sigma[E_{\\pi^0}]/\\mu[E_{\\pi^0}]$")
+plt.xlabel("$p_{\\pi^0}$ [GeV]")
+
+fnc=lambda E,a: a/np.sqrt(E)
+#pvals, resvals, dresvals
+coeff, var_matrix = curve_fit(fnc, pvals, resvals, p0=(1,),
+ sigma=dresvals)
+xx=np.linspace(55, 200, 100)
+plt.plot(xx, fnc(xx, *coeff), label=f'fit: $\\frac{{{coeff[0]:.2f}\\%}}{{\\sqrt{{E}}}}$')
+plt.legend()
+plt.ylim(0)
+plt.sca(axs[2])
+plt.errorbar(pvals, scalevals, dscalevals, ls='', marker='o')
+plt.ylim(0.8, 1.2)
+plt.ylabel("$\\mu[E_{\\pi^0}]/E_{\\pi^0}$")
+plt.xlabel("$p_{\\pi^0}$ [GeV]")
+plt.axhline(1, ls='--', alpha=0.7, color='0.5')
+plt.tight_layout()
+plt.savefig(outdir+"/pi0_energy_res.pdf")
+
+
+fig,axs=plt.subplots(1,2, figsize=(16, 8))
+pvals=[]
+resvals=[]
+dresvals=[]
+for p in momenta:
+ selection=[len(arrays_sim[p]["HcalFarForwardZDCClusters.energy"][i])==2 for i in range(len(arrays_sim[p]))]
+ x=arrays_sim[p][selection]["HcalFarForwardZDCClusters.position.x"]
+ y=arrays_sim[p][selection]["HcalFarForwardZDCClusters.position.y"]
+ z=arrays_sim[p][selection]["HcalFarForwardZDCClusters.position.z"]
+ E=arrays_sim[p][selection]["HcalFarForwardZDCClusters.energy"]
+ r=np.sqrt(x**2+y**2+z**2)
+ px=np.sum(E*x/r, axis=-1)
+ py=np.sum(E*y/r, axis=-1)
+ pz=np.sum(E*z/r, axis=-1)
+
+ theta_recon=np.arctan2(np.hypot(px*np.cos(-.025)-pz*np.sin(-.025), py), pz*np.cos(-.025)+px*np.sin(-.025))
+ if len(theta_recon)<25:
+ continue
+ px=arrays_sim[p][selection]["MCParticles.momentum.x"][::,2]
+ py=arrays_sim[p][selection]["MCParticles.momentum.y"][::,2]
+ pz=arrays_sim[p][selection]["MCParticles.momentum.z"][::,2]
+
+ theta_truth=np.arctan2(np.hypot(px*np.cos(-.025)-pz*np.sin(-.025), py), pz*np.cos(-.025)+px*np.sin(-.025))
+
+ Etot=np.sum(E, axis=-1)
+ #print(p, res, mrecon)
+ if p==100:
+ plt.sca(axs[0])
+ y, x, _=plt.hist(1000*(theta_recon-theta_truth), bins=100, range=(-0.5, 0.5), histtype='step')
+ plt.ylabel("events")
+ plt.title(f"$p_{{\\pi^0}}$={p} GeV")
+ plt.xlabel("$\\theta^{\\pi^0}_{recon}$ [mrad]")
+ else:
+ y, x = np.histogram(1000*(theta_recon-theta_truth), bins=100, range=(-0.5, 0.5))
+
+ bc=(x[1:]+x[:-1])/2
+ from scipy.optimize import curve_fit
+ slc=abs(bc)<0.2#1.5*np.std(1000*(theta_recon-theta_truth))
+ fnc=gauss
+ p0=[100, 0, 0.1]
+ #print(list(y), list(x))
+ coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,
+ sigma=list(np.sqrt(y[slc])+(y[slc]==0)))
+ if p==100:
+ xx=np.linspace(-0.5,0.5, 100)
+ plt.plot(xx, fnc(xx,*coeff))
+ pvals.append(p)
+ resvals.append(np.abs(coeff[2]))
+ dresvals.append(np.sqrt(var_matrix[2][2]))
+
+plt.sca(axs[1])
+plt.errorbar(pvals, resvals, dresvals, ls='', marker='o')
+#print(dresvals)
+
+fnc=lambda E,a: a/np.sqrt(E)
+#pvals, resvals, dresvals
+coeff, var_matrix = curve_fit(fnc, pvals, resvals, p0=(1,),
+ sigma=dresvals)
+
+xx=np.linspace(55, 200, 100)
+
+plt.plot(xx, fnc(xx, *coeff), label=f'fit: $\\frac{{{coeff[0]:.2f}}}{{\\sqrt{{E}}}}$ mrad')
+
+plt.ylabel("$\\sigma[\\theta_{\\pi^0}]$ [mrad]")
+plt.xlabel("$p_{\\pi^0}$ [GeV]")
+
+plt.ylim(0, 0.1)
+plt.legend()
+plt.tight_layout()
+plt.savefig(outdir+"/pi0_theta_res.pdf")
+
+fig,axs=plt.subplots(1,2, figsize=(16, 8))
+pvals=[]
+resvals=[]
+dresvals=[]
+for p in momenta:
+ selection=[len(arrays_sim[p]["HcalFarForwardZDCClusters.energy"][i])==2 for i in range(len(arrays_sim[p]))]
+ E=arrays_sim[p][selection]["HcalFarForwardZDCClusters.energy"]
+ cx=arrays_sim[p][selection]["HcalFarForwardZDCClusters.position.x"]
+ cy=arrays_sim[p][selection]["HcalFarForwardZDCClusters.position.y"]
+ cz=arrays_sim[p][selection]["HcalFarForwardZDCClusters.position.z"]
+ r=np.sqrt(cx**2+cy**2+cz**2)
+ px=E*cx/r
+ py=E*cy/r
+ pz=E*cz/r
+
+ cos_opening_angle=(cx/r)[::,0]*(cx/r)[::,1]+(cy/r)[::,0]*(cy/r)[::,1]+(cz/r)[::,0]*(cz/r)[::,1]
+ mrecon=np.sqrt(2*E[::,0]*E[::,1]*(1-cos_opening_angle))
+
+ if len(mrecon)<25:
+ continue
+
+ #print(p, res, mrecon)
+ if p==100:
+ plt.sca(axs[0])
+ y, x, _=plt.hist(mrecon, bins=100, range=(0, 0.2), histtype='step')
+ plt.ylabel("events")
+ plt.title(f"$p_{{\pi^0}}$={p} GeV")
+ plt.xlabel("$m^{\\pi^{0}}_{recon}$ [GeV]")
+ else:
+ #y, x, _=plt.hist(mrecon, bins=100, range=(0, 0.2), histtype='step')#y, x =np.histogram(mrecon, bins=100, range=(0, 0.2))
+ y, x = np.histogram(mrecon, bins=100, range=(0, 0.2))
+
+ bc=(x[1:]+x[:-1])/2
+ from scipy.optimize import curve_fit
+ slc=abs(bc-.135)<.1
+ fnc=gauss
+ p0=[100, .135, 0.2]
+ #print(list(y), list(x))
+ coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,
+ sigma=list(np.sqrt(y[slc])+(y[slc]==0)))
+ if p==100:
+ xx=np.linspace(0,0.2)
+ plt.plot(xx, fnc(xx,*coeff))
+ pvals.append(p)
+ resvals.append(np.abs(coeff[2]))
+ dresvals.append(np.sqrt(var_matrix[2][2]))
+
+plt.sca(axs[1])
+plt.errorbar(pvals, resvals, dresvals, ls='', marker='o')
+plt.ylim(0)
+plt.ylabel("$\\sigma[m_{\\pi^0}]$ [GeV]")
+plt.xlabel("$p_{\\pi^0}$ [GeV]")
+
+fnc=lambda E,a,b: a+b*E
+#pvals, resvals, dresvals
+coeff, var_matrix = curve_fit(fnc, pvals, resvals, p0=(1,1),
+ sigma=dresvals)
+xx=np.linspace(55, 200, 100)
+#plt.plot(xx, fnc(xx, *coeff), label=f'fit: ${coeff[0]*1000:.1f}+{coeff[1]*1000:.4f}\\times E$ MeV')
+plt.plot(xx, fnc(xx, *coeff), label=f'fit: $({coeff[0]*1000:.1f}+{coeff[1]*1000:.4f}\\times [E\,in\,GeV])$ MeV')
+plt.legend()
+
+
+plt.tight_layout()
+plt.savefig(outdir+"/pi0_mass_res.pdf")
diff --git a/benchmarks/zdc_pi0/config.yml b/benchmarks/zdc_pi0/config.yml
new file mode 100644
index 00000000..4059d236
--- /dev/null
+++ b/benchmarks/zdc_pi0/config.yml
@@ -0,0 +1,34 @@
+sim:zdc_pi0:
+ stage: simulate
+ extends: .det_benchmark
+ needs: ["common:detector"]
+ parallel:
+ matrix:
+ - P: 60
+ - P: 80
+ - P: 100
+ - P: 130
+ - P: 160
+ timeout: 1 hours
+ script:
+ - snakemake --cores 1 sim_output/zdc_pi0/epic_zdc_sipm_on_tile_only_rec_zdc_pi0_${P}GeV.edm4hep.root
+ retry:
+ max: 2
+ when:
+ - runner_system_failure
+
+bench:zdc_pi0:
+ stage: benchmarks
+
+ extends: .det_benchmark
+ needs: ["sim:zdc_pi0"]
+ script:
+ - mkdir -p results/epic_zdc_sipm_on_tile_only
+ - python benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py results/epic_zdc_sipm_on_tile_only/zdc_pi0
+
+collect_results:zdc_pi0:
+ stage: collect
+ extends: .det_benchmark
+ needs: ["bench:zdc_pi0"]
+ script:
+ - ls -al
--
GitLab