File indexing completed on 2025-07-03 07:55:53
0001
0002
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
0045 m_irt_det_coll = irt_det_coll;
0046
0047
0048 debug() << m_cfg << endmsg;
0049
0050
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
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
0075 m_cell_mask = m_irt_det->GetReadoutCellMask();
0076 debug("readout cellMask = {:#X}", m_cell_mask);
0077
0078
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
0085
0086 }
0087
0088
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
0098 for (auto [rad_name, irt_rad] : m_pid_radiators) {
0099
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
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
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
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
0151 std::map<std::string, edm4eic::CherenkovParticleIDCollection*> out_cherenkov_pids{
0152 {"Aerogel", out_aerogel_particleIDs}, {"Gas", out_gas_particleIDs}};
0153
0154
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
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
0178 auto irt_particle = std::make_unique<ChargedParticle>();
0179
0180
0181
0182
0183 std::lock_guard<std::mutex> lock(m_pid_radiators_mutex);
0184 for (auto [rad_name, irt_rad] : m_pid_radiators) {
0185
0186
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
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
0203
0204
0205 auto* irt_rad_history = new RadiatorHistory();
0206 irt_particle->StartRadiatorHistory({irt_rad, irt_rad_history});
0207
0208
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
0220 trace("{:#<70}", "### SENSOR HITS ");
0221 for (const auto& raw_hit : *in_raw_hits) {
0222
0223
0224
0225
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
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
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);
0262 if (mc_rad != irt_rad) {
0263 continue;
0264 }
0265 trace(Tools::PrintTVector3(
0266 fmt::format("cheat: radiator '{}' determined from photon vertex", rad_name), vtx));
0267 }
0268
0269
0270
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
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
0289 auto* irt_sensor = m_irt_det->m_PhotonDetectors[0];
0290 auto* irt_photon =
0291 new OpticalPhoton();
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
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
0304
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
0320
0321
0322
0323 irt_rad_history->AddOpticalPhoton(irt_photon);
0324
0325
0326
0327
0328
0329 }
0330
0331 }
0332
0333
0334
0335
0336 trace("{:+^70}", " PARTICLE IDENTIFICATION ");
0337 CherenkovPID irt_pid;
0338 std::unordered_map<int, MassHypothesis*> pdg_to_hyp;
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
0345 irt_particle->PIDReconstruction(irt_pid);
0346 trace("{:-^70}", " IRT RESULTS ");
0347
0348
0349 for (auto [rad_name, irt_rad] : m_pid_radiators) {
0350 trace("-> {} Radiator (ID={}):", rad_name, irt_rad->m_ID);
0351
0352
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
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
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
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
0386 auto phot_theta = irt_photon->_m_PDF[irt_rad].GetAverage();
0387 auto phot_phi = irt_photon->m_Phi[irt_rad];
0388
0389
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 }
0398
0399
0400 if (npe > 0) {
0401 rindex_ave /= npe;
0402 energy_ave /= npe;
0403 }
0404
0405
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
0418 for (auto [pdg, mass] : m_pdg_mass) {
0419
0420
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
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
0435 out_cherenkov_pid.addToHypotheses(out_hypothesis);
0436
0437 }
0438
0439
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);
0452 trace("{:{}} {:>16}: {:>10.8}", "", indent, "<rindex>", pid.getRefractiveIndex());
0453 trace("{:{}} {:>16}: {:>10.8} eV", "", indent, "<energy>",
0454 pid.getPhotonEnergy() * 1e9);
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
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
0476 for (const auto& hit_assoc : *in_hit_assocs) {
0477 out_cherenkov_pid.addToRawHitAssociations(hit_assoc);
0478 }
0479
0480 }
0481
0482
0483
0484
0485
0486
0487 }
0488 }
0489
0490 }