rich/src/ReconstructedParticleAnalysis.cc

0001 // Copyright 2023, Christopher Dilks
0002 // Subject to the terms in the LICENSE file found in the top-level directory.
0003
0004 #include "ReconstructedParticleAnalysis.h"
0005
0006 namespace benchmarks {
0007
0008   // AlgorithmInit
0009   //---------------------------------------------------------------------------
0010   void ReconstructedParticleAnalysis::AlgorithmInit(std::shared_ptr<spdlog::logger>& logger) {
0011     m_log = logger;
0012
0013     // initialize plots
0014     m_correct_vs_p_transient = new TH2D("correct_vs_p_transient", "Fraction with Correct PDG vs. Thrown Momentum;p [GeV];fraction",
0015         Tools::momentum_bins, 0, Tools::momentum_max,
0016         2,                    0, 2
0017         );
0018     m_correct_vs_p = new TGraphErrors();
0019     m_correct_vs_p->SetName("correct_vs_p");
0020     m_correct_vs_p->SetTitle(m_correct_vs_p_transient->GetTitle());
0021     m_correct_vs_p->GetXaxis()->SetTitle(m_correct_vs_p_transient->GetXaxis()->GetTitle());
0022     m_correct_vs_p->GetYaxis()->SetTitle(m_correct_vs_p_transient->GetYaxis()->GetTitle());
0023     m_correct_vs_p->SetMarkerStyle(kFullCircle);
0024   }
0025
0026
0027   // AlgorithmProcess
0028   //---------------------------------------------------------------------------
0029   void ReconstructedParticleAnalysis::AlgorithmProcess(
0030       const edm4eic::MCRecoParticleAssociationCollection& in_assocs
0031       )
0032   {
0033     if(!in_assocs.isValid()) { m_log->error("invalid input collection 'in_assocs'"); return; }
0034
0035     // loop over input associations
0036     for(const auto& assoc : in_assocs) {
0037
0038       // get particles
0039       m_log->trace("{:-^50}"," Particle ");
0040       auto& simpart = assoc.getSim();
0041       auto& recpart = assoc.getRec();
0042       if(!recpart.isAvailable()) { m_log->warn("reconstructed particle not available for this association"); continue; }
0043       if(!simpart.isAvailable()) { m_log->warn("simulated particle not available for this association");     continue; }
0044
0045       // get PDG values
0046       auto simpart_pdg      = simpart.getPDG(); // from MC truth
0047       auto recpart_pdg      = recpart.getPDG(); // PDG member of ReconstructedParticle
0048       auto recpart_pid_used = recpart.getParticleIDUsed();
0049       auto pid_pdg          = recpart_pid_used.isAvailable() ? recpart_pid_used.getPDG() : 0;
0050       if(!recpart_pid_used.isAvailable())
0051         m_log->error("reconstructed particle does not have particleIDUsed relation");
0052       m_log->trace("PDGs: sim.PDG | rec.PDG | PDG from PID  =  {:^6} | {:^6} | {:^6}", simpart_pdg, recpart_pdg, pid_pdg);
0053
0054       // check goodness of PID: if near zero, assume PID was not used for this particle and skip
0055       if(recpart.getGoodnessOfPID()<0.01) {
0056         m_log->warn("skip particle with low goodnessOfPID (likely has no PID)");
0057         continue;
0058       }
0059
0060       // get momenta
0061       auto simpart_p = edm4hep::utils::p(simpart);
0062       // auto recpart_p = edm4hep::utils::p(recpart); // unused
0063
0064       // fill plots
0065       m_correct_vs_p_transient->Fill(
0066           simpart_p,
0067           simpart_pdg == pid_pdg ? 1 : 0
0068           );
0069
0070     } // end loop over input associations
0071
0072   }
0073
0074
0075   // AlgorithmFinish
0076   //---------------------------------------------------------------------------
0077   void ReconstructedParticleAnalysis::AlgorithmFinish() {
0078
0079     // compute fraction of correct PDG for each x-axis bin
0080     for(int bx=1; bx<=m_correct_vs_p_transient->GetNbinsX(); bx++) {
0081       auto n_wrong  = m_correct_vs_p_transient->GetBinContent(bx,1);
0082       auto n_right  = m_correct_vs_p_transient->GetBinContent(bx,2);
0083       auto n_total  = n_wrong + n_right;
0084       auto momentum = m_correct_vs_p_transient->GetXaxis()->GetBinCenter(bx);
0085       if(n_total>0) {
0086         auto frac = n_right / n_total;
0087         m_correct_vs_p->AddPoint(momentum, frac);
0088       }
0089     }
0090
0091     // delete transient histograms, so they don't get written
0092     delete m_correct_vs_p_transient;
0093
0094     // write objects which are not automatically saved
0095     m_correct_vs_p->Write();
0096   }
0097
0098 }