Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-07-03 07:55:53

0001 // Copyright 2023, Christopher Dilks, adapted from Alexander Kiselev's Juggler implementation `IRTAlgorithm`
0002 // Subject to the terms in the LICENSE file found in the top-level directory.
0003 
0004 #include "IrtCherenkovParticleID.h"
0005 
0006 #include <IRT/ChargedParticle.h>
0007 #include <IRT/CherenkovPID.h>
0008 #include <IRT/OpticalPhoton.h>
0009 #include <IRT/RadiatorHistory.h>
0010 #include <IRT/SinglePDF.h>
0011 #include <TString.h>
0012 #include <TVector3.h>
0013 #include <algorithms/logger.h>
0014 #include <edm4eic/CherenkovParticleIDHypothesis.h>
0015 #include <edm4eic/EDM4eicVersion.h>
0016 #include <edm4eic/TrackPoint.h>
0017 #include <edm4hep/EDM4hepVersion.h>
0018 #include <edm4hep/MCParticleCollection.h>
0019 #include <edm4hep/SimTrackerHitCollection.h>
0020 #include <edm4hep/Vector2f.h>
0021 #include <fmt/core.h>
0022 #include <fmt/format.h>
0023 #include <cmath>
0024 
0025 #include <algorithm>
0026 #include <cstddef>
0027 #include <functional>
0028 #include <gsl/pointers>
0029 #include <iterator>
0030 #include <memory>
0031 #include <podio/ObjectID.h>
0032 #include <podio/RelationRange.h>
0033 #include <set>
0034 #include <stdexcept>
0035 #include <utility>
0036 #include <vector>
0037 
0038 #include "algorithms/pid/IrtCherenkovParticleIDConfig.h"
0039 #include "algorithms/pid/Tools.h"
0040 
0041 namespace eicrecon {
0042 
0043 void IrtCherenkovParticleID::init(CherenkovDetectorCollection* irt_det_coll) {
0044   // members
0045   m_irt_det_coll = irt_det_coll;
0046 
0047   // print the configuration parameters
0048   debug() << m_cfg << endmsg;
0049 
0050   // inform the user if a cheat mode is enabled
0051   auto print_cheat = [this](auto name, bool val, auto desc) {
0052     if (val) {
0053       warning("CHEAT MODE '{}' ENABLED: {}", name, desc);
0054     }
0055   };
0056   print_cheat("cheatPhotonVertex", m_cfg.cheatPhotonVertex,
0057               "use MC photon vertex, wavelength, refractive index");
0058   print_cheat("cheatTrueRadiator", m_cfg.cheatTrueRadiator, "use MC truth to obtain true radiator");
0059 
0060   // extract the the relevant `CherenkovDetector`, set to `m_irt_det`
0061   const auto& detectors = m_irt_det_coll->GetDetectors();
0062   if (detectors.empty()) {
0063     throw std::runtime_error("No CherenkovDetectors found in input collection `irt_det_coll`");
0064   }
0065   if (detectors.size() > 1) {
0066     warning("IrtCherenkovParticleID currently only supports 1 CherenkovDetector at a time; "
0067             "taking the first");
0068   }
0069   auto this_detector = *detectors.begin();
0070   m_det_name         = this_detector.first;
0071   m_irt_det          = this_detector.second;
0072   debug("Initializing IrtCherenkovParticleID algorithm for CherenkovDetector '{}'", m_det_name);
0073 
0074   // readout decoding
0075   m_cell_mask = m_irt_det->GetReadoutCellMask();
0076   debug("readout cellMask = {:#X}", m_cell_mask);
0077 
0078   // rebin refractive index tables to have `m_cfg.numRIndexBins` bins
0079   trace("Rebinning refractive index tables to have {} bins", m_cfg.numRIndexBins);
0080   for (auto [rad_name, irt_rad] : m_irt_det->Radiators()) {
0081     auto ri_lookup_table_orig = irt_rad->m_ri_lookup_table;
0082     irt_rad->m_ri_lookup_table.clear();
0083     irt_rad->m_ri_lookup_table = Tools::ApplyFineBinning(ri_lookup_table_orig, m_cfg.numRIndexBins);
0084     // trace("- {}", rad_name);
0085     // for(auto [energy,rindex] : irt_rad->m_ri_lookup_table) trace("  {:>5} eV   {:<}", energy, rindex);
0086   }
0087 
0088   // build `m_pid_radiators`, the list of radiators to use for PID
0089   debug("Obtain List of Radiators:");
0090   for (auto [rad_name, irt_rad] : m_irt_det->Radiators()) {
0091     if (rad_name != "Filter") {
0092       m_pid_radiators.insert({std::string(rad_name), irt_rad});
0093       debug("- {}", rad_name.Data());
0094     }
0095   }
0096 
0097   // check radiators' configuration, and pass it to `m_irt_det`'s radiators
0098   for (auto [rad_name, irt_rad] : m_pid_radiators) {
0099     // find `cfg_rad`, the associated `IrtCherenkovParticleIDConfig` radiator
0100     auto cfg_rad_it = m_cfg.radiators.find(rad_name);
0101     if (cfg_rad_it != m_cfg.radiators.end()) {
0102       auto cfg_rad = cfg_rad_it->second;
0103       // pass `cfg_rad` params to `irt_rad`, the IRT radiator
0104       irt_rad->m_ID                     = Tools::GetRadiatorID(std::string(rad_name));
0105       irt_rad->m_AverageRefractiveIndex = cfg_rad.referenceRIndex;
0106       irt_rad->SetReferenceRefractiveIndex(cfg_rad.referenceRIndex);
0107       if (cfg_rad.attenuation > 0) {
0108         irt_rad->SetReferenceAttenuationLength(cfg_rad.attenuation);
0109       }
0110       if (cfg_rad.smearing > 0) {
0111         if (cfg_rad.smearingMode == "uniform") {
0112           irt_rad->SetUniformSmearing(cfg_rad.smearing);
0113         } else if (cfg_rad.smearingMode == "gaussian") {
0114           irt_rad->SetGaussianSmearing(cfg_rad.smearing);
0115         } else {
0116           error("Unknown smearing mode '{}' for {} radiator", cfg_rad.smearingMode, rad_name);
0117         }
0118       }
0119     } else {
0120       error("Cannot find radiator '{}' in IrtCherenkovParticleIDConfig instance", rad_name);
0121     }
0122   }
0123 
0124   // get PDG info for the particles we want to identify in PID
0125   debug("List of particles for PID:");
0126   for (auto pdg : m_cfg.pdgList) {
0127     auto mass = m_particleSvc.particle(pdg).mass;
0128     m_pdg_mass.insert({pdg, mass});
0129     debug("  {:>8}  M={} GeV", pdg, mass);
0130   }
0131 }
0132 
0133 void IrtCherenkovParticleID::process(const IrtCherenkovParticleID::Input& input,
0134                                      const IrtCherenkovParticleID::Output& output) const {
0135   const auto [in_aerogel_tracks, in_gas_tracks, in_merged_tracks, in_raw_hits, in_hit_assocs] =
0136       input;
0137   auto [out_aerogel_particleIDs, out_gas_particleIDs] = output;
0138 
0139   // logging
0140   trace("{:=^70}", " call IrtCherenkovParticleID::AlgorithmProcess ");
0141   trace("number of raw sensor hits: {}", in_raw_hits->size());
0142   trace("number of raw sensor hit with associated photons: {}", in_hit_assocs->size());
0143 
0144   std::map<std::string, const edm4eic::TrackSegmentCollection*> in_charged_particles{
0145       {"Aerogel", in_aerogel_tracks},
0146       {"Gas", in_gas_tracks},
0147       {"Merged", in_merged_tracks},
0148   };
0149 
0150   // start output collections
0151   std::map<std::string, edm4eic::CherenkovParticleIDCollection*> out_cherenkov_pids{
0152       {"Aerogel", out_aerogel_particleIDs}, {"Gas", out_gas_particleIDs}};
0153 
0154   // check `in_charged_particles`: each radiator should have the same number of TrackSegments
0155   std::unordered_map<std::size_t, std::size_t> in_charged_particle_size_distribution;
0156   for (const auto& [rad_name, in_charged_particle] : in_charged_particles) {
0157     ++in_charged_particle_size_distribution[in_charged_particle->size()];
0158   }
0159   if (in_charged_particle_size_distribution.size() != 1) {
0160     std::vector<std::size_t> in_charged_particle_sizes;
0161     std::transform(
0162         in_charged_particles.begin(), in_charged_particles.end(),
0163         std::back_inserter(in_charged_particle_sizes),
0164         [](const auto& in_charged_particle) { return in_charged_particle.second->size(); });
0165     error("radiators have differing numbers of TrackSegments {}",
0166           fmt::join(in_charged_particle_sizes, ", "));
0167     return;
0168   }
0169 
0170   // loop over charged particles ********************************************
0171   trace("{:#<70}", "### CHARGED PARTICLES ");
0172   std::size_t num_charged_particles = in_charged_particle_size_distribution.begin()->first;
0173   for (std::size_t i_charged_particle = 0; i_charged_particle < num_charged_particles;
0174        i_charged_particle++) {
0175     trace("{:-<70}", fmt::format("--- charged particle #{} ", i_charged_particle));
0176 
0177     // start an `irt_particle`, for `IRT`
0178     auto irt_particle = std::make_unique<ChargedParticle>();
0179 
0180     // loop over radiators
0181     // note: this must run exclusively since irt_rad points to shared IRT objects that are
0182     // owned by the RichGeo_service; it holds state (e.g. irt_rad->ResetLocation())
0183     std::lock_guard<std::mutex> lock(m_pid_radiators_mutex);
0184     for (auto [rad_name, irt_rad] : m_pid_radiators) {
0185 
0186       // get the `charged_particle` for this radiator
0187       auto charged_particle_list_it = in_charged_particles.find(rad_name);
0188       if (charged_particle_list_it == in_charged_particles.end()) {
0189         error("Cannot find radiator '{}' in `in_charged_particles`", rad_name);
0190         continue;
0191       }
0192       const auto* charged_particle_list = charged_particle_list_it->second;
0193       auto charged_particle             = charged_particle_list->at(i_charged_particle);
0194 
0195       // set number of bins for this radiator and charged particle
0196       if (charged_particle.points_size() == 0) {
0197         trace("No propagated track points in radiator '{}'", rad_name);
0198         continue;
0199       }
0200       irt_rad->SetTrajectoryBinCount(charged_particle.points_size() - 1);
0201 
0202       // start a new IRT `RadiatorHistory`
0203       // - must be a raw pointer for `irt` compatibility
0204       // - it will be destroyed when `irt_particle` is destroyed
0205       auto* irt_rad_history = new RadiatorHistory();
0206       irt_particle->StartRadiatorHistory({irt_rad, irt_rad_history});
0207 
0208       // loop over `TrackPoint`s of this `charged_particle`, adding each to the IRT radiator
0209       irt_rad->ResetLocations();
0210       trace("TrackPoints in '{}' radiator:", rad_name);
0211       for (const auto& point : charged_particle.getPoints()) {
0212         TVector3 position = Tools::PodioVector3_to_TVector3(point.position);
0213         TVector3 momentum = Tools::PodioVector3_to_TVector3(point.momentum);
0214         irt_rad->AddLocation(position, momentum);
0215         trace(Tools::PrintTVector3(" point: x", position));
0216         trace(Tools::PrintTVector3("        p", momentum));
0217       }
0218 
0219       // loop over raw hits ***************************************************
0220       trace("{:#<70}", "### SENSOR HITS ");
0221       for (const auto& raw_hit : *in_raw_hits) {
0222 
0223         // get MC photon(s), typically only used by cheat modes or trace logging
0224         // - loop over `in_hit_assocs`, searching for the matching hit association
0225         // - will not exist for noise hits
0226         edm4hep::MCParticle mc_photon;
0227         bool mc_photon_found = false;
0228         if (m_cfg.cheatPhotonVertex || m_cfg.cheatTrueRadiator) {
0229           for (const auto& hit_assoc : *in_hit_assocs) {
0230             if (hit_assoc.getRawHit().isAvailable()) {
0231               if (hit_assoc.getRawHit().id() == raw_hit.id()) {
0232 #if EDM4EIC_VERSION_MAJOR >= 6
0233 #if EDM4HEP_BUILD_VERSION >= EDM4HEP_VERSION(0, 99, 0)
0234                 mc_photon = hit_assoc.getSimHit().getParticle();
0235 #else
0236                 mc_photon = hit_assoc.getSimHit().getMCParticle();
0237 #endif
0238 #else
0239                 // hit association found, get the MC photon and break the loop
0240                 if (hit_assoc.simHits_size() > 0) {
0241                   mc_photon = hit_assoc.getSimHits(0).getMCParticle();
0242 #endif
0243                 mc_photon_found = true;
0244                 if (mc_photon.getPDG() != -22) {
0245                   warning("non-opticalphoton hit: PDG = {}", mc_photon.getPDG());
0246                 }
0247 #if EDM4EIC_VERSION_MAJOR >= 6
0248 #else
0249                 } else if (m_cfg.CheatModeEnabled())
0250                   error("cheat mode enabled, but no MC photons provided");
0251 #endif
0252                 break;
0253               }
0254             }
0255           }
0256         }
0257 
0258         // cheat mode, for testing only: use MC photon to get the actual radiator
0259         if (m_cfg.cheatTrueRadiator && mc_photon_found) {
0260           auto vtx     = Tools::PodioVector3_to_TVector3(mc_photon.getVertex());
0261           auto* mc_rad = m_irt_det->GuessRadiator(vtx, vtx); // assume IP is at (0,0,0)
0262           if (mc_rad != irt_rad) {
0263             continue; // skip this photon, if not from radiator `irt_rad`
0264           }
0265           trace(Tools::PrintTVector3(
0266               fmt::format("cheat: radiator '{}' determined from photon vertex", rad_name), vtx));
0267         }
0268 
0269         // get sensor and pixel info
0270         // FIXME: signal and timing cuts (ADC, TDC, ToT, ...)
0271         auto cell_id       = raw_hit.getCellID();
0272         uint64_t sensor_id = cell_id & m_cell_mask;
0273         TVector3 pixel_pos = m_irt_det->m_ReadoutIDToPosition(cell_id);
0274 
0275         // trace logging
0276         if (level() <= algorithms::LogLevel::kTrace) {
0277           trace("cell_id={:#X}  sensor_id={:#X}", cell_id, sensor_id);
0278           trace(Tools::PrintTVector3("pixel position", pixel_pos));
0279           if (mc_photon_found) {
0280             TVector3 mc_endpoint = Tools::PodioVector3_to_TVector3(mc_photon.getEndpoint());
0281             trace(Tools::PrintTVector3("photon endpoint", mc_endpoint));
0282             trace("{:>30} = {}", "dist( pixel,  photon )", (pixel_pos - mc_endpoint).Mag());
0283           } else {
0284             trace("  no MC photon found; probably a noise hit");
0285           }
0286         }
0287 
0288         // start new IRT photon
0289         auto* irt_sensor = m_irt_det->m_PhotonDetectors[0]; // NOTE: assumes one sensor type
0290         auto* irt_photon =
0291             new OpticalPhoton(); // new raw pointer; it will also be destroyed when `irt_particle` is destroyed
0292         irt_photon->SetVolumeCopy(sensor_id);
0293         irt_photon->SetDetectionPosition(pixel_pos);
0294         irt_photon->SetPhotonDetector(irt_sensor);
0295         irt_photon->SetDetected(true);
0296 
0297         // cheat mode: get photon vertex info from MC truth
0298         if ((m_cfg.cheatPhotonVertex || m_cfg.cheatTrueRadiator) && mc_photon_found) {
0299           irt_photon->SetVertexPosition(Tools::PodioVector3_to_TVector3(mc_photon.getVertex()));
0300           irt_photon->SetVertexMomentum(Tools::PodioVector3_to_TVector3(mc_photon.getMomentum()));
0301         }
0302 
0303         // cheat mode: retrieve a refractive index estimate; it is not exactly the one, which
0304         // was used in GEANT, but should be very close
0305         if (m_cfg.cheatPhotonVertex) {
0306           double ri   = NAN;
0307           auto mom    = 1e9 * irt_photon->GetVertexMomentum().Mag();
0308           auto ri_set = Tools::GetFinelyBinnedTableEntry(irt_rad->m_ri_lookup_table, mom, &ri);
0309           if (ri_set) {
0310             irt_photon->SetVertexRefractiveIndex(ri);
0311             trace("{:>30} = {}", "refractive index", ri);
0312           } else {
0313             warning("Tools::GetFinelyBinnedTableEntry failed to lookup refractive index for "
0314                     "momentum {} eV",
0315                     mom);
0316           }
0317         }
0318 
0319         // add each `irt_photon` to the radiator history
0320         // - unless cheating, we don't know which photon goes with which
0321         // radiator, thus we add them all to each radiator; the radiators'
0322         // photons are mixed in `ChargedParticle::PIDReconstruction`
0323         irt_rad_history->AddOpticalPhoton(irt_photon);
0324         /* FIXME: this considers ALL of the `irt_photon`s... we can limit this
0325          * once we add the ability to get a fiducial volume for each track, i.e.,
0326          * a region of sensors where we expect to see this `irt_particle`'s
0327          * Cherenkov photons; this should also combat sensor noise
0328          */
0329       } // end `in_hit_assocs` loop
0330 
0331     } // end radiator loop
0332 
0333     // particle identification +++++++++++++++++++++++++++++++++++++++++++++++++++++
0334 
0335     // define a mass hypothesis for each particle we want to check
0336     trace("{:+^70}", " PARTICLE IDENTIFICATION ");
0337     CherenkovPID irt_pid;
0338     std::unordered_map<int, MassHypothesis*> pdg_to_hyp; // `pdg` -> hypothesis
0339     for (auto [pdg, mass] : m_pdg_mass) {
0340       irt_pid.AddMassHypothesis(mass);
0341       pdg_to_hyp.insert({pdg, irt_pid.GetHypothesis(irt_pid.GetHypothesesCount() - 1)});
0342     }
0343 
0344     // run IRT PID
0345     irt_particle->PIDReconstruction(irt_pid);
0346     trace("{:-^70}", " IRT RESULTS ");
0347 
0348     // loop over radiators
0349     for (auto [rad_name, irt_rad] : m_pid_radiators) {
0350       trace("-> {} Radiator (ID={}):", rad_name, irt_rad->m_ID);
0351 
0352       // Cherenkov angle (theta) estimate
0353       unsigned npe      = 0;
0354       double rindex_ave = 0.0;
0355       double energy_ave = 0.0;
0356       std::vector<std::pair<double, double>> phot_theta_phi;
0357 
0358       // loop over this radiator's photons, and decide which to include in the theta estimate
0359       auto* irt_rad_history = irt_particle->FindRadiatorHistory(irt_rad);
0360       if (irt_rad_history == nullptr) {
0361         trace("  No radiator history; skip");
0362         continue;
0363       }
0364       trace("  Photoelectrons:");
0365       for (auto* irt_photon : irt_rad_history->Photons()) {
0366 
0367         // check whether this photon was selected by at least one mass hypothesis
0368         bool photon_selected = false;
0369         for (auto irt_photon_sel : irt_photon->_m_Selected) {
0370           if (irt_photon_sel.second == irt_rad) {
0371             photon_selected = true;
0372             break;
0373           }
0374         }
0375         if (!photon_selected) {
0376           continue;
0377         }
0378 
0379         // trace logging
0380         trace(
0381             Tools::PrintTVector3(fmt::format("- sensor_id={:#X}: hit", irt_photon->GetVolumeCopy()),
0382                                  irt_photon->GetDetectionPosition()));
0383         trace(Tools::PrintTVector3("photon vertex", irt_photon->GetVertexPosition()));
0384 
0385         // get this photon's theta and phi estimates
0386         auto phot_theta = irt_photon->_m_PDF[irt_rad].GetAverage();
0387         auto phot_phi   = irt_photon->m_Phi[irt_rad];
0388 
0389         // add to the total
0390         npe++;
0391         phot_theta_phi.emplace_back(phot_theta, phot_phi);
0392         if (m_cfg.cheatPhotonVertex) {
0393           rindex_ave += irt_photon->GetVertexRefractiveIndex();
0394           energy_ave += irt_photon->GetVertexMomentum().Mag();
0395         }
0396 
0397       } // end loop over this radiator's photons
0398 
0399       // compute averages
0400       if (npe > 0) {
0401         rindex_ave /= npe;
0402         energy_ave /= npe;
0403       }
0404 
0405       // fill photon info
0406       auto out_cherenkov_pid = out_cherenkov_pids.at(rad_name)->create();
0407       out_cherenkov_pid.setNpe(static_cast<decltype(edm4eic::CherenkovParticleIDData::npe)>(npe));
0408       out_cherenkov_pid.setRefractiveIndex(
0409           static_cast<decltype(edm4eic::CherenkovParticleIDData::refractiveIndex)>(rindex_ave));
0410       out_cherenkov_pid.setPhotonEnergy(
0411           static_cast<decltype(edm4eic::CherenkovParticleIDData::photonEnergy)>(energy_ave));
0412       for (auto [phot_theta, phot_phi] : phot_theta_phi) {
0413         out_cherenkov_pid.addToThetaPhiPhotons(
0414             edm4hep::Vector2f{static_cast<float>(phot_theta), static_cast<float>(phot_phi)});
0415       }
0416 
0417       // relate mass hypotheses
0418       for (auto [pdg, mass] : m_pdg_mass) {
0419 
0420         // get hypothesis results
0421         auto* irt_hypothesis = pdg_to_hyp.at(pdg);
0422         auto hyp_weight      = irt_hypothesis->GetWeight(irt_rad);
0423         auto hyp_npe         = irt_hypothesis->GetNpe(irt_rad);
0424 
0425         // fill `ParticleID` output collection
0426         edm4eic::CherenkovParticleIDHypothesis out_hypothesis;
0427         out_hypothesis.PDG =
0428             static_cast<decltype(edm4eic::CherenkovParticleIDHypothesis::PDG)>(pdg);
0429         out_hypothesis.weight =
0430             static_cast<decltype(edm4eic::CherenkovParticleIDHypothesis::weight)>(hyp_weight);
0431         out_hypothesis.npe =
0432             static_cast<decltype(edm4eic::CherenkovParticleIDHypothesis::npe)>(hyp_npe);
0433 
0434         // relate
0435         out_cherenkov_pid.addToHypotheses(out_hypothesis);
0436 
0437       } // end hypothesis loop
0438 
0439       // logging: Cherenkov angle estimate
0440       auto PrintCherenkovEstimate = [this](edm4eic::CherenkovParticleID pid,
0441                                            bool printHypotheses = true, int indent = 2) {
0442         double thetaAve = 0;
0443         if (pid.getNpe() > 0) {
0444           for (const auto& [theta, phi] : pid.getThetaPhiPhotons()) {
0445             thetaAve += theta / pid.getNpe();
0446           }
0447         }
0448         trace("{:{}}Cherenkov Angle Estimate:", "", indent);
0449         trace("{:{}}  {:>16}:  {:>10}", "", indent, "NPE", pid.getNpe());
0450         trace("{:{}}  {:>16}:  {:>10.8} mrad", "", indent, "<theta>",
0451               thetaAve * 1e3); // [rad] -> [mrad]
0452         trace("{:{}}  {:>16}:  {:>10.8}", "", indent, "<rindex>", pid.getRefractiveIndex());
0453         trace("{:{}}  {:>16}:  {:>10.8} eV", "", indent, "<energy>",
0454               pid.getPhotonEnergy() * 1e9); // [GeV] -> [eV]
0455         if (printHypotheses) {
0456           trace("{:{}}Mass Hypotheses:", "", indent);
0457           trace("{}", Tools::HypothesisTableHead(indent + 2));
0458           for (const auto& hyp : pid.getHypotheses()) {
0459             trace("{}", Tools::HypothesisTableLine(hyp, indent + 2));
0460           }
0461         }
0462       };
0463       PrintCherenkovEstimate(out_cherenkov_pid);
0464 
0465       // relate charged particle projection
0466       auto charged_particle_list_it = in_charged_particles.find("Merged");
0467       if (charged_particle_list_it != in_charged_particles.end()) {
0468         const auto* charged_particle_list = charged_particle_list_it->second;
0469         auto charged_particle             = charged_particle_list->at(i_charged_particle);
0470         out_cherenkov_pid.setChargedParticle(charged_particle);
0471       } else {
0472         error("Cannot find radiator 'Merged' in `in_charged_particles`");
0473       }
0474 
0475       // relate hit associations
0476       for (const auto& hit_assoc : *in_hit_assocs) {
0477         out_cherenkov_pid.addToRawHitAssociations(hit_assoc);
0478       }
0479 
0480     } // end radiator loop
0481 
0482     /* NOTE: `unique_ptr irt_particle` goes out of scope and will now be destroyed, and along with it:
0483      * - raw pointer `irt_rad_history` for each radiator
0484      * - all `irt_photon` raw pointers
0485      */
0486 
0487   } // end `in_charged_particles` loop
0488 }
0489 
0490 } // namespace eicrecon