add "rich" benchmark from reconstruction_benchmarks (#172)

.gitlab-ci.yml

@@ -134,6 +134,7 @@ include:
   - local: 'benchmarks/zdc_pi0/config.yml'
   - local: 'benchmarks/material_scan/config.yml'
   - local: 'benchmarks/pid/config.yml'
+  - local: 'benchmarks/rich/config.yml'
   - local: 'benchmarks/timing/config.yml'
   - local: 'benchmarks/b0_tracker/config.yml'
   - local: 'benchmarks/insert_muon/config.yml'
@@ -157,6 +158,7 @@ deploy_results:
     - "collect_results:lfhcal"
     - "collect_results:material_scan"
     - "collect_results:pid"
+    - "collect_results:rich"
     - "collect_results:tracking_performance"
     - "collect_results:tracking_performance_campaigns"
     - "collect_results:zdc_sigma"

benchmarks/rich/CMakeLists.txt

+cmake_minimum_required(VERSION 3.0.0 FATAL_ERROR)
+
+get_filename_component(det_name ${CMAKE_CURRENT_LIST_DIR} NAME)
+project(${CMAKE_PROJECT_NAME}_${det_name} 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)
+
+# 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_lib}
+  PUBLIC
+  ROOT::Core
+  ROOT::Hist
+  podio::podio
+  podio::podioRootIO
+  EDM4EIC::edm4eic
+  EDM4HEP::edm4hep
+  spdlog::spdlog
+  )
+target_link_libraries(${algo_exe}
+  PRIVATE
+  ${algo_lib}
+  )
+
+# installation
+install(FILES ${algo_headers} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/${CMAKE_PROJECT_NAME}/${det_name})
+install(TARGETS ${algo_exe} ${algo_lib})

benchmarks/rich/README.md

+# 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

benchmarks/rich/config.yml

+bench:drich_fixed_eta:
+  extends: .det_benchmark
+  stage: benchmarks
+  script:
+    - cmake -S benchmarks/rich -B build -DCMAKE_INSTALL_PREFIX=.local
+    - cmake --build build -j2
+    - cmake --install build
+    - |
+      for mode in fixedEtaIdeal fixedEtaMin fixedEtaMax ; do
+        ruby benchmarks/rich/run_benchmark.rb --ci -s $mode --num-events 50
+      done
+
+collect_results:rich:
+  extends: .det_benchmark
+  stage: collect
+  needs:
+    - "bench:drich_fixed_eta"
+  script:
+    - echo "Done collecting artifacts."

benchmarks/rich/draw_benchmark.py

+#!/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()

benchmarks/rich/include/ChargedParticle.h

+// 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;
+
+  };
+
+}

benchmarks/rich/include/CherenkovPIDAnalysis.h

+// 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;
+
+  };
+
+}

benchmarks/rich/include/RawHitAnalysis.h

+// 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;
+
+  };
+
+}

benchmarks/rich/include/ReconstructedParticleAnalysis.h

+// 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;
+
+  };
+
+}

benchmarks/rich/include/SimHitAnalysis.h

+// 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;
+
+  };
+
+}

benchmarks/rich/include/Tools.h

+// 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/EDM4hepVersion.h>
+#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) {
+#if EDM4HEP_BUILD_VERSION >= EDM4HEP_VERSION(0, 99, 0)
+        auto phot = hit.getParticle();
+#else
+        auto phot = hit.getMCParticle();
+#endif
+        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

benchmarks/rich/run_benchmark.rb

