algorithms/pid/IrtInterface.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2024, 2025, 2026, Alexander Kiselev
0003
0004 #if IRT2_VERSION_MAJOR > 2 || (IRT2_VERSION_MAJOR == 2 && IRT2_VERSION_MINOR >= 2)
0005
0006 #include <DD4hep/Detector.h>
0007 #include <DD4hep/Objects.h>
0008 #include <Evaluator/DD4hepUnits.h>
0009 #include <IRT2/ChargedParticle.h>
0010 #include <IRT2/ChargedParticleStep.h>
0011 #include <IRT2/CherenkovDetector.h>
0012 #include <IRT2/CherenkovDetectorCollection.h>
0013 #include <IRT2/CherenkovEvent.h>
0014 #include <IRT2/CherenkovPhotonDetector.h>
0015 #include <IRT2/CherenkovRadiator.h>
0016 #include <IRT2/DataInterpolation.h>
0017 #include <IRT2/OpticalPhoton.h>
0018 #include <IRT2/RadiatorHistory.h>
0019 #include <IRT2/ReconstructionFactory.h>
0020 #include <TGDMLMatrix.h>
0021 #include <TString.h>
0022 #include <TVector3.h>
0023 #include <edm4eic/TrackCollection.h>
0024 #include <edm4eic/TrackPoint.h>
0025 #include <edm4eic/unit_system.h>
0026 #include <edm4hep/SimTrackerHitCollection.h>
0027 #include <edm4hep/Vector3d.h>
0028 #include <edm4hep/Vector3f.h>
0029 #include <edm4hep/utils/vector_utils.h>
0030 #include <podio/LinkNavigator.h>
0031 #include <podio/ObjectID.h>
0032 #include <podio/RelationRange.h>
0033 #include <stdint.h>
0034 #include <gsl/pointers>
0035 #include <map>
0036 #include <mutex>
0037 #include <random>
0038 #include <set>
0039 #include <tuple>
0040 #include <utility>
0041 #include <vector>
0042
0043 #include "algorithms/pid/IrtInterfaceConfig.h"
0044 #include "algorithms/pid/Tools.h"
0045
0046 using namespace IRT2;
0047
0048 #include "IrtInterface.h"
0049
0050 // -------------------------------------------------------------------------------------
0051
0052 namespace eicrecon {
0053
0054 std::mutex IrtInterface::m_irt_geom_mutex;
0055
0056 IrtInterface::~IrtInterface() {
0057   if (!m_irt_detector)
0058     return;
0059   std::lock_guard<std::mutex> lock(m_irt_geom_mutex);
0060   for (auto [name, rad] : m_irt_detector->Radiators()) {
0061     if (rad && rad->m_RefractiveIndex) {
0062       delete rad->m_RefractiveIndex;
0063       rad->m_RefractiveIndex = nullptr;
0064     }
0065   }
0066 }
0067
0068 // -------------------------------------------------------------------------------------
0069
0070 void IrtInterface::init() {
0071   // Cannot fail (see RICH-IRT.cc);
0072   m_irt_geometry = IRT2::CherenkovDetectorCollection::Instance();
0073   m_irt_detector = m_irt_geometry->GetDetector(m_cfg.m_detector_name.c_str());
0074
0075   m_Event = std::make_unique<IRT2::CherenkovEvent>();
0076
0077   const dd4hep::Detector* det = m_geo.detector();
0078
0079   std::lock_guard<std::mutex> lock(m_irt_geom_mutex);
0080
0081 #if IRT2_VERSION_MAJOR == 2 && IRT2_VERSION_MINOR == 2 && IRT2_VERSION_PATCH < 1
0082   // A workaround for when missing https://github.com/eic/irt/pull/61
0083   {
0084     static std::set<IRT2::CherenkovDetector*> zeroed;
0085     if (zeroed.insert(m_irt_detector).second) {
0086       for (auto [name, rad] : m_irt_detector->Radiators())
0087         rad->m_RefractiveIndex = nullptr;
0088     }
0089   }
0090 #endif
0091
0092   m_ReconstructionFactory = std::make_unique<IRT2::ReconstructionFactory>(
0093       m_irt_geometry, m_irt_detector, m_Event.get(), m_cfg.m_json_config_file_name.c_str());
0094   // JANA2 prints out event progress; the rest is kind of irrelevant;
0095   m_ReconstructionFactory->SetQuietMode();
0096
0097   for (auto [name, rad] : m_irt_detector->Radiators()) {
0098     if (rad->m_RefractiveIndex)
0099       continue;
0100     const auto* rindex_matrix = det->material(rad->GetAlternativeMaterialName()).property("RINDEX");
0101     if (rindex_matrix) {
0102       const unsigned dim = rindex_matrix->GetRows();
0103       std::unique_ptr<double[]> e(new double[dim]);
0104       std::unique_ptr<double[]> ri(new double[dim]);
0105
0106       for (unsigned row = 0; row < rindex_matrix->GetRows(); row++) {
0107         e[row]  = rindex_matrix->Get(row, 0) / dd4hep::eV;
0108         ri[row] = rindex_matrix->Get(row, 1);
0109       } //for row
0110
0111       auto ptr = rad->m_RefractiveIndex = new DataInterpolation(e.get(), ri.get(), dim);
0112       // FIXME: 100 hardcoded;
0113       ptr->CreateLookupTable(100);
0114     } //if
0115   } //for radiators
0116 } // IrtInterface::init()
0117
0118 // -------------------------------------------------------------------------------------
0119
0120 void IrtInterface::process(const IrtInterface::Input& input,
0121                            const IrtInterface::Output& output) const {
0122   // Reset the event structure;
0123   m_Event->Reset();
0124
0125   // Intermediate variables, for less typing;
0126   const auto [headers, in_mc_particles, in_tracks, in_track_associations, in_track_projections,
0127               in_sim_hits]                        = input;
0128   auto [out_irt_radiator_info, out_irt_particles] = output;
0129
0130   // local random generator
0131   auto seed = m_uid.getUniqueID(*headers, name());
0132   std::default_random_engine generator(seed);
0133   std::uniform_real_distribution<double> uniform(0.0, 1.0);
0134
0135   // Build fast MC->reco lookup using podio::LinkNavigator;
0136   const podio::LinkNavigator<edm4eic::MCRecoTrackParticleLinkCollection> link_nav(
0137       *in_track_associations);
0138
0139   // Then track -> track projection lookup table; FIXME: other radiators;
0140   std::map<podio::ObjectID, edm4eic::TrackSegment> Track_to_TrackSegment_lut;
0141   for (auto segment : *in_track_projections) {
0142     auto track                                     = segment.getTrack();
0143     Track_to_TrackSegment_lut[track.getObjectID()] = segment;
0144   } //for particle
0145
0146   // Help optical photons to find their parents;
0147   std::map<unsigned, ChargedParticle*> MCParticle_to_ChargedParticle;
0148
0149   // Create event structure a la standalone pfRICH/IRT code; use MC particles (in_mc_particles)
0150   // in this first iteration (and only select primary ones); later on should do it probably
0151   // the same way as in ATHENA IRT codes, where one had an option to build this event
0152   // structure using reconstructed tracks (in_tracks);
0153   for (const auto& mcparticle : *in_mc_particles) {
0154     // Deal only with charged primary ones, for the time being; FIXME: low momentum cutoff?;
0155     if (mcparticle.isCreatedInSimulation() || !mcparticle.getCharge())
0156       continue;
0157
0158     unsigned mcid = mcparticle.id().index;
0159
0160     // Now check that MC->reco association exists; for now ignore cases where more than one
0161     // reconstructed track is associated with a given MC particle;
0162     const auto rctracks = link_nav.getLinkedFrom(mcparticle);
0163     if (rctracks.empty() || rctracks.size() > 1)
0164       continue;
0165     const auto rctrack_id = rctracks[0].o.getObjectID();
0166     unsigned rctrack      = rctrack_id.index;
0167
0168     // Do not want to deal with particles outside of the nominal acceptance; FIXME: do it better later;
0169     double eta = edm4hep::utils::eta(mcparticle.getMomentum());
0170     if (eta < m_cfg.m_eta_min || eta > m_cfg.m_eta_max)
0171       continue;
0172
0173     const auto seg_it = Track_to_TrackSegment_lut.find(rctrack_id);
0174     if (seg_it == Track_to_TrackSegment_lut.end() || seg_it->second.getPoints().empty()) {
0175       debug("MC particle {} (eta {:.2f}): skipping, track {} has {}", mcid, eta, rctrack,
0176             seg_it == Track_to_TrackSegment_lut.end() ? "no track segment"
0177                                                       : "a track segment without points");
0178       continue;
0179     }
0180     const auto& segment = seg_it->second;
0181
0182     // Now add a charged particle to the event structure; 'true': primary;
0183     auto particle = new ChargedParticle(mcparticle.getPDG(), true);
0184
0185     // For now, just record a reference to EICrecon track; the actual loop with assignments
0186     // will be at the end of processing;
0187     particle->SetEICreconParticleID(rctrack);
0188
0189     // FIXME: check units;
0190     particle->SetVertexPosition(Tools::PodioVector3_to_TVector3(mcparticle.getVertex()));
0191     particle->SetVertexMomentum(Tools::PodioVector3_to_TVector3(mcparticle.getMomentum()));
0192     particle->SetVertexTime(mcparticle.getTime());
0193
0194     m_Event->AddChargedParticle(particle);
0195     MCParticle_to_ChargedParticle[mcid] = particle;
0196
0197     // Create history records for all known radiators; FIXME: may want to optimize a bit;
0198     for (auto [name, rad] : m_irt_detector->Radiators()) {
0199       auto history = new RadiatorHistory();
0200       particle->StartRadiatorHistory(std::make_pair(rad, history));
0201     } //for radiator
0202
0203     for (const auto& point : segment.getPoints()) {
0204       TVector3 position = Tools::PodioVector3_to_TVector3(point.position);
0205       TVector3 momentum = Tools::PodioVector3_to_TVector3(point.momentum);
0206
0207       // FIXME: this call is relatively CPU-intensive; may want to create a lookup
0208       // table, since in principle it is known which projection point corresponds to
0209       // which radiator; however, this way would not be exactly clean because of
0210       // spherical boundaries; leave as it is for now and optimize later;
0211       auto radiator = m_irt_detector->GuessRadiator(position, momentum.Unit());
0212       if (radiator) {
0213         auto history = particle->FindRadiatorHistory(radiator);
0214
0215         // FIXME: this check is redundant?;
0216         if (history) {
0217           auto step = new ChargedParticleStep(position, momentum);
0218           history->AddStep(step);
0219         } //if
0220       } //if
0221     } //for point
0222   } //for mcparticle
0223
0224   // Now loop through simulated hits;
0225   for (auto mchit : *in_sim_hits) {
0226     auto cell_id      = mchit.getCellID();
0227     uint64_t sensorID = cell_id & m_irt_detector->GetReadoutCellMask();
0228
0229     // Get photon which created this hit; filter out charged particles;
0230     auto const& mcparticle = mchit.getParticle();
0231     if (mcparticle.getPDG() != -22)
0232       continue;
0233
0234     // Create an optical photon class instance and populate it; units: [mm], [ns], [eV];
0235     auto photon = new IRT2::OpticalPhoton();
0236
0237     // Information provided by the hit itself: detection position and time;
0238     photon->SetDetectionPosition(Tools::PodioVector3_to_TVector3(mchit.getPosition()));
0239     photon->SetDetectionTime(mchit.getTime());
0240
0241     // Information inherited from photon MCParticle;
0242     photon->SetVertexPosition(Tools::PodioVector3_to_TVector3(mcparticle.getVertex()));
0243     photon->SetVertexMomentum((1 / edm4eic::unit::eV) *
0244                               Tools::PodioVector3_to_TVector3(mcparticle.getMomentum()));
0245     photon->SetVertexTime(mcparticle.getTime());
0246
0247     auto parents = mcparticle.getParents();
0248     if (parents.size() != 1)
0249       continue;
0250     unsigned parent_id = mcparticle.parents_begin()->id().index;
0251     if (MCParticle_to_ChargedParticle.find(parent_id) == MCParticle_to_ChargedParticle.end())
0252       continue;
0253     auto parent = MCParticle_to_ChargedParticle[parent_id];
0254
0255     TVector3 vtx = Tools::PodioVector3_to_TVector3(mcparticle.getVertex());
0256     // FIXME: may want to use the very first projection rather than the IP info?;
0257     auto radiator = m_irt_detector->GuessRadiator(vtx, parent->GetVertexMomentum().Unit());
0258
0259     if (radiator) {
0260       {
0261         double e = photon->GetVertexMomentum().Mag();
0262         double ri =
0263             radiator->m_RefractiveIndex->GetInterpolatedValue(e, DataInterpolation::FirstOrder);
0264         photon->SetVertexRefractiveIndex(ri);
0265
0266         // Will be stored in a (fixed) order in which radiators were defined for this Cherenkov detector;
0267         for (auto [name, rad] : m_irt_detector->Radiators())
0268           photon->StoreRefractiveIndex(
0269               rad->m_RefractiveIndex->GetInterpolatedValue(e, DataInterpolation::FirstOrder));
0270       }
0271
0272       auto history = parent->FindRadiatorHistory(radiator);
0273       // FIXME: this check is redundant?;
0274       if (history) {
0275         history->AddOpticalPhoton(photon);
0276       } //if
0277     } else {
0278       parent->AddOrphanPhoton(photon);
0279     } //if
0280
0281     // FIXME: this is kind of a hack (assume a single type of photodetectors); should be fine
0282     // for ePIC, though a standalone GEANT code has amore generic implementation;
0283     if (m_irt_detector->m_PhotonDetectors.size() != 1)
0284       continue;
0285     auto pd = m_irt_detector->m_PhotonDetectors[0];
0286
0287     if (pd) {
0288       photon->SetPhotonDetector(pd);
0289       // Would be a VolumeCopy in a standalone GEANT code, but is an encoded sensor ID (with
0290       // a blanked out 'sector' field for dRICH) in ePIC;
0291       photon->SetVolumeCopy(sensorID);
0292
0293       if (pd->GetQE() == nullptr) {
0294         static std::once_flag qe_warning_flag;
0295         std::call_once(qe_warning_flag, [this]() {
0296           error("QE is not set - skipping photon selection. Likely, IRT2 was built without "
0297                 "JSON_IMPORT_EXPORT.");
0298         });
0299       } else {
0300         // The logic behind this multiplication and division by the same number is
0301         // to select calibration photons, which originate from the same QE(lambda)
0302         // parent distribution, but do not pass the overall efficiency test;
0303         double e  = photon->GetVertexMomentum().Mag();
0304         double qe = pd->GetQE()->WithinRange(e)
0305                         ? pd->GetQE()->GetInterpolatedValue(e, DataInterpolation::FirstOrder)
0306                         : 0.0;
0307
0308         if (qe * pd->GetScaleFactor() > uniform(generator)) {
0309           if (pd->GetGeometricEfficiency() / pd->GetScaleFactor() > uniform(generator))
0310             photon->SetDetected(true);
0311           else
0312             photon->SetCalibrationFlag();
0313         } //if
0314       }
0315     } //if
0316
0317     // FIXME: should be added? if (!info->Parent()) m_EventPtr->AddOrphanPhoton(photon);
0318   } //for mchit
0319
0320   // FIXME: this is a hack to the moment; also, should one check for
0321   // m_Event->ChargedParticles().size() before mchit loop?;
0322   if (m_ReconstructionFactory && m_Event->ChargedParticles().size())
0323     m_ReconstructionFactory->GetEvent(0, false);
0324
0325   // And eventually, populate PODIO output tables;
0326   {
0327     for (auto particle : m_Event->ChargedParticles()) {
0328       if (!particle->IsPrimary())
0329         continue;
0330
0331       unsigned npe_per_track = 0, nhits_per_track = 0;
0332
0333       edm4eic::MutableIrtParticle irtParticle = out_irt_particles->create();
0334       irtParticle.setTrack((*in_tracks)[particle->m_EICreconParticleID]);
0335
0336       for (auto rhptr : particle->GetRadiatorHistory()) {
0337         auto radiator = particle->GetRadiator(rhptr);
0338         if (!radiator->UsedInRingImaging())
0339           continue;
0340
0341         edm4eic::MutableIrtRadiatorInfo irtRadiator = out_irt_radiator_info->create();
0342         unsigned npe_per_radiator = 0, nhits_per_radiator = 0;
0343
0344         nhits_per_radiator = particle->GetRecoCherenkovPhotonCount(radiator);
0345         irtRadiator.setNhits(nhits_per_radiator);
0346         nhits_per_track += nhits_per_radiator;
0347
0348         for (auto photon : particle->GetHistory(rhptr)->Photons())
0349           if (photon->WasDetected())
0350             npe_per_radiator++;
0351
0352         irtRadiator.setNpe(npe_per_radiator);
0353         npe_per_track += npe_per_radiator;
0354
0355         irtRadiator.setAngle(1000 * particle->GetRecoCherenkovAverageTheta(radiator));
0356
0357         irtParticle.addToRadiators(irtRadiator);
0358       } //for rhistory
0359
0360       irtParticle.setPDG(particle->GetRecoPdgCode());
0361       irtParticle.setNpe(npe_per_track);
0362       irtParticle.setNhits(nhits_per_track);
0363     } //for particle
0364   }
0365 } // IrtInterface::process()
0366 } // namespace eicrecon
0367
0368 #endif