rich/src/SimHitAnalysis.cc

0001 // Copyright 2023, Christopher Dilks
0002 // Subject to the terms in the LICENSE file found in the top-level directory.
0003
0004 #include "SimHitAnalysis.h"
0005
0006 namespace benchmarks {
0007
0008   // AlgorithmInit
0009   //---------------------------------------------------------------------------
0010   void SimHitAnalysis::AlgorithmInit(
0011       std::shared_ptr<spdlog::logger>& logger
0012       )
0013   {
0014     m_log = logger;
0015
0016     // initialize histograms
0017     m_nphot_dist = new TH1D("nphot_dist", "Number of incident photons;N_{photons}",
0018         Tools::nphot_max, 0, Tools::nphot_max
0019         );
0020     m_nphot_vs_p = new TH2D("nphot_vs_p", "Number of incident photons vs. Thrown Momentum;p [GeV];N_{photons}",
0021         Tools::momentum_bins, 0, Tools::momentum_max,
0022         Tools::nphot_max,     0, Tools::nphot_max
0023         );
0024     m_nphot_vs_p__transient = new TH1D("nphot_vs_p__transient", "",
0025         m_nphot_vs_p->GetNbinsX(),
0026         m_nphot_vs_p->GetXaxis()->GetXmin(),
0027         m_nphot_vs_p->GetXaxis()->GetXmax()
0028         );
0029     m_phot_spectrum = new TH1D(
0030         "phot_spectrum_sim",
0031         "Incident photon wavelength;#lambda [nm], at emission point",
0032         Tools::wl_bins, Tools::wl_min, Tools::wl_max
0033         );
0034     m_phot_spectrum->SetLineColor(kViolet+2);
0035     m_phot_spectrum->SetFillColor(kViolet+2);
0036     m_nphot_dist->SetLineColor(kCyan+2);
0037     m_nphot_dist->SetFillColor(kCyan+2);
0038   }
0039
0040
0041   // AlgorithmProcess
0042   //---------------------------------------------------------------------------
0043   void SimHitAnalysis::AlgorithmProcess(
0044       const edm4hep::MCParticleCollection&    mc_parts,
0045       const edm4hep::SimTrackerHitCollection& sim_hits
0046       )
0047   {
0048
0049     // get thrown momentum from a single MCParticle
0050     auto part = std::make_unique<ChargedParticle>(m_log);
0051     part->SetSingleParticle(mc_parts);
0052     auto thrown_momentum = part->GetMomentum();
0053
0054     // fill nphot distribution
0055     m_nphot_dist->Fill(sim_hits.size());
0056
0057     // loop over hits
0058     for(const auto& sim_hit : sim_hits) {
0059
0060       // - fill 1D thrown momentum histogram `m_nphot_vs_p__transient` for each (pre-digitized) sensor hit
0061       // FIXME: get thrown momentum from multi-track events; but in this attempt
0062       // below the parent momentum is not consistent; instead for now we use
0063       // `thrown_momentum` from truth above
0064       /*
0065       auto photon = sim_hit.getMCParticle();
0066       if(photon.parents_size()>0) thrown_momentum = edm4hep::utils::p(photon.getParents(0));
0067       */
0068       m_nphot_vs_p__transient->Fill(thrown_momentum);
0069
0070       // fill photon spectrum
0071       auto wavelength = Tools::GetPhotonWavelength(sim_hit, m_log);
0072       if(wavelength>=0)
0073         m_phot_spectrum->Fill(wavelength);
0074     }
0075
0076     // use `m_nphot_vs_p__transient` results to fill 2D hist `m_nphot_vs_p` for this event
0077     for(int b=1; b<=m_nphot_vs_p__transient->GetNbinsX(); b++) {
0078       auto nphot = m_nphot_vs_p__transient->GetBinContent(b);
0079       if(nphot>0) {
0080         auto momentum = m_nphot_vs_p__transient->GetBinCenter(b);
0081         m_nphot_vs_p->Fill(momentum,nphot);
0082       }
0083     }
0084
0085     // clear `m_nphot_vs_p__transient` to be ready for the next event
0086     m_nphot_vs_p__transient->Reset();
0087
0088   }
0089
0090
0091   // AlgorithmFinish
0092   //---------------------------------------------------------------------------
0093   void SimHitAnalysis::AlgorithmFinish() {
0094     // delete transient histograms, so they don't get written
0095     delete m_nphot_vs_p__transient;
0096   }
0097
0098 }