From 6c895b6d8f12203c96be5695c382b22a7e8270f0 Mon Sep 17 00:00:00 2001
From: Sebouh Paul <sebouh.paul@gmail.com>
Date: Tue, 5 Nov 2024 10:13:39 -0500
Subject: [PATCH] Insert tau (#95)
Co-authored-by: Dmitry Kalinkin <dmitry.kalinkin@gmail.com>
---
.gitlab-ci.yml | 2 +
Snakefile | 3 +-
benchmarks/insert_tau/Snakefile | 68 ++++++++
.../insert_tau/analysis/gen_particles.cxx | 127 +++++++++++++++
benchmarks/insert_tau/analysis/tau_plots.py | 154 ++++++++++++++++++
benchmarks/insert_tau/config.yml | 36 ++++
6 files changed, 389 insertions(+), 1 deletion(-)
create mode 100644 benchmarks/insert_tau/Snakefile
create mode 100644 benchmarks/insert_tau/analysis/gen_particles.cxx
create mode 100644 benchmarks/insert_tau/analysis/tau_plots.py
create mode 100644 benchmarks/insert_tau/config.yml
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 47d126fd..1df546a9 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -143,6 +143,7 @@ include:
- local: 'benchmarks/timing/config.yml'
- local: 'benchmarks/b0_tracker/config.yml'
- local: 'benchmarks/insert_muon/config.yml'
+ - local: 'benchmarks/insert_tau/config.yml'
- local: 'benchmarks/zdc_sigma/config.yml'
- local: 'benchmarks/zdc_lambda/config.yml'
- local: 'benchmarks/insert_neutron/config.yml'
@@ -170,6 +171,7 @@ deploy_results:
- "collect_results:zdc"
- "collect_results:zdc_lyso"
- "collect_results:insert_muon"
+ - "collect_results:insert_tau"
- "collect_results:zdc_photon"
- "collect_results:zdc_pi0"
- "collect_results:femc_electron"
diff --git a/Snakefile b/Snakefile
index 1390e867..9c125483 100644
--- a/Snakefile
+++ b/Snakefile
@@ -8,13 +8,14 @@ include: "benchmarks/material_scan/Snakefile"
include: "benchmarks/tracking_performances/Snakefile"
include: "benchmarks/tracking_performances_dis/Snakefile"
include: "benchmarks/lfhcal/Snakefile"
+include: "benchmarks/zdc_lyso/Snakefile"
include: "benchmarks/insert_muon/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"
+include: "benchmarks/insert_tau/Snakefile"
include: "benchmarks/femc_electron/Snakefile"
include: "benchmarks/femc_photon/Snakefile"
include: "benchmarks/femc_pi0/Snakefile"
diff --git a/benchmarks/insert_tau/Snakefile b/benchmarks/insert_tau/Snakefile
new file mode 100644
index 00000000..57e4390e
--- /dev/null
+++ b/benchmarks/insert_tau/Snakefile
@@ -0,0 +1,68 @@
+def get_n_events(wildcards):
+ energy = float(wildcards.P)
+ n_events = 1000
+ n_events = int(n_events // ((energy / 20) ** 0.5))
+ return n_events
+
+rule insert_tau_generate:
+ input:
+ script="benchmarks/insert_tau/analysis/gen_particles.cxx",
+ params:
+ N_EVENTS=get_n_events,
+ th_max=7.0,
+ th_min=1.7,
+ output:
+ GEN_FILE="sim_output/insert_tau/tau-_{P}GeV.hepmc",
+ shell:
+ """
+root -l -b -q '{input.script}({params.N_EVENTS},"{output.GEN_FILE}", "tau-", {params.th_min}, {params.th_max}, 0., 360., {wildcards.P})'
+"""
+
+rule insert_tau_simulate:
+ input:
+ GEN_FILE="sim_output/insert_tau/tau-_{P}GeV.hepmc",
+ warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root",
+ params:
+ N_EVENTS=get_n_events,
+ PHYSICS_LIST="FTFP_BERT",
+ output:
+ SIM_FILE="sim_output/insert_tau/{DETECTOR_CONFIG}_sim_tau-_{P}GeV_{INDEX}.edm4hep.root"
+ shell:
+ """
+# Running simulation
+npsim \
+ --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \
+ --numberOfEvents {params.N_EVENTS} \
+ --skipNEvents $(( {params.N_EVENTS} * {wildcards.INDEX} )) \
+ --physicsList {params.PHYSICS_LIST} \
+ --inputFiles {input.GEN_FILE} \
+ --outputFile {output.SIM_FILE}
+"""
+
+rule insert_tau_recon:
+ input:
+ SIM_FILE="sim_output/insert_tau/{DETECTOR_CONFIG}_sim_tau-_{P}GeV_{INDEX}.edm4hep.root",
+ params:
+ N_EVENTS=get_n_events,
+ output:
+ REC_FILE="sim_output/insert_tau/{DETECTOR_CONFIG}_rec_tau-_{P}GeV_{INDEX}.edm4eic.root",
+ shell:
+ """
+eicrecon {input.SIM_FILE} -Ppodio:output_file={output.REC_FILE} -Pdd4hep:xml_files=$DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml -Ppodio:output_collections=MCParticles,HcalEndcapPInsertRecHits,HcalEndcapPInsertClusters,HcalEndcapPInsertSubcellHits,EcalEndcapPInsertRecHits,EcalEndcapPInsertClusters,EcalEndcapPClusters,LFHCALClusters -Pjana:nevents={params.N_EVENTS}
+"""
+
+rule insert_tau_analysis:
+ input:
+ expand("sim_output/insert_tau/{DETECTOR_CONFIG}_rec_tau-_{P}GeV_{INDEX}.edm4eic.root",
+ P=[20, 30, 40, 50, 60, 80, 100],
+ DETECTOR_CONFIG=["{DETECTOR_CONFIG}"],
+ INDEX=range(5),
+ ),
+ script="benchmarks/insert_tau/analysis/tau_plots.py",
+ output:
+ results_dir=directory("results/{DETECTOR_CONFIG}/insert_tau"),
+ shell:
+ """
+mkdir -p {output.results_dir}
+python {input.script} {output.results_dir}
+"""
diff --git a/benchmarks/insert_tau/analysis/gen_particles.cxx b/benchmarks/insert_tau/analysis/gen_particles.cxx
new file mode 100644
index 00000000..0877a852
--- /dev/null
+++ b/benchmarks/insert_tau/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/insert_tau/analysis/tau_plots.py b/benchmarks/insert_tau/analysis/tau_plots.py
new file mode 100644
index 00000000..77c2dd32
--- /dev/null
+++ b/benchmarks/insert_tau/analysis/tau_plots.py
@@ -0,0 +1,154 @@
+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'
+
+plt.rcParams["figure.figsize"] = (7, 7)
+
+config=sys.argv[1].split("/")[1] #results/{config}/{benchmark_name}
+outdir=sys.argv[1]+"/"
+try:
+ import os
+ os.mkdir(outdir[:-1])
+except:
+ pass
+
+import uproot as ur
+arrays_sim={}
+momenta=20, 30, 40, 50, 60,80,100
+for p in momenta:
+ arrays_sim[p] = ur.concatenate({
+ f'sim_output/insert_tau/{config}_rec_tau-_{p}GeV_{index}.edm4eic.root': 'events'
+ for index in range(5)
+ })
+
+
+for a in arrays_sim.values():
+ #recon
+ Etot=0
+ px=0
+ py=0
+ pz=0
+ for det in "HcalEndcapPInsert", "EcalEndcapPInsert", "EcalEndcapP", "LFHCAL":
+
+ E=a[f'{det}Clusters.energy']
+
+ #todo apply corrections depending on whether this is an electromagnetic or hadronic shower.
+
+ x=a[f'{det}Clusters.position.x']
+ y=a[f'{det}Clusters.position.y']
+ z=a[f'{det}Clusters.position.z']
+ r=np.sqrt(x**2+y**2+z**2)
+ Etot=Etot+np.sum(E, axis=-1)
+ px=px+np.sum(E*x/r,axis=-1)
+ py=py+np.sum(E*y/r,axis=-1)
+ pz=pz+np.sum(E*z/r,axis=-1)
+
+ a['jet_p_recon']=np.sqrt(px**2+py**2+pz**2)
+ a['jet_E_recon']=Etot
+
+ a['jet_theta_recon']=np.arctan2(np.hypot(px*np.cos(-.025)-pz*np.sin(-.025),py),
+ pz*np.cos(-.025)+px*np.sin(-.025))
+
+ #truth
+ m=a['MCParticles.mass'][::,2:]
+ px=a['MCParticles.momentum.x'][::,2:]
+ py=a['MCParticles.momentum.y'][::,2:]
+ pz=a['MCParticles.momentum.z'][::,2:]
+ E=np.sqrt(m**2+px**2+py**2+pz**2)
+ status=a['MCParticles.simulatorStatus'][::,2:]
+ PDG=a['MCParticles.PDG'][::,2:]
+
+ #find the hadronic part: initial-state tau - all leptons
+ selection=1*(PDG==15)-1*(np.abs(PDG)==16)
+ is_hadronic=1*(np.sum((PDG==-14)+(PDG==-12), axis=-1)==0)
+
+ px_hfs, py_hfs, pz_hfs= np.sum(px*selection,axis=-1)*is_hadronic,np.sum(py*selection,axis=-1)*is_hadronic, np.sum(pz*selection,axis=-1)*is_hadronic
+
+ a['hfs_p_truth']=np.sqrt(px_hfs**2+py_hfs**2+pz_hfs**2)
+ a['hfs_E_truth']=np.sum(E*selection,axis=-1)*is_hadronic
+
+
+ a['hfs_theta_truth']=np.arctan2(np.hypot(px_hfs*np.cos(-.025)-pz_hfs*np.sin(-.025),py_hfs),
+ pz_hfs*np.cos(-.025)+px_hfs*np.sin(-.025))
+ a['hfs_eta_truth']=-np.log(np.tan(a['hfs_theta_truth']/2))
+ a['n_mu']=np.sum(np.abs(PDG)==13, axis=-1)
+ a['n_e']=np.sum(np.abs(PDG)==13, axis=-1)
+ a['hfs_mass_truth']=np.sqrt(a['hfs_E_truth']**2-a['hfs_p_truth']**2)
+
+for a in arrays_sim.values():
+ selection=(a['hfs_eta_truth']>3.1) & (a['hfs_eta_truth']<3.8)\
+ &(a['n_mu']==0)&(a['n_e']==0)&(a['hfs_mass_truth']>.140)&(a['jet_E_recon']>0)
+
+Etruth=[]
+Erecon=[]
+
+theta_truth=[]
+theta_recon=[]
+
+eta_max=3.7
+eta_min=3.3
+for a in arrays_sim.values():
+ selection=(a['hfs_eta_truth']>eta_min) & (a['hfs_eta_truth']<eta_max)\
+ &(a['n_mu']==0)&(a['n_e']==0)&(a['hfs_mass_truth']>.140)&(a['jet_E_recon']>1)
+ theta_truth=np.concatenate((theta_truth,a['hfs_theta_truth'][selection]))
+ theta_recon=np.concatenate((theta_recon,a['jet_theta_recon'][selection]))
+ Etruth=np.concatenate((Etruth,a['hfs_E_truth'][selection]))
+ Erecon=np.concatenate((Erecon,a['jet_E_recon'][selection]))
+
+plt.figure()
+plt.scatter(theta_truth, theta_recon, 1)
+plt.xlabel("$\\theta^{hfs}_{\\rm truth}$ [rad]")
+plt.ylabel("$\\theta^{hfs}_{\\rm rec}$ [rad]")
+plt.title(f"$E_{{\\tau}}$=20-100 GeV, ${eta_min}<\\eta_{{hfs}}<{eta_max}$")
+plt.plot([0.04,0.1], [0.04, 0.1], color='tab:orange')
+plt.ylim(0, 0.15)
+plt.savefig(outdir+"/theta_scatter.pdf")
+
+plt.figure()
+plt.scatter(Etruth, Erecon, 1)
+plt.xlabel("$E^{hfs}_{\\rm truth}$ [GeV]")
+plt.ylabel("$E^{hfs}_{\\rm rec}$ [GeV]")
+plt.title(f"$E_{{\\tau}}$=20-100 GeV, ${eta_min}<\\eta_{{hfs}}<{eta_max}$")
+plt.plot((0,100), (0, 100), color='tab:orange')
+plt.savefig(outdir+"/energy_scatter.pdf")
+
+def gauss(x, A,mu, sigma):
+ return A * np.exp(-(x-mu)**2/(2*sigma**2))
+from scipy.optimize import curve_fit
+
+res=[]
+dres=[]
+emid=[]
+ew=[]
+partitions=(20,30, 40, 60,80,100)
+for emin, emax in zip(partitions[:-1], partitions[1:]):
+
+ y,x = np.histogram((theta_recon-theta_truth)[(Etruth>emin)&(Etruth<emax)], bins=100, range=(-0.03,0.03))
+ bc=(x[1:]+x[:-1])/2
+ slc=abs(bc)<0.5
+ # try:
+ p0=(100, 0, 0.15)
+ fnc=gauss
+ sigma=np.sqrt(y[slc])+(y[slc]==0)
+
+ coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000)
+ res.append(abs(coeff[2]))
+ dres.append(np.sqrt(var_matrix[2][2]))
+ emid.append((emin+emax)/2)
+ ew.append((emax-emin)/2)
+plt.errorbar(emid, 1000*np.array(res),1000*np.array(dres), ew, ls='', label=f'{eta_min}<$\\eta_{{hfs}}$<{eta_max}')
+plt.xlabel('$E_{hfs}$ [GeV]')
+plt.ylabel('$\\theta$ resolution [mrad]')
+plt.ylim(0)
+
+fnc=lambda E,B:B/E
+p0=[1,]
+coeff, var_matrix = curve_fit(fnc, emid, res, p0=p0, sigma=list(dres), maxfev=10000)
+xx=np.linspace(10, 100, 100)
+plt.plot(xx, 1000*fnc(xx, *coeff), label=f"fit: ${coeff[0]:.2f}/E$ mrad")
+plt.legend()
+plt.savefig(outdir+"/theta_res.pdf")
diff --git a/benchmarks/insert_tau/config.yml b/benchmarks/insert_tau/config.yml
new file mode 100644
index 00000000..1c43d0e2
--- /dev/null
+++ b/benchmarks/insert_tau/config.yml
@@ -0,0 +1,36 @@
+sim:insert_tau:
+ stage: simulate
+ extends: .det_benchmark
+ parallel:
+ matrix:
+ - P: 20
+ - P: 30
+ - P: 40
+ - P: 50
+ - P: 60
+ - P: 80
+ - P: 100
+ timeout: 1 hours
+ script:
+ - snakemake $SNAKEMAKE_FLAGS --cores 5 sim_output/insert_tau/epic_craterlake_sim_tau-_${P}GeV_{0,1,2,3,4}.edm4hep.root
+ retry:
+ max: 2
+ when:
+ - runner_system_failure
+
+bench:insert_tau:
+ stage: benchmarks
+ extends: .det_benchmark
+ needs: ["sim:insert_tau"]
+ script:
+ - snakemake $SNAKEMAKE_FLAGS --cores 1 results/epic_craterlake/insert_tau
+
+collect_results:insert_tau:
+ extends: .det_benchmark
+ stage: collect
+ needs: ["bench:insert_tau"]
+ script:
+ - ls -al
+ - mv results{,_save}/ # move results directory out of the way to preserve it
+ - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output results/epic_craterlake/insert_tau
+ - mv results{_save,}/
--
GitLab