Merge branch 'scint_calo_digi' into 'master'

First full working version of calo/scint digi

See merge request faser/calypso!208

Calorimeter/CaloDigiAlgs/CMakeLists.txt

@@ -9,7 +9,7 @@ atlas_subdir( CaloDigiAlgs )
 atlas_add_component( CaloDigiAlgs
                      src/*.cxx src/*.h
                      src/components/*.cxx
-                     LINK_LIBRARIES AthenaBaseComps Identifier StoreGateLib WaveRawEvent FaserCaloSimEvent WaveDigiToolsLib)
+                     LINK_LIBRARIES AthenaBaseComps Identifier FaserCaloIdentifier  StoreGateLib WaveRawEvent FaserCaloSimEvent WaveDigiToolsLib)
 
 atlas_install_python_modules( python/*.py )
 

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");
 

Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h

@@ -19,11 +19,15 @@
 #include "GaudiKernel/ServiceHandle.h"
 #include "GaudiKernel/ToolHandle.h"
 
+// Helpers
+#include "FaserCaloIdentifier/EcalID.h"
+
 // ROOT
 #include "TF1.h"
 
 // STL
 #include <string>
+#include <vector>
 
 class CaloWaveformDigiAlg : public AthReentrantAlgorithm {
 
@@ -48,6 +52,8 @@ class CaloWaveformDigiAlg : public AthReentrantAlgorithm {
   CaloWaveformDigiAlg &operator=(const CaloWaveformDigiAlg&) = delete;
   //@}
 
+  /** @name Steerable pameters for crystal ball and baseline **/
+  //@{
   Gaudi::Property<double> m_CB_alpha {this, "CB_alpha", 0, "Alpha of the crystal ball function"};
   Gaudi::Property<double> m_CB_n {this, "CB_n", 0, "n of the crystal ball function"};
   Gaudi::Property<double> m_CB_mean {this, "CB_mean", 0, "Mean of the crystal ball function"};
@@ -56,10 +62,16 @@ class CaloWaveformDigiAlg : public AthReentrantAlgorithm {
 
   Gaudi::Property<double> m_base_mean {this, "base_mean", 0, "Mean of the baseline"};
   Gaudi::Property<double> m_base_rms {this, "base_rms", 0, "RMS of the baseline"};
+  //@}
 
-  /// Kernel PDF
-  TF1* m_kernel;
+  /** Kernel PDF and evaluated values **/
+  //@{
+  TF1*                            m_kernel;
+  std::vector<float>              m_timekernel;
+  //@}
 
+  /// Detector ID helper
+  const EcalID* m_ecalID{nullptr};
 
   /**
    * @name Digitisation tool
@@ -89,4 +101,6 @@ class CaloWaveformDigiAlg : public AthReentrantAlgorithm {
 
 };
 
+
+
 #endif // CALODIGIALGS_CALODIGIALG_H

Calorimeter/FaserCaloSimEvent/CMakeLists.txt

@@ -18,11 +18,11 @@ atlas_add_library( FaserCaloSimEvent
                    PRIVATE_INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${GEANT4_INCLUDE_DIRS}
                    DEFINITIONS ${CLHEP_DEFINITIONS}
                    LINK_LIBRARIES ${CLHEP_LIBRARIES} AthAllocators AthenaKernel CxxUtils GeneratorObjects HitManagement StoreGateLib SGtests
-                   PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES} FaserCaloIdentifier )
+                   PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES} FaserCaloIdentifier Identifier)
 
 atlas_add_dictionary( FaserCaloSimEventDict
                       FaserCaloSimEvent/CaloSimEventDict.h
                       FaserCaloSimEvent/selection.xml
                       INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS} ${GEANT4_INCLUDE_DIRS}
-                      LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} AthAllocators CxxUtils GeneratorObjects HitManagement StoreGateLib SGtests FaserCaloIdentifier FaserCaloSimEvent )
+                      LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} AthAllocators CxxUtils GeneratorObjects HitManagement StoreGateLib SGtests FaserCaloIdentifier FaserCaloSimEvent Identifier)
 

Calorimeter/FaserCaloSimEvent/FaserCaloSimEvent/CaloHitIdHelper.h

@@ -23,6 +23,11 @@
 
 // This class is singleton and static method and variable are used.
 #include "CxxUtils/checker_macros.h"
+
+#include "Identifier/Identifier.h"
+
+#include "FaserCaloIdentifier/EcalID.h"
+
 ATLAS_NO_CHECK_FILE_THREAD_SAFETY;
 
 class CaloHitIdHelper : HitIdHelper {
@@ -43,6 +48,8 @@ class CaloHitIdHelper : HitIdHelper {
   // Left or Right
   int getModule(const int& hid) const;
 
+  Identifier getIdentifier(const int& hid) const;
+
   //
   // Info packing:
   int buildHitId(const int, const int) const;
@@ -54,6 +61,9 @@ class CaloHitIdHelper : HitIdHelper {
   //
   // Initialize the helper, only called by the constructor
   void Initialize();
+
+  /// Detector ID helper
+  const EcalID* m_ecalID{nullptr};
 };
 
 #endif // CALOSIMEVENT_CALOHITIDHELPER

Calorimeter/FaserCaloSimEvent/src/CaloHitIdHelper.cxx

@@ -6,8 +6,6 @@
 
 #include "FaserCaloSimEvent/CaloHitIdHelper.h"
 #include "StoreGate/StoreGateSvc.h"
-#include "StoreGate/StoreGateSvc.h"
-#include "FaserCaloIdentifier/EcalID.h"
 
 #include "G4Types.hh"
 #ifdef G4MULTITHREADED
@@ -42,10 +40,9 @@ void CaloHitIdHelper::Initialize() {
   // determine whether hits were created with an SLHC dictionary
   // in which case eta module field is expanded.
   // Need to lock this thread-unsafe retrieval
-  const EcalID* pix;
   ServiceHandle<StoreGateSvc> detStore ("DetectorStore", "CaloHitIdHelper");
   if (detStore.retrieve().isSuccess()) {
-    if (detStore->retrieve(pix, "EcalID").isFailure()) { pix = 0; }
+    if (detStore->retrieve(m_ecalID, "EcalID").isFailure()) { m_ecalID = 0; }
   }
 
   InitializeField("Row", 0, 2);
@@ -64,6 +61,14 @@ int CaloHitIdHelper::getModule(const int& hid) const
   return this->GetFieldValue("Module", hid);
 }
 
+// identifier
+Identifier CaloHitIdHelper::getIdentifier(const int& hid) const
+{
+  return m_ecalID->pmt_id(getRow(hid), getModule(hid), 0);
+
+}
+
+
 //
 // Info packing:
 int CaloHitIdHelper::buildHitId( const int row, 

Scintillator/ScintDigiAlgs/CMakeLists.txt

@@ -9,7 +9,7 @@ atlas_subdir( ScintDigiAlgs )
 atlas_add_component( ScintDigiAlgs
                      src/*.cxx src/*.h
                      src/components/*.cxx
-                     LINK_LIBRARIES AthenaBaseComps Identifier StoreGateLib WaveRawEvent ScintSimEvent WaveDigiToolsLib)
+                     LINK_LIBRARIES AthenaBaseComps Identifier ScintIdentifier StoreGateLib WaveRawEvent ScintSimEvent WaveDigiToolsLib)
 
 atlas_install_python_modules( python/*.py )
 

Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py

@@ -15,6 +15,13 @@ dict_CB_param["Timing"]=dict(CB_alpha=-0.32, CB_n=65, CB_mean=846, CB_sigma=5.3,
 dict_CB_param["Veto"]=dict(CB_alpha=-0.38, CB_n=25, CB_mean=815, CB_sigma=7.7, CB_norm = 1000) # copy from Trigger; Veto was not in TestBeam, but in sim "Veto" is the TestBeam Trigger component
 dict_CB_param["Preshower"]=dict(CB_alpha=-0.32, CB_n=65, CB_mean=846, CB_sigma=5.3, CB_norm = 500)
 
+dict_baseline_params = {
+    "Trigger"   : {"mean" : 15000, "rms" : 3},
+    "Timing"    : {"mean" : 15000, "rms" : 3},
+    "Veto"      : {"mean" : 15000, "rms" : 3},
+    "Preshower" : {"mean" : 15000, "rms" : 3},        
+    }
+
 # One stop shopping for normal FASER data
 def ScintWaveformDigitizationCfg(flags):
     """ Return all algorithms and tools for Waveform digitization """
@@ -48,8 +55,8 @@ def ScintWaveformDigiCfg(flags, name="ScintWaveformDigiAlg", source="", **kwargs
     digiAlg.CB_sigma = dict_CB_param[source]["CB_sigma"]
     digiAlg.CB_norm = dict_CB_param[source]["CB_norm"]    
 
-    digiAlg.base_mean = 15000
-    digiAlg.base_rms = 3
+    digiAlg.base_mean = dict_baseline_params[source]["mean"]
+    digiAlg.base_rms = dict_baseline_params[source]["rms"]
      
     kwargs.setdefault("WaveformDigitisationTool", tool)
 

Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx

-
 #include "ScintWaveformDigiAlg.h"
 
-#include "Identifier/Identifier.h"
+#include "ScintSimEvent/ScintHitIdHelper.h"
 
-#include <vector>
 #include <map>
-#include <utility>
 
 
 ScintWaveformDigiAlg::ScintWaveformDigiAlg(const std::string& name, 
@@ -21,25 +18,28 @@ ScintWaveformDigiAlg::initialize() {
   // Initalize tools
   ATH_CHECK( m_digiTool.retrieve() );
 
-
   // Set key to read waveform from
   ATH_CHECK( m_scintHitContainerKey.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: Better parameters
- 
+  // Set up helpers
+  ATH_CHECK(detStore()->retrieve(m_vetoID, "VetoID"));
+  ATH_CHECK(detStore()->retrieve(m_triggerID, "TriggerID"));
+  ATH_CHECK(detStore()->retrieve(m_preshowerID, "PreshowerID"));
 
-  
+  // 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->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,11 +55,9 @@ ScintWaveformDigiAlg::finalize() {
 StatusCode 
 ScintWaveformDigiAlg::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
+  // Find the input HITS collection
   SG::ReadHandle<ScintHitCollection> scintHitHandle(m_scintHitContainerKey, ctx);
 
   ATH_CHECK( scintHitHandle.isValid() );
@@ -70,16 +68,59 @@ ScintWaveformDigiAlg::execute(const EventContext& ctx) const {
   ATH_CHECK( waveformContainerHandle.record( std::make_unique<RawWaveformContainer>()) );
 
   ATH_MSG_DEBUG("WaveformsContainer '" << waveformContainerHandle.name() << "' initialized");
-
+  
   if (scintHitHandle->size() == 0) {
     ATH_MSG_DEBUG("ScintHitCollection found with zero length!");
     return StatusCode::SUCCESS;
   }
+  
+  // Create structure to store pulse for each channel
+  std::map<Identifier, std::vector<uint16_t>> waveforms;
+
+  auto first = *scintHitHandle->begin();
+  if (first.isVeto()) {
+    waveforms = m_digiTool->create_waveform_map(m_vetoID);
+  } else if (first.isTrigger()) {
+    waveforms = m_digiTool->create_waveform_map(m_triggerID);
+  } else if (first.isPreshower()) {
+    waveforms = m_digiTool->create_waveform_map(m_preshowerID);
+  } 
+
+  // 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 : *scintHitHandle) {            
+      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 = ScintHitIdHelper::GetHelper()->getIdentifier(c.first);
+      waveforms[id].push_back(value);
+    }
+  }
 
-  // Digitise the hits
-  CHECK( m_digiTool->digitise<ScintHitCollection>(scintHitHandle.ptr(),
-						  waveformContainerHandle.ptr(), m_kernel,
-						  std::pair<float, float>(m_base_mean, m_base_rms)) );
+  // Loop over wavefrom vectors to make and store raw 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);
+  }
 
   ATH_MSG_DEBUG("WaveformsHitContainer " << waveformContainerHandle.name() << "' filled with "<< waveformContainerHandle->size() <<" items");
 

Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h

@@ -19,12 +19,21 @@
 #include "GaudiKernel/ServiceHandle.h"
 #include "GaudiKernel/ToolHandle.h"
 
+//Helpers
+#include "ScintIdentifier/VetoID.h"
+#include "ScintIdentifier/TriggerID.h"
+#include "ScintIdentifier/PreshowerID.h"
+#include "ScintSimEvent/ScintHitIdHelper.h"
+#include "Identifier/Identifier.h"
+
+
 // ROOT
 #include "TF1.h"
 
 
 // STL
 #include <string>
+#include <vector>
 
 class ScintWaveformDigiAlg : public AthReentrantAlgorithm {
 
@@ -50,6 +59,10 @@ class ScintWaveformDigiAlg : public AthReentrantAlgorithm {
   ScintWaveformDigiAlg &operator=(const ScintWaveformDigiAlg&) = delete;
   //@}
 
+  /// 
+
+  /** @name Steerable pameters for crystal ball and baseline **/
+  //@{
   Gaudi::Property<double> m_CB_alpha {this, "CB_alpha", 0, "Alpha of the crystal ball function"};
   Gaudi::Property<double> m_CB_n {this, "CB_n", 0, "n of the crystal ball function"};
   Gaudi::Property<double> m_CB_mean {this, "CB_mean", 0, "Mean of the crystal ball function"};  
@@ -58,10 +71,20 @@ class ScintWaveformDigiAlg : public AthReentrantAlgorithm {
 
   Gaudi::Property<double> m_base_mean {this, "base_mean", 0, "Mean of the baseline"};
   Gaudi::Property<double> m_base_rms {this, "base_rms", 0, "RMS of the baseline"};
+  //@}
 
-  /// Kernel PDF
-  TF1* m_kernel;
 
+  /** Kernel PDF and evaluated values **/
+  //@{
+  TF1*                            m_kernel;
+  std::vector<float>              m_timekernel;
+  //@}
+  
+
+  /// Detector ID helpers
+  const VetoID* m_vetoID{nullptr};
+  const TriggerID* m_triggerID{nullptr};
+  const PreshowerID* m_preshowerID{nullptr};
 
   /**
    * @name Digitisation tool
@@ -91,4 +114,5 @@ class ScintWaveformDigiAlg : public AthReentrantAlgorithm {
 
 };
 
+
 #endif // SCINTDIGIALGS_SCINTDIGIALG_H

Scintillator/ScintSimEvent/CMakeLists.txt

@@ -18,11 +18,11 @@ atlas_add_library( ScintSimEvent
                    PRIVATE_INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${GEANT4_INCLUDE_DIRS}
                    DEFINITIONS ${CLHEP_DEFINITIONS}
                    LINK_LIBRARIES ${CLHEP_LIBRARIES} AthAllocators AthenaKernel CxxUtils GeneratorObjects HitManagement StoreGateLib SGtests
-                   PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES} ScintIdentifier )
+                   PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES} ScintIdentifier Identifier)
 
 atlas_add_dictionary( ScintSimEventDict
                       ScintSimEvent/ScintSimEventDict.h
                       ScintSimEvent/selection.xml
                       INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS} ${GEANT4_INCLUDE_DIRS}
-                      LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} AthAllocators CxxUtils GeneratorObjects HitManagement StoreGateLib SGtests ScintIdentifier ScintSimEvent )
+                      LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} AthAllocators CxxUtils GeneratorObjects HitManagement StoreGateLib SGtests ScintIdentifier ScintSimEvent Identifier)
 

Scintillator/ScintSimEvent/ScintSimEvent/ScintHitIdHelper.h

@@ -23,6 +23,15 @@
 
 // This class is singleton and static method and variable are used.
 #include "CxxUtils/checker_macros.h"
+
+//Helpers
+#include "ScintIdentifier/VetoID.h"
+#include "ScintIdentifier/TriggerID.h"
+#include "ScintIdentifier/PreshowerID.h"
+
+#include "Identifier/Identifier.h"
+
+
 ATLAS_NO_CHECK_FILE_THREAD_SAFETY;
 
 class ScintHitIdHelper : HitIdHelper {
@@ -43,6 +52,9 @@ class ScintHitIdHelper : HitIdHelper {
   // Layer/Disk
   int getPlate(const int& hid) const;
 
+  // identifier
+  Identifier getIdentifier(const int& hid) const;
+
   //
   // Info packing:
   int buildHitId(const int, const int, const int) const;
@@ -54,6 +66,12 @@ class ScintHitIdHelper : HitIdHelper {
   //
   // Initialize the helper, only called by the constructor
   void Initialize();
+
+  /// Detector ID helpers
+  const VetoID* m_vetoID{nullptr};
+  const TriggerID* m_triggerID{nullptr};
+  const PreshowerID* m_preshowerID{nullptr};
+
 };
 
 #endif // SCINTSIMEVENT_SCINTHITIDHELPER

Scintillator/ScintSimEvent/src/ScintHitIdHelper.cxx

@@ -6,7 +6,6 @@
 
 #include "ScintSimEvent/ScintHitIdHelper.h"
 #include "StoreGate/StoreGateSvc.h"
-#include "ScintIdentifier/VetoID.h"
 
 #include "G4Types.hh"
 #ifdef G4MULTITHREADED
@@ -44,7 +43,9 @@ void ScintHitIdHelper::Initialize() {
   const VetoID* pix;
   ServiceHandle<StoreGateSvc> detStore ("DetectorStore", "ScitHitIdHelper");
   if (detStore.retrieve().isSuccess()) {
-    if (detStore->retrieve(pix, "VetoID").isFailure()) { pix = 0; }
+    if (detStore->retrieve(m_vetoID, "VetoID").isFailure()) { m_vetoID = 0; }
+    if (detStore->retrieve(m_triggerID, "TriggerID").isFailure()) { m_triggerID = 0; }
+    if (detStore->retrieve(m_preshowerID, "PreshowerID").isFailure()) { m_preshowerID = 0; }
   }
 
   InitializeField("VetoTriggerPreshower", 0, 2);
@@ -87,6 +88,20 @@ int ScintHitIdHelper::getPlate(const int& hid) const
   return this->GetFieldValue("Plate", hid);
 }
 
+
+// identifier
+Identifier ScintHitIdHelper::getIdentifier(const int& hid) const
+{
+  if (isVeto(hid)) {
+    return m_vetoID->pmt_id(getStation(hid), getPlate(hid), 0);
+  } else if (isTrigger(hid)) {
+    return m_triggerID->pmt_id(getStation(hid), getPlate(hid), 0);
+  } else if (isPreshower(hid)) {
+    return m_preshowerID->pmt_id(getStation(hid), getPlate(hid), 0);
+  }
+  return Identifier();
+}
+
 //
 // Info packing:
 int ScintHitIdHelper::buildHitId(const int veto_trigger_preshower,

Waveform/WaveDigiTools/CMakeLists.txt

@@ -13,13 +13,12 @@ atlas_add_library( WaveDigiToolsLib
                    WaveDigiTools/*.h src/*.cxx src/*.h
                    PUBLIC_HEADERS WaveDigiTools
                    PRIVATE_INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
-                   LINK_LIBRARIES AthenaBaseComps AthenaKernel GeoPrimitives WaveRawEvent
+                   LINK_LIBRARIES AthenaBaseComps AthenaKernel GeoPrimitives WaveRawEvent Identifier 
                    PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES}
 		   )
 
 atlas_add_component( WaveDigiTools
 		     src/components/*.cxx 
 		     INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
-                     LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps GaudiKernel WaveDigiToolsLib )
-
+                     LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps GaudiKernel WaveDigiToolsLib)
 

Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h

@@ -23,10 +23,15 @@
 #include "WaveRawEvent/RawWaveformContainer.h"
 #include "WaveRawEvent/RawWaveform.h"
 
+#include "Identifier/Identifier.h"
+
 #include "TF1.h"
 #include "TRandom3.h"
 
 #include <utility>
+#include <map>
+#include <vector>
+
 
 ///Interface for waveform digitisation tools
 class IWaveformDigitisationTool : virtual public IAlgTool 
@@ -42,19 +47,25 @@ public:
 
   virtual ~IWaveformDigitisationTool() = default;
 
-  // Digitise HITS to Raw waveform
-  template<class CONT>
-  StatusCode digitise(const CONT* hitCollection, 
-		      RawWaveformContainer* waveContainer, 
-		      TF1* kernel, std::pair<float, float> base
-		      ) const;
+  /// Evaluate time kernel over time samples  
+  virtual std::vector<float> evaluate_timekernel(TF1* kernel) const = 0;
+
+  /// Generate random baseline 
+  virtual unsigned int generate_baseline(int mean, int rms) const = 0;
+ 
+  /// Create structure to store pulse for each channel
+  template <class T> 
+  std::map<Identifier, std::vector<uint16_t>> create_waveform_map(const T* idHelper) const;
+
+  /// Number of time samples
+  unsigned int nsamples() const { return m_nsamples; }
 
 private:
   ServiceHandle<IMessageSvc>      m_msgSvc;
 
 protected:
   TRandom3*                       m_random;
-
+  unsigned int                    m_nsamples;
 };
 
 #include "WaveDigiTools/IWaveformDigitisationTool.icc"

Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.icc

-#include <vector>
-#include <map>
+#include "Identifier/Identifier.h"
+#include "Identifier/ExpandedIdentifier.h"
 
-template<class CONT>
-StatusCode IWaveformDigitisationTool::digitise(const CONT* hitCollection,
-					       RawWaveformContainer* container, TF1* kernel,
-					       std::pair<float, float> base) const {
+template <class ID> 
+std::map<Identifier, std::vector<uint16_t>>  IWaveformDigitisationTool::create_waveform_map(const ID* idHelper) const {
 
+  std::map<Identifier, std::vector<uint16_t>> waveforms;
 
-  // Check the container
-  if (!container) {
-    MsgStream log(&(*m_msgSvc), name());
-    log << MSG::ERROR << "HitCollection passed to digitise() is null!" << endmsg;
-    return StatusCode::FAILURE;
+  for (auto itr = idHelper->pmt_begin(); itr != idHelper->pmt_end(); ++itr) {
+    const ExpandedIdentifier& extId = *itr;
+    Identifier id = idHelper->pmt_id(extId);
+    waveforms[id] = std::vector<uint16_t>();
+    waveforms[id].reserve(m_nsamples);
   }
 
-  unsigned int size = 600;  // TODO: how know the correct number of time samples?
-  std::vector<float> time(size);  
-  for (unsigned int i=0; i<size; i++) time[i] = 2.*i;
-
-  std::map<unsigned int, std::vector<uint16_t>> waveforms;
-  //unsigned int baseline = 800; 
-
-  // TODO: varying time to centre of edge of bin (odd = centre, even = edge)
-
-  // Loop over time samples
-  for (const auto& t : time) {
-    std::map<unsigned int, float> counts;
-   
-    unsigned int baseline = m_random->Gaus(base.first, base.second); 
-
-    // Convolve hit energy with kernel and sum for each ID (i.e. channel)
-    for (const auto& hit : *hitCollection) { 
-      counts[hit.identify()] += kernel->Eval(t) * hit.energyLoss();
-      //std::cout << "HIT " << hit.identify() << " @ " << t << ": " <<  kernel->Eval(t) << " " << hit.energyLoss() << " -> " << counts[hit.identify()] << std::endl;
-    }
-
-    // Add count to correct waveform vec
-    for (const auto& c : counts) {
-      int value = baseline - c.second;
-      if (value < 0) {
-	MsgStream log(&(*m_msgSvc), name());
-	log << MSG::WARNING << "Found pulse " << c.second << " larger than baseline " << c.first << endmsg;
-	value = 0; // Protect against scaling signal above baseline
-      }
-      waveforms[c.first].push_back(value);
-      //std::cout << "ADC " << c.first << " @ " << t << ": " << baseline - c.second << std::endl;
-    }
-  }
-
-  // Loop over wavefrom vecs to make and store waveform
-  for (const auto& w : waveforms) {
-    RawWaveform* wfm = new RawWaveform();
-    wfm->setWaveform(0, w.second);
-    wfm->setIdentifier(Identifier(w.first));
-    wfm->setSamples(size);
-    container->push_back(wfm);
-  }
-
-
-  return StatusCode::SUCCESS;
+  return waveforms;
 }

Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx

@@ -20,8 +20,28 @@ WaveformDigitisationTool::WaveformDigitisationTool(const std::string& type, cons
 StatusCode
 WaveformDigitisationTool::initialize() {
   ATH_MSG_INFO( name() << "::initalize()" );
+
+  m_nsamples = 600;
   m_random = new TRandom3();
+
   return StatusCode::SUCCESS;
 }
 
 
+std::vector<float>
+WaveformDigitisationTool::evaluate_timekernel(TF1* kernel) const {
+  
+  std::vector<float> timekernel;
+  timekernel.reserve(m_nsamples);
+  
+  for (unsigned int i=0; i<m_nsamples; i++) {
+    timekernel.push_back(kernel->Eval(2.*i));  
+  }
+  
+  return timekernel;
+}
+
+unsigned int 
+WaveformDigitisationTool::generate_baseline(int mean, int rms) const {
+  return m_random->Gaus(mean, rms);
+}

Waveform/WaveDigiTools/src/WaveformDigitisationTool.h

@@ -28,6 +28,13 @@ class WaveformDigitisationTool: public extends<AthAlgTool, IWaveformDigitisation
   /// Retrieve the necessary services in initialize
   StatusCode initialize();
 
+  /// Evaluate time kernel over samples  
+  std::vector<float> evaluate_timekernel(TF1* kernel) const;
+
+  /// Generate random baseline 
+  unsigned int generate_baseline(int mean, int rms) const;
+
+
  private:
   // None
 

Waveform/WaveRecAlgs/src/RawWaveformRecAlg.cxx

@@ -52,6 +52,13 @@ RawWaveformRecAlg::execute(const EventContext& ctx) const {
   ATH_CHECK( waveformHandle.isValid() );
   ATH_MSG_DEBUG("Found ReadHandle for RawWaveformContainer " << m_waveformContainerKey);
 
+  // Find the output waveform container
+  SG::WriteHandle<xAOD::WaveformHitContainer> hitContainerHandle(m_waveformHitContainerKey, ctx);
+  ATH_CHECK( hitContainerHandle.record( std::make_unique<xAOD::WaveformHitContainer>(),
+					std::make_unique<xAOD::WaveformHitAuxContainer>() ) );
+
+  ATH_MSG_DEBUG("WaveformsHitContainer '" << hitContainerHandle.name() << "' initialized");
+
   if (waveformHandle->size() == 0) {
     ATH_MSG_DEBUG("Waveform container found with zero length!");
     return StatusCode::SUCCESS;
@@ -69,12 +76,6 @@ RawWaveformRecAlg::execute(const EventContext& ctx) const {
     ATH_MSG_WARNING("Didn't find ReadHandle for WaveformClock!");
   }
 
-  // Find the output waveform container
-  SG::WriteHandle<xAOD::WaveformHitContainer> hitContainerHandle(m_waveformHitContainerKey, ctx);
-  ATH_CHECK( hitContainerHandle.record( std::make_unique<xAOD::WaveformHitContainer>(),
-					std::make_unique<xAOD::WaveformHitAuxContainer>() ) );
-
-  ATH_MSG_DEBUG("WaveformsHitContainer '" << hitContainerHandle.name() << "' initialized");
 
   // Reconstruct all waveforms
   CHECK( m_recoTool->reconstructAll(*waveformHandle, clockptr, hitContainerHandle.ptr()) );