Add TS DQM for digi diagnostics, have energy of the hit deducted from PE instead of Sim info (#1825)

Author: tvami

.github/validation_samples/target_pn_lyso/config.py

@@ -75,6 +75,7 @@
     dqm.TrigScintSimDQM('TargetSimHits','TargetSimHits','target'),
     dqm.TrigScintDigiDQM('TargetDigis','TargetDigis','target'),
     dqm.TrigScintClusterDQM('TargetClusters','TargetClusters','target'),
+    dqm.TrigScintDigiVerifierDQM('TrigScintDigiVerifier','TargetSimHits','TargetDigis'),
     ]
 
 

DQM/include/DQM/TrigScintDigiVerifier.h

@@ -0,0 +1,54 @@
+#ifndef DQM_TSDIGIVERIFIER_H
+#define DQM_TSDIGIVERIFIER_H
+
+#include <algorithm>
+
+// LDMX Framework
+#include "DetDescr/TrigScintID.h"
+#include "Framework/Configure/Parameters.h"
+#include "Framework/EventProcessor.h"
+#include "SimCore/Event/SimCalorimeterHit.h"
+#include "TrigScint/Event/TrigScintHit.h"
+
+namespace dqm {
+
+/**
+ * @class TrigScintDigiVerifier
+ * @brief Generate histograms to check digi pipeline performance
+ */
+class TrigScintDigiVerifier : public framework::Analyzer {
+ public:
+  /**
+   * Constructor
+   *
+   * Blank Analyzer constructor
+   */
+  TrigScintDigiVerifier(const std::string& name, framework::Process& process)
+      : framework::Analyzer(name, process) {}
+
+  /**
+   * Input python configuration parameters
+   */
+  virtual void configure(framework::config::Parameters& ps);
+
+  /**
+   * Fills histograms
+   */
+  virtual void analyze(const framework::Event& event);
+
+ private:
+  /// Collection Name for SimHits
+  std::string ts_simhit_coll_;
+
+  /// Pass Name for SimHits
+  std::string ts_simhit_pass_;
+
+  /// Collection Name for digis
+  std::string ts_digi_coll_;
+
+  /// Pass Name for digis
+  std::string ts_digi_pass_;
+};
+}  // namespace dqm
+
+#endif /* DQM_TSDIGIVERIFIER_H */

DQM/python/dqm.py

@@ -904,6 +904,24 @@ def __init__(self,name='TrigScintDigiUp',hit_coll='trigScintDigisUp',pad='up') :
         self.pad = pad
         self.trig_scint_passname = ''
 
