Skip to content
Snippets Groups Projects
Commit 4150e622 authored by Christopher Dilks's avatar Christopher Dilks
Browse files

feat: add photon spectra; plot aesthetics

parent dd838b0e
No related branches found
No related tags found
2 merge requests!309Irt algo,!293feat: dRICH benchmarks
...@@ -10,6 +10,9 @@ ...@@ -10,6 +10,9 @@
#include <TMath.h> #include <TMath.h>
#include <edm4eic/RawTrackerHitCollection.h> #include <edm4eic/RawTrackerHitCollection.h>
#include <edm4eic/MCRecoTrackerHitAssociationCollection.h>
#include "Tools.h"
namespace benchmarks { namespace benchmarks {
...@@ -22,7 +25,10 @@ namespace benchmarks { ...@@ -22,7 +25,10 @@ namespace benchmarks {
// algorithm methods // algorithm methods
void AlgorithmInit(std::shared_ptr<spdlog::logger>& logger); 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(); void AlgorithmFinish();
private: private:
...@@ -35,6 +41,7 @@ namespace benchmarks { ...@@ -35,6 +41,7 @@ namespace benchmarks {
TH1D *m_adc_dist; TH1D *m_adc_dist;
TH1D *m_tdc_dist; TH1D *m_tdc_dist;
TH2D *m_tdc_vs_adc; TH2D *m_tdc_vs_adc;
TH1D *m_phot_spectrum;
// logging // logging
std::shared_ptr<spdlog::logger> m_log; std::shared_ptr<spdlog::logger> m_log;
......
...@@ -35,6 +35,7 @@ namespace benchmarks { ...@@ -35,6 +35,7 @@ namespace benchmarks {
// histograms // histograms
TH2D *m_nphot_vs_p; TH2D *m_nphot_vs_p;
TH1D *m_nphot_vs_p__transient; // transient (not written) TH1D *m_nphot_vs_p__transient; // transient (not written)
TH1D *m_phot_spectrum;
// logger // logger
std::shared_ptr<spdlog::logger> m_log; std::shared_ptr<spdlog::logger> m_log;
......
...@@ -17,8 +17,8 @@ ...@@ -17,8 +17,8 @@
#include <Evaluator/DD4hepUnits.h> #include <Evaluator/DD4hepUnits.h>
// data model // data model
#include <edm4eic/CherenkovParticleIDCollection.h> #include <edm4hep/SimTrackerHitCollection.h>
#include <edm4hep/ParticleIDCollection.h> #include <edm4hep/utils/kinematics.h>
namespace benchmarks { namespace benchmarks {
...@@ -33,6 +33,9 @@ namespace benchmarks { ...@@ -33,6 +33,9 @@ namespace benchmarks {
// obtained from `dd4hep::h_Planck * dd4hep::c_light / (dd4hep::GeV * dd4hep::nm)` // obtained from `dd4hep::h_Planck * dd4hep::c_light / (dd4hep::GeV * dd4hep::nm)`
static constexpr double HC = 1.239841875e-06; static constexpr double HC = 1.239841875e-06;
// opticalphoton PDG
static constexpr int PHOTON_PDG = -22;
// ------------------------------------------------------------------------------------- // -------------------------------------------------------------------------------------
// PDG mass lookup // PDG mass lookup
...@@ -60,6 +63,24 @@ namespace benchmarks { ...@@ -60,6 +63,24 @@ namespace benchmarks {
static int GetNumPDGs() { return GetPDGMasses().size(); }; 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 // binning
static constexpr int n_bins = 100; static constexpr int n_bins = 100;
...@@ -72,6 +93,9 @@ namespace benchmarks { ...@@ -72,6 +93,9 @@ namespace benchmarks {
static constexpr double theta_max = 300; static constexpr double theta_max = 300;
static constexpr double thetaResid_max = 100; static constexpr double thetaResid_max = 100;
static constexpr int phi_bins = 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 }; // class Tools
} // namespace benchmarks } // namespace benchmarks
...@@ -90,6 +90,18 @@ namespace benchmarks { ...@@ -90,6 +90,18 @@ namespace benchmarks {
Tools::momentum_bins, 0, Tools::momentum_max, Tools::momentum_bins, 0, Tools::momentum_max,
Tools::GetNumPDGs()+1, 0, Tools::GetNumPDGs()+1 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);
} }
......
...@@ -20,6 +20,11 @@ namespace benchmarks { ...@@ -20,6 +20,11 @@ namespace benchmarks {
adc_max/2, 0, adc_max/2, adc_max/2, 0, adc_max/2,
tdc_max/2, 0, tdc_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 // format histograms
auto format1D = [] (auto h) { auto format1D = [] (auto h) {
...@@ -28,16 +33,20 @@ namespace benchmarks { ...@@ -28,16 +33,20 @@ namespace benchmarks {
}; };
format1D(m_adc_dist); format1D(m_adc_dist);
format1D(m_tdc_dist); format1D(m_tdc_dist);
m_phot_spectrum->SetLineColor(kViolet+2);
m_phot_spectrum->SetFillColor(kViolet+2);
} }
// AlgorithmProcess // AlgorithmProcess
//--------------------------------------------------------------------------- //---------------------------------------------------------------------------
void RawHitAnalysis::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) { for(const auto& raw_hit : raw_hits) {
auto adc = raw_hit.getCharge(); auto adc = raw_hit.getCharge();
auto tdc = raw_hit.getTimeStamp(); auto tdc = raw_hit.getTimeStamp();
...@@ -45,6 +54,15 @@ namespace benchmarks { ...@@ -45,6 +54,15 @@ namespace benchmarks {
m_tdc_dist->Fill(tdc); m_tdc_dist->Fill(tdc);
m_tdc_vs_adc->Fill(adc,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);
}
}
} }
......
...@@ -23,6 +23,13 @@ namespace benchmarks { ...@@ -23,6 +23,13 @@ namespace benchmarks {
m_nphot_vs_p->GetXaxis()->GetXmin(), m_nphot_vs_p->GetXaxis()->GetXmin(),
m_nphot_vs_p->GetXaxis()->GetXmax() 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 { ...@@ -39,8 +46,10 @@ namespace benchmarks {
part->SetSingleParticle(mc_parts); part->SetSingleParticle(mc_parts);
auto thrown_momentum = part->GetMomentum(); 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) { 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 // FIXME: get thrown momentum from multi-track events; but in this attempt
// below the parent momentum is not consistent; instead for now we use // below the parent momentum is not consistent; instead for now we use
// `thrown_momentum` from truth above // `thrown_momentum` from truth above
...@@ -49,6 +58,11 @@ namespace benchmarks { ...@@ -49,6 +58,11 @@ namespace benchmarks {
if(photon.parents_size()>0) thrown_momentum = edm4hep::utils::p(photon.getParents(0)); if(photon.parents_size()>0) thrown_momentum = edm4hep::utils::p(photon.getParents(0));
*/ */
m_nphot_vs_p__transient->Fill(thrown_momentum); 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 // use `m_nphot_vs_p__transient` results to fill 2D hist `m_nphot_vs_p` for this event
......
...@@ -195,7 +195,8 @@ int main(int argc, char** argv) { ...@@ -195,7 +195,8 @@ int main(int argc, char** argv) {
digi_algo->AlgorithmInit(algo.log); digi_algo->AlgorithmInit(algo.log);
algo.process = [digi_algo] (podio::Frame& frame) { algo.process = [digi_algo] (podio::Frame& frame) {
digi_algo->AlgorithmProcess( 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(); }; algo.finish = [digi_algo] () { digi_algo->AlgorithmFinish(); };
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment