From f40ec6d7594c68ab781d5ba616e2998e3cfd3d3a Mon Sep 17 00:00:00 2001
From: Christopher Dilks <c-dilks@users.noreply.github.com>
Date: Wed, 7 Jun 2023 20:10:05 -0400
Subject: [PATCH] feat: photon and hit multiplicities

---
 benchmarks/rich/draw_benchmark.py        |  8 ++++++--
 benchmarks/rich/include/RawHitAnalysis.h |  1 +
 benchmarks/rich/include/SimHitAnalysis.h |  1 +
 benchmarks/rich/src/RawHitAnalysis.cc    | 10 +++++++++-
 benchmarks/rich/src/SimHitAnalysis.cc    | 12 ++++++++++--
 5 files changed, 27 insertions(+), 5 deletions(-)

diff --git a/benchmarks/rich/draw_benchmark.py b/benchmarks/rich/draw_benchmark.py
index 8d1711b3..9041184e 100755
--- a/benchmarks/rich/draw_benchmark.py
+++ b/benchmarks/rich/draw_benchmark.py
@@ -61,14 +61,18 @@ canv_dict = {
 
 # draw photon spectra
 canv = canv_dict["photon_spectra"]
-canv.Divide(1,2)
-for i in range(2):
+canv.Divide(2,2)
+for i in range(4):
     canv.GetPad(i+1).SetGrid(1,1)
     canv.GetPad(i+1).SetLogy()
 canv.cd(1)
 ana_file.Get("phot/phot_spectrum_sim").Draw()
 canv.cd(2)
+ana_file.Get("phot/nphot_dist").Draw()
+canv.cd(3)
 ana_file.Get("digi/phot_spectrum_rec").Draw()
+canv.cd(4)
+ana_file.Get("digi/nhits_dist").Draw()
 
 # draw digitization
 canv = canv_dict["digitization"]
diff --git a/benchmarks/rich/include/RawHitAnalysis.h b/benchmarks/rich/include/RawHitAnalysis.h
index 047d8efa..56b0b839 100644
--- a/benchmarks/rich/include/RawHitAnalysis.h
+++ b/benchmarks/rich/include/RawHitAnalysis.h
@@ -42,6 +42,7 @@ namespace benchmarks {
       TH1D *m_tdc_dist;
       TH2D *m_tdc_vs_adc;
       TH1D *m_phot_spectrum;
+      TH1D *m_nhits_dist;
 
       // logging
       std::shared_ptr<spdlog::logger> m_log;
diff --git a/benchmarks/rich/include/SimHitAnalysis.h b/benchmarks/rich/include/SimHitAnalysis.h
index 7e21f161..836a6ad2 100644
--- a/benchmarks/rich/include/SimHitAnalysis.h
+++ b/benchmarks/rich/include/SimHitAnalysis.h
@@ -33,6 +33,7 @@ namespace benchmarks {
     private:
 
       // histograms
+      TH1D *m_nphot_dist;
       TH2D *m_nphot_vs_p;
       TH1D *m_nphot_vs_p__transient; // transient (not written)
       TH1D *m_phot_spectrum;
diff --git a/benchmarks/rich/src/RawHitAnalysis.cc b/benchmarks/rich/src/RawHitAnalysis.cc
index 91bfde6e..d55ac99d 100644
--- a/benchmarks/rich/src/RawHitAnalysis.cc
+++ b/benchmarks/rich/src/RawHitAnalysis.cc
@@ -22,9 +22,12 @@ namespace benchmarks {
         );
     m_phot_spectrum = new TH1D(
         "phot_spectrum_rec",
-        "Photon wavelength at emission point -- only photons from digitized hits;#lambda [nm]",
+        "Photon wavelength for digitized hits;#lambda [nm], at emission point",
         Tools::wl_bins, Tools::wl_min, Tools::wl_max
         );
+    m_nhits_dist = new TH1D("nhits_dist", "Number of digitized hits;N_{hits}",
+        Tools::npe_max, 0, Tools::npe_max
+        );
 
     // format histograms
     auto format1D = [] (auto h) {
@@ -35,6 +38,8 @@ namespace benchmarks {
     format1D(m_tdc_dist);
     m_phot_spectrum->SetLineColor(kViolet+2);
     m_phot_spectrum->SetFillColor(kViolet+2);
+    m_nhits_dist->SetLineColor(kCyan+2);
+    m_nhits_dist->SetFillColor(kCyan+2);
   }
 
 
@@ -46,6 +51,9 @@ namespace benchmarks {
       )
   {
 
+    // fill nhits
+    m_nhits_dist->Fill(raw_hits.size());
+
     // loop over all raw hits (including noise)
     for(const auto& raw_hit : raw_hits) {
       auto adc = raw_hit.getCharge();
diff --git a/benchmarks/rich/src/SimHitAnalysis.cc b/benchmarks/rich/src/SimHitAnalysis.cc
index f39db5ab..f204b512 100644
--- a/benchmarks/rich/src/SimHitAnalysis.cc
+++ b/benchmarks/rich/src/SimHitAnalysis.cc
@@ -14,7 +14,10 @@ namespace benchmarks {
     m_log = logger;
 
     // initialize histograms
-    m_nphot_vs_p = new TH2D("nphot_vs_p", "N_{photons} vs. Thrown Momentum;p [GeV];N_{photons}",
+    m_nphot_dist = new TH1D("nphot_dist", "Number of incident photons;N_{photons}",
+        Tools::nphot_max, 0, Tools::nphot_max
+        );
+    m_nphot_vs_p = new TH2D("nphot_vs_p", "Number of incident photons vs. Thrown Momentum;p [GeV];N_{photons}",
         Tools::momentum_bins, 0, Tools::momentum_max,
         Tools::nphot_max,     0, Tools::nphot_max
         );
@@ -25,11 +28,13 @@ namespace benchmarks {
         );
     m_phot_spectrum = new TH1D(
         "phot_spectrum_sim",
-        "Photon wavelength at emission point -- all photons incident on sensors;#lambda [nm]",
+        "Incident photon wavelength;#lambda [nm], at emission point",
         Tools::wl_bins, Tools::wl_min, Tools::wl_max
         );
     m_phot_spectrum->SetLineColor(kViolet+2);
     m_phot_spectrum->SetFillColor(kViolet+2);
+    m_nphot_dist->SetLineColor(kCyan+2);
+    m_nphot_dist->SetFillColor(kCyan+2);
   }
 
 
@@ -46,6 +51,9 @@ namespace benchmarks {
     part->SetSingleParticle(mc_parts);
     auto thrown_momentum = part->GetMomentum();
 
+    // fill nphot distribution
+    m_nphot_dist->Fill(sim_hits.size());
+
     // loop over hits
     for(const auto& sim_hit : sim_hits) {
 
-- 
GitLab