+#!/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'
+require 'fileutils'
+
+###################################
+# constants:
+EtaTestValues = {
+  :ideal => 2.0,
+  :min   => 1.6,
+  :max   => 3.5,
+}
+IdealEnergy   = 12.0 # [GeV]
+IdealParticle = 'pi+'
+###################################
+
+# setup
+# ---------------------------------------------------
+
+# default opt
+opt = OpenStruct.new
+opt.sim_mode      = ''
+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
+opt.using_ci      = false
+
+# available simulation modes
+avail_sim_modes = [
+  'fixedEtaIdeal',
+  'fixedEtaMin',
+  'fixedEtaMax',
+]
+
+# 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("-s", "--sim-mode [SIMULATION_MODE]", "Run the simulation, reconstruction, and analysis",
+       "[SIMULATION_MODE] must be one of:") do |a|
+    unless avail_sim_modes.include? a
+      $stderr.puts "ERROR: unknown simulation mode '#{a}'"
+      exit 1
+    end
+    opt.sim_mode = a
+    required_set = true
+  end
+  avail_sim_modes.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.on("--ci", "output plots to ./results, for CI artifact collection") { |a| opt.using_ci=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"
+
+
+# simulation command generators
+# ---------------------------------------------------
+def theta2xyz(theta)
+  [ Math.sin(theta), 0.0, Math.cos(theta) ]
+end
+
+def eta2theta(eta)
+  2 * Math.atan(Math.exp(-eta))
+end
+
+
+# fixed angle particle gun
+simulate_fixed_angle = Proc.new do |theta, energy, particle|
+  [
+    "npsim",
+    "--runType batch",
+    "--compactFile #{compact_file}",
+    "--outputFile #{opt.sim_file}",
+    "--part.userParticleHandler=''",  # allow opticalphotons in MC particles
+    "--enableGun",
+    "--numberOfEvents #{opt.num_events}",
+    "--gun.particle #{particle}",
+    "--gun.energy #{energy}*GeV",
+    "--gun.direction '(#{theta2xyz(theta).join ", "})'",
+  ]
+end
+
+
+# define simulation command
+# ---------------------------------------------------
+sim_cmd = Array.new
+case opt.sim_mode
+when /^fixedEta/
+  key       = opt.sim_mode.sub('fixedEta','').downcase.to_sym
+  fixed_eta = EtaTestValues[key]
+  if fixed_eta.nil?
+    $stderr.puts "ERROR: EtaTestValues[#{key}] is not defined"
+    exit 1
+  end
+  fixed_theta = eta2theta fixed_eta
+  puts """Simulating fixed-eta events:
+  - eta      = #{fixed_eta}
+  - theta    = #{180.0 * fixed_theta / Math::PI} degrees
+  - energy   = #{IdealEnergy} GeV
+  - particle = #{IdealParticle}
+  """
+  sim_cmd = simulate_fixed_angle.call fixed_theta, 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_collections=\"#{output_collections.join ','}\"",
+  '-Pjana:nevents="0"',
+  '-Pjana:debug_plugin_loading="1"',
+  '-Pacts:MaterialMap="calibrations/materials-map.cbor"',
+  "-Ppodio:output_file=\"#{opt.rec_file}\"",
+  '-PDRICH:DRICHIrtCherenkovParticleID:cheatPhotonVertex=true',  # allow knowledge of true photons, for accurate residual determination
+  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.using_ci ? "results/#{opt.sim_mode}" : 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]
+    case step
+    when :sim
+      FileUtils.mkdir_p File.dirname(opt.sim_file)
+    when :rec
+      FileUtils.mkdir_p File.dirname(opt.rec_file)
+    when :ana
+      FileUtils.mkdir_p File.dirname(opt.ana_file)
+    end
+    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 sim_cmd,      'simulation',     :sim
+exe.call recon_cmd,    'reconstruction', :rec
+exe.call analysis_cmd, 'analysis',       :ana
+exe.call draw_cmd,     'draw',           :ana

benchmarks/rich/scripts/rich_data_gen.py