+class TrigScintDigiVerifierDQM(ldmxcfg.Analyzer) :
+    def __init__(self, name = 'TrigScintDigiVerifier', ts_simhit_coll = 'TriggerPadUpSimHits', ts_digi_coll = 'trigScintDigisUp') :
+        super().__init__(name,'dqm::TrigScintDigiVerifier','DQM')
+
+        self.ts_simhit_coll = ts_simhit_coll
+        self.ts_simhit_pass = ''
+        self.ts_digi_coll = ts_digi_coll
+        self.ts_digi_pass = ''
+
+        self.build2DHistogram( "sim_edep:rec_amplitude" ,
+                "Simulated Energy [MeV]" , 1000 , 0. , 50. ,
+                "Reconstructed Amplitude [MeV]" , 1000 , 0. , 50. )
+
+        self.build2DHistogram( "sim_edep:rec_energy" ,
+                "Simulated Energy [MeV]" , 1000 , 0. , 50. ,
+                "Reconstructed Energy [MeV]" , 1000 , 0. , 50. )
+
+
 class TrigScintClusterDQM(ldmxcfg.Analyzer) :
     """Configured TrigScintClusterDQM python object
     

DQM/src/DQM/TrigScintDigiVerifier.cxx

@@ -0,0 +1,79 @@
+#include "DQM/TrigScintDigiVerifier.h"
+
+namespace dqm {
+
+void TrigScintDigiVerifier::configure(framework::config::Parameters &ps) {
+  ts_simhit_coll_ = ps.get<std::string>("ts_simhit_coll");
+  ts_simhit_pass_ = ps.get<std::string>("ts_simhit_pass");
+  ts_digi_coll_ = ps.get<std::string>("ts_digi_coll");
+  ts_digi_pass_ = ps.get<std::string>("ts_digi_pass");
+
+  return;
+}
+
+void TrigScintDigiVerifier::analyze(const framework::Event &event) {
+  // get truth information sorted into an ID based map
+  auto ts_simhits = event.getCollection<ldmx::SimCalorimeterHit>(
+      ts_simhit_coll_, ts_simhit_pass_);
+
+  // sort sim hits by ID
+  std::sort(ts_simhits.begin(), ts_simhits.end(),
+            [](const ldmx::SimCalorimeterHit &lhs,
+               const ldmx::SimCalorimeterHit &rhs) {
+              return lhs.getID() < rhs.getID();
+            });
+
+  auto ts_digis{
+      event.getCollection<ldmx::TrigScintHit>(ts_digi_coll_, ts_digi_pass_)};
+
+  // sort digi hits by ID
+  std::sort(ts_digis.begin(), ts_digis.end(),
+            [](const ldmx::TrigScintHit &lhs, const ldmx::TrigScintHit &rhs) {
+              return lhs.getID() < rhs.getID();
+            });
+
+  // Loop on the ts rechits
+  ldmx_log(info) << "There are " << ts_digis.size()
+                 << " ts digis in this event";
+  for (const auto &ts_digi : ts_digis) {
+    // skip anything that digi flagged as noise
+    if (ts_digi.isNoise()) {
+      ldmx_log(debug) << "Digi with raw ID " << ts_digi.getID()
+                      << " and bar ID " << ts_digi.getBarID()
+                      << " is flagged as noise, skipping";
+      continue;
+    }
+    int raw_id = ts_digi.getID();
+    ldmx_log(debug) << "Digi with raw ID " << raw_id << " and bar ID "
+                    << ts_digi.getBarID() << " has energy "
+                    << ts_digi.getEnergy() << " and amplitude "
+                    << ts_digi.getAmplitude();
+
+    // get information for this hit
+
+    double total_sim_energy_dep = 0.;
+    for (const auto &ts_simhit : ts_simhits) {
+      if (raw_id == ts_simhit.getID()) {
+        total_sim_energy_dep += ts_simhit.getEdep();
+      } else if (raw_id < ts_simhit.getID()) {
+        // later sim hits - all done
+        break;
+      }
+    }
+
+    ldmx_log(info) << " There are " << ts_simhits.size()
+                   << " sim hits in this event, adding up to a total energy of "
+                   << total_sim_energy_dep;
+
+    histograms_.fill("sim_edep:rec_amplitude", total_sim_energy_dep,
+                     ts_digi.getAmplitude());
+    histograms_.fill("sim_edep:rec_energy", total_sim_energy_dep,
+                     ts_digi.getEnergy());
+  }  // end of loop on ts digis
+
+  return;
+}
+
+}  // namespace dqm
+
+DECLARE_ANALYZER(dqm::TrigScintDigiVerifier);

TrigScint/python/trigScint.py

@@ -27,9 +27,6 @@ def __init__(self,name) :
         self.input_collection="TriggerPad3SimHits"
         self.input_pass_name="" #take any pass
         self.output_collection="trigScintDigisPad3"
-        import time
-        self.randomSeed = int(time.time())
-        self.verbose = False
 
         self.sim_particles_passname = ""
 
@@ -59,6 +56,8 @@ def target() :
         digi = TrigScintDigiProducer( 'TargetDigis' )
         digi.input_collection = 'TargetSimHits'
         digi.output_collection= 'TargetDigis'
+        digi.mev_per_mip = 1.
+        digi.pe_per_mip = 300.
         return digi
 
 

TrigScint/src/TrigScint/TrigScintDigiProducer.cxx

@@ -6,19 +6,17 @@ TrigScintDigiProducer::TrigScintDigiProducer(const std::string &name,
                                              framework::Process &process)
     : Producer(name, process) {}
 
-void TrigScintDigiProducer::configure(
-    framework::config::Parameters &parameters) {
+void TrigScintDigiProducer::configure(framework::config::Parameters &ps) {
   // Configure this instance of the producer
-  strips_per_array_ = parameters.get<int>("number_of_strips");
-  number_of_arrays_ = parameters.get<int>("number_of_arrays");
-  mean_noise_ = parameters.get<double>("mean_noise");
-  mev_per_mip_ = parameters.get<double>("mev_per_mip");
-  pe_per_mip_ = parameters.get<double>("pe_per_mip");
-  input_collection_ = parameters.get<std::string>("input_collection");
-  input_pass_name_ = parameters.get<std::string>("input_pass_name");
-  output_collection_ = parameters.get<std::string>("output_collection");
-  sim_particles_passname_ =
-      parameters.get<std::string>("sim_particles_passname");
+  strips_per_array_ = ps.get<int>("number_of_strips");
+  number_of_arrays_ = ps.get<int>("number_of_arrays");
+  mean_noise_ = ps.get<double>("mean_noise");
+  mev_per_mip_ = ps.get<double>("mev_per_mip");
+  pe_per_mip_ = ps.get<double>("pe_per_mip");
+  input_collection_ = ps.get<std::string>("input_collection");
+  input_pass_name_ = ps.get<std::string>("input_pass_name");
+  output_collection_ = ps.get<std::string>("output_collection");
+  sim_particles_passname_ = ps.get<std::string>("sim_particles_passname");
 }
 
 void TrigScintDigiProducer::onNewRun(const ldmx::RunHeader &) {
@@ -46,7 +44,7 @@ ldmx::TrigScintID TrigScintDigiProducer::generateRandomID(int module) {
 }
 
 void TrigScintDigiProducer::produce(framework::Event &event) {
-  std::map<ldmx::TrigScintID, int> cell_p_es, cell_min_p_es;
+  std::map<ldmx::TrigScintID, int> cell_pes, cell_min_p_es;
   std::map<ldmx::TrigScintID, float> xpos, ypos, zpos, edep, time, beam_frac;
   std::set<ldmx::TrigScintID> noise_hit_i_ds;
 
@@ -73,7 +71,7 @@ void TrigScintDigiProducer::produce(framework::Event &event) {
     for (int i = 0; i < sim_hit.getNumberOfContribs(); i++) {
       auto contrib = sim_hit.getContrib(i);
 
-      ldmx_log(trace) << "contrib " << i << " trackID: " << contrib.trackID
+      ldmx_log(trace) << "Contrib " << i << " trackID: " << contrib.trackID
                       << " pdgID: " << contrib.pdgCode
                       << " edep: " << contrib.edep;
       ldmx_log(trace) << "\t particle id: "
@@ -128,21 +126,26 @@ void TrigScintDigiProducer::produce(framework::Event &event) {
     xpos[id] = xpos[id] / edep[id];
     ypos[id] = ypos[id] / edep[id];
     zpos[id] = zpos[id] / edep[id];
+    //  mean number of photoelectrons produced for the given deposited energy
     double mean_pe = dep_energy / mev_per_mip_ * pe_per_mip_;
     std::poisson_distribution<int> poisson_dist(mean_pe + mean_noise_);
-    cell_p_es[id] = poisson_dist(rng_);
+    cell_pes[id] = poisson_dist(rng_);
+    // energy corresponding to the number of PEs observed
+    // the minimum number of PEs is the mean number of PEs minus the noise
+    double energy_per_pe = mev_per_mip_ / pe_per_mip_;
+    double cell_energy = energy_per_pe * cell_pes[id];
 
     // If a cell has a PE count above threshold, persit the hit.
     // Thresholds are introduced (and configurable) in clustering.
     // the cell PE >=1 suppresses artifical noise that is below one light
     // quantum in the SiPM and unphysical.
-    if (cell_p_es[id] >= 1) {
+    if (cell_pes[id] >= 1) {
       ldmx::TrigScintHit hit;
       hit.setID(id.raw());
-      hit.setPE(cell_p_es[id]);
+      hit.setPE(cell_pes[id]);
       hit.setMinPE(cell_min_p_es[id]);
-      hit.setAmplitude(cell_p_es[id]);
-      hit.setEnergy(dep_energy);
+      hit.setAmplitude(cell_pes[id]);
+      hit.setEnergy(cell_energy);
       hit.setTime(time[id]);
       hit.setXPos(xpos[id]);
       hit.setYPos(ypos[id]);
@@ -155,11 +158,11 @@ void TrigScintDigiProducer::produce(framework::Event &event) {
       trig_scint_hits.push_back(hit);
     }
 
-    ldmx_log(trace) << " ID = " << id.raw() << " Edep: " << edep[id]
-                    << " numPEs: " << cell_p_es[id] << " time: " << time[id]
+    ldmx_log(debug) << " ID = " << id.raw() << " Edep: " << edep[id]
+                    << " numPEs: " << cell_pes[id] << " time: " << time[id]
                     << " z: " << zpos[id] << "\t X: " << xpos[id]
                     << " Y: " << ypos[id] << " Z: " << zpos[id];
-  }
+  }  // end of loop over detIDs
 
   // ------------------------------- Noise simulation -----------------------//
   // ------------------------------------------------------------------------//