From 40f66c6edabd38da22d24be6802a7dbb30943688 Mon Sep 17 00:00:00 2001
From: Carl Gwilliam <gwilliam@hep.ph.liv.ac.uk>
Date: Fri, 27 Jan 2023 12:11:45 +0000
Subject: [PATCH] Speed up scint/calo digi by factorising e and time loops

---
 .../CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx  | 16 +++-
 .../src/ScintWaveformDigiAlg.cxx              | 76 +++++++++++++------
 2 files changed, 65 insertions(+), 27 deletions(-)

diff --git a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx
index 9f62fd57..8ea1d4a3 100644
--- a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx
+++ b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx
@@ -99,12 +99,24 @@ CaloWaveformDigiAlg::execute(const EventContext& ctx) const {
   // Create structure to store pulse for each channel
   std::map<Identifier, std::vector<uint16_t>> waveforms = m_digiTool->create_waveform_map(m_ecalID);
 
+  // Sum energy for each channel (i.e. identifier)
+  std::map<unsigned int, float> esum;  
+  for (const auto& hit : *caloHitHandle) { 
+    esum[hit.identify()] += hit.energyLoss();
+  }
+
+  // Loop over time samples
   for (const auto& tk : m_timekernel) {
     std::map<unsigned int, float> counts;    
 
     // Convolve hit energy with evaluated kernel and sum for each hit id (i.e. channel)
-    for (const auto& hit : *caloHitHandle) { 
-      counts[hit.identify()] += tk * hit.energyLoss();
+    //for (const auto& hit : *caloHitHandle) { 
+    //  counts[hit.identify()] += tk * hit.energyLoss();
+    //}
+
+    // Convolve summed energy with evaluated kernel for each hit id (i.e. channel)
+    for (const auto& e : esum) {
+      counts[e.first] = tk * e.second;
     }
 
     // Subtract count from basleine and add result to correct waveform vector
diff --git a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
index ea03a77f..3cb70517 100644
--- a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
+++ b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
@@ -110,35 +110,61 @@ ScintWaveformDigiAlg::execute(const EventContext& ctx) const {
     waveforms = m_digiTool->create_waveform_map(m_preshowerID);
   } 
 
+  // Sum energy for each channel (i.e. identifier)
+  std::map<Identifier, float> esum;  
+
+  Identifier id;
+  for (const auto& hit : *scintHitHandle) { 
+    if (first.isTrigger()) {
+      Identifier plate_id = m_triggerID->plate_id(hit.getIdentifier());      
+      // Need to do something better here, but for now just fill both
+      id = m_triggerID->pmt_id(plate_id, 0); // ID for PMT 0
+      esum[id] += hit.energyLoss();
+      id = m_triggerID->pmt_id(plate_id, 1); // ID for PMT 1
+      esum[id] += hit.energyLoss();
+    } else {
+      // All others have only 1 PMT
+      // Use detector id (not hit id) throughout
+      id = hit.getIdentifier();
+      esum[id] += hit.energyLoss();
+    }
+  }
+
   // Loop over time samples
-  // Should really do sums first, as repeating these in the loop is a waste
   for (const auto& tk : m_timekernel) {
     std::map<Identifier, float> counts;
    
-    // Convolve hit energy with evaluated kernel and sum for each hit id (i.e. channel)
-    Identifier id;
-    for (const auto& hit : *scintHitHandle) {            
-
-      // Special handling for trigger scintillator
-      if (first.isTrigger()) {
-	// Hits are created per plate
-	// Need to split between 2 PMTs
-
-	// Identifier from hit is plate ID
-	Identifier plate_id = m_triggerID->plate_id(hit.getIdentifier());
-
-	// Need to do something better here, but for now just fill both
-	id = m_triggerID->pmt_id(plate_id, 0); // ID for PMT 0
-	counts[id] += tk * hit.energyLoss();
-	id = m_triggerID->pmt_id(plate_id, 1); // ID for PMT 1
-	counts[id] += tk * hit.energyLoss();
-
-      } else {
-	// All others have only 1 PMT
-	// Use detector id (not hit id) throughout
-	id = hit.getIdentifier();
-	counts[id] += tk * hit.energyLoss();
-      }
+//     // Convolve hit energy with evaluated kernel and sum for each hit id (i.e. channel)
+//     Identifier id;
+//     for (const auto& hit : *scintHitHandle) {            
+// 
+//       // Special handling for trigger scintillator
+//       if (first.isTrigger()) {
+// 	// Hits are created per plate
+// 	// Need to split between 2 PMTs
+// 
+// 	// Identifier from hit is plate ID
+// 	Identifier plate_id = m_triggerID->plate_id(hit.getIdentifier());
+// 
+// 	// Need to do something better here, but for now just fill both
+// 	id = m_triggerID->pmt_id(plate_id, 0); // ID for PMT 0
+// 	counts[id] += tk * hit.energyLoss();
+// 	id = m_triggerID->pmt_id(plate_id, 1); // ID for PMT 1
+// 	counts[id] += tk * hit.energyLoss();
+// 
+//       } else {
+// 	// All others have only 1 PMT
+// 	// Use detector id (not hit id) throughout
+// 	id = hit.getIdentifier();
+// 	counts[id] += tk * hit.energyLoss();
+//       }
+//     }
+
+    // Convolve summed energy with evaluated kernel for each hit id (i.e. channel)
+    for (const auto& e : esum) {
+      // Convert hit id to Identifier and store 
+      id = e.first;
+      counts[id] = tk * e.second;
     }
 
     // Subtract count from basleine and add result to correct waveform vector
-- 
GitLab