File indexing completed on 2024-09-27 07:02:59
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 <edm4eic/CherenkovParticleIDHypothesis.h>
0014 #include <edm4eic/EDM4eicVersion.h>
0015 #include <edm4eic/TrackPoint.h>
0016 #include <edm4hep/MCParticleCollection.h>
0017 #include <edm4hep/SimTrackerHitCollection.h>
0018 #include <edm4hep/Vector2f.h>
0019 #include <fmt/core.h>
0020 #include <fmt/format.h>
0021 #include <podio/ObjectID.h>
0022 #include <podio/RelationRange.h>
0023 #include <spdlog/common.h>
0024 #include <algorithm>
0025 #include <cstddef>
0026 #include <functional>
0027 #include <gsl/pointers>
0028 #include <iterator>
0029 #include <set>
0030 #include <stdexcept>
0031 #include <utility>
0032 #include <vector>
0033
0034 #include "algorithms/pid/IrtCherenkovParticleIDConfig.h"
0035 #include "algorithms/pid/Tools.h"
0036
0037 namespace eicrecon {
0038
0039 void IrtCherenkovParticleID::init(
0040 CherenkovDetectorCollection* irt_det_coll,
0041 std::shared_ptr<spdlog::logger>& logger
0042 )
0043 {
0044
0045 m_irt_det_coll = irt_det_coll;
0046 m_log = logger;
0047
0048
0049 m_cfg.Print(m_log, spdlog::level::debug);
0050
0051
0052 m_cfg.PrintCheats(m_log);
0053
0054
0055 auto& detectors = m_irt_det_coll->GetDetectors();
0056 if(detectors.size() == 0)
0057 throw std::runtime_error("No CherenkovDetectors found in input collection `irt_det_coll`");
0058 if(detectors.size() > 1)
0059 m_log->warn("IrtCherenkovParticleID currently only supports 1 CherenkovDetector at a time; taking the first");
0060 auto this_detector = *detectors.begin();
0061 m_det_name = this_detector.first;
0062 m_irt_det = this_detector.second;
0063 m_log->debug("Initializing IrtCherenkovParticleID algorithm for CherenkovDetector '{}'", m_det_name);
0064
0065
0066 m_cell_mask = m_irt_det->GetReadoutCellMask();
0067 m_log->debug("readout cellMask = {:#X}", m_cell_mask);
0068
0069
0070 m_log->trace("Rebinning refractive index tables to have {} bins",m_cfg.numRIndexBins);
0071 for(auto [rad_name,irt_rad] : m_irt_det->Radiators()) {
0072 auto ri_lookup_table_orig = irt_rad->m_ri_lookup_table;
0073 irt_rad->m_ri_lookup_table.clear();
0074 irt_rad->m_ri_lookup_table = Tools::ApplyFineBinning( ri_lookup_table_orig, m_cfg.numRIndexBins );
0075
0076
0077 }
0078
0079
0080 m_log->debug("Obtain List of Radiators:");
0081 for(auto [rad_name,irt_rad] : m_irt_det->Radiators()) {
0082 if(rad_name!="Filter") {
0083 m_pid_radiators.insert({ std::string(rad_name), irt_rad });
0084 m_log->debug("- {}", rad_name.Data());
0085 }
0086 }
0087
0088
0089 for(auto [rad_name,irt_rad] : m_pid_radiators) {
0090
0091 auto cfg_rad_it = m_cfg.radiators.find(rad_name);
0092 if(cfg_rad_it != m_cfg.radiators.end()) {
0093 auto cfg_rad = cfg_rad_it->second;
0094
0095 irt_rad->m_ID = Tools::GetRadiatorID(std::string(rad_name));
0096 irt_rad->m_AverageRefractiveIndex = cfg_rad.referenceRIndex;
0097 irt_rad->SetReferenceRefractiveIndex(cfg_rad.referenceRIndex);
0098 if(cfg_rad.attenuation>0)
0099 irt_rad->SetReferenceAttenuationLength(cfg_rad.attenuation);
0100 if(cfg_rad.smearing>0) {
0101 if(cfg_rad.smearingMode=="uniform")
0102 irt_rad->SetUniformSmearing(cfg_rad.smearing);
0103 else if(cfg_rad.smearingMode=="gaussian")
0104 irt_rad->SetGaussianSmearing(cfg_rad.smearing);
0105 else
0106 m_log->error("Unknown smearing mode '{}' for {} radiator", cfg_rad.smearingMode, rad_name);
0107 }
0108 }
0109 else
0110 m_log->error("Cannot find radiator '{}' in IrtCherenkovParticleIDConfig instance", rad_name);
0111 }
0112
0113
0114 m_log->debug("List of particles for PID:");
0115 for(auto pdg : m_cfg.pdgList) {
0116 auto mass = m_particleSvc.particle(pdg).mass;
0117 m_pdg_mass.insert({ pdg, mass });
0118 m_log->debug(" {:>8} M={} GeV", pdg, mass);
0119 }
0120
0121 }
0122
0123 void IrtCherenkovParticleID::process(
0124 const IrtCherenkovParticleID::Input& input,
0125 const IrtCherenkovParticleID::Output& output) const
0126 {
0127 const auto [in_aerogel_tracks, in_gas_tracks, in_merged_tracks, in_raw_hits, in_hit_assocs] = input;
0128 auto [out_aerogel_particleIDs, out_gas_particleIDs] = output;
0129
0130
0131 m_log->trace("{:=^70}"," call IrtCherenkovParticleID::AlgorithmProcess ");
0132 m_log->trace("number of raw sensor hits: {}", in_raw_hits->size());
0133 m_log->trace("number of raw sensor hit with associated photons: {}", in_hit_assocs->size());
0134
0135 std::map<std::string, const edm4eic::TrackSegmentCollection*> in_charged_particles{
0136 {"Aerogel", in_aerogel_tracks},
0137 {"Gas", in_gas_tracks},
0138 {"Merged", in_merged_tracks},
0139 };
0140
0141
0142 std::map<std::string, edm4eic::CherenkovParticleIDCollection*> out_cherenkov_pids{
0143 {"Aerogel", out_aerogel_particleIDs},
0144 {"Gas", out_gas_particleIDs}
0145 };
0146
0147
0148 std::unordered_map<std::size_t, std::size_t> in_charged_particle_size_distribution;
0149 for(const auto& [rad_name, in_charged_particle] : in_charged_particles) {
0150 ++in_charged_particle_size_distribution[in_charged_particle->size()];
0151 }
0152 if (in_charged_particle_size_distribution.size() != 1) {
0153 std::vector<size_t> in_charged_particle_sizes;
0154 std::transform(in_charged_particles.begin(), in_charged_particles.end(),
0155 std::back_inserter(in_charged_particle_sizes),
0156 [](const auto& in_charged_particle) { return in_charged_particle.second->size(); });
0157 m_log->error("radiators have differing numbers of TrackSegments {}", fmt::join(in_charged_particle_sizes, ", "));
0158 return;
0159 }
0160
0161
0162 m_log->trace("{:#<70}","### CHARGED PARTICLES ");
0163 std::size_t num_charged_particles = in_charged_particle_size_distribution.begin()->first;
0164 for(long i_charged_particle=0; i_charged_particle<num_charged_particles; i_charged_particle++) {
0165 m_log->trace("{:-<70}", fmt::format("--- charged particle #{} ", i_charged_particle));
0166
0167
0168 auto irt_particle = std::make_unique<ChargedParticle>();
0169
0170
0171 for(auto [rad_name,irt_rad] : m_pid_radiators) {
0172
0173
0174 auto charged_particle_list_it = in_charged_particles.find(rad_name);
0175 if(charged_particle_list_it == in_charged_particles.end()) {
0176 m_log->error("Cannot find radiator '{}' in `in_charged_particles`", rad_name);
0177 continue;
0178 }
0179 const auto *charged_particle_list = charged_particle_list_it->second;
0180 auto charged_particle = charged_particle_list->at(i_charged_particle);
0181
0182
0183 if(charged_particle.points_size()==0) {
0184 m_log->trace("No propagated track points in radiator '{}'", rad_name);
0185 continue;
0186 }
0187 irt_rad->SetTrajectoryBinCount(charged_particle.points_size() - 1);
0188
0189
0190
0191
0192 auto *irt_rad_history = new RadiatorHistory();
0193 irt_particle->StartRadiatorHistory({ irt_rad, irt_rad_history });
0194
0195
0196 irt_rad->ResetLocations();
0197 m_log->trace("TrackPoints in '{}' radiator:", rad_name);
0198 for(const auto& point : charged_particle.getPoints()) {
0199 TVector3 position = Tools::PodioVector3_to_TVector3(point.position);
0200 TVector3 momentum = Tools::PodioVector3_to_TVector3(point.momentum);
0201 irt_rad->AddLocation(position, momentum);
0202 Tools::PrintTVector3(m_log, " point: x", position);
0203 Tools::PrintTVector3(m_log, " p", momentum);
0204 }
0205
0206
0207
0208 m_log->trace("{:#<70}","### SENSOR HITS ");
0209 for(const auto& raw_hit : *in_raw_hits) {
0210
0211
0212
0213
0214 edm4hep::MCParticle mc_photon;
0215 bool mc_photon_found = false;
0216 if(m_cfg.cheatPhotonVertex || m_cfg.cheatTrueRadiator) {
0217 for(const auto& hit_assoc : *in_hit_assocs) {
0218 if(hit_assoc.getRawHit().isAvailable()) {
0219 if(hit_assoc.getRawHit().id() == raw_hit.id()) {
0220 #if EDM4EIC_VERSION_MAJOR >= 6
0221 mc_photon = hit_assoc.getSimHit().getMCParticle();
0222 #else
0223
0224 if(hit_assoc.simHits_size() > 0) {
0225 mc_photon = hit_assoc.getSimHits(0).getMCParticle();
0226 #endif
0227 mc_photon_found = true;
0228 if(mc_photon.getPDG() != -22)
0229 m_log->warn("non-opticalphoton hit: PDG = {}",mc_photon.getPDG());
0230 #if EDM4EIC_VERSION_MAJOR >= 6
0231 #else
0232 }
0233 else if(m_cfg.CheatModeEnabled())
0234 m_log->error("cheat mode enabled, but no MC photons provided");
0235 #endif
0236 break;
0237 }
0238 }
0239 }
0240 }
0241
0242
0243 if(m_cfg.cheatTrueRadiator && mc_photon_found) {
0244 auto vtx = Tools::PodioVector3_to_TVector3(mc_photon.getVertex());
0245 auto mc_rad = m_irt_det->GuessRadiator(vtx, vtx);
0246 if(mc_rad != irt_rad) continue;
0247 Tools::PrintTVector3(m_log, fmt::format("cheat: radiator '{}' determined from photon vertex", rad_name), vtx);
0248 }
0249
0250
0251
0252 auto cell_id = raw_hit.getCellID();
0253 uint64_t sensor_id = cell_id & m_cell_mask;
0254 TVector3 pixel_pos = m_irt_det->m_ReadoutIDToPosition(cell_id);
0255
0256
0257 if(m_log->level() <= spdlog::level::trace) {
0258 m_log->trace("cell_id={:#X} sensor_id={:#X}", cell_id, sensor_id);
0259 Tools::PrintTVector3(m_log, "pixel position", pixel_pos);
0260 if(mc_photon_found) {
0261 TVector3 mc_endpoint = Tools::PodioVector3_to_TVector3(mc_photon.getEndpoint());
0262 Tools::PrintTVector3(m_log, "photon endpoint", mc_endpoint);
0263 m_log->trace("{:>30} = {}", "dist( pixel, photon )", (pixel_pos - mc_endpoint).Mag());
0264 }
0265 else m_log->trace(" no MC photon found; probably a noise hit");
0266 }
0267
0268
0269 auto *irt_sensor = m_irt_det->m_PhotonDetectors[0];
0270 auto *irt_photon = new OpticalPhoton();
0271 irt_photon->SetVolumeCopy(sensor_id);
0272 irt_photon->SetDetectionPosition(pixel_pos);
0273 irt_photon->SetPhotonDetector(irt_sensor);
0274 irt_photon->SetDetected(true);
0275
0276
0277 if((m_cfg.cheatPhotonVertex || m_cfg.cheatTrueRadiator) && mc_photon_found) {
0278 irt_photon->SetVertexPosition(Tools::PodioVector3_to_TVector3(mc_photon.getVertex()));
0279 irt_photon->SetVertexMomentum(Tools::PodioVector3_to_TVector3(mc_photon.getMomentum()));
0280 }
0281
0282
0283
0284 if(m_cfg.cheatPhotonVertex) {
0285 double ri;
0286 auto mom = 1e9 * irt_photon->GetVertexMomentum().Mag();
0287 auto ri_set = Tools::GetFinelyBinnedTableEntry(irt_rad->m_ri_lookup_table, mom, &ri);
0288 if(ri_set) {
0289 irt_photon->SetVertexRefractiveIndex(ri);
0290 m_log->trace("{:>30} = {}", "refractive index", ri);
0291 }
0292 else
0293 m_log->warn("Tools::GetFinelyBinnedTableEntry failed to lookup refractive index for momentum {} eV", mom);
0294 }
0295
0296
0297
0298
0299
0300 irt_rad_history->AddOpticalPhoton(irt_photon);
0301
0302
0303
0304
0305
0306 }
0307
0308 }
0309
0310
0311
0312
0313
0314
0315 m_log->trace("{:+^70}"," PARTICLE IDENTIFICATION ");
0316 CherenkovPID irt_pid;
0317 std::unordered_map<int,MassHypothesis*> pdg_to_hyp;
0318 for(auto [pdg,mass] : m_pdg_mass) {
0319 irt_pid.AddMassHypothesis(mass);
0320 pdg_to_hyp.insert({ pdg, irt_pid.GetHypothesis(irt_pid.GetHypothesesCount()-1) });
0321 }
0322
0323
0324 irt_particle->PIDReconstruction(irt_pid);
0325 m_log->trace("{:-^70}"," IRT RESULTS ");
0326
0327
0328 for(auto [rad_name,irt_rad] : m_pid_radiators) {
0329 m_log->trace("-> {} Radiator (ID={}):", rad_name, irt_rad->m_ID);
0330
0331
0332 unsigned npe = 0;
0333 double rindex_ave = 0.0;
0334 double energy_ave = 0.0;
0335 std::vector<std::pair<double,double>> phot_theta_phi;
0336
0337
0338 auto *irt_rad_history = irt_particle->FindRadiatorHistory(irt_rad);
0339 if(irt_rad_history==nullptr) {
0340 m_log->trace(" No radiator history; skip");
0341 continue;
0342 }
0343 m_log->trace(" Photoelectrons:");
0344 for(auto *irt_photon : irt_rad_history->Photons()) {
0345
0346
0347 bool photon_selected = false;
0348 for(auto irt_photon_sel : irt_photon->_m_Selected) {
0349 if(irt_photon_sel.second == irt_rad) {
0350 photon_selected = true;
0351 break;
0352 }
0353 }
0354 if(!photon_selected) continue;
0355
0356
0357 Tools::PrintTVector3(
0358 m_log,
0359 fmt::format("- sensor_id={:#X}: hit",irt_photon->GetVolumeCopy()),
0360 irt_photon->GetDetectionPosition()
0361 );
0362 Tools::PrintTVector3(m_log, "photon vertex", irt_photon->GetVertexPosition());
0363
0364
0365 auto phot_theta = irt_photon->_m_PDF[irt_rad].GetAverage();
0366 auto phot_phi = irt_photon->m_Phi[irt_rad];
0367
0368
0369 npe++;
0370 phot_theta_phi.emplace_back( phot_theta, phot_phi );
0371 if(m_cfg.cheatPhotonVertex) {
0372 rindex_ave += irt_photon->GetVertexRefractiveIndex();
0373 energy_ave += irt_photon->GetVertexMomentum().Mag();
0374 }
0375
0376 }
0377
0378
0379 if(npe>0) {
0380 rindex_ave /= npe;
0381 energy_ave /= npe;
0382 }
0383
0384
0385 auto out_cherenkov_pid = out_cherenkov_pids.at(rad_name)->create();
0386 out_cherenkov_pid.setNpe(static_cast<decltype(edm4eic::CherenkovParticleIDData::npe)>(npe));
0387 out_cherenkov_pid.setRefractiveIndex(static_cast<decltype(edm4eic::CherenkovParticleIDData::refractiveIndex)>(rindex_ave));
0388 out_cherenkov_pid.setPhotonEnergy(static_cast<decltype(edm4eic::CherenkovParticleIDData::photonEnergy)>(energy_ave));
0389 for(auto [phot_theta,phot_phi] : phot_theta_phi)
0390 out_cherenkov_pid.addToThetaPhiPhotons(edm4hep::Vector2f{
0391 static_cast<float>(phot_theta),
0392 static_cast<float>(phot_phi)
0393 });
0394
0395
0396 for(auto [pdg,mass] : m_pdg_mass) {
0397
0398
0399 auto *irt_hypothesis = pdg_to_hyp.at(pdg);
0400 auto hyp_weight = irt_hypothesis->GetWeight(irt_rad);
0401 auto hyp_npe = irt_hypothesis->GetNpe(irt_rad);
0402
0403
0404 edm4eic::CherenkovParticleIDHypothesis out_hypothesis;
0405 out_hypothesis.PDG = static_cast<decltype(edm4eic::CherenkovParticleIDHypothesis::PDG)>(pdg);
0406 out_hypothesis.weight = static_cast<decltype(edm4eic::CherenkovParticleIDHypothesis::weight)>(hyp_weight);
0407 out_hypothesis.npe = static_cast<decltype(edm4eic::CherenkovParticleIDHypothesis::npe)>(hyp_npe);
0408
0409
0410 out_cherenkov_pid.addToHypotheses(out_hypothesis);
0411
0412 }
0413
0414
0415 Tools::PrintCherenkovEstimate(m_log, out_cherenkov_pid);
0416
0417
0418 auto charged_particle_list_it = in_charged_particles.find("Merged");
0419 if(charged_particle_list_it != in_charged_particles.end()) {
0420 const auto *charged_particle_list = charged_particle_list_it->second;
0421 auto charged_particle = charged_particle_list->at(i_charged_particle);
0422 out_cherenkov_pid.setChargedParticle(charged_particle);
0423 }
0424 else
0425 m_log->error("Cannot find radiator 'Merged' in `in_charged_particles`");
0426
0427
0428 for(const auto& hit_assoc : *in_hit_assocs)
0429 out_cherenkov_pid.addToRawHitAssociations(hit_assoc);
0430
0431 }
0432
0433
0434
0435
0436
0437
0438 }
0439 }
0440
0441 }