rich/src/CherenkovPIDAnalysis.cc

0001 // Copyright 2023, Christopher Dilks
0002 // Subject to the terms in the LICENSE file found in the top-level directory.
0003
0004 #include "CherenkovPIDAnalysis.h"
0005
0006 namespace benchmarks {
0007
0008   // AlgorithmInit
0009   //---------------------------------------------------------------------------
0010   void CherenkovPIDAnalysis::AlgorithmInit(
0011       std::string                      rad_name,
0012       std::shared_ptr<spdlog::logger>& logger
0013       )
0014   {
0015     m_rad_name = rad_name;
0016     m_log      = logger;
0017
0018     // set radiator title
0019     m_rad_title = TString(m_rad_name);
0020     if(m_rad_name=="Merged") m_rad_title += " Aerogel+Gas";
0021
0022     // truncate `m_rad_name` (for object names)
0023     m_rad_name_trun = TString(m_rad_name);
0024     m_rad_name_trun(TRegexp(" .*")) = "";
0025
0026     // distributions
0027     m_npe_dist = new TH1D(
0028         "npe_dist_"+m_rad_name_trun,
0029         "Overall NPE for "+m_rad_title+";NPE",
0030         Tools::npe_max, 0, Tools::npe_max
0031         );
0032     m_theta_dist = new TH1D(
0033         "theta_dist_"+m_rad_name_trun,
0034         "Estimated Cherenkov Angle for "+m_rad_title+";#theta [mrad]",
0035         Tools::theta_bins, 0, Tools::theta_max
0036         );
0037     m_thetaResid_dist = new TH1D(
0038         "thetaResid_dist_"+m_rad_name_trun,
0039         "Estimated Cherenkov Angle Residual for "+m_rad_title+";#Delta#theta [mrad]",
0040         Tools::theta_bins, -Tools::thetaResid_max, Tools::thetaResid_max
0041         );
0042     m_photonTheta_vs_photonPhi = new TH2D(
0043         "photonTheta_vs_photonPhi_"+m_rad_name_trun,
0044         "Estimated Photon #theta vs #phi for "+m_rad_title+";#phi [rad];#theta [mrad]",
0045         Tools::phi_bins, -TMath::Pi(), TMath::Pi(),
0046         Tools::theta_bins, 0, Tools::theta_max
0047         );
0048
0049     // truth distributions
0050     m_mcWavelength_dist = new TH1D(
0051         "mcWavelength_"+m_rad_name_trun,
0052         "MC Photon Wavelength for "+m_rad_title+";#lambda [nm]",
0053         Tools::n_bins, 0, 1000
0054         );
0055     m_mcRindex_dist = new TH1D(
0056         "mcRindex_"+m_rad_name_trun,
0057         "MC Refractive Index for "+m_rad_title+";n",
0058         10*Tools::n_bins, 0.99, 1.03
0059         );
0060
0061     // PID
0062     m_highestWeight_dist = new TH1D(
0063         "highestWeight_dist_"+m_rad_name_trun,
0064         "Highest PDG Weight for "+m_rad_title+";PDG",
0065         Tools::GetNumPDGs()+1, 0, Tools::GetNumPDGs()+1
0066         );
0067
0068     // momentum scans
0069     m_npe_vs_p = new TH2D(
0070         "npe_vs_p_"+m_rad_name_trun,
0071         "Overall NPE vs. Particle Momentum for "+m_rad_title+";p [GeV];NPE",
0072         Tools::momentum_bins, 0, Tools::momentum_max,
0073         Tools::npe_max, 0, Tools::npe_max
0074         );
0075     m_theta_vs_p = new TH2D(
0076         "theta_vs_p_"+m_rad_name_trun,
0077         "Estimated Cherenkov Angle vs. Particle Momentum for "+m_rad_title+";p [GeV];#theta [mrad]",
0078         Tools::momentum_bins, 0, Tools::momentum_max,
0079         Tools::theta_bins, 0, Tools::theta_max
0080         );
0081     m_thetaResid_vs_p = new TH2D(
0082         "thetaResid_vs_p_"+m_rad_name_trun,
0083         "Estimated Cherenkov Angle Residual vs. Particle Momentum for "+m_rad_title+";p [GeV];#Delta#theta [mrad]",
0084         Tools::momentum_bins, 0, Tools::momentum_max,
0085         Tools::theta_bins, -Tools::thetaResid_max, Tools::thetaResid_max
0086         );
0087     m_highestWeight_vs_p = new TH2D(
0088         "highestWeight_vs_p_"+m_rad_name_trun,
0089         "Highest PDG Weight vs. Particle Momentum for "+m_rad_title+";p [GeV]",
0090         Tools::momentum_bins, 0, Tools::momentum_max,
0091         Tools::GetNumPDGs()+1, 0, Tools::GetNumPDGs()+1
0092         );
0093
0094     // pseudorapidity scans
0095     m_npe_vs_eta = new TH2D(
0096         "npe_vs_eta_"+m_rad_name_trun,
0097         "Overall NPE vs. Pseudorapidity for "+m_rad_title+";#eta;NPE",
0098         Tools::eta_bins, Tools::eta_min, Tools::eta_max,
0099         Tools::npe_max, 0, Tools::npe_max
0100         );
0101     m_theta_vs_eta = new TH2D(
0102         "theta_vs_eta_"+m_rad_name_trun,
0103         "Estimated Cherenkov Angle vs. Pseudorapidity for "+m_rad_title+";#eta;#theta [mrad]",
0104         Tools::eta_bins, Tools::eta_min, Tools::eta_max,
0105         Tools::theta_bins, 0, Tools::theta_max
0106         );
0107     m_thetaResid_vs_eta = new TH2D(
0108         "thetaResid_vs_eta_"+m_rad_name_trun,
0109         "Estimated Cherenkov Angle Residual vs. Pseudorapidity for "+m_rad_title+";#eta;#Delta#theta [mrad]",
0110         Tools::eta_bins, Tools::eta_min, Tools::eta_max,
0111         Tools::theta_bins, -Tools::thetaResid_max, Tools::thetaResid_max
0112         );
0113     m_highestWeight_vs_eta = new TH2D(
0114         "highestWeight_vs_eta_"+m_rad_name_trun,
0115         "Highest PDG Weight vs. Pseudorapidity for "+m_rad_title+";#eta",
0116         Tools::eta_bins, Tools::eta_min, Tools::eta_max,
0117         Tools::GetNumPDGs()+1, 0, Tools::GetNumPDGs()+1
0118         );
0119
0120     // format histograms
0121     auto format1D = [] (auto h) {
0122       h->SetLineColor(kAzure);
0123       h->SetFillColor(kAzure);
0124     };
0125     format1D(m_npe_dist);
0126     format1D(m_theta_dist);
0127     format1D(m_thetaResid_dist);
0128     format1D(m_mcWavelength_dist);
0129     format1D(m_mcRindex_dist);
0130     format1D(m_highestWeight_dist);
0131   }
0132
0133
0134   // AlgorithmProcess
0135   //---------------------------------------------------------------------------
0136   void CherenkovPIDAnalysis::AlgorithmProcess(
0137       const edm4hep::MCParticleCollection&          mc_parts,
0138       const edm4eic::CherenkovParticleIDCollection& cherenkov_pids
0139       )
0140   {
0141     m_log->trace("-> {} Radiator:", m_rad_name);
0142
0143     // get thrown momentum from a single MCParticle
0144     auto part = std::make_unique<ChargedParticle>(m_log);
0145     part->SetSingleParticle(mc_parts);
0146     auto thrown_momentum = part->GetMomentum();
0147     auto thrown_eta      = part->GetEta();
0148     auto thrown_pdg      = part->GetPDG();
0149
0150     // loop over `CherenkovParticleID` objects
0151     for(const auto& cherenkov_pid : cherenkov_pids) {
0152
0153       // skip if NPE==0
0154       if(cherenkov_pid.getNpe() == 0) {
0155         m_log->warn("Event found with NPE=0");
0156         continue;
0157       }
0158
0159       // estimate the charged particle momentum using the momentum of the first TrackPoint at this radiator's entrance
0160       /*
0161       auto charged_particle = cherenkov_pid.getChargedParticle();
0162       if(!charged_particle.isAvailable())   { m_log->warn("Charged particle not available in this radiator");      continue; }
0163       if(charged_particle.points_size()==0) { m_log->warn("Charged particle has no TrackPoints in this radiator"); continue; }
0164       auto charged_particle_momentum = edm4hep::utils::magnitude( charged_particle.getPoints(0).momentum );
0165       m_log->trace("  Charged Particle p = {} GeV at radiator entrance", charged_particle_momentum);
0166       m_log->trace("  If it is a pion, E = {} GeV", std::hypot(charged_particle_momentum, Tools::GetPDGMass(211)));
0167       */
0168       // alternatively: use `thrown_momentum` (FIXME: will not work for multi-track events)
0169       auto charged_particle_momentum = thrown_momentum;
0170       auto charged_particle_eta      = thrown_eta;
0171       auto charged_particle_pdg      = thrown_pdg;
0172       auto charged_particle_mass     = Tools::GetPDGMass(charged_particle_pdg);
0173       m_log->trace("  Charged Particle PDG={}, mass={} GeV, p={} GeV from truth",
0174           charged_particle_pdg, charged_particle_mass, charged_particle_momentum);
0175
0176       // get average Cherenkov angle: `theta_rec`
0177       double theta_rec = 0.0;
0178       for(const auto& [theta,phi] : cherenkov_pid.getThetaPhiPhotons())
0179         theta_rec += theta;
0180       theta_rec /= cherenkov_pid.getNpe();
0181       auto theta_rec_mrad = theta_rec * 1e3; // [rad] -> [mrad]
0182
0183       // calculate expected Cherenkov angle `theta_exp` and residual `theta_resid`,
0184       // using refractive index from MC truth
0185       auto mc_rindex = cherenkov_pid.getRefractiveIndex(); // average refractive index for photons used in this `cherenkov_pid`
0186       auto theta_exp = TMath::ACos(
0187           TMath::Hypot( charged_particle_momentum, charged_particle_mass ) /
0188           ( mc_rindex * charged_particle_momentum )
0189           );
0190       auto theta_resid = theta_rec - theta_exp;
0191       auto theta_resid_mrad = theta_resid * 1e3; // [rad] -> [mrad]
0192
0193       // fill PID histograms
0194       m_npe_dist->Fill(   cherenkov_pid.getNpe());
0195       m_npe_vs_p->Fill(   charged_particle_momentum, cherenkov_pid.getNpe());
0196       m_npe_vs_eta->Fill( charged_particle_eta,      cherenkov_pid.getNpe());
0197       m_theta_dist->Fill(   theta_rec_mrad);
0198       m_theta_vs_p->Fill(   charged_particle_momentum, theta_rec_mrad);
0199       m_theta_vs_eta->Fill( charged_particle_eta,      theta_rec_mrad);
0200       m_thetaResid_dist->Fill(   theta_resid_mrad);
0201       m_thetaResid_vs_p->Fill(   charged_particle_momentum, theta_resid_mrad);
0202       m_thetaResid_vs_eta->Fill( charged_particle_eta,      theta_resid_mrad);
0203       for(const auto& [theta,phi] : cherenkov_pid.getThetaPhiPhotons())
0204         m_photonTheta_vs_photonPhi->Fill( phi, theta*1e3); // [rad] -> [mrad]
0205       m_mcWavelength_dist->Fill(  Tools::HC / cherenkov_pid.getPhotonEnergy() ); // energy [GeV] -> wavelength [nm]
0206       m_mcRindex_dist->Fill( mc_rindex);
0207
0208       // find the PDG hypothesis with the highest weight
0209       float max_weight     = -1000;
0210       int   pdg_max_weight = 0;
0211       for(const auto& hyp : cherenkov_pid.getHypotheses()) {
0212         if(hyp.weight > max_weight) {
0213           max_weight     = hyp.weight;
0214           pdg_max_weight = hyp.PDG;
0215         }
0216       }
0217       std::string pdg_max_weight_str = "UNKNOWN";
0218       if(pdg_max_weight!=0 && !std::isnan(pdg_max_weight))
0219         pdg_max_weight_str = std::to_string(pdg_max_weight);
0220       m_log->trace("  Highest weight is {} for PDG {} (string='{}')", max_weight, pdg_max_weight, pdg_max_weight_str);
0221       m_highestWeight_dist->Fill(   pdg_max_weight_str.c_str(), 1);
0222       m_highestWeight_vs_p->Fill(   charged_particle_momentum,  pdg_max_weight_str.c_str(), 1);
0223       m_highestWeight_vs_eta->Fill( charged_particle_eta,       pdg_max_weight_str.c_str(), 1);
0224
0225     } // end loop over `CherenkovParticleID` objects
0226   }
0227
0228
0229   // AlgorithmFinish
0230   //---------------------------------------------------------------------------
0231   void CherenkovPIDAnalysis::AlgorithmFinish() {
0232   }
0233
0234 }