diff --git a/benchmarks/rich/include/RawHitAnalysis.h b/benchmarks/rich/include/RawHitAnalysis.h index 94d375c58b404fad350ff4e6f838a936f2ac3adc..047d8efa0a8ac2969ab358eeebcf051d03d696e4 100644 --- a/benchmarks/rich/include/RawHitAnalysis.h +++ b/benchmarks/rich/include/RawHitAnalysis.h @@ -10,6 +10,9 @@ #include <TMath.h> #include <edm4eic/RawTrackerHitCollection.h> +#include <edm4eic/MCRecoTrackerHitAssociationCollection.h> + +#include "Tools.h" namespace benchmarks { @@ -22,7 +25,10 @@ namespace benchmarks { // algorithm methods void AlgorithmInit(std::shared_ptr<spdlog::logger>& logger); - void AlgorithmProcess(const edm4eic::RawTrackerHitCollection& raw_hits); + void AlgorithmProcess( + const edm4eic::RawTrackerHitCollection& raw_hits, + const edm4eic::MCRecoTrackerHitAssociationCollection& assocs + ); void AlgorithmFinish(); private: @@ -35,6 +41,7 @@ namespace benchmarks { TH1D *m_adc_dist; TH1D *m_tdc_dist; TH2D *m_tdc_vs_adc; + TH1D *m_phot_spectrum; // logging std::shared_ptr<spdlog::logger> m_log; diff --git a/benchmarks/rich/include/SimHitAnalysis.h b/benchmarks/rich/include/SimHitAnalysis.h index d440d65ddec5396fe3d7f8bf953f202408ce36a0..7e21f1611b27ea1646bee13f1335050b68a1ecae 100644 --- a/benchmarks/rich/include/SimHitAnalysis.h +++ b/benchmarks/rich/include/SimHitAnalysis.h @@ -35,6 +35,7 @@ namespace benchmarks { // histograms TH2D *m_nphot_vs_p; TH1D *m_nphot_vs_p__transient; // transient (not written) + TH1D *m_phot_spectrum; // logger std::shared_ptr<spdlog::logger> m_log; diff --git a/benchmarks/rich/include/Tools.h b/benchmarks/rich/include/Tools.h index e62ec4deb07db0d03bb6f46f6eae472d0618f8e1..9b596e54a256d6f16084ed12fc4a782fafa801ab 100644 --- a/benchmarks/rich/include/Tools.h +++ b/benchmarks/rich/include/Tools.h @@ -17,8 +17,8 @@ #include <Evaluator/DD4hepUnits.h> // data model -#include <edm4eic/CherenkovParticleIDCollection.h> -#include <edm4hep/ParticleIDCollection.h> +#include <edm4hep/SimTrackerHitCollection.h> +#include <edm4hep/utils/kinematics.h> namespace benchmarks { @@ -33,6 +33,9 @@ namespace benchmarks { // obtained from `dd4hep::h_Planck * dd4hep::c_light / (dd4hep::GeV * dd4hep::nm)` static constexpr double HC = 1.239841875e-06; + // opticalphoton PDG + static constexpr int PHOTON_PDG = -22; + // ------------------------------------------------------------------------------------- // PDG mass lookup @@ -60,6 +63,24 @@ namespace benchmarks { static int GetNumPDGs() { return GetPDGMasses().size(); }; + + // ------------------------------------------------------------------------------------- + // get photon wavelength + static double GetPhotonWavelength(const edm4hep::SimTrackerHit& hit, std::shared_ptr<spdlog::logger> log) { + auto phot = hit.getMCParticle(); + if(!phot.isAvailable()) { + log->error("no MCParticle in hit"); + return -1.0; + } + if(phot.getPDG() != PHOTON_PDG) { + log->warn("hit MCParticle is not an opticalphoton; PDG is {}; ignoring", phot.getPDG()); + return -1.0; + } + auto momentum = edm4hep::utils::p(phot); + auto wavelength = momentum > 0 ? HC / momentum : 0.0; + return wavelength; + } + // ------------------------------------------------------------------------------------- // binning static constexpr int n_bins = 100; @@ -72,6 +93,9 @@ namespace benchmarks { static constexpr double theta_max = 300; static constexpr double thetaResid_max = 100; static constexpr int phi_bins = 100; + static constexpr int wl_bins = 120; + static constexpr double wl_min = 0; + static constexpr double wl_max = 1200; }; // class Tools } // namespace benchmarks diff --git a/benchmarks/rich/src/CherenkovPIDAnalysis.cc b/benchmarks/rich/src/CherenkovPIDAnalysis.cc index 9c1c5a3fcda145da55f8ca509d0d3f5807055b31..84b424c73016aa587b9fc64e66bc5174688db577 100644 --- a/benchmarks/rich/src/CherenkovPIDAnalysis.cc +++ b/benchmarks/rich/src/CherenkovPIDAnalysis.cc @@ -90,6 +90,18 @@ namespace benchmarks { Tools::momentum_bins, 0, Tools::momentum_max, Tools::GetNumPDGs()+1, 0, Tools::GetNumPDGs()+1 ); + + // format histograms + auto format1D = [] (auto h) { + h->SetLineColor(kAzure); + h->SetFillColor(kAzure); + }; + format1D(m_npe_dist); + format1D(m_theta_dist); + format1D(m_thetaResid_dist); + format1D(m_mcWavelength_dist); + format1D(m_mcRindex_dist); + format1D(m_highestWeight_dist); } diff --git a/benchmarks/rich/src/RawHitAnalysis.cc b/benchmarks/rich/src/RawHitAnalysis.cc index 3e05d855710c8d76912e708a73312809333e5f74..91bfde6e180d4b60a9a3e16603c4514cb9489a9a 100644 --- a/benchmarks/rich/src/RawHitAnalysis.cc +++ b/benchmarks/rich/src/RawHitAnalysis.cc @@ -20,6 +20,11 @@ namespace benchmarks { adc_max/2, 0, adc_max/2, tdc_max/2, 0, tdc_max/2 ); + m_phot_spectrum = new TH1D( + "phot_spectrum_rec", + "Photon wavelength at emission point -- only photons from digitized hits;#lambda [nm]", + Tools::wl_bins, Tools::wl_min, Tools::wl_max + ); // format histograms auto format1D = [] (auto h) { @@ -28,16 +33,20 @@ namespace benchmarks { }; format1D(m_adc_dist); format1D(m_tdc_dist); + m_phot_spectrum->SetLineColor(kViolet+2); + m_phot_spectrum->SetFillColor(kViolet+2); } // AlgorithmProcess //--------------------------------------------------------------------------- void RawHitAnalysis::AlgorithmProcess( - const edm4eic::RawTrackerHitCollection& raw_hits + const edm4eic::RawTrackerHitCollection& raw_hits, + const edm4eic::MCRecoTrackerHitAssociationCollection& assocs ) { - // loop over digitized hits + + // loop over all raw hits (including noise) for(const auto& raw_hit : raw_hits) { auto adc = raw_hit.getCharge(); auto tdc = raw_hit.getTimeStamp(); @@ -45,6 +54,15 @@ namespace benchmarks { m_tdc_dist->Fill(tdc); m_tdc_vs_adc->Fill(adc,tdc); } + + // loop over hits with associations (no noise) + for(const auto& assoc : assocs) { + for(const auto& hit : assoc.getSimHits()) { + auto wavelength = Tools::GetPhotonWavelength(hit, m_log); + if(wavelength>=0) + m_phot_spectrum->Fill(wavelength); + } + } } diff --git a/benchmarks/rich/src/SimHitAnalysis.cc b/benchmarks/rich/src/SimHitAnalysis.cc index d16424defd7e1aae8968b95004d7a94801cf0691..f39db5ab133e5c3f95f37bffba738abfaa556893 100644 --- a/benchmarks/rich/src/SimHitAnalysis.cc +++ b/benchmarks/rich/src/SimHitAnalysis.cc @@ -23,6 +23,13 @@ namespace benchmarks { m_nphot_vs_p->GetXaxis()->GetXmin(), m_nphot_vs_p->GetXaxis()->GetXmax() ); + m_phot_spectrum = new TH1D( + "phot_spectrum_sim", + "Photon wavelength at emission point -- all photons incident on sensors;#lambda [nm]", + Tools::wl_bins, Tools::wl_min, Tools::wl_max + ); + m_phot_spectrum->SetLineColor(kViolet+2); + m_phot_spectrum->SetFillColor(kViolet+2); } @@ -39,8 +46,10 @@ namespace benchmarks { part->SetSingleParticle(mc_parts); auto thrown_momentum = part->GetMomentum(); - // fill 1D thrown momentum histogram `m_nphot_vs_p__transient` for each (pre-digitized) sensor hit + // loop over hits for(const auto& sim_hit : sim_hits) { + + // - fill 1D thrown momentum histogram `m_nphot_vs_p__transient` for each (pre-digitized) sensor hit // FIXME: get thrown momentum from multi-track events; but in this attempt // below the parent momentum is not consistent; instead for now we use // `thrown_momentum` from truth above @@ -49,6 +58,11 @@ namespace benchmarks { if(photon.parents_size()>0) thrown_momentum = edm4hep::utils::p(photon.getParents(0)); */ m_nphot_vs_p__transient->Fill(thrown_momentum); + + // fill photon spectrum + auto wavelength = Tools::GetPhotonWavelength(sim_hit, m_log); + if(wavelength>=0) + m_phot_spectrum->Fill(wavelength); } // use `m_nphot_vs_p__transient` results to fill 2D hist `m_nphot_vs_p` for this event diff --git a/benchmarks/rich/src/benchmark.cc b/benchmarks/rich/src/benchmark.cc index 36af285c8f09e640a05adca45e728ab3bcad58b0..d988601f1727dd71d7e6a13548e4aa021b470359 100644 --- a/benchmarks/rich/src/benchmark.cc +++ b/benchmarks/rich/src/benchmark.cc @@ -195,7 +195,8 @@ int main(int argc, char** argv) { digi_algo->AlgorithmInit(algo.log); algo.process = [digi_algo] (podio::Frame& frame) { digi_algo->AlgorithmProcess( - GetCollection<edm4eic::RawTrackerHitCollection>(frame,"DRICHRawHits") + GetCollection<edm4eic::RawTrackerHitCollection>(frame,"DRICHRawHits"), + GetCollection<edm4eic::MCRecoTrackerHitAssociationCollection>(frame,"DRICHRawHitsAssociations") ); }; algo.finish = [digi_algo] () { digi_algo->AlgorithmFinish(); };