+import os
+from pyHepMC3 import HepMC3 as hm
+import numpy as np
+import argparse
+
+
+PARTICLES = [
+#    (111, 0.134.9766),      # pi0
+    (211, 0.13957018),      # pi+
+#    (-211, 0.13957018),     # pi-
+#    (311, 0.497648),        # K0
+    (321, 0.493677),        # K+
+#    (-321, 0.493677),       # K-
+    (2212, 0.938272),       # proton
+#    (2112, 0.939565),       # neutron
+#    (11, 0.51099895e-3),    # electron
+#    (-11, 0.51099895e-3),   # positron
+]
+
+
+# p in GeV, angle in degree, vertex in mm
+def gen_event(prange=(8, 100), arange=(0, 20)):
+    evt = hm.GenEvent(hm.Units.MomentumUnit.GEV, hm.Units.LengthUnit.MM)
+    pid, mass = PARTICLES[np.random.randint(len(PARTICLES))]
+
+    # final state
+    state = 1
+
+    # momentum, angles, energy
+    p = np.random.uniform(*prange)
+    theta = np.random.uniform(*arange)*np.pi/180.
+    phi = np.random.uniform(0., 2.*np.pi)
+    e0 = np.sqrt(p*p + mass*mass)
+
+    px = np.cos(phi)*np.sin(theta)
+    py = np.sin(phi)*np.sin(theta)
+    pz = np.cos(theta)
+
+    # beam
+    pbeam = hm.GenParticle(hm.FourVector(0, 0, 0, 0.938272), 2212, 4)
+    ebeam = hm.GenParticle(hm.FourVector(0, 0, e0, np.sqrt(e0*e0 + 0.511e-3*0.511e-3)), -11, 4)
+
+    hout = hm.GenParticle(hm.FourVector(px*p, py*p, pz*p, e0), pid, state)
+    # evt.add_particle(part)
+    vert = hm.GenVertex()
+    vert.add_particle_in(ebeam)
+    vert.add_particle_in(pbeam)
+    vert.add_particle_out(hout)
+    evt.add_vertex(vert)
+    return evt
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser('RICH dataset generator')
+
+    parser.add_argument('output', help='path to the output file')
+    parser.add_argument('-n', type=int, default=1000, dest='nev', help='number of events to generate')
+    parser.add_argument('--pmin', type=float, default=8.0, dest='pmin', help='minimum momentum in GeV')
+    parser.add_argument('--pmax', type=float, default=100.0, dest='pmax', help='maximum momentum in GeV')
+    parser.add_argument('--angmin', type=float, default=0.0, dest='angmin', help='minimum angle in degree')
+    parser.add_argument('--angmax', type=float, default=20.0, dest='angmax', help='maximum angle in degree')
+    parser.add_argument('-s', type=int, default=-1, dest='seed', help='seed for random generator')
+
+    args = parser.parse_args()
+
+    # random seed
+    if args.seed < 0:
+        args.seed = os.environ.get('SEED', int.from_bytes(os.urandom(4), byteorder='big', signed=False))
+    np.random.seed(args.seed)
+
+    output = hm.WriterAscii(args.output);
+    if output.failed():
+        print("Cannot open file \"{}\"".format(args.output))
+        sys.exit(2)
+
+    count = 0
+    while count < args.nev:
+        if (count % 1000 == 0):
+            print("Generated {} events".format(count), end='\r')
+        evt = gen_event((args.pmin, args.pmax), (args.angmin, args.angmax))
+        output.write_event(evt)
+        evt.clear()
+        count += 1
+
+    print("Generated {} events".format(args.nev))
+    output.close()
+

benchmarks/rich/src/ChargedParticle.cc

+// 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);
+  }
+
+}

benchmarks/rich/src/CherenkovPIDAnalysis.cc

+// 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() {
+  }
+
+}

benchmarks/rich/src/RawHitAnalysis.cc

+// Copyright 2023, Christopher Dilks
+// Subject to the terms in the LICENSE file found in the top-level directory.
+
+#include "RawHitAnalysis.h"
+#include <edm4eic/EDM4eicVersion.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) {
+#if EDM4EIC_VERSION_MAJOR >= 6
+      auto hit = assoc.getSimHit();
+#else
+    for(const auto& hit : assoc.getSimHits()) {
+#endif
+      auto wavelength = Tools::GetPhotonWavelength(hit, m_log);
+      if(wavelength>=0)
+        m_phot_spectrum->Fill(wavelength);
+#if EDM4EIC_VERSION_MAJOR >= 6
+#else
+    }
+#endif
+    }
+  }
+
+
+  // AlgorithmFinish
+  //---------------------------------------------------------------------------
+  void RawHitAnalysis::AlgorithmFinish() {
+  }
+
+}

benchmarks/rich/src/ReconstructedParticleAnalysis.cc

+// 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();
+  }
+
+}

benchmarks/rich/src/SimHitAnalysis.cc

+// 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;
+  }
+
+}

benchmarks/rich/src/benchmark.cc

+// 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/podioVersion.h>
+#if podio_VERSION >= PODIO_VERSION(0, 99, 0)
+#include <podio/ROOTReader.h>
+#else
+#include <podio/ROOTFrameReader.h>
+#endif
+#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
+#if podio_VERSION >= PODIO_VERSION(0, 99, 0)
+  podio::ROOTReader podioReader;
+#else
+  podio::ROOTFrameReader podioReader;
+#endif
+  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);
+}