diff --git a/.gitignore b/.gitignore
index 8fdeae7104dc3f3bc8bdc0586e2f8cab338df867..370a1989db9c164205c6b3144b09f234db1bdb93 100644
--- a/.gitignore
+++ b/.gitignore
@@ -18,6 +18,8 @@ BUILD*/*
RELEASE*/*
TEST*/*
.local/*
+build/*
+install/*
# cmake
CMakeCache.txt
diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..df59b35cc09bf3967a092c4415e6d1cbb0f13242
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,22 @@
+cmake_minimum_required(VERSION 3.0.0 FATAL_ERROR)
+
+project(reconstruction_benchmarks LANGUAGES CXX)
+
+cmake_policy(SET CMP0074 NEW) # use `<PackageName>_ROOT` variables
+include(GNUInstallDirs)
+
+# ROOT
+find_package(ROOT 6 REQUIRED COMPONENTS Core RIO Hist)
+include(${ROOT_USE_FILE})
+
+# data model
+find_package(podio REQUIRED)
+find_package(EDM4HEP REQUIRED)
+find_package(EDM4EIC REQUIRED)
+
+# logging
+find_package(spdlog REQUIRED)
+add_compile_definitions(SPDLOG_FMT_EXTERNAL)
+
+# benchmarks
+add_subdirectory(benchmarks/rich)
diff --git a/benchmarks/rich/CMakeLists.txt b/benchmarks/rich/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..9afbcef280e6c37a6097880ac6fde901f0f894cf
--- /dev/null
+++ b/benchmarks/rich/CMakeLists.txt
@@ -0,0 +1,34 @@
+get_filename_component(det_name ${CMAKE_CURRENT_LIST_DIR} NAME)
+project(${CMAKE_PROJECT_NAME}_${det_name})
+
+# sources
+set(algo_exe_source ${CMAKE_CURRENT_SOURCE_DIR}/src/benchmark.cc)
+file(GLOB algo_headers CONFIGURE_DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/include/*.h)
+file(GLOB algo_sources CONFIGURE_DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/src/*.cc)
+list(REMOVE_ITEM algo_sources ${algo_exe_source})
+
+# library
+set(algo_lib ${PROJECT_NAME})
+add_library(${algo_lib} SHARED ${algo_sources})
+target_include_directories(${algo_lib} PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/include)
+target_compile_options(${algo_lib} PRIVATE -Wall -Wno-misleading-indentation)
+
+# executable
+set(algo_exe benchmark_${det_name}_reconstruction)
+add_executable(${algo_exe} ${algo_exe_source})
+
+# linking
+target_link_libraries(${algo_exe}
+ ${algo_lib}
+ ROOT::Core
+ ROOT::Hist
+ podio::podio
+ podio::podioRootIO
+ EDM4EIC::edm4eic
+ EDM4HEP::edm4hep
+ spdlog::spdlog
+ )
+
+# installation
+install(FILES ${algo_headers} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/${CMAKE_PROJECT_NAME}/${det_name})
+install(TARGETS ${algo_exe} ${algo_lib})
diff --git a/benchmarks/rich/README.md b/benchmarks/rich/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..3f2a53f74034149bcadb7a7f1778924dbd65414e
--- /dev/null
+++ b/benchmarks/rich/README.md
@@ -0,0 +1,16 @@
+# RICH Benchmarks
+
+## Setup
+- Build `reconstruction_benchmarks` with `cmake`
+ - Builds code in `src/` and `include/` to a library
+ - Benchmarks are designed as independent algorithms, similar to reconstruction algorithms
+ - Algorithm inputs are PODIO collections produced from reconstruction, and the output
+ is typically a set of `ROOT` plots
+ - Builds executable `benchmark_rich_reconstruction`, which runs any or all of the benchmark
+ algorithms and streams any plots to an output `ROOT` file
+
+## Running
+- Run the benchmark executable using the wrapper script `run_benchmark.rb`
+ - Run with no arguments for usage guide
+ - This script can run simulation, reconstruction and/or benchmarks
+ - Run `benchmark_rich_reconstruction` to just run the benchmark algorithms
diff --git a/benchmarks/rich/draw_benchmark.py b/benchmarks/rich/draw_benchmark.py
new file mode 100755
index 0000000000000000000000000000000000000000..df220f6e1dcc56fc48958e6ae1efc85019dc655a
--- /dev/null
+++ b/benchmarks/rich/draw_benchmark.py
@@ -0,0 +1,241 @@
+#!/usr/bin/env python
+# Copyright 2023, Christopher Dilks
+# Subject to the terms in the LICENSE file found in the top-level directory.
+
+import ROOT as r
+import sys, getopt, pathlib, math
+
+# suppress graphics
+r.gROOT.SetBatch(True)
+
+# GLOBAL SETTINGS
+################################################################
+RESID_MAX = 10 # cherenkov angle |residual| maximum
+RADIATORS = {
+ 'Aerogel': { 'p_max': 30, 'p_rebin': 12, 'rindex_ref': 1.0190},
+ 'Gas': { 'p_max': 70, 'p_rebin': 30, 'rindex_ref': 1.00076},
+ 'Merged': { 'p_max': 70, 'p_rebin': 30, },
+ }
+PARTICLE_MASS = {
+ 'e': 0.00051,
+ 'pi': 0.13957,
+ 'K': 0.49368,
+ 'p': 0.93827,
+ }
+
+# ARGUMENTS
+################################################################
+
+ana_file_name = 'out/ana.edm4hep.root'
+output_dir = 'out/ana.plots'
+
+helpStr = f'''
+{sys.argv[0]} [OPTIONS]
+
+ -i <input file>: specify an input file, e.g., hepmc
+ default: {ana_file_name}
+
+ -o <output dir>: specify an output directory
+ default: {output_dir}
+
+ -h: show this usage guide
+
+ '''
+
+try:
+ opts, args = getopt.getopt(sys.argv[1:], 'i:o:h')
+except getopt.GetoptError:
+ print('\n\nERROR: invalid argument\n', helpStr, file=sys.stderr)
+ sys.exit(2)
+for opt, arg in opts:
+ if(opt == '-i'): ana_file_name = arg.lstrip()
+ if(opt == '-o'): output_dir = arg.lstrip()
+ if(opt == '-h'):
+ print(helpStr)
+ sys.exit(2)
+print(f'''
+ana_file_name = {ana_file_name}
+output_dir = {output_dir}
+''')
+
+
+# PLOTTING
+################################################################
+
+# make canvases
+ana_file = r.TFile.Open(ana_file_name, "READ")
+def make_canv(name, nx, ny):
+ canv = r.TCanvas(name, name, nx*600, ny*400)
+ canv.Divide(nx,ny)
+ return canv
+canv_dict = {
+ "photon_spectra": make_canv("photon_spectra", 2, 2),
+ "digitization": make_canv("digitization", 2, 2),
+ "pidAerogel": make_canv("pidAerogel", 5, 3),
+ "pidGas": make_canv("pidGas", 5, 3),
+ "pidMerged": make_canv("pidMerged", 2, 3),
+}
+
+# helper functions
+### minimum momentum for Cherenkov radiation
+def calculate_mom_min(mass,n):
+ return mass / math.sqrt(n**2-1)
+### Cherenkov angle at saturation
+def calculate_theta_max(n):
+ return 1000 * math.acos(1/n)
+### execute `op` for every pad of `canv`
+def loop_pads(canv, op):
+ for pad in canv.GetListOfPrimitives():
+ if pad.InheritsFrom(r.TVirtualPad.Class()):
+ op(pad)
+### draw a profile on 2D histogram `hist`
+def draw_profile(hist, rad, style="BOX"):
+ prof = hist.ProfileX("_pfx", 1, -1, "i")
+ prof.SetMarkerStyle(r.kFullCircle)
+ prof.SetMarkerColor(r.kRed)
+ prof.SetMarkerSize(2)
+ hist_rebinned = hist.Clone(f'{hist.GetName()}_rebin')
+ if 'vs_p' in hist.GetName():
+ hist_rebinned.RebinX(RADIATORS[rad]['p_rebin'])
+ else:
+ hist_rebinned.RebinX(4)
+ hist_rebinned.SetLineColor(r.kBlue)
+ hist_rebinned.SetFillColor(r.kBlue)
+ for p in [hist_rebinned, prof]:
+ if 'vs_p' in hist.GetName():
+ p.GetXaxis().SetRangeUser(0, RADIATORS[rad]['p_max'])
+ if 'rindex_ref' in RADIATORS[rad] and 'theta_vs_p' in hist.GetName():
+ p.GetYaxis().SetRangeUser(0, 1.5*calculate_theta_max(RADIATORS[rad]['rindex_ref']))
+ if 'thetaResid_vs' in hist.GetName():
+ p.GetYaxis().SetRangeUser(-RESID_MAX, RESID_MAX)
+ hist_rebinned.Draw(style)
+ prof.Draw("SAME")
+
+# Cherenkov angle vs. p functions
+for rad, radHash in RADIATORS.items():
+ if rad == "Merged":
+ continue
+ radHash['cherenkov_curves'] = []
+ n = radHash['rindex_ref']
+ for part, mass in PARTICLE_MASS.items():
+ ftn = r.TF1(
+ f'ftn_theta_{part}_{rad}',
+ f'1000*TMath::ACos(TMath::Sqrt(x^2+{mass}^2)/({n}*x))',
+ calculate_mom_min(mass, n),
+ radHash['p_max']
+ )
+ ftn.SetLineColor(r.kBlack)
+ ftn.SetLineWidth(1)
+ radHash['cherenkov_curves'].append(ftn)
+
+# draw photon spectra
+canv = canv_dict["photon_spectra"]
+def canv_op(pad):
+ pad.SetGrid(1,1)
+ pad.SetLogy()
+loop_pads(canv, canv_op)
+canv.cd(1)
+ana_file.Get("phot/phot_spectrum_sim").Draw()
+canv.cd(2)
+ana_file.Get("phot/nphot_dist").Draw()
+canv.cd(3)
+ana_file.Get("digi/phot_spectrum_rec").Draw()
+canv.cd(4)
+ana_file.Get("digi/nhits_dist").Draw()
+
+# draw digitization
+canv = canv_dict["digitization"]
+def canv_op(pad):
+ pad.SetGrid(1,1)
+ if pad.GetNumber() < 3:
+ pad.SetLogy()
+ else:
+ pad.SetLogz()
+loop_pads(canv, canv_op)
+canv.cd(1)
+ana_file.Get("digi/adc_dist").Draw()
+canv.cd(2)
+ana_file.Get("digi/tdc_dist").Draw()
+canv.cd(3)
+ana_file.Get("digi/tdc_vs_adc").Draw("COLZ")
+
+# draw CherenkovPID for each radiator
+for rad in ["Aerogel", "Gas"]:
+ pid_name = f'pid{rad}'
+ canv = canv_dict[pid_name]
+ def canv_op(pad):
+ pad.SetGrid(1,1)
+ loop_pads(canv, canv_op)
+ canv.cd(1)
+ ana_file.Get(f'{pid_name}/npe_dist_{rad}').Draw()
+ canv.cd(2)
+ ana_file.Get(f'{pid_name}/theta_dist_{rad}').Draw()
+ canv.cd(3)
+ hist = ana_file.Get(f'{pid_name}/thetaResid_dist_{rad}')
+ hist.DrawCopy()
+ canv.cd(4)
+ hist.SetTitle(hist.GetTitle() + " - ZOOM")
+ hist.GetXaxis().SetRangeUser(-RESID_MAX, RESID_MAX)
+ hist.Draw()
+ resid_mode = hist.GetBinCenter(hist.GetMaximumBin())
+ resid_dev = hist.GetStdDev()
+ nsigma = 3
+ hist.Fit("gaus", "", "", resid_mode - nsigma*resid_dev, resid_mode + nsigma*resid_dev)
+ if hist.GetFunction("gaus"):
+ hist.GetFunction("gaus").SetLineWidth(5)
+ canv.cd(5)
+ ana_file.Get(f'{pid_name}/highestWeight_dist_{rad}').Draw()
+ canv.cd(6)
+ hist = ana_file.Get(f'{pid_name}/npe_vs_p_{rad}')
+ draw_profile(hist, rad)
+ canv.cd(7)
+ hist = ana_file.Get(f'{pid_name}/theta_vs_p_{rad}')
+ draw_profile(hist, rad)
+ if 'cherenkov_curves' in RADIATORS[rad]:
+ for ftn in RADIATORS[rad]['cherenkov_curves']:
+ ftn.Draw("SAME")
+ canv.cd(8)
+ hist = ana_file.Get(f'{pid_name}/thetaResid_vs_p_{rad}')
+ hist.GetYaxis().SetRangeUser(-RESID_MAX, RESID_MAX)
+ draw_profile(hist, rad)
+ canv.cd(9)
+ ana_file.Get(f'{pid_name}/photonTheta_vs_photonPhi_{rad}').Draw("COLZ")
+ canv.cd(10)
+ ana_file.Get(f'{pid_name}/mcRindex_{rad}').Draw()
+ canv.cd(11)
+ hist = ana_file.Get(f'{pid_name}/npe_vs_eta_{rad}')
+ draw_profile(hist, rad)
+ canv.cd(12)
+ hist = ana_file.Get(f'{pid_name}/theta_vs_eta_{rad}')
+ draw_profile(hist, rad)
+ canv.cd(13)
+ hist = ana_file.Get(f'{pid_name}/thetaResid_vs_eta_{rad}')
+ hist.GetYaxis().SetRangeUser(-RESID_MAX, RESID_MAX)
+ draw_profile(hist, rad)
+
+# draw CherenkovPID for merged radiators
+pid_name = f'pidMerged'
+canv = canv_dict[pid_name]
+def canv_op(pad):
+ pad.SetGrid(1,1)
+loop_pads(canv, canv_op)
+canv.cd(1)
+ana_file.Get(f'{pid_name}/npe_dist_Merged').Draw()
+canv.cd(2)
+ana_file.Get(f'{pid_name}/highestWeight_dist_Merged').Draw()
+canv.cd(3)
+hist = ana_file.Get(f'{pid_name}/npe_vs_p_Merged')
+draw_profile(hist, 'Merged')
+canv.cd(4)
+ana_file.Get(f'{pid_name}/photonTheta_vs_photonPhi_Merged').Draw("COLZ")
+canv.cd(5)
+hist = ana_file.Get(f'{pid_name}/npe_vs_eta_Merged')
+draw_profile(hist, 'Merged')
+
+# FINISH
+################################################################
+
+pathlib.Path(output_dir).mkdir(parents=True, exist_ok=True)
+for name, canvas in canv_dict.items():
+ canvas.SaveAs(f'{output_dir}/{name}.png')
+ana_file.Close()
diff --git a/benchmarks/rich/forward_hadrons.sh b/benchmarks/rich/forward_hadrons.sh
deleted file mode 100644
index aebf4ea89974a50507854fade7216405e18fce26..0000000000000000000000000000000000000000
--- a/benchmarks/rich/forward_hadrons.sh
+++ /dev/null
@@ -1,89 +0,0 @@
-#!/bin/bash
-
-./util/print_env.sh
-
-## To run the reconstruction, we need the following global variables:
-## - DETECTOR: the detector package we want to use for this benchmark
-## - DETECTOR_VERSION: the detector package we want to use for this benchmark
-## - DETECTOR_PATH: full path to the detector definitions
-##
-## You can ready options/env.sh for more in-depth explanations of the variables
-## and how they can be controlled.
-source options/env.sh
-export DETECTOR_PATH=${DETECTOR_PATH}
-
-if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then
- export JUGGLER_N_EVENTS=100
-fi
-
-export JUGGLER_FILE_NAME_TAG="rich_forward_hadrons"
-export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
-
-export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.edm4hep.root"
-export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
-
-echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
-echo "DETECTOR = ${DETECTOR}"
-
-
-# generate the input events
-# note datasets is now only used to develop datasets.
-#git clone https://eicweb.phy.anl.gov/EIC/datasets.git datasets
-python benchmarks/rich/scripts/rich_data_gen.py \
- ${JUGGLER_FILE_NAME_TAG}.hepmc \
- -n ${JUGGLER_N_EVENTS} \
- --pmin 5.0 \
- --pmax 100.0 \
- --angmin 3.0 \
- --angmax 8.0
-if [[ "$?" -ne "0" ]] ; then
- echo "ERROR running script"
- exit 1
-fi
-
-# This script appears to be missing
-## run geant4 simulations
-#python benchmakrs/options/ForwardRICH/simu.py \
-# --compact=${DETECTOR_PATH}/${DETECTOR_CONFIG}.xml \
-# -i ${JUGGLER_FILE_NAME_TAG}.hepmc \
-# -o ${JUGGLER_SIM_FILE}
-# -n ${JUGGLER_N_EVENTS}
-#if [[ "$?" -ne "0" ]] ; then
-# echo "ERROR running script"
-# exit 1
-#fi
-
-#### skipping
-exit 0
-
-
-ddsim --runType batch \
- --part.minimalKineticEnergy 1000*GeV \
- --filter.tracker edep0 \
- -v WARNING \
- --numberOfEvents ${JUGGLER_N_EVENTS} \
- --compactFile ${DETECTOR_PATH}/${DETECTOR_CONFIG}.xml \
- --inputFiles ${JUGGLER_FILE_NAME_TAG}.hepmc \
- --outputFile ${JUGGLER_SIM_FILE}
-if [[ "$?" -ne "0" ]] ; then
- echo "ERROR running npdet"
- exit 1
-fi
-
-# @TODO changeable simulation file name and detector xml file name
-gaudirun.py benchmarks/rich/options/rich_reco.py
-status=$?
-if [[ "$status" -ne "0" && "$status" -ne "4" ]] ; then
- echo "ERROR running juggler, got $status"
- exit 1
-fi
-
-# @TODO add analysis scripts
-#root_filesize=$(stat --format=%s "${DETECTOR}/${JUGGLER_REC_FILE}")
-#if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then
-# # file must be less than 10 MB to upload
-# if [[ "${root_filesize}" -lt "10000000" ]] ; then
-# cp ${DETECTOR}/${JUGGLER_REC_FILE} results/.
-# fi
-#fi
-
diff --git a/benchmarks/rich/include/ChargedParticle.h b/benchmarks/rich/include/ChargedParticle.h
new file mode 100644
index 0000000000000000000000000000000000000000..254819a2a8145e9332da0194013e5b3ed80b869c
--- /dev/null
+++ b/benchmarks/rich/include/ChargedParticle.h
@@ -0,0 +1,42 @@
+// Copyright 2023, Christopher Dilks
+// Subject to the terms in the LICENSE file found in the top-level directory.
+
+#pragma once
+
+#include <spdlog/spdlog.h>
+
+#include <edm4hep/MCParticleCollection.h>
+#include <edm4hep/utils/vector_utils.h>
+
+namespace benchmarks {
+
+ class ChargedParticle {
+
+ public:
+
+ ChargedParticle(std::shared_ptr<spdlog::logger> log) : m_log(log) {};
+ ~ChargedParticle() {}
+
+ // accessors
+ double GetMomentum() { return m_momentum; };
+ double GetEta() { return m_eta; };
+ int GetPDG() { return m_pdg; };
+
+ // set info from a single MCParticle
+ void SetSingleParticle(const edm4hep::MCParticle& mc_part);
+ // set info from a (thrown) MCParticle of a collection
+ void SetSingleParticle(const edm4hep::MCParticleCollection& mc_parts);
+
+ private:
+
+ // members
+ double m_momentum;
+ double m_eta;
+ int m_pdg = 0;
+
+ // logger
+ std::shared_ptr<spdlog::logger> m_log;
+
+ };
+
+}
diff --git a/benchmarks/rich/include/CherenkovPIDAnalysis.h b/benchmarks/rich/include/CherenkovPIDAnalysis.h
new file mode 100644
index 0000000000000000000000000000000000000000..76a310e73d028856d5068c5984059deea75863dc
--- /dev/null
+++ b/benchmarks/rich/include/CherenkovPIDAnalysis.h
@@ -0,0 +1,72 @@
+// Copyright 2023, Christopher Dilks
+// Subject to the terms in the LICENSE file found in the top-level directory.
+
+#pragma once
+
+#include <spdlog/spdlog.h>
+
+#include <TH1D.h>
+#include <TH2D.h>
+#include <TMath.h>
+#include <TRegexp.h>
+
+#include <edm4eic/CherenkovParticleIDCollection.h>
+#include <edm4hep/MCParticleCollection.h>
+#include <edm4hep/utils/vector_utils.h>
+#include <edm4hep/utils/kinematics.h>
+
+#include "Tools.h"
+#include "ChargedParticle.h"
+
+namespace benchmarks {
+
+ class CherenkovPIDAnalysis {
+
+ public:
+ CherenkovPIDAnalysis() = default;
+ ~CherenkovPIDAnalysis() {}
+
+ // algorithm methods
+ void AlgorithmInit(
+ std::string rad_name,
+ std::shared_ptr<spdlog::logger>& logger
+ );
+ void AlgorithmProcess(
+ const edm4hep::MCParticleCollection& mc_parts,
+ const edm4eic::CherenkovParticleIDCollection& cherenkov_pids
+ );
+ void AlgorithmFinish();
+
+ private:
+
+ // radiator name
+ std::string m_rad_name;
+ TString m_rad_name_trun;
+ TString m_rad_title;
+
+ // histograms
+ // - distributions
+ TH1D *m_npe_dist;
+ TH1D *m_theta_dist;
+ TH1D *m_thetaResid_dist;
+ TH1D *m_mcWavelength_dist;
+ TH1D *m_mcRindex_dist;
+ TH1D *m_highestWeight_dist;
+ TH2D *m_photonTheta_vs_photonPhi;
+ // - momentum scans
+ TH2D *m_npe_vs_p;
+ TH2D *m_theta_vs_p;
+ TH2D *m_thetaResid_vs_p;
+ TH2D *m_highestWeight_vs_p;
+ // - pseudorapidity scans
+ TH2D *m_npe_vs_eta;
+ TH2D *m_theta_vs_eta;
+ TH2D *m_thetaResid_vs_eta;
+ TH2D *m_highestWeight_vs_eta;
+
+ // logger
+ std::shared_ptr<spdlog::logger> m_log;
+
+ };
+
+}
diff --git a/benchmarks/rich/include/RawHitAnalysis.h b/benchmarks/rich/include/RawHitAnalysis.h
new file mode 100644
index 0000000000000000000000000000000000000000..56b0b839b0e3cac38b1368ef8a816fbb81fe5d9e
--- /dev/null
+++ b/benchmarks/rich/include/RawHitAnalysis.h
@@ -0,0 +1,52 @@
+// Copyright 2023, Christopher Dilks
+// Subject to the terms in the LICENSE file found in the top-level directory.
+
+#pragma once
+
+#include <spdlog/spdlog.h>
+
+#include <TH1D.h>
+#include <TH2D.h>
+#include <TMath.h>
+
+#include <edm4eic/RawTrackerHitCollection.h>
+#include <edm4eic/MCRecoTrackerHitAssociationCollection.h>
+
+#include "Tools.h"
+
+namespace benchmarks {
+
+ class RawHitAnalysis {
+
+ public:
+
+ RawHitAnalysis() = default;
+ ~RawHitAnalysis() {}
+
+ // algorithm methods
+ void AlgorithmInit(std::shared_ptr<spdlog::logger>& logger);
+ void AlgorithmProcess(
+ const edm4eic::RawTrackerHitCollection& raw_hits,
+ const edm4eic::MCRecoTrackerHitAssociationCollection& assocs
+ );
+ void AlgorithmFinish();
+
+ private:
+
+ // binning
+ const int adc_max = std::pow(2,10);
+ const int tdc_max = std::pow(2,10);
+
+ // histograms
+ TH1D *m_adc_dist;
+ TH1D *m_tdc_dist;
+ TH2D *m_tdc_vs_adc;
+ TH1D *m_phot_spectrum;
+ TH1D *m_nhits_dist;
+
+ // logging
+ std::shared_ptr<spdlog::logger> m_log;
+
+ };
+
+}
diff --git a/benchmarks/rich/include/ReconstructedParticleAnalysis.h b/benchmarks/rich/include/ReconstructedParticleAnalysis.h
new file mode 100644
index 0000000000000000000000000000000000000000..5ca8f6945d0db1b22234dfae9c1569dc56265f83
--- /dev/null
+++ b/benchmarks/rich/include/ReconstructedParticleAnalysis.h
@@ -0,0 +1,42 @@
+// Copyright 2023, Christopher Dilks
+// Subject to the terms in the LICENSE file found in the top-level directory.
+
+#pragma once
+
+#include <spdlog/spdlog.h>
+
+#include <TH2D.h>
+#include <TGraphErrors.h>
+
+#include <edm4eic/MCRecoParticleAssociationCollection.h>
+#include <edm4hep/utils/kinematics.h>
+
+#include "Tools.h"
+
+namespace benchmarks {
+
+ class ReconstructedParticleAnalysis {
+
+ public:
+ ReconstructedParticleAnalysis() = default;
+ ~ReconstructedParticleAnalysis() {}
+
+ // algorithm methods
+ void AlgorithmInit(std::shared_ptr<spdlog::logger>& logger);
+ void AlgorithmProcess(
+ const edm4eic::MCRecoParticleAssociationCollection& in_assocs
+ );
+ void AlgorithmFinish();
+
+ private:
+
+ // plots
+ TGraphErrors *m_correct_vs_p;
+ TH2D *m_correct_vs_p_transient; // (not written)
+
+ // additional objects
+ std::shared_ptr<spdlog::logger> m_log;
+
+ };
+
+}
diff --git a/benchmarks/rich/include/SimHitAnalysis.h b/benchmarks/rich/include/SimHitAnalysis.h
new file mode 100644
index 0000000000000000000000000000000000000000..836a6ad2999c4d25a21e967d5cccaee1acab4ff2
--- /dev/null
+++ b/benchmarks/rich/include/SimHitAnalysis.h
@@ -0,0 +1,46 @@
+// Copyright 2023, Christopher Dilks
+// Subject to the terms in the LICENSE file found in the top-level directory.
+
+#pragma once
+
+#include <spdlog/spdlog.h>
+
+#include <TH1D.h>
+#include <TH2D.h>
+
+#include <edm4hep/SimTrackerHitCollection.h>
+#include <edm4hep/MCParticleCollection.h>
+
+#include "Tools.h"
+#include "ChargedParticle.h"
+
+namespace benchmarks {
+
+ class SimHitAnalysis {
+
+ public:
+ SimHitAnalysis() = default;
+ ~SimHitAnalysis() {}
+
+ // algorithm methods
+ void AlgorithmInit(std::shared_ptr<spdlog::logger>& logger);
+ void AlgorithmProcess(
+ const edm4hep::MCParticleCollection& mc_parts,
+ const edm4hep::SimTrackerHitCollection& sim_hits
+ );
+ void AlgorithmFinish();
+
+ private:
+
+ // histograms
+ TH1D *m_nphot_dist;
+ TH2D *m_nphot_vs_p;
+ TH1D *m_nphot_vs_p__transient; // transient (not written)
+ TH1D *m_phot_spectrum;
+
+ // logger
+ std::shared_ptr<spdlog::logger> m_log;
+
+ };
+
+}
diff --git a/benchmarks/rich/include/Tools.h b/benchmarks/rich/include/Tools.h
new file mode 100644
index 0000000000000000000000000000000000000000..c34a142a9708ce45730ed6113f3ae77733de05ca
--- /dev/null
+++ b/benchmarks/rich/include/Tools.h
@@ -0,0 +1,103 @@
+// Copyright 2023, Christopher Dilks
+// Subject to the terms in the LICENSE file found in the top-level directory.
+
+#pragma once
+
+// general
+#include <map>
+#include <spdlog/spdlog.h>
+
+// ROOT
+#include <TVector3.h>
+
+// IRT
+#include <IRT/ParametricSurface.h>
+
+// DD4hep
+#include <Evaluator/DD4hepUnits.h>
+
+// data model
+#include <edm4hep/SimTrackerHitCollection.h>
+#include <edm4hep/utils/kinematics.h>
+
+namespace benchmarks {
+
+ // Tools class, filled with miscellanous helper functions
+ class Tools {
+ public:
+
+ // -------------------------------------------------------------------------------------
+ // physics constants
+
+ // h*c constant, for wavelength <=> energy conversion [GeV*nm]
+ // obtained from `dd4hep::h_Planck * dd4hep::c_light / (dd4hep::GeV * dd4hep::nm)`
+ static constexpr double HC = 1.239841875e-06;
+
+ // opticalphoton PDG
+ static constexpr int PHOTON_PDG = -22;
+
+ // -------------------------------------------------------------------------------------
+ // PDG mass lookup
+
+ // local PDG mass database
+ // FIXME: cannot use `TDatabasePDG` since it is not thread safe; until we
+ // have a proper PDG database service, we hard-code the masses we need;
+ // use Tools::GetPDGMass for access
+ static std::unordered_map<int,double> GetPDGMasses() {
+ return std::unordered_map<int,double>{
+ { 11, 0.000510999 },
+ { 211, 0.13957 },
+ { 321, 0.493677 },
+ { 2212, 0.938272 }
+ };
+ }
+
+ static double GetPDGMass(int pdg) {
+ double mass;
+ try { mass = GetPDGMasses().at(std::abs(pdg)); }
+ catch(const std::out_of_range& e) {
+ throw std::runtime_error(fmt::format("RUNTIME ERROR: unknown PDG={} in algorithms/pid/Tools::GetPDGMass",pdg));
+ }
+ return mass;
+ }
+
+ static int GetNumPDGs() { return GetPDGMasses().size(); };
+
+
+ // -------------------------------------------------------------------------------------
+ // get photon wavelength
+ static double GetPhotonWavelength(const edm4hep::SimTrackerHit& hit, std::shared_ptr<spdlog::logger> log) {
+ auto phot = hit.getMCParticle();
+ if(!phot.isAvailable()) {
+ log->error("no MCParticle in hit");
+ return -1.0;
+ }
+ if(phot.getPDG() != PHOTON_PDG) {
+ log->warn("hit MCParticle is not an opticalphoton; PDG is {}; ignoring", phot.getPDG());
+ return -1.0;
+ }
+ auto momentum = edm4hep::utils::p(phot);
+ auto wavelength = momentum > 0 ? HC / momentum : 0.0;
+ return wavelength;
+ }
+
+ // -------------------------------------------------------------------------------------
+ // binning
+ static constexpr int n_bins = 100;
+ static constexpr int momentum_bins = 500;
+ static constexpr double momentum_max = 70;
+ static constexpr int eta_bins = 50;
+ static constexpr double eta_min = 0;
+ static constexpr double eta_max = 5;
+ static constexpr int npe_max = 50;
+ static constexpr int nphot_max = 3000;
+ static constexpr int theta_bins = 1500;
+ static constexpr double theta_max = 300;
+ static constexpr double thetaResid_max = 100;
+ static constexpr int phi_bins = 100;
+ static constexpr int wl_bins = 120;
+ static constexpr double wl_min = 0;
+ static constexpr double wl_max = 1200;
+
+ }; // class Tools
+} // namespace benchmarks
diff --git a/benchmarks/rich/options/rich_reco.py b/benchmarks/rich/options/rich_reco.py
deleted file mode 100644
index f27f62a534302aa0f266e1d10f046e867ef7d283..0000000000000000000000000000000000000000
--- a/benchmarks/rich/options/rich_reco.py
+++ /dev/null
@@ -1,43 +0,0 @@
-from Gaudi.Configuration import *
-from GaudiKernel import SystemOfUnits as units
-
-from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
-
-detector_path = str(os.environ.get("DETECTOR_PATH", "."))
-detector_name = str(os.environ.get("DETECTOR_CONFIG", "epic"))
-detector_config = str(os.environ.get("DETECTOR_CONFIG", detector_name))
-detector_version = str(os.environ.get("DETECTOR_VERSION", "main"))
-
-# todo add checks
-input_sim_file = str(os.environ["JUGGLER_SIM_FILE"])
-output_rec_file = str(os.environ["JUGGLER_REC_FILE"])
-n_events = str(os.environ["JUGGLER_N_EVENTS"])
-
-geo_service = GeoSvc("GeoSvc", detectors=["{}/{}.xml".format(detector_path, detector_name)])
-podioevent = EICDataSvc("EventDataSvc", inputs=[input_sim_file], OutputLevel=DEBUG)
-
-from Configurables import PodioInput
-from Configurables import Jug__Digi__PhotoMultiplierDigi as PhotoMultiplierDigi
-from Configurables import Jug__Reco__PhotoMultiplierReco as PhotoMultiplierReco
-from Configurables import Jug__Reco__PhotoRingClusters as PhotoRingClusters
-
-qe_data = [(1.0, 0.25), (7.5, 0.25),]
-podioinput = PodioInput("PodioReader", collections=["MCParticles", "ForwardRICHHits"], OutputLevel=DEBUG)
-
-pmtdigi = PhotoMultiplierDigi(inputHitCollection="ForwardRICHHits", outputHitCollection="DigiForwardRICHHits",
- quantumEfficiency=[(a*units.eV, b) for a, b in qe_data])
-pmtreco = PhotoMultiplierReco(inputHitCollection="DigiForwardRICHHits", outputHitCollection="RecoForwardRICHHits")
-richcluster = PhotoRingClusters(inputHitCollection="RecoForwardRICHHits", #inputTrackCollection="MCParticles",
- outputClusterCollection="RICHClusters")
-
-out = PodioOutput("out", filename=output_rec_file)
-out.outputCommands = ["keep *"]
-
-ApplicationMgr(
- TopAlg = [podioinput, pmtdigi, pmtreco, richcluster, out],
- EvtSel = 'NONE',
- EvtMax = n_events,
- ExtSvc = [podioevent],
- OutputLevel=DEBUG
- )
-
diff --git a/benchmarks/rich/run_benchmark.rb b/benchmarks/rich/run_benchmark.rb
new file mode 100755
index 0000000000000000000000000000000000000000..0a1a748e3643529aafad8756a412f52090c7cee3
--- /dev/null
+++ b/benchmarks/rich/run_benchmark.rb
@@ -0,0 +1,240 @@
+#!/usr/bin/env ruby
+# Copyright 2023, Christopher Dilks
+# Subject to the terms in the LICENSE file found in the top-level directory.
+
+require 'optparse'
+require 'ostruct'
+
+###################################
+# constants:
+IdealTheta = 23.5 # [degrees]
+IdealEnergy = 20.0 # [GeV]
+IdealParticle = 'pi+'
+###################################
+
+# setup
+# ---------------------------------------------------
+
+# default opt
+opt = OpenStruct.new
+opt.evgen = ''
+opt.sim_file = 'out/sim.edm4hep.root'
+opt.rec_file = 'out/rec.edm4hep.root'
+opt.ana_file = 'out/ana.edm4hep.root'
+opt.run_sim = true
+opt.run_rec = true
+opt.run_ana = true
+opt.run_rec_down = false
+opt.run_ana_down = false
+opt.num_events = 10
+opt.benchmark_exe = 'benchmark_rich_reconstruction'
+opt.algos = Array.new
+opt.plots = Array.new
+opt.verbosity = 0
+opt.dry_run = false
+
+# available event generation modes
+avail_evgens = [
+ 'idealAngle',
+ 'minTheta',
+ 'maxTheta',
+]
+
+# parse options
+required_set = false
+OptionParser.new do |o|
+ o.banner = "USAGE: #{$0} [OPTIONS]..."
+ o.separator ''
+ o.separator 'required options, one of either:'.upcase
+ o.separator ''
+ o.on("-r", "--rec-only", "Run only the reconstruction, then the analysis benchmark") do |a|
+ opt.run_rec_down = true
+ required_set = true
+ end
+ o.separator ''
+ o.on("-b", "--ana-only", "Run only the analysis benchmark") do |a|
+ opt.run_ana_down = true
+ required_set = true
+ end
+ o.separator ''
+ o.on("-e", "--evgen [EVGEN_MODE]", "Run the event generation, reconstruction, and analysis",
+ "[EVGEN_MODE] must be one of:") do |a|
+ unless avail_evgens.include? a
+ $stderr.puts "ERROR: unknown event generation mode '#{a}'"
+ exit 1
+ end
+ opt.evgen = a
+ required_set = true
+ end
+ avail_evgens.each{ |it| o.separator ' '*40+it }
+ o.separator ''
+ o.separator 'optional options:'.upcase
+ o.separator ''
+ o.on("--sim-file [FILE]", "simulation file name", "default = #{opt.sim_file}") { |a| opt.sim_file=a }
+ o.on("--rec-file [FILE]", "reconstruction file name", "default = #{opt.rec_file}") { |a| opt.rec_file=a }
+ o.on("--ana-file [FILE]", "analysis file name", "default = #{opt.ana_file}") { |a| opt.ana_file=a }
+ o.separator ''
+ o.on("--[no-]run-sim", "simulation on/off", "default = #{opt.run_sim ? 'on' : 'off'}") { |a| opt.run_sim=a }
+ o.on("--[no-]run-rec", "reconstruction on/off", "default = #{opt.run_rec ? 'on' : 'off'}") { |a| opt.run_rec=a }
+ o.on("--[no-]run-ana", "analysis benchmark on/off", "default = #{opt.run_ana ? 'on' : 'off'}") { |a| opt.run_ana=a }
+ o.separator ''
+ o.on("-n", "--num-events [NUM_EVENTS]", Integer, "Number of events", "default = #{opt.num_events}") { |a| opt.num_events=a }
+ o.on("-a", "--algos [ALGORITHMS]...", Array, "List of analysis algorithms to run; default = all",
+ "delimit by commas, no spaces",
+ "for more info, run: #{opt.benchmark_exe}") { |a| opt.algos=a }
+ o.on("-p", "--plots [PLOTS]...", Array, "List of plots to draw",
+ "delimit by commas, no spaces",
+ "default = all") { |a| opt.plots=a }
+ o.on("-x", "--benchmark-exe [EXECUTABLE]", "benchmark executable",
+ "default = #{opt.benchmark_exe} (from $PATH)") { |a| opt.benchmark_exe=a}
+ o.separator ''
+ o.on("-v", "--verbose", "Increase verbosity (-vv for more verbose)") { |a| opt.verbosity+=1 }
+ o.on("-d", "--dry-run", "Dry run (just print commands)") { |a| opt.dry_run=true }
+ o.separator ''
+ o.on_tail("-h", "--help", "Show this message") do
+ puts o
+ exit 2
+ end
+end.parse!(ARGV.length>0 ? ARGV : ['--help'])
+
+# print options
+puts 'settings: {'.upcase
+opt.each_pair { |k,v| puts "#{k.to_s.rjust(20)} => #{v}" }
+puts '}'
+
+# check for required options
+unless required_set
+ $stderr.puts "ERROR: required options have not been set"
+ $stderr.puts "run '#{$0} --help' for guidance"
+ exit 1
+end
+
+# figure out which steps to run
+run_step = { :sim=>false, :rec=>false, :ana=>false, }
+if opt.run_ana_down
+ run_step[:ana] = true
+elsif opt.run_rec_down
+ run_step[:rec] = true
+ run_step[:ana] = true
+else
+ run_step[:sim] = opt.run_sim
+ run_step[:rec] = opt.run_rec
+ run_step[:ana] = opt.run_ana
+end
+puts "steps to run: #{run_step}"
+
+# get compact file
+if ENV['DETECTOR_PATH'].nil? or ENV['DETECTOR_CONFIG'].nil?
+ $stderr.puts "ERROR: unknown DETECTOR_PATH or DETECTOR_CONFIG"
+ exit 1
+end
+compact_file = "#{ENV['DETECTOR_PATH']}/#{ENV['DETECTOR_CONFIG']}.xml"
+
+
+# event generator command generators
+# ---------------------------------------------------
+def theta2xyz(theta)
+ rad = theta * Math::PI / 180
+ [ Math.sin(rad), 0.0, Math.cos(rad) ]
+end
+
+# fixed angle particle gun
+evgen_fixed_angle = Proc.new do |theta, energy, particle|
+ [
+ "npsim",
+ "--runType batch",
+ "--compactFile #{compact_file}",
+ "--outputFile #{opt.sim_file}",
+ "--part.userParticleHandler=''",
+ "--enableGun",
+ "--numberOfEvents #{opt.num_events}",
+ "--gun.particle #{particle}",
+ "--gun.energy #{energy}*GeV",
+ "--gun.direction (#{theta2xyz(theta).join ", "})",
+ ]
+end
+
+
+# define event generator command
+# ---------------------------------------------------
+evgen_cmd = Array.new
+case opt.evgen
+when 'idealAngle'
+ evgen_cmd = evgen_fixed_angle.call IdealTheta, IdealEnergy, IdealParticle
+else
+ exit 1 if run_step[:sim]
+end
+
+
+# define reconstruction command
+# ---------------------------------------------------
+output_collections = [
+ "DRICHHits",
+ "MCParticles",
+ "DRICHRawHits",
+ "DRICHRawHitsAssociations",
+ "DRICHAerogelTracks",
+ "DRICHGasTracks",
+ "DRICHAerogelIrtCherenkovParticleID",
+ "DRICHGasIrtCherenkovParticleID",
+ "DRICHMergedIrtCherenkovParticleID",
+ "ReconstructedChargedParticleAssociationsWithDRICHPID",
+]
+recon_cmd = [
+ 'eicrecon',
+ "-Ppodio:output_include_collections=\"#{output_collections.join ','}\"",
+ '-Pjana:nevents="0"',
+ '-Pjana:debug_plugin_loading="1"',
+ '-Pacts:MaterialMap="calibrations/materials-map.cbor"',
+ "-Ppodio:output_file=\"#{opt.rec_file}\"",
+ opt.sim_file,
+]
+
+
+# define analysis benchmark command
+# ---------------------------------------------------
+analysis_cmd = [
+ opt.benchmark_exe,
+ "-i #{opt.rec_file}",
+ "-o #{opt.ana_file}",
+]
+analysis_cmd.append "-a #{opt.algos.join ' '}" if opt.algos.size > 0
+analysis_cmd.append '-' + 'v'*opt.verbosity if opt.verbosity > 0
+
+# define analysis draw command
+# ---------------------------------------------------
+draw_cmd = [
+ "#{__dir__}/draw_benchmark.py",
+ "-i #{opt.ana_file}",
+ "-o #{opt.ana_file.gsub(/edm4hep.root$/,"plots")}"
+]
+
+# execute commands
+# ---------------------------------------------------
+
+# proc: execute a command; raise runtime exception if failure (exit nonzero)
+exe = Proc.new do |cmd_args, name, step|
+ if run_step[step]
+ cmd = cmd_args.join ' '
+ puts "benchmark #{name} command:".upcase
+ cmd_args.each_with_index do |arg,i|
+ line = i==0 ? '' : ' '
+ line += arg
+ line += ' \\' unless i+1==cmd_args.size
+ puts line
+ end
+ unless opt.dry_run
+ puts '-'*50
+ puts "#{name} execution:".upcase
+ system cmd or raise "benchmark #{name} failed!".upcase
+ end
+ puts '-'*50
+ end
+end
+puts '-'*50
+
+# execute the commands
+exe.call evgen_cmd, 'event generation', :sim
+exe.call recon_cmd, 'reconstruction', :rec
+exe.call analysis_cmd, 'analysis', :ana
+exe.call draw_cmd, 'draw', :ana
diff --git a/benchmarks/rich/scripts/rich_samples.py b/benchmarks/rich/scripts/rich_samples.py
deleted file mode 100644
index 0dcc2bdbcd5e73e31bb70fb4e0028dc37cf1112f..0000000000000000000000000000000000000000
--- a/benchmarks/rich/scripts/rich_samples.py
+++ /dev/null
@@ -1,44 +0,0 @@
-import ROOT
-import numpy as np
-from ROOT import gROOT
-import argparse
-from matplotlib import pyplot as plt
-from matplotlib.patches import Ellipse
-import matplotlib.transforms as transforms
-
-parser = argparse.ArgumentParser(description='rich reconstruction analysis')
-parser.add_argument('nevents', type=int, help='number of events')
-args = parser.parse_args()
-
-gROOT.Macro('rootlogon.C')
-f = ROOT.TFile("rich_test_reco.root")
-events = f.events
-
-nev = 0
-while nev < args.nevents:
- iev = np.random.randint(0, events.GetEntries())
- events.GetEntry(iev)
- if events.RICHClusters.size() == 0:
- continue
-
- nev += 1
- x, y = [], []
- for hit in events.RecoForwardRICHHits:
- x.append(hit.position.x)
- y.append(hit.position.y)
-
- fig, ax = plt.subplots(figsize=(9, 9), dpi=120)
- colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
- ax.scatter(x, y, s=9, marker='s', color=colors[0])
-
- for cluster in events.RICHClusters:
- ellipse = Ellipse((cluster.position.x, cluster.position.y), width=cluster.radius * 2, height=cluster.radius * 2,
- edgecolor='red', fill=False, lw=2, alpha=0.8)
- ellipse.set_transform(ax.transData)
- ax.add_patch(ellipse)
- ax.set_xlim(-60, 60)
- ax.set_ylim(-60, 60)
- ax.set_xlabel('X (mm)')
- ax.set_ylabel('Y (mm)')
- fig.savefig('rich_cluster_{}.png'.format(iev))
-
diff --git a/benchmarks/rich/src/ChargedParticle.cc b/benchmarks/rich/src/ChargedParticle.cc
new file mode 100644
index 0000000000000000000000000000000000000000..60fbd938b56481ee974ea4a2c5a5487bf5bd4597
--- /dev/null
+++ b/benchmarks/rich/src/ChargedParticle.cc
@@ -0,0 +1,28 @@
+// Copyright 2023, Christopher Dilks
+// Subject to the terms in the LICENSE file found in the top-level directory.
+
+#include "ChargedParticle.h"
+
+namespace benchmarks {
+
+ // set info from a single MCParticle
+ void ChargedParticle::SetSingleParticle(const edm4hep::MCParticle& mc_part) {
+ m_momentum = edm4hep::utils::magnitude(mc_part.getMomentum());
+ m_eta = edm4hep::utils::eta(mc_part.getMomentum());
+ m_pdg = mc_part.getPDG();
+ }
+
+ // set info from a (thrown) MCParticle of a collection
+ void ChargedParticle::SetSingleParticle(const edm4hep::MCParticleCollection& mc_parts) {
+ int n_thrown = 0;
+ for(const auto& mc_part : mc_parts) {
+ if(mc_part.getGeneratorStatus() == 1) {
+ this->SetSingleParticle(mc_part);
+ n_thrown++;
+ }
+ }
+ if(n_thrown == 0) m_log->warn("ChargedParticle::SetSingleParticle: no particle found");
+ if(n_thrown > 1) m_log->warn("ChargedParticle::SetSingleParticle: found {} particles (more than one)", n_thrown);
+ }
+
+}
diff --git a/benchmarks/rich/src/CherenkovPIDAnalysis.cc b/benchmarks/rich/src/CherenkovPIDAnalysis.cc
new file mode 100644
index 0000000000000000000000000000000000000000..795dbe814d9a7640b253fddec5f7570063a4ef47
--- /dev/null
+++ b/benchmarks/rich/src/CherenkovPIDAnalysis.cc
@@ -0,0 +1,234 @@
+// Copyright 2023, Christopher Dilks
+// Subject to the terms in the LICENSE file found in the top-level directory.
+
+#include "CherenkovPIDAnalysis.h"
+
+namespace benchmarks {
+
+ // AlgorithmInit
+ //---------------------------------------------------------------------------
+ void CherenkovPIDAnalysis::AlgorithmInit(
+ std::string rad_name,
+ std::shared_ptr<spdlog::logger>& logger
+ )
+ {
+ m_rad_name = rad_name;
+ m_log = logger;
+
+ // set radiator title
+ m_rad_title = TString(m_rad_name);
+ if(m_rad_name=="Merged") m_rad_title += " Aerogel+Gas";
+
+ // truncate `m_rad_name` (for object names)
+ m_rad_name_trun = TString(m_rad_name);
+ m_rad_name_trun(TRegexp(" .*")) = "";
+
+ // distributions
+ m_npe_dist = new TH1D(
+ "npe_dist_"+m_rad_name_trun,
+ "Overall NPE for "+m_rad_title+";NPE",
+ Tools::npe_max, 0, Tools::npe_max
+ );
+ m_theta_dist = new TH1D(
+ "theta_dist_"+m_rad_name_trun,
+ "Estimated Cherenkov Angle for "+m_rad_title+";#theta [mrad]",
+ Tools::theta_bins, 0, Tools::theta_max
+ );
+ m_thetaResid_dist = new TH1D(
+ "thetaResid_dist_"+m_rad_name_trun,
+ "Estimated Cherenkov Angle Residual for "+m_rad_title+";#Delta#theta [mrad]",
+ Tools::theta_bins, -Tools::thetaResid_max, Tools::thetaResid_max
+ );
+ m_photonTheta_vs_photonPhi = new TH2D(
+ "photonTheta_vs_photonPhi_"+m_rad_name_trun,
+ "Estimated Photon #theta vs #phi for "+m_rad_title+";#phi [rad];#theta [mrad]",
+ Tools::phi_bins, -TMath::Pi(), TMath::Pi(),
+ Tools::theta_bins, 0, Tools::theta_max
+ );
+
+ // truth distributions
+ m_mcWavelength_dist = new TH1D(
+ "mcWavelength_"+m_rad_name_trun,
+ "MC Photon Wavelength for "+m_rad_title+";#lambda [nm]",
+ Tools::n_bins, 0, 1000
+ );
+ m_mcRindex_dist = new TH1D(
+ "mcRindex_"+m_rad_name_trun,
+ "MC Refractive Index for "+m_rad_title+";n",
+ 10*Tools::n_bins, 0.99, 1.03
+ );
+
+ // PID
+ m_highestWeight_dist = new TH1D(
+ "highestWeight_dist_"+m_rad_name_trun,
+ "Highest PDG Weight for "+m_rad_title+";PDG",
+ Tools::GetNumPDGs()+1, 0, Tools::GetNumPDGs()+1
+ );
+
+ // momentum scans
+ m_npe_vs_p = new TH2D(
+ "npe_vs_p_"+m_rad_name_trun,
+ "Overall NPE vs. Particle Momentum for "+m_rad_title+";p [GeV];NPE",
+ Tools::momentum_bins, 0, Tools::momentum_max,
+ Tools::npe_max, 0, Tools::npe_max
+ );
+ m_theta_vs_p = new TH2D(
+ "theta_vs_p_"+m_rad_name_trun,
+ "Estimated Cherenkov Angle vs. Particle Momentum for "+m_rad_title+";p [GeV];#theta [mrad]",
+ Tools::momentum_bins, 0, Tools::momentum_max,
+ Tools::theta_bins, 0, Tools::theta_max
+ );
+ m_thetaResid_vs_p = new TH2D(
+ "thetaResid_vs_p_"+m_rad_name_trun,
+ "Estimated Cherenkov Angle Residual vs. Particle Momentum for "+m_rad_title+";p [GeV];#Delta#theta [mrad]",
+ Tools::momentum_bins, 0, Tools::momentum_max,
+ Tools::theta_bins, -Tools::thetaResid_max, Tools::thetaResid_max
+ );
+ m_highestWeight_vs_p = new TH2D(
+ "highestWeight_vs_p_"+m_rad_name_trun,
+ "Highest PDG Weight vs. Particle Momentum for "+m_rad_title+";p [GeV]",
+ Tools::momentum_bins, 0, Tools::momentum_max,
+ Tools::GetNumPDGs()+1, 0, Tools::GetNumPDGs()+1
+ );
+
+ // pseudorapidity scans
+ m_npe_vs_eta = new TH2D(
+ "npe_vs_eta_"+m_rad_name_trun,
+ "Overall NPE vs. Pseudorapidity for "+m_rad_title+";#eta;NPE",
+ Tools::eta_bins, Tools::eta_min, Tools::eta_max,
+ Tools::npe_max, 0, Tools::npe_max
+ );
+ m_theta_vs_eta = new TH2D(
+ "theta_vs_eta_"+m_rad_name_trun,
+ "Estimated Cherenkov Angle vs. Pseudorapidity for "+m_rad_title+";#eta;#theta [mrad]",
+ Tools::eta_bins, Tools::eta_min, Tools::eta_max,
+ Tools::theta_bins, 0, Tools::theta_max
+ );
+ m_thetaResid_vs_eta = new TH2D(
+ "thetaResid_vs_eta_"+m_rad_name_trun,
+ "Estimated Cherenkov Angle Residual vs. Pseudorapidity for "+m_rad_title+";#eta;#Delta#theta [mrad]",
+ Tools::eta_bins, Tools::eta_min, Tools::eta_max,
+ Tools::theta_bins, -Tools::thetaResid_max, Tools::thetaResid_max
+ );
+ m_highestWeight_vs_eta = new TH2D(
+ "highestWeight_vs_eta_"+m_rad_name_trun,
+ "Highest PDG Weight vs. Pseudorapidity for "+m_rad_title+";#eta",
+ Tools::eta_bins, Tools::eta_min, Tools::eta_max,
+ Tools::GetNumPDGs()+1, 0, Tools::GetNumPDGs()+1
+ );
+
+ // format histograms
+ auto format1D = [] (auto h) {
+ h->SetLineColor(kAzure);
+ h->SetFillColor(kAzure);
+ };
+ format1D(m_npe_dist);
+ format1D(m_theta_dist);
+ format1D(m_thetaResid_dist);
+ format1D(m_mcWavelength_dist);
+ format1D(m_mcRindex_dist);
+ format1D(m_highestWeight_dist);
+ }
+
+
+ // AlgorithmProcess
+ //---------------------------------------------------------------------------
+ void CherenkovPIDAnalysis::AlgorithmProcess(
+ const edm4hep::MCParticleCollection& mc_parts,
+ const edm4eic::CherenkovParticleIDCollection& cherenkov_pids
+ )
+ {
+ m_log->trace("-> {} Radiator:", m_rad_name);
+
+ // get thrown momentum from a single MCParticle
+ auto part = std::make_unique<ChargedParticle>(m_log);
+ part->SetSingleParticle(mc_parts);
+ auto thrown_momentum = part->GetMomentum();
+ auto thrown_eta = part->GetEta();
+ auto thrown_pdg = part->GetPDG();
+
+ // loop over `CherenkovParticleID` objects
+ for(const auto& cherenkov_pid : cherenkov_pids) {
+
+ // skip if NPE==0
+ if(cherenkov_pid.getNpe() == 0) {
+ m_log->warn("Event found with NPE=0");
+ continue;
+ }
+
+ // estimate the charged particle momentum using the momentum of the first TrackPoint at this radiator's entrance
+ /*
+ auto charged_particle = cherenkov_pid.getChargedParticle();
+ if(!charged_particle.isAvailable()) { m_log->warn("Charged particle not available in this radiator"); continue; }
+ if(charged_particle.points_size()==0) { m_log->warn("Charged particle has no TrackPoints in this radiator"); continue; }
+ auto charged_particle_momentum = edm4hep::utils::magnitude( charged_particle.getPoints(0).momentum );
+ m_log->trace(" Charged Particle p = {} GeV at radiator entrance", charged_particle_momentum);
+ m_log->trace(" If it is a pion, E = {} GeV", std::hypot(charged_particle_momentum, Tools::GetPDGMass(211)));
+ */
+ // alternatively: use `thrown_momentum` (FIXME: will not work for multi-track events)
+ auto charged_particle_momentum = thrown_momentum;
+ auto charged_particle_eta = thrown_eta;
+ auto charged_particle_pdg = thrown_pdg;
+ auto charged_particle_mass = Tools::GetPDGMass(charged_particle_pdg);
+ m_log->trace(" Charged Particle PDG={}, mass={} GeV, p={} GeV from truth",
+ charged_particle_pdg, charged_particle_mass, charged_particle_momentum);
+
+ // get average Cherenkov angle: `theta_rec`
+ double theta_rec = 0.0;
+ for(const auto& [theta,phi] : cherenkov_pid.getThetaPhiPhotons())
+ theta_rec += theta;
+ theta_rec /= cherenkov_pid.getNpe();
+ auto theta_rec_mrad = theta_rec * 1e3; // [rad] -> [mrad]
+
+ // calculate expected Cherenkov angle `theta_exp` and residual `theta_resid`,
+ // using refractive index from MC truth
+ auto mc_rindex = cherenkov_pid.getRefractiveIndex(); // average refractive index for photons used in this `cherenkov_pid`
+ auto theta_exp = TMath::ACos(
+ TMath::Hypot( charged_particle_momentum, charged_particle_mass ) /
+ ( mc_rindex * charged_particle_momentum )
+ );
+ auto theta_resid = theta_rec - theta_exp;
+ auto theta_resid_mrad = theta_resid * 1e3; // [rad] -> [mrad]
+
+ // fill PID histograms
+ m_npe_dist->Fill( cherenkov_pid.getNpe());
+ m_npe_vs_p->Fill( charged_particle_momentum, cherenkov_pid.getNpe());
+ m_npe_vs_eta->Fill( charged_particle_eta, cherenkov_pid.getNpe());
+ m_theta_dist->Fill( theta_rec_mrad);
+ m_theta_vs_p->Fill( charged_particle_momentum, theta_rec_mrad);
+ m_theta_vs_eta->Fill( charged_particle_eta, theta_rec_mrad);
+ m_thetaResid_dist->Fill( theta_resid_mrad);
+ m_thetaResid_vs_p->Fill( charged_particle_momentum, theta_resid_mrad);
+ m_thetaResid_vs_eta->Fill( charged_particle_eta, theta_resid_mrad);
+ for(const auto& [theta,phi] : cherenkov_pid.getThetaPhiPhotons())
+ m_photonTheta_vs_photonPhi->Fill( phi, theta*1e3); // [rad] -> [mrad]
+ m_mcWavelength_dist->Fill( Tools::HC / cherenkov_pid.getPhotonEnergy() ); // energy [GeV] -> wavelength [nm]
+ m_mcRindex_dist->Fill( mc_rindex);
+
+ // find the PDG hypothesis with the highest weight
+ float max_weight = -1000;
+ int pdg_max_weight = 0;
+ for(const auto& hyp : cherenkov_pid.getHypotheses()) {
+ if(hyp.weight > max_weight) {
+ max_weight = hyp.weight;
+ pdg_max_weight = hyp.PDG;
+ }
+ }
+ std::string pdg_max_weight_str = "UNKNOWN";
+ if(pdg_max_weight!=0 && !std::isnan(pdg_max_weight))
+ pdg_max_weight_str = std::to_string(pdg_max_weight);
+ m_log->trace(" Highest weight is {} for PDG {} (string='{}')", max_weight, pdg_max_weight, pdg_max_weight_str);
+ m_highestWeight_dist->Fill( pdg_max_weight_str.c_str(), 1);
+ m_highestWeight_vs_p->Fill( charged_particle_momentum, pdg_max_weight_str.c_str(), 1);
+ m_highestWeight_vs_eta->Fill( charged_particle_eta, pdg_max_weight_str.c_str(), 1);
+
+ } // end loop over `CherenkovParticleID` objects
+ }
+
+
+ // AlgorithmFinish
+ //---------------------------------------------------------------------------
+ void CherenkovPIDAnalysis::AlgorithmFinish() {
+ }
+
+}
diff --git a/benchmarks/rich/src/RawHitAnalysis.cc b/benchmarks/rich/src/RawHitAnalysis.cc
new file mode 100644
index 0000000000000000000000000000000000000000..d55ac99d7819489c6a408e537e15d2ca7455cd5f
--- /dev/null
+++ b/benchmarks/rich/src/RawHitAnalysis.cc
@@ -0,0 +1,82 @@
+// Copyright 2023, Christopher Dilks
+// Subject to the terms in the LICENSE file found in the top-level directory.
+
+#include "RawHitAnalysis.h"
+
+namespace benchmarks {
+
+ // AlgorithmInit
+ //---------------------------------------------------------------------------
+ void RawHitAnalysis::AlgorithmInit(
+ std::shared_ptr<spdlog::logger>& logger
+ )
+ {
+ m_log = logger;
+
+ // define histograms
+ m_adc_dist = new TH1D("adc_dist", "ADC Distribution;counts", adc_max, 0, adc_max);
+ m_tdc_dist = new TH1D("tdc_dist", "TDC Distribution;counts", tdc_max, 0, tdc_max);
+ m_tdc_vs_adc = new TH2D("tdc_vs_adc", "TDC vs. ADC;ADC counts;TDC counts",
+ adc_max/2, 0, adc_max/2,
+ tdc_max/2, 0, tdc_max/2
+ );
+ m_phot_spectrum = new TH1D(
+ "phot_spectrum_rec",
+ "Photon wavelength for digitized hits;#lambda [nm], at emission point",
+ Tools::wl_bins, Tools::wl_min, Tools::wl_max
+ );
+ m_nhits_dist = new TH1D("nhits_dist", "Number of digitized hits;N_{hits}",
+ Tools::npe_max, 0, Tools::npe_max
+ );
+
+ // format histograms
+ auto format1D = [] (auto h) {
+ h->SetLineColor(kBlack);
+ h->SetFillColor(kBlack);
+ };
+ format1D(m_adc_dist);
+ format1D(m_tdc_dist);
+ m_phot_spectrum->SetLineColor(kViolet+2);
+ m_phot_spectrum->SetFillColor(kViolet+2);
+ m_nhits_dist->SetLineColor(kCyan+2);
+ m_nhits_dist->SetFillColor(kCyan+2);
+ }
+
+
+ // AlgorithmProcess
+ //---------------------------------------------------------------------------
+ void RawHitAnalysis::AlgorithmProcess(
+ const edm4eic::RawTrackerHitCollection& raw_hits,
+ const edm4eic::MCRecoTrackerHitAssociationCollection& assocs
+ )
+ {
+
+ // fill nhits
+ m_nhits_dist->Fill(raw_hits.size());
+
+ // loop over all raw hits (including noise)
+ for(const auto& raw_hit : raw_hits) {
+ auto adc = raw_hit.getCharge();
+ auto tdc = raw_hit.getTimeStamp();
+ m_adc_dist->Fill(adc);
+ m_tdc_dist->Fill(tdc);
+ m_tdc_vs_adc->Fill(adc,tdc);
+ }
+
+ // loop over hits with associations (no noise)
+ for(const auto& assoc : assocs) {
+ for(const auto& hit : assoc.getSimHits()) {
+ auto wavelength = Tools::GetPhotonWavelength(hit, m_log);
+ if(wavelength>=0)
+ m_phot_spectrum->Fill(wavelength);
+ }
+ }
+ }
+
+
+ // AlgorithmFinish
+ //---------------------------------------------------------------------------
+ void RawHitAnalysis::AlgorithmFinish() {
+ }
+
+}
diff --git a/benchmarks/rich/src/ReconstructedParticleAnalysis.cc b/benchmarks/rich/src/ReconstructedParticleAnalysis.cc
new file mode 100644
index 0000000000000000000000000000000000000000..974551701976825d134d6766dae8c73c46a227f3
--- /dev/null
+++ b/benchmarks/rich/src/ReconstructedParticleAnalysis.cc
@@ -0,0 +1,98 @@
+// Copyright 2023, Christopher Dilks
+// Subject to the terms in the LICENSE file found in the top-level directory.
+
+#include "ReconstructedParticleAnalysis.h"
+
+namespace benchmarks {
+
+ // AlgorithmInit
+ //---------------------------------------------------------------------------
+ void ReconstructedParticleAnalysis::AlgorithmInit(std::shared_ptr<spdlog::logger>& logger) {
+ m_log = logger;
+
+ // initialize plots
+ m_correct_vs_p_transient = new TH2D("correct_vs_p_transient", "Fraction with Correct PDG vs. Thrown Momentum;p [GeV];fraction",
+ Tools::momentum_bins, 0, Tools::momentum_max,
+ 2, 0, 2
+ );
+ m_correct_vs_p = new TGraphErrors();
+ m_correct_vs_p->SetName("correct_vs_p");
+ m_correct_vs_p->SetTitle(m_correct_vs_p_transient->GetTitle());
+ m_correct_vs_p->GetXaxis()->SetTitle(m_correct_vs_p_transient->GetXaxis()->GetTitle());
+ m_correct_vs_p->GetYaxis()->SetTitle(m_correct_vs_p_transient->GetYaxis()->GetTitle());
+ m_correct_vs_p->SetMarkerStyle(kFullCircle);
+ }
+
+
+ // AlgorithmProcess
+ //---------------------------------------------------------------------------
+ void ReconstructedParticleAnalysis::AlgorithmProcess(
+ const edm4eic::MCRecoParticleAssociationCollection& in_assocs
+ )
+ {
+ if(!in_assocs.isValid()) { m_log->error("invalid input collection 'in_assocs'"); return; }
+
+ // loop over input associations
+ for(const auto& assoc : in_assocs) {
+
+ // get particles
+ m_log->trace("{:-^50}"," Particle ");
+ auto& simpart = assoc.getSim();
+ auto& recpart = assoc.getRec();
+ if(!recpart.isAvailable()) { m_log->warn("reconstructed particle not available for this association"); continue; }
+ if(!simpart.isAvailable()) { m_log->warn("simulated particle not available for this association"); continue; }
+
+ // get PDG values
+ auto simpart_pdg = simpart.getPDG(); // from MC truth
+ auto recpart_pdg = recpart.getPDG(); // PDG member of ReconstructedParticle
+ auto recpart_pid_used = recpart.getParticleIDUsed();
+ auto pid_pdg = recpart_pid_used.isAvailable() ? recpart_pid_used.getPDG() : 0;
+ if(!recpart_pid_used.isAvailable())
+ m_log->error("reconstructed particle does not have particleIDUsed relation");
+ m_log->trace("PDGs: sim.PDG | rec.PDG | PDG from PID = {:^6} | {:^6} | {:^6}", simpart_pdg, recpart_pdg, pid_pdg);
+
+ // check goodness of PID: if near zero, assume PID was not used for this particle and skip
+ if(recpart.getGoodnessOfPID()<0.01) {
+ m_log->warn("skip particle with low goodnessOfPID (likely has no PID)");
+ continue;
+ }
+
+ // get momenta
+ auto simpart_p = edm4hep::utils::p(simpart);
+ // auto recpart_p = edm4hep::utils::p(recpart); // unused
+
+ // fill plots
+ m_correct_vs_p_transient->Fill(
+ simpart_p,
+ simpart_pdg == pid_pdg ? 1 : 0
+ );
+
+ } // end loop over input associations
+
+ }
+
+
+ // AlgorithmFinish
+ //---------------------------------------------------------------------------
+ void ReconstructedParticleAnalysis::AlgorithmFinish() {
+
+ // compute fraction of correct PDG for each x-axis bin
+ for(int bx=1; bx<=m_correct_vs_p_transient->GetNbinsX(); bx++) {
+ auto n_wrong = m_correct_vs_p_transient->GetBinContent(bx,1);
+ auto n_right = m_correct_vs_p_transient->GetBinContent(bx,2);
+ auto n_total = n_wrong + n_right;
+ auto momentum = m_correct_vs_p_transient->GetXaxis()->GetBinCenter(bx);
+ if(n_total>0) {
+ auto frac = n_right / n_total;
+ m_correct_vs_p->AddPoint(momentum, frac);
+ }
+ }
+
+ // delete transient histograms, so they don't get written
+ delete m_correct_vs_p_transient;
+
+ // write objects which are not automatically saved
+ m_correct_vs_p->Write();
+ }
+
+}
diff --git a/benchmarks/rich/src/SimHitAnalysis.cc b/benchmarks/rich/src/SimHitAnalysis.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f204b5125073a3d5ae3297e49638897979662513
--- /dev/null
+++ b/benchmarks/rich/src/SimHitAnalysis.cc
@@ -0,0 +1,98 @@
+// Copyright 2023, Christopher Dilks
+// Subject to the terms in the LICENSE file found in the top-level directory.
+
+#include "SimHitAnalysis.h"
+
+namespace benchmarks {
+
+ // AlgorithmInit
+ //---------------------------------------------------------------------------
+ void SimHitAnalysis::AlgorithmInit(
+ std::shared_ptr<spdlog::logger>& logger
+ )
+ {
+ m_log = logger;
+
+ // initialize histograms
+ m_nphot_dist = new TH1D("nphot_dist", "Number of incident photons;N_{photons}",
+ Tools::nphot_max, 0, Tools::nphot_max
+ );
+ m_nphot_vs_p = new TH2D("nphot_vs_p", "Number of incident photons vs. Thrown Momentum;p [GeV];N_{photons}",
+ Tools::momentum_bins, 0, Tools::momentum_max,
+ Tools::nphot_max, 0, Tools::nphot_max
+ );
+ m_nphot_vs_p__transient = new TH1D("nphot_vs_p__transient", "",
+ m_nphot_vs_p->GetNbinsX(),
+ m_nphot_vs_p->GetXaxis()->GetXmin(),
+ m_nphot_vs_p->GetXaxis()->GetXmax()
+ );
+ m_phot_spectrum = new TH1D(
+ "phot_spectrum_sim",
+ "Incident photon wavelength;#lambda [nm], at emission point",
+ Tools::wl_bins, Tools::wl_min, Tools::wl_max
+ );
+ m_phot_spectrum->SetLineColor(kViolet+2);
+ m_phot_spectrum->SetFillColor(kViolet+2);
+ m_nphot_dist->SetLineColor(kCyan+2);
+ m_nphot_dist->SetFillColor(kCyan+2);
+ }
+
+
+ // AlgorithmProcess
+ //---------------------------------------------------------------------------
+ void SimHitAnalysis::AlgorithmProcess(
+ const edm4hep::MCParticleCollection& mc_parts,
+ const edm4hep::SimTrackerHitCollection& sim_hits
+ )
+ {
+
+ // get thrown momentum from a single MCParticle
+ auto part = std::make_unique<ChargedParticle>(m_log);
+ part->SetSingleParticle(mc_parts);
+ auto thrown_momentum = part->GetMomentum();
+
+ // fill nphot distribution
+ m_nphot_dist->Fill(sim_hits.size());
+
+ // loop over hits
+ for(const auto& sim_hit : sim_hits) {
+
+ // - fill 1D thrown momentum histogram `m_nphot_vs_p__transient` for each (pre-digitized) sensor hit
+ // FIXME: get thrown momentum from multi-track events; but in this attempt
+ // below the parent momentum is not consistent; instead for now we use
+ // `thrown_momentum` from truth above
+ /*
+ auto photon = sim_hit.getMCParticle();
+ if(photon.parents_size()>0) thrown_momentum = edm4hep::utils::p(photon.getParents(0));
+ */
+ m_nphot_vs_p__transient->Fill(thrown_momentum);
+
+ // fill photon spectrum
+ auto wavelength = Tools::GetPhotonWavelength(sim_hit, m_log);
+ if(wavelength>=0)
+ m_phot_spectrum->Fill(wavelength);
+ }
+
+ // use `m_nphot_vs_p__transient` results to fill 2D hist `m_nphot_vs_p` for this event
+ for(int b=1; b<=m_nphot_vs_p__transient->GetNbinsX(); b++) {
+ auto nphot = m_nphot_vs_p__transient->GetBinContent(b);
+ if(nphot>0) {
+ auto momentum = m_nphot_vs_p__transient->GetBinCenter(b);
+ m_nphot_vs_p->Fill(momentum,nphot);
+ }
+ }
+
+ // clear `m_nphot_vs_p__transient` to be ready for the next event
+ m_nphot_vs_p__transient->Reset();
+
+ }
+
+
+ // AlgorithmFinish
+ //---------------------------------------------------------------------------
+ void SimHitAnalysis::AlgorithmFinish() {
+ // delete transient histograms, so they don't get written
+ delete m_nphot_vs_p__transient;
+ }
+
+}
diff --git a/benchmarks/rich/src/benchmark.cc b/benchmarks/rich/src/benchmark.cc
new file mode 100644
index 0000000000000000000000000000000000000000..d988601f1727dd71d7e6a13548e4aa021b470359
--- /dev/null
+++ b/benchmarks/rich/src/benchmark.cc
@@ -0,0 +1,299 @@
+// Copyright 2023, Christopher Dilks
+// Subject to the terms in the LICENSE file found in the top-level directory.
+
+#include <functional>
+#include <iostream>
+#include <unistd.h>
+#include <spdlog/spdlog.h>
+
+#include <TFile.h>
+
+#include <podio/ROOTFrameReader.h>
+#include <podio/Frame.h>
+
+#include "SimHitAnalysis.h"
+#include "RawHitAnalysis.h"
+#include "CherenkovPIDAnalysis.h"
+#include "ReconstructedParticleAnalysis.h"
+
+using namespace benchmarks;
+
+// -------------------------------------------------------------
+
+// main logger
+std::shared_ptr<spdlog::logger> m_log;
+
+// get a collection, and check its validity
+template<class C>
+const C& GetCollection(podio::Frame& frame, std::string name) {
+ const auto& coll = frame.get<C>(name);
+ if(!coll.isValid())
+ m_log->error("invalid collection '{}'", name);
+ return coll;
+}
+
+// -------------------------------------------------------------
+// MAIN
+// -------------------------------------------------------------
+int main(int argc, char** argv) {
+
+ // -------------------------------------------------------------
+ // setup and CLI option parsing
+ // -------------------------------------------------------------
+
+ // loggers
+ m_log = spdlog::default_logger()->clone("benchmark_rich");
+ auto log_level = spdlog::level::info;
+ m_log->set_level(log_level);
+
+ // default options
+ int verbosity = 0;
+ long long num_events = 0;
+ std::string ana_file_name = "out_rich.root";
+ std::vector<std::string> rec_files;
+ std::vector<std::string> analysis_algorithms = {
+ "SimHit",
+ "RawHit",
+ "CherenkovPID",
+ "ReconstructedParticle"
+ };
+ auto list_algorithms = [&analysis_algorithms] () {
+ for(auto analysis_algorithm : analysis_algorithms)
+ std::cout << " " << analysis_algorithm << std::endl;
+ std::cout << std::endl;
+ };
+
+ // usage guide
+ auto PrintUsage = [&argv, &list_algorithms, ana_file_name] () {
+ std::cout << "\nUSAGE: " << argv[0] << " -i [INPUT FILES]... [OPTIONS]..." << std::endl
+ << "\nREQUIRED OPTIONS:\n"
+ << " -i [INPUT FILES]...\n"
+ << " input file from reconstruction output;\n"
+ << " can specify more than one (delimited by space)\n"
+ << "\nOPTIONAL OPTIONS:\n"
+ << " -o [OUTPUT FILE]\n"
+ << " output analysis file (default: " << ana_file_name << ")\n"
+ << " -n [NUMBER OF EVENTS]\n"
+ << " number of events to process (default: process all)\n"
+ << " -v: verbose output\n"
+ << " -vv: more verbose output\n"
+ << " -a [ALGORITHMS]...\n"
+ << " list of analysis algorithms to run (delimited by spaces)\n"
+ << " default: run all of them\n"
+ << " available algorithms:\n";
+ list_algorithms();
+ };
+ if(argc==1) {
+ PrintUsage();
+ return 2;
+ }
+
+ // parse options
+ int opt;
+ while( (opt=getopt(argc, argv, "i:o:n:v|a:")) != -1) {
+ switch(opt) {
+ case 'i':
+ optind--;
+ for(; optind<argc && *argv[optind]!='-'; optind++)
+ rec_files.push_back(std::string(argv[optind]));
+ break;
+ case 'o':
+ ana_file_name = std::string(optarg);
+ break;
+ case 'n':
+ num_events = std::atoll(optarg);
+ break;
+ case 'v':
+ verbosity++;
+ break;
+ case 'a':
+ analysis_algorithms.clear();
+ optind--;
+ for(; optind<argc && *argv[optind]!='-'; optind++)
+ analysis_algorithms.push_back(std::string(argv[optind]));
+ break;
+ default:
+ m_log->error("Unknown option '{}'", opt);
+ PrintUsage();
+ return 1;
+ }
+ }
+
+ // print options
+ m_log->info("{:-^50}", " User options ");
+ m_log->info(" input files:");
+ for(auto rec_file : rec_files)
+ m_log->info(" {}", rec_file);
+ m_log->info(" output file: {}", ana_file_name);
+ m_log->info(" algorithms:");
+ for(auto analysis_algorithm : analysis_algorithms)
+ m_log->info(" {}", analysis_algorithm);
+ m_log->info(" num_events: {}", num_events);
+ m_log->info(" verbosity: {}", verbosity);
+ m_log->info("{:-^50}", "");
+
+ // check options
+ std::vector<std::string> bad_options;
+ if(rec_files.size()==0) bad_options.push_back("no [INPUT FILES] have been specified");
+ if(analysis_algorithms.size()==0) bad_options.push_back("no [ALGORITHMS] have been specified");
+ for(auto bad_option : bad_options) m_log->error(bad_option);
+ if(bad_options.size()>0) {
+ PrintUsage();
+ return 1;
+ }
+
+ // set log level
+ switch(verbosity) {
+ case 0: log_level = spdlog::level::info; break;
+ case 1: log_level = spdlog::level::debug; break;
+ default: log_level = spdlog::level::trace;
+ }
+ m_log->set_level(log_level);
+
+ // start the output file
+ auto ana_file = new TFile(TString(ana_file_name), "RECREATE");
+
+ // struct to hold an algorithm execution processor
+ struct AnalysisAlgorithm {
+ std::shared_ptr<spdlog::logger> log; // algorithm-specific logger
+ std::function<void(podio::Frame&)> process; // call AlgorithmProcess
+ std::function<void()> finish; // call AlgorithmFinish
+ };
+ std::vector<AnalysisAlgorithm> algorithm_processors;
+
+ // loop over algorithms to run, and define how to run them
+ for(auto analysis_algorithm : analysis_algorithms) {
+ AnalysisAlgorithm algo;
+
+ // algorithm logger
+ algo.log = spdlog::default_logger()->clone(analysis_algorithm);
+ algo.log->set_level(log_level);
+ algo.log->debug("{:-<50}","INITIALIZE ");
+
+ // --------------------------------------------------------------
+ // define how to run each algorithm
+ // --------------------------------------------------------------
+
+ // truth hits (photons) .........................................
+ if(analysis_algorithm == "SimHit") {
+ auto sim_algo = std::make_shared<SimHitAnalysis>();
+ ana_file->mkdir("phot")->cd();
+ sim_algo->AlgorithmInit(algo.log);
+ algo.process = [sim_algo] (podio::Frame& frame) {
+ sim_algo->AlgorithmProcess(
+ GetCollection<edm4hep::MCParticleCollection>(frame,"MCParticles"),
+ GetCollection<edm4hep::SimTrackerHitCollection>(frame,"DRICHHits")
+ );
+ };
+ algo.finish = [sim_algo] () { sim_algo->AlgorithmFinish(); };
+ }
+
+ // digitizer ....................................................
+ else if(analysis_algorithm == "RawHit") {
+ auto digi_algo = std::make_shared<RawHitAnalysis>();
+ ana_file->mkdir("digi")->cd();
+ digi_algo->AlgorithmInit(algo.log);
+ algo.process = [digi_algo] (podio::Frame& frame) {
+ digi_algo->AlgorithmProcess(
+ GetCollection<edm4eic::RawTrackerHitCollection>(frame,"DRICHRawHits"),
+ GetCollection<edm4eic::MCRecoTrackerHitAssociationCollection>(frame,"DRICHRawHitsAssociations")
+ );
+ };
+ algo.finish = [digi_algo] () { digi_algo->AlgorithmFinish(); };
+ }
+
+ // IRT ..........................................................
+ else if(analysis_algorithm == "CherenkovPID") {
+ std::vector<std::string> radiators = { "Aerogel", "Gas", "Merged" };
+ std::map<std::string, std::shared_ptr<CherenkovPIDAnalysis>> irt_algos;
+ for(auto radiator : radiators)
+ irt_algos.insert({radiator, std::make_shared<CherenkovPIDAnalysis>()});
+
+ for(auto& [radiator, irt_algo] : irt_algos) {
+ ana_file->mkdir(("pid"+radiator).c_str())->cd();
+ irt_algo->AlgorithmInit(radiator, algo.log);
+ }
+
+ algo.process = [irt_algos] (podio::Frame& frame) {
+ const auto& mc_parts = GetCollection<edm4hep::MCParticleCollection>(frame,"MCParticles");
+ for(auto& [radiator, irt_algo] : irt_algos) {
+ const auto& cherenkov_pids = GetCollection<edm4eic::CherenkovParticleIDCollection>(
+ frame,
+ "DRICH" + radiator + "IrtCherenkovParticleID"
+ );
+ irt_algo->AlgorithmProcess(mc_parts, cherenkov_pids);
+ }
+ };
+
+ algo.finish = [irt_algos] () {
+ for(auto& [radiator, irt_algo] : irt_algos)
+ irt_algo->AlgorithmFinish();
+ };
+
+ }
+
+ // linking to reconstructed particles ...........................
+ else if(analysis_algorithm == "ReconstructedParticle") {
+ auto link_algo = std::make_shared<ReconstructedParticleAnalysis>();
+ ana_file->mkdir("link")->cd();
+ link_algo->AlgorithmInit(algo.log);
+ algo.process = [link_algo] (podio::Frame& frame) {
+ link_algo->AlgorithmProcess(
+ GetCollection<edm4eic::MCRecoParticleAssociationCollection>(
+ frame,
+ "ReconstructedChargedParticleAssociations"
+ )
+ );
+ };
+ algo.finish = [link_algo] () { link_algo->AlgorithmFinish(); };
+ }
+
+ // unknown algorithm ............................................
+ else {
+ m_log->error("Unknown [ALGORITHM] '{}'; ignoring", analysis_algorithm);
+ continue;
+ }
+
+ // --------------------------------------------------------------
+
+ algorithm_processors.push_back(algo);
+ }
+
+ // -------------------------------------------------------------
+ // read the input file and run all algorithms
+ // -------------------------------------------------------------
+
+ // open the input files
+ podio::ROOTFrameReader podioReader;
+ m_log->warn("podio::ROOTFrameReader cannot yet support multiple files; reading only the first");
+ // podioReader.openFiles(rec_files);
+ podioReader.openFile(rec_files.front());
+
+ // get the number of events to process
+ const std::string tree_name = "events";
+ long long num_entries = podioReader.getEntries(tree_name);
+ if(num_events>0) num_entries = std::min(num_events, num_entries);
+
+ // event loop
+ m_log->info("BEGIN EVENT LOOP");
+ for(long long e=0; e<num_entries; e++) {
+ if(e%100==0) m_log->info("Progress: {:.3f}%", 100.0 * e/num_entries);
+ auto podioFrame = podio::Frame(podioReader.readNextEntry(tree_name));
+ for(auto& algo : algorithm_processors) {
+ algo.log->debug("{:-<50}","PROCESS EVENT ");
+ algo.process(podioFrame);
+ }
+ }
+
+ // finish
+ for(auto& algo : algorithm_processors) {
+ algo.log->debug("{:-<50}","FINISH ");
+ algo.finish();
+ }
+
+ // write output
+ ana_file->Write();
+ ana_file->Close();
+ m_log->info("Done. Wrote analysis file:");
+ m_log->info(" {}", ana_file_name);
+}