Merge branch 'scint_calo_digi' into 'master'

Add digitized waveform reconstruction for simulated calo and scintillator hits

See merge request faser/calypso!195

Calorimeter/CaloDigiAlgs/CMakeLists.txt

+################################################################################
+# Package: CaloDigiAlgs
+################################################################################
+
+# Declare the package name:
+atlas_subdir( CaloDigiAlgs )
+
+# Component(s) in the package:
+atlas_add_component( CaloDigiAlgs
+                     src/*.cxx src/*.h
+                     src/components/*.cxx
+                     LINK_LIBRARIES AthenaBaseComps Identifier StoreGateLib WaveRawEvent FaserCaloSimEvent WaveDigiToolsLib)
+
+atlas_install_python_modules( python/*.py )
+

Calorimeter/CaloDigiAlgs/python/CaloDigiAlgsConfig.py

+# Copyright (C) 2020-2021 CERN for the benefit of the FASER collaboration
+
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+
+from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+
+# One stop shopping for normal FASER data
+def CaloWaveformDigitizationCfg(flags):
+    """ Return all algorithms and tools for Waveform digitization """
+    acc = ComponentAccumulator()
+
+    if not flags.Input.isMC:
+        return
+
+    acc.merge(CaloWaveformDigiCfg(flags, "CaloWaveformDigiAlg"))
+    return acc
+
+# Return configured digitization algorithm from SIM hits
+def CaloWaveformDigiCfg(flags, name="CaloWaveformDigiAlg", **kwargs):
+
+    acc = ComponentAccumulator()
+
+    tool = CompFactory.WaveformDigitisationTool(name="CaloWaveformDigtisationTool", **kwargs)
+    
+    kwargs.setdefault("CaloHitContainerKey", "EcalHits")
+    kwargs.setdefault("WaveformContainerKey", "CaloWaveforms")
+
+    digiAlg = CompFactory.CaloWaveformDigiAlg(name, **kwargs)
+    kwargs.setdefault("WaveformDigitisationTool", tool)
+    
+    digiAlg.CB_alpha = -0.9
+    digiAlg.CB_n = 10
+    digiAlg.CB_sigma = 4
+    digiAlg.CB_mean = 820
+
+
+    acc.addEventAlgo(digiAlg)
+
+    return acc
+
+def CaloWaveformDigitizationOutputCfg(flags, **kwargs):
+    """ Return ComponentAccumulator with output for Waveform Digi"""
+    acc = ComponentAccumulator()
+    ItemList = [
+        "RawWaveformContainer#*"
+    ]
+    acc.merge(OutputStreamCfg(flags, "RDO", ItemList))
+    # ostream = acc.getEventAlgo("OutputStreamRDO")
+    # ostream.TakeItemsFromInput = True # Don't know what this does
+    return acc

Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx

+#include "CaloWaveformDigiAlg.h"
+
+#include "Identifier/Identifier.h"
+
+#include <vector>
+#include <map>
+
+CaloWaveformDigiAlg::CaloWaveformDigiAlg(const std::string& name, 
+					 ISvcLocator* pSvcLocator)
+  : AthReentrantAlgorithm(name, pSvcLocator) { 
+
+}
+
+StatusCode 
+CaloWaveformDigiAlg::initialize() {
+  ATH_MSG_INFO(name() << "::initalize()" );
+
+  // 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);
+  
+  m_kernel = new TF1("PDF", "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);
+
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode 
+CaloWaveformDigiAlg::finalize() {
+  ATH_MSG_INFO(name() << "::finalize()");
+
+  delete m_kernel;
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode 
+CaloWaveformDigiAlg::execute(const EventContext& ctx) const {
+  ATH_MSG_DEBUG("Executing");
+
+  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);
+
+  ATH_CHECK( caloHitHandle.isValid() );
+  ATH_MSG_DEBUG("Found ReadHandle for CaloHitCollection " << m_caloHitContainerKey);
+
+  if (caloHitHandle->size() == 0) {
+    ATH_MSG_DEBUG("CaloHitCollection found with zero length!");
+    return StatusCode::SUCCESS;
+  }
+
+  // Find the output waveform container
+  SG::WriteHandle<RawWaveformContainer> waveformContainerHandle(m_waveformContainerKey, ctx);
+  ATH_CHECK( waveformContainerHandle.record( std::make_unique<RawWaveformContainer>()) );
+
+  ATH_MSG_DEBUG("WaveformsContainer '" << waveformContainerHandle.name() << "' initialized");
+
+  // Digitise the hits
+  CHECK( m_digiTool->digitise<CaloHitCollection>(caloHitHandle.ptr(),
+						  waveformContainerHandle.ptr(), m_kernel) );
+
+  ATH_MSG_DEBUG("WaveformsHitContainer " << waveformContainerHandle.name() << "' filled with "<< waveformContainerHandle->size() <<" items");
+
+  return StatusCode::SUCCESS;
+}

Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h

+#ifndef CALODIGIALGS_CALOWAVEFORMDIGIALG_H
+#define CALODIGIALGS_CALOWAVEFORMDIGIALG_H
+
+// Base class
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+
+// Data classes
+#include "WaveRawEvent/RawWaveformContainer.h"
+#include "FaserCaloSimEvent/CaloHitCollection.h"
+
+// Tool classes
+#include "WaveDigiTools/IWaveformDigitisationTool.h"
+
+// Handles
+#include "StoreGate/ReadHandleKey.h"
+#include "StoreGate/WriteHandleKey.h"
+
+// Gaudi
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/ToolHandle.h"
+
+// ROOT
+#include "TF1.h"
+
+// STL
+#include <string>
+
+class CaloWaveformDigiAlg : public AthReentrantAlgorithm {
+
+ public:
+  // Constructor
+  CaloWaveformDigiAlg(const std::string& name, ISvcLocator* pSvcLocator);
+  virtual ~CaloWaveformDigiAlg() = default;
+
+  /** @name Usual algorithm methods */
+  //@{
+  virtual StatusCode initialize() override;
+  virtual StatusCode execute(const EventContext& ctx) const override;
+  virtual StatusCode finalize() override;
+  //@}
+
+ private:
+
+  /** @name Disallow default instantiation, copy, assignment */
+  //@{
+  CaloWaveformDigiAlg() = delete;
+  CaloWaveformDigiAlg(const CaloWaveformDigiAlg&) = delete;
+  CaloWaveformDigiAlg &operator=(const CaloWaveformDigiAlg&) = delete;
+  //@}
+
+  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"};
+  Gaudi::Property<double> m_CB_sigma {this, "CB_sigma", 0, "Sigma of the crystal ball function"};
+
+
+
+
+
+  /// Kernel PDF
+  TF1* m_kernel;
+
+
+  /**
+   * @name Digitisation tool
+   */
+  ToolHandle<IWaveformDigitisationTool> m_digiTool
+    {this, "WaveformDigitisationTool", "WaveformDigitisationTool"};
+
+
+  /**
+   * @name Input HITS using SG::ReadHandleKey
+   */
+  //@{
+
+  SG::ReadHandleKey<CaloHitCollection> m_caloHitContainerKey 
+  {this, "CaloHitContainerKey", ""};
+
+  //@}
+
+
+  /**
+   * @name Output data using SG::WriteHandleKey
+   */
+  //@{
+  SG::WriteHandleKey<RawWaveformContainer> m_waveformContainerKey
+    {this, "WaveformContainerKey", ""};
+  //@}
+
+};
+
+#endif // CALODIGIALGS_CALODIGIALG_H

Calorimeter/CaloDigiAlgs/src/components/CaloDigiAlgs_entries.cxx

+#include "../CaloWaveformDigiAlg.h"
+DECLARE_COMPONENT( CaloWaveformDigiAlg )

Scintillator/ScintDigiAlgs/CMakeLists.txt

+################################################################################
+# Package: ScintDigiAlgs
+################################################################################
+
+# Declare the package name:
+atlas_subdir( ScintDigiAlgs )
+
+# Component(s) in the package:
+atlas_add_component( ScintDigiAlgs
+                     src/*.cxx src/*.h
+                     src/components/*.cxx
+                     LINK_LIBRARIES AthenaBaseComps Identifier StoreGateLib WaveRawEvent ScintSimEvent WaveDigiToolsLib)
+
+atlas_install_python_modules( python/*.py )
+

Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py

+# Copyright (C) 2020-2021 CERN for the benefit of the FASER collaboration
+
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+
+from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+# Crystallball parameter dictionary used in simulated digitized wave reconstruction.
+# Crystalball function Parameters estimated from Deion's slides uploaded at
+# https://indico.cern.ch/event/1099652/contributions/4626975/attachments/2352595/4013927/Faser-Physics-run3933-plots.pdf  (20/01/2022)
+# Parameters are per scintillator source, but not per channel.
+dict_CB_param = {}
+dict_CB_param["Trigger"]=dict(CB_alpha=-0.38, CB_n=25, CB_mean=815, CB_sigma=7.7)
+dict_CB_param["Timing"]=dict(CB_alpha=-0.32, CB_n=65, CB_mean=846, CB_sigma=5.3) # copy from Preshower; Timing was not in TestBeam
+dict_CB_param["Veto"]=dict(CB_alpha=-0.38, CB_n=25, CB_mean=815, CB_sigma=7.7) # 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)
+
+# One stop shopping for normal FASER data
+def ScintWaveformDigitizationCfg(flags):
+    """ Return all algorithms and tools for Waveform digitization """
+    acc = ComponentAccumulator()
+
+    if not flags.Input.isMC:
+        return
+
+    if "TB" not in flags.GeoModel.FaserVersion:
+        acc.merge(ScintWaveformDigiCfg(flags, "TimingWaveformDigiAlg", "Trigger"))
+    acc.merge(ScintWaveformDigiCfg(flags, "VetoWaveformDigiAlg", "Veto"))
+    acc.merge(ScintWaveformDigiCfg(flags, "PreshowerWaveformDigiAlg", "Preshower"))
+    return acc
+
+# Return configured digitization algorithm from SIM hits
+# Specify data source (Veto, Trigger, Preshower)
+def ScintWaveformDigiCfg(flags, name="ScintWaveformDigiAlg", source="", **kwargs):
+
+    acc = ComponentAccumulator()
+
+    tool = CompFactory.WaveformDigitisationTool(name=source+"WaveformDigtisationTool", **kwargs)
+    
+    kwargs.setdefault("ScintHitContainerKey", source+"Hits")
+    kwargs.setdefault("WaveformContainerKey", source+"Waveforms")
+
+    digiAlg = CompFactory.ScintWaveformDigiAlg(name, **kwargs)
+    digiAlg.CB_alpha = dict_CB_param[source]["CB_alpha"]
+    digiAlg.CB_n = dict_CB_param[source]["CB_n"]
+    digiAlg.CB_mean = dict_CB_param[source]["CB_mean"]
+    digiAlg.CB_sigma = dict_CB_param[source]["CB_sigma"]
+    
+    kwargs.setdefault("WaveformDigitisationTool", tool)
+
+    acc.addEventAlgo(digiAlg)
+
+    return acc
+
+def ScintWaveformDigitizationOutputCfg(flags, **kwargs):
+    """ Return ComponentAccumulator with output for Waveform Digi"""
+    acc = ComponentAccumulator()
+    ItemList = [
+        "RawWaveformContainer#*"
+    ]
+    acc.merge(OutputStreamCfg(flags, "RDO", ItemList))
+    # ostream = acc.getEventAlgo("OutputStreamRDO")
+    # ostream.TakeItemsFromInput = True # Don't know what this does
+    return acc

Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx

+
+#include "ScintWaveformDigiAlg.h"
+
+#include "Identifier/Identifier.h"
+
+#include <vector>
+#include <map>
+
+
+ScintWaveformDigiAlg::ScintWaveformDigiAlg(const std::string& name, 
+					 ISvcLocator* pSvcLocator)
+  : AthReentrantAlgorithm(name, pSvcLocator) { 
+
+}
+
+StatusCode 
+ScintWaveformDigiAlg::initialize() {
+  ATH_MSG_INFO(name() << "::initalize()" );
+
+  // 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
+ 
+
+  
+  m_kernel = new TF1("PDF", "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);
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode 
+ScintWaveformDigiAlg::finalize() {
+  ATH_MSG_INFO(name() << "::finalize()");
+
+  delete m_kernel;
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode 
+ScintWaveformDigiAlg::execute(const EventContext& ctx) const {
+  ATH_MSG_DEBUG("Executing");
+
+  ATH_MSG_DEBUG("Run: " << ctx.eventID().run_number() 
+		<< " Event: " << ctx.eventID().event_number());
+
+  // Find the input HIT collection
+  SG::ReadHandle<ScintHitCollection> scintHitHandle(m_scintHitContainerKey, ctx);
+
+  ATH_CHECK( scintHitHandle.isValid() );
+  ATH_MSG_DEBUG("Found ReadHandle for ScintHitCollection " << m_scintHitContainerKey);
+
+  if (scintHitHandle->size() == 0) {
+    ATH_MSG_DEBUG("ScintHitCollection found with zero length!");
+    return StatusCode::SUCCESS;
+  }
+
+  // Find the output waveform container
+  SG::WriteHandle<RawWaveformContainer> waveformContainerHandle(m_waveformContainerKey, ctx);
+  ATH_CHECK( waveformContainerHandle.record( std::make_unique<RawWaveformContainer>()) );
+
+  ATH_MSG_DEBUG("WaveformsContainer '" << waveformContainerHandle.name() << "' initialized");
+
+  // Digitise the hits
+  CHECK( m_digiTool->digitise<ScintHitCollection>(scintHitHandle.ptr(),
+						  waveformContainerHandle.ptr(), m_kernel) );
+
+  ATH_MSG_DEBUG("WaveformsHitContainer " << waveformContainerHandle.name() << "' filled with "<< waveformContainerHandle->size() <<" items");
+
+  return StatusCode::SUCCESS;
+}

Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h

+#ifndef SCINTDIGIALGS_SCINTWAVEFORMDIGIALG_H
+#define SCINTDIGIALGS_SCINTWAVEFORMDIGIALG_H
+
+// Base class
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+
+// Data classes
+#include "WaveRawEvent/RawWaveformContainer.h"
+#include "ScintSimEvent/ScintHitCollection.h"
+
+// Tool classes
+#include "WaveDigiTools/IWaveformDigitisationTool.h"
+
+// Handles
+#include "StoreGate/ReadHandleKey.h"
+#include "StoreGate/WriteHandleKey.h"
+
+// Gaudi
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/ToolHandle.h"
+
+// ROOT
+#include "TF1.h"
+
+
+// STL
+#include <string>
+
+class ScintWaveformDigiAlg : public AthReentrantAlgorithm {
+
+ public:
+  // Constructor
+  ScintWaveformDigiAlg(const std::string& name, ISvcLocator* pSvcLocator);
+  virtual ~ScintWaveformDigiAlg() = default;
+
+  /** @name Usual algorithm methods */
+  //@{
+  virtual StatusCode initialize() override;
+  virtual StatusCode execute(const EventContext& ctx) const override;
+  virtual StatusCode finalize() override;
+  //@}
+
+
+ private:
+
+  /** @name Disallow default instantiation, copy, assignment */
+  //@{
+  ScintWaveformDigiAlg() = delete;
+  ScintWaveformDigiAlg(const ScintWaveformDigiAlg&) = delete;
+  ScintWaveformDigiAlg &operator=(const ScintWaveformDigiAlg&) = delete;
+  //@}
+
+  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"};  
+  Gaudi::Property<double> m_CB_sigma {this, "CB_sigma", 0, "Sigma of the crystal ball function"};
+
+
+  /// Kernel PDF
+  TF1* m_kernel;
+
+
+  /**
+   * @name Digitisation tool
+   */
+  ToolHandle<IWaveformDigitisationTool> m_digiTool
+    {this, "WaveformDigitisationTool", "WaveformDigitisationTool"};
+
+
+  /**
+   * @name Input HITS using SG::ReadHandleKey
+   */
+  //@{
+
+  SG::ReadHandleKey<ScintHitCollection> m_scintHitContainerKey 
+  {this, "ScintHitContainerKey", ""};
+
+  //@}
+
+
+  /**
+   * @name Output data using SG::WriteHandleKey
+   */
+  //@{
+  SG::WriteHandleKey<RawWaveformContainer> m_waveformContainerKey
+    {this, "WaveformContainerKey", ""};
+  //@}
+
+};
+
+#endif // SCINTDIGIALGS_SCINTDIGIALG_H

Scintillator/ScintDigiAlgs/src/components/ScintDigiAlgs_entries.cxx

+#include "../ScintWaveformDigiAlg.h"
+DECLARE_COMPONENT( ScintWaveformDigiAlg )

Waveform/WaveDigiTools/CMakeLists.txt

+################################################################################
+# Package: WaveDigiTools
+################################################################################
+
+# Declare the package name:
+atlas_subdir( WaveDigiTools )
+
+# External dependencies:
+find_package( ROOT )
+
+# Component(s) in the package:
+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
+                   PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES}
+		   )
+
+atlas_add_component( WaveDigiTools
+		     src/components/*.cxx 
+		     INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
+                     LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps GaudiKernel WaveDigiToolsLib )
+
+

Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h

+/*
+  Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+*/
+
+/**
+ * @file IWaveformDigitisationTool.h
+ * Header file for the IWaveformDigitisationTool class
+ * @author Carl Gwilliam, 2021
+ */
+
+
+#ifndef WAVEDIGITOOLS_IWAVEFORMDIGITISATIONTOOL_H
+#define WAVEDIGITOOLS_IWAVEFORMDIGITISATIONTOOL_H
+
+// Base class
+#include "GaudiKernel/IAlgTool.h"
+#include "GaudiKernel/ToolHandle.h"
+
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/IMessageSvc.h"
+#include "GaudiKernel/MsgStream.h"
+
+#include "WaveRawEvent/RawWaveformContainer.h"
+#include "WaveRawEvent/RawWaveform.h"
+
+#include "TF1.h"
+
+///Interface for waveform digitisation tools
+class IWaveformDigitisationTool : virtual public IAlgTool 
+{
+public:
+
+  // InterfaceID
+  DeclareInterfaceID(IWaveformDigitisationTool, 1, 0);
+
+  IWaveformDigitisationTool():
+    m_msgSvc         ( "MessageSvc",   "ITrkEventCnvTool" )
+  {}
+
+  virtual ~IWaveformDigitisationTool() = default;
+
+  // Digitise HITS to Raw waveform
+  template<class CONT>
+  StatusCode digitise(const CONT* hitCollection, 
+		      RawWaveformContainer* waveContainer, TF1* kernel) const;
+
+private:
+  ServiceHandle<IMessageSvc>      m_msgSvc;
+
+};
+
+#include "WaveDigiTools/IWaveformDigitisationTool.icc"
+
+
+#endif //WAVEDIGITOOLS_IWAVEFORMDIGITISATIONTOOL_H

Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.icc

+#include <vector>
+#include <map>
+
+template<class CONT>
+StatusCode IWaveformDigitisationTool::digitise(const CONT* hitCollection,
+					       RawWaveformContainer* container, TF1* kernel) const {
+
+
+  // Check the container
+  if (!container) {
+    MsgStream log(&(*m_msgSvc), name());
+    log << MSG::ERROR << "HitCollection passed to digitise() is null!" << endmsg;
+    return StatusCode::FAILURE;
+  }
+
+  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 = 8000; // TODO: vary this + add noise
+
+  // Loop over time samples
+  for (const auto& t : time) {
+    std::map<unsigned int, float> counts;
+   
+    // 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) {
+      waveforms[c.first].push_back(baseline - c.second);
+      //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;
+}

Waveform/WaveDigiTools/share/WaveformDigiAndRecoExample_jobOptions.py

+#!/usr/bin/env python
+
+import sys
+
+if __name__ == "__main__":
+
+    fileroot = "output"
+    if len(sys.argv) > 1:
+       filename = sys.argv[1]
+
+    doRDO = False
+    if len(sys.argv) > 2:
+       filename = sys.argv[2]
+
+    from AthenaCommon.Logging import log, logging
+    from AthenaCommon.Constants import DEBUG, VERBOSE, INFO
+    from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+    from AthenaCommon.Configurable import Configurable
+        
+    Configurable.configurableRun3Behavior = True
+
+    log.setLevel(VERBOSE)
+    
+    ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02"             # Always needed; must match FaserVersion
+    ConfigFlags.IOVDb.DatabaseInstance = "OFLP200"               # Use MC conditions for now
+    ConfigFlags.Input.ProjectName = "mc21"                       # Needed to bypass autoconfig
+    ConfigFlags.Input.isMC = True                                # Needed to bypass autoconfig
+    ConfigFlags.GeoModel.FaserVersion     = "FASER-TB00"         # FASER geometry
+    ConfigFlags.Common.isOnline = False
+    ConfigFlags.GeoModel.Align.Dynamic = False
+    
+    ConfigFlags.Input.Files = [
+        "/eos/project-f/faser-commissioning/Simulation/Test/TB.Elec.200GeV.SIM.root"
+        ]
+
+    if doRDO:
+        ConfigFlags.Output.RDOFileName = f"{fileroot}.RDO.root"
+    else:
+        ConfigFlags.addFlag("Output.xAODFileName", f"{fileroot}.xAOD.root")
+        ConfigFlags.Output.ESDFileName = f"{fileroot}.ESD.root"
+
+    ConfigFlags.lock()
+
+    from CalypsoConfiguration.MainServicesConfig import MainServicesCfg
+    acc = MainServicesCfg(ConfigFlags)
+
+    from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+    from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg
+
+    acc.merge(PoolReadCfg(ConfigFlags))
+    acc.merge(PoolWriteCfg(ConfigFlags))
+
+    if doRDO:
+         from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+         itemList = [
+             "RawWaveformContainer#*"
+             ]
+         acc.merge(OutputStreamCfg(ConfigFlags, "RDO", itemList,disableEventTag=True))
+
+
+    else:
+        from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+        itemList = [
+            "xAOD::EventInfo#*",
+            "xAOD::WaveformHitContainer#*",
+            "xAOD::WaveformHitAuxContainer#*",
+            ]
+
+        acc.merge(OutputStreamCfg(ConfigFlags, "xAOD", itemList, disableEventTag=True))
+
+    from ScintDigiAlgs.ScintDigiAlgsConfig import ScintWaveformDigitizationCfg
+    acc.merge(ScintWaveformDigitizationCfg(ConfigFlags))
+
+    from CaloDigiAlgs.CaloDigiAlgsConfig import CaloWaveformDigitizationCfg
+    acc.merge(CaloWaveformDigitizationCfg(ConfigFlags))
+
+    if not doRDO:
+        from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionCfg
+        acc.merge(WaveformReconstructionCfg(ConfigFlags))
+
+    #acc.foreach_component("*").OutputLevel = VERBOSE
+
+    # Execute and finish
+    sc = acc.run(maxEvents=100)
+
+    # Success should be 0
+    sys.exit(not sc.isSuccess())

Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx

+/*
+  Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+*/
+
+/**
+ * @file WaveformDigitisationTool.cxx
+ * Implementation file for the WaveformDigitisationTool class
+ * @ author C. Gwilliam, 2021
+ **/
+
+#include "WaveformDigitisationTool.h"
+
+// Constructor
+WaveformDigitisationTool::WaveformDigitisationTool(const std::string& type, const std::string& name, const IInterface* parent) :
+  base_class(type, name, parent)
+{
+}
+
+// Initialization
+StatusCode
+WaveformDigitisationTool::initialize() {
+  ATH_MSG_INFO( name() << "::initalize()" );
+  return StatusCode::SUCCESS;
+}
+
+

Waveform/WaveDigiTools/src/WaveformDigitisationTool.h

+/*
+   Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+*/
+
+/** @file WaveformDigitisationTool.h
+ *  Header file for WaveformDigitisationTool.h
+ *
+ */
+#ifndef WAVEDIGITOOLS_WAVEFORMDIGITISATIONTOOL_H
+#define WAVEDIGITOOLS_WAVEFORMDIGITISATIONTOOL_H
+
+//Athena
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "WaveDigiTools/IWaveformDigitisationTool.h"
+
+//Gaudi
+#include "GaudiKernel/ToolHandle.h"
+
+//STL
+
+class WaveformDigitisationTool: public extends<AthAlgTool, IWaveformDigitisationTool> {
+ public:
+
+  /// Normal constructor for an AlgTool; 'properties' are also declared here
+ WaveformDigitisationTool(const std::string& type, 
+			  const std::string& name, const IInterface* parent);
+
+  /// Retrieve the necessary services in initialize
+  StatusCode initialize();
+
+ private:
+  // None
+
+};
+
+#endif // WAVEDIGITOOLS_WAVEFORMDIGITISATIONTOOL_H

Waveform/WaveDigiTools/src/components/WaveDigiTools_entries.cxx

+#include "../WaveformDigitisationTool.h"
+
+DECLARE_COMPONENT( WaveformDigitisationTool )

Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py

@@ -11,24 +11,27 @@ WaveformReconstructionTool = CompFactory.WaveformReconstructionTool
 ClockReconstructionTool = CompFactory.ClockReconstructionTool
 
 # One stop shopping for normal FASER data
-def WaveformReconstructionCfg(flags, naive = True):
+def WaveformReconstructionCfg(flags, naive = False):
     """ Return all algorithms and tools for Waveform reconstruction """
     acc = ComponentAccumulator()
 
+    if not flags.Input.isMC:
+        acc.merge(WaveformClockRecCfg(flags, "ClockRecAlg"))
+
     if flags.Input.isMC  and naive:
-        if not "TB" in flags.GeoModel.FaserVersion:
+        if "TB" not in flags.GeoModel.FaserVersion:
             acc.merge(PseudoScintHitToWaveformRecCfg(flags, "PseudoTimingHitWaveformRecAlg", "Trigger"))
         acc.merge(PseudoScintHitToWaveformRecCfg(flags, "PseudoVetoHitToWaveformRecAlg", "Veto"))
         acc.merge(PseudoScintHitToWaveformRecCfg(flags, "PseudoPresehowerHitWaveformRecAlg", "Preshower"))
         acc.merge(PseudoCaloHitToWaveformRecCfg(flags, "PseudoCaloHitWaveformRecAlg"))                
         return acc
 
-    acc.merge(WaveformClockRecCfg(flags, "ClockRecAlg"))
-
+    if "TB" not in flags.GeoModel.FaserVersion:
+        acc.merge(WaveformHitRecCfg(flags, "TimingWaveformRecAlg", "Trigger"))
     acc.merge(WaveformHitRecCfg(flags, "VetoWaveformRecAlg", "Veto"))
-    acc.merge(WaveformHitRecCfg(flags, "TimingWaveformRecAlg", "Trigger"))
     acc.merge(WaveformHitRecCfg(flags, "PreshowerWaveformRecAlg", "Preshower"))
     acc.merge(WaveformHitRecCfg(flags, "CaloWaveformRecAlg", "Calo"))
+
     return acc
 
 # Return configured WaveformClock reconstruction algorithm
@@ -52,7 +55,13 @@ def WaveformHitRecCfg(flags, name="WaveformRecAlg", source="", **kwargs):
 
     acc = ComponentAccumulator()
 
+    if flags.Input.isMC:
+        kwargs.setdefault("PeakThreshold", 5)
+
     tool = WaveformReconstructionTool(name=source+"WaveformRecTool", **kwargs)
+
+    if "PeakThreshold" in kwargs: kwargs.pop("PeakThreshold")
+
     kwargs.setdefault("WaveformContainerKey", source+"Waveforms")
     kwargs.setdefault("WaveformHitContainerKey", source+"WaveformHits")
     kwargs.setdefault("WaveformReconstructionTool", tool)