First full working version of calo/scint digi

Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx

@@ -2,13 +2,15 @@
 
 #include "Identifier/Identifier.h"
 
-#include <vector>
+#include "FaserCaloSimEvent/CaloHitIdHelper.h"
+
 #include <map>
 #include <utility>
 
 CaloWaveformDigiAlg::CaloWaveformDigiAlg(const std::string& name, 
 					 ISvcLocator* pSvcLocator)
-  : AthReentrantAlgorithm(name, pSvcLocator) { 
+  : AthReentrantAlgorithm(name, pSvcLocator)
+{ 
 
 }
 
@@ -19,26 +21,25 @@ CaloWaveformDigiAlg::initialize() {
   // Initalize tools
   ATH_CHECK( m_digiTool.retrieve() );
 
-
   // Set key to read waveform from
   ATH_CHECK( m_caloHitContainerKey.initialize() );
 
   // Set key to write container
   ATH_CHECK( m_waveformContainerKey.initialize() );
 
-  // Will eventually depend on the type of detector
-  // TODO: Vary time at which centre it?
-  // TODO: Change params compared to scint
-  //  m_kernel = new TF1("PDF", " ROOT::Math::crystalball_pdf(x, -0.9, 10, 4, 900)", 0, 1200);
-  
+  // Set up  helper
+  ATH_CHECK(detStore()->retrieve(m_ecalID, "EcalID"));
+
+  // Create CB time kernel and pre-evaluate for number of samples
   m_kernel = new TF1("PDF", "[4] * ROOT::Math::crystalball_pdf(x, [0],[1],[2],[3])", 0, 1200);
-  //m_kernel->SetParameters(-0.25,10,4,900);                                                      
   m_kernel->SetParameter(0, m_CB_alpha);
   m_kernel->SetParameter(1, m_CB_n);
   m_kernel->SetParameter(2, m_CB_sigma);
   m_kernel->SetParameter(3, m_CB_mean);
   m_kernel->SetParameter(4, m_CB_norm);
 
+  // Pre-evaluate time kernel for each bin
+  m_timekernel = m_digiTool->evaluate_timekernel(m_kernel);
 
   return StatusCode::SUCCESS;
 }
@@ -55,12 +56,10 @@ CaloWaveformDigiAlg::finalize() {
 StatusCode 
 CaloWaveformDigiAlg::execute(const EventContext& ctx) const {
   ATH_MSG_DEBUG("Executing");
-
-  ATH_MSG_DEBUG("Run: " << ctx.eventID().run_number() 
-		<< " Event: " << ctx.eventID().event_number());
+  ATH_MSG_DEBUG("Run: " << ctx.eventID().run_number() << " Event: " << ctx.eventID().event_number());
 
   // Find the input HIT collection
-  SG::ReadHandle<CaloHitCollection> caloHitHandle(m_caloHitContainerKey, ctx);
+   SG::ReadHandle<CaloHitCollection> caloHitHandle(m_caloHitContainerKey, ctx);
 
   ATH_CHECK( caloHitHandle.isValid() );
   ATH_MSG_DEBUG("Found ReadHandle for CaloHitCollection " << m_caloHitContainerKey);
@@ -75,11 +74,50 @@ CaloWaveformDigiAlg::execute(const EventContext& ctx) const {
     ATH_MSG_DEBUG("CaloHitCollection found with zero length!");
     return StatusCode::SUCCESS;
   }
-  
-  // Digitise the hits
-  CHECK( m_digiTool->digitise<CaloHitCollection>(caloHitHandle.ptr(),
-						 waveformContainerHandle.ptr(), m_kernel,
-						 std::pair<float, float>(m_base_mean, m_base_rms)) );
+
+  // Create structure to store pulse for each channel
+  std::map<Identifier, std::vector<uint16_t>> waveforms = m_digiTool->create_waveform_map(m_ecalID);
+
+  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();
+    }
+
+    // Subtract count from basleine and add result to correct waveform vector
+    for (const auto& c : counts) {
+
+      unsigned int baseline = m_digiTool->generate_baseline(m_base_mean, m_base_rms);
+      int value = baseline - c.second;
+
+      if (value < 0) {
+	ATH_MSG_WARNING("Found pulse " << c.second << " larger than baseline " << c.first);
+	value = 0; // Protect against scaling signal above baseline
+      }
+
+      // Convert hit id to Identifier and store 
+      Identifier id = CaloHitIdHelper::GetHelper()->getIdentifier(c.first);
+      waveforms[id].push_back(value);
+    }
+  }
+
+
+  //m_chrono->chronoStop("Digit");
+  //m_chrono->chronoStart("Write");
+
+  // Loop over wavefrom vectors to make and store waveform
+  unsigned int nsamples = m_digiTool->nsamples();
+  for (const auto& w : waveforms) {
+    RawWaveform* wfm = new RawWaveform();
+    wfm->setWaveform(0, w.second);
+    wfm->setIdentifier(w.first);
+    wfm->setSamples(nsamples);
+    waveformContainerHandle->push_back(wfm);
+  }
+
+  //m_chrono->chronoStop("Write");
 
   ATH_MSG_DEBUG("WaveformsHitContainer " << waveformContainerHandle.name() << "' filled with "<< waveformContainerHandle->size() <<" items");