First commit of scintillator reconstruction code.

Code doesn't actually do anything, but this compiles and give the basic framework
for adding real algorithims.

Scintillator/ScintEventCnv/ScintByteStream/src/ScintWaveformDecoderTool.cxx

@@ -95,6 +95,7 @@ ScintWaveformDecoderTool::convert(const DAQFormats::EventFull* re,
     ATH_MSG_DEBUG("Fragment:\n" << *frag);
 
     digitizer = new DigitizerDataFragment(frag->payload<const uint32_t*>(), frag->payload_size()); 
+
     break;
   }
 
@@ -103,7 +104,12 @@ ScintWaveformDecoderTool::convert(const DAQFormats::EventFull* re,
     return StatusCode::FAILURE;
   }
 
-  // ATH_MSG_DEBUG("Digitizer Fragment:\n" << *digitizer);
+  // Check validity here
+  if (!digitizer->valid()) {
+    ATH_MSG_WARNING("Found invalid digitizer fragment:\n" << *digitizer);
+  } else {
+    ATH_MSG_DEBUG("Found valid digitizer fragment");
+  }
 
   std::vector<unsigned int>* channelList;
 
@@ -127,13 +133,18 @@ ScintWaveformDecoderTool::convert(const DAQFormats::EventFull* re,
   for (unsigned int channel: *channelList) {
     ATH_MSG_DEBUG("Converting channel "+std::to_string(channel)+" for "+key);
 
+    // Check if this has data
+    if (!digitizer->channel_has_data(channel)) {
+      ATH_MSG_INFO("Channel " << channel << " has no data - skipping!");
+      continue;
+    } 
+
     ScintWaveform* wfm = new ScintWaveform();
 
-    // Finally, set values in Waveform object from Digitizer fragment
     try {
       wfm->setWaveform( channel, digitizer->channel_adc_counts( channel ) );
     } catch ( DigitizerDataException& e ) {
-      ATH_MSG_INFO("ScintWaveformDecoderTool:\n"
+      ATH_MSG_WARNING("ScintWaveformDecoderTool:\n"
 		   <<e.what()
 		   << "\nChannel "
 		   << channel
@@ -151,12 +162,19 @@ ScintWaveformDecoderTool::convert(const DAQFormats::EventFull* re,
     }
 
     container->push_back(wfm);    
+
+    // Sanity check
+    if (wfm->adc_counts().size() != wfm->n_samples()) {
+      ATH_MSG_WARNING("Created waveform channel" << channel << "with length " << wfm->adc_counts().size() << " but header reports n_samples = " << wfm->n_samples() << "!");
+      ATH_MSG_WARNING(*wfm);
+    }
+
   }
 
   // Don't spring a leak
   delete digitizer;
 
-ATH_MSG_DEBUG( "ScintWaveformDecoderTool created container " << key 
-<< " with size=" << container->size());
+  ATH_MSG_DEBUG( "ScintWaveformDecoderTool created container " << key 
+		 << " with size=" << container->size());
   return StatusCode::SUCCESS; 
 }

Scintillator/ScintRawEvent/src/ScintWaveform.cxx

@@ -52,14 +52,15 @@ ScintWaveform::setWaveform(unsigned int channel, const std::vector<uint16_t> wav
 
 std::ostream 
 &operator<<(std::ostream &out, const ScintWaveform &wfm) {
-  out << "Waveform data:" << std::endl
-  << std::setw(30) << " board_id:             "<<std::setfill(' ')<<std::setw(32)<<std::dec<<wfm.board_id()<<std::setfill(' ')<<std::endl
-  << std::setw(30) << " board_fail_flag:      "<<std::setfill(' ')<<std::setw(32)<<std::dec<<wfm.board_fail_flag()<<std::setfill(' ')<<std::endl
-  << std::setw(30) << " board_trig_options:   "<<std::setfill(' ')<<std::setw(32)<<std::dec<<wfm.pattern_trig_options()<<std::setfill(' ')<<std::endl
-  << std::setw(30) << " channel_mask:         "<<std::setfill(' ')<<std::setw(32)<<std::dec<<wfm.channel_mask()<<std::setfill(' ')<<std::endl
-  << std::setw(30) << " event_counter:        "<<std::setfill(' ')<<std::setw(32)<<std::dec<<wfm.event_counter()<<std::setfill(' ')<<std::endl
-  << std::setw(30) << " trigger_time_tag:     "<<std::setfill(' ')<<std::setw(32)<<std::dec<<wfm.trigger_time_tag()<<std::setfill(' ')<<std::endl
-      << std::setw(30) << " n_samples:            "<<std::setfill(' ')<<std::setw(32)<<std::dec<<wfm.n_samples()<<std::setfill(' ')<<std::endl;
+  out << "Waveform data channel:" << std::dec << wfm.channel() << std::endl;
+  out << std::setw(30) << " board_id:             "<<std::setfill(' ')<<std::setw(32)<<std::dec<<wfm.board_id()<<std::setfill(' ')<<std::endl;
+  out << std::setw(30) << " board_fail_flag:      "<<std::setfill(' ')<<std::setw(32)<<std::dec<<wfm.board_fail_flag()<<std::setfill(' ')<<std::endl;
+  out << std::setw(30) << " board_trig_options:   "<<std::setfill(' ')<<std::setw(32)<<std::dec<<wfm.pattern_trig_options()<<std::setfill(' ')<<std::endl;
+  out << std::setw(30) << " channel_mask:         "<<std::setfill(' ')<<std::setw(32)<<std::hex<<wfm.channel_mask()<<std::setfill(' ')<<std::endl;
+  out << std::setw(30) << " event_counter:        "<<std::setfill(' ')<<std::setw(32)<<std::dec<<wfm.event_counter()<<std::setfill(' ')<<std::endl;
+  out << std::setw(30) << " trigger_time_tag:     "<<std::setfill(' ')<<std::setw(32)<<std::dec<<wfm.trigger_time_tag()<<std::setfill(' ')<<std::endl;
+  out << std::setw(30) << " n_samples:            "<<std::setfill(' ')<<std::setw(32)<<std::dec<<wfm.n_samples()<<std::setfill(' ')<<std::endl;
+  out << std::setw(30) << " adc_counts.size():    "<<std::setfill(' ')<<std::setw(32)<<std::dec<<wfm.adc_counts().size()<<std::setfill(' ')<<std::endl;
 
   return out;
 }

Scintillator/ScintRecAlgs/CMakeLists.txt

+################################################################################
+# Package: ScintRecAlgs
+################################################################################
+
+# Declare the package name:
+atlas_subdir( ScintRecAlgs )
+
+# Component(s) in the package:
+atlas_add_component( ScintRecAlgs
+                     src/*.cxx src/*.h
+                     src/components/*.cxx
+                     LINK_LIBRARIES AthenaBaseComps StoreGateLib ScintRawEvent xAODFaserWaveform ScintRecToolsLib )
+
+atlas_install_python_modules( python/*.py )
+

Scintillator/ScintRecAlgs/python/ScintRecAlgsConfig.py

+""" Define methods used to instantiate configured Waveform reconstruction tools and algorithms
+
+Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+"""
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+
+from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
+
+WaveformReconstructionTool = CompFactory.WaveformReconstructionTool
+ClockReconstructionTool = CompFactory.ClockReconstructionTool
+
+# One stop shopping for normal FASER data
+def WaveformReconstructionCfg(flags):
+    """ Return all algorithms and tools for Waveform reconstruction """
+    acc = ComponentAccumulator()
+
+    acc.merge(WaveformClockRecCfg(flags, "ClockRecAlg"))
+
+    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
+def WaveformClockRecCfg(flags, name="ClockRecAlg", **kwargs):
+
+    acc = ComponentAccumulator()
+
+    tool = ClockReconstructionTool(name="ClockReconstructionTool")
+
+    kwargs.setdefault("ClockReconstructionTool", tool)
+
+    recoAlg = CompFactory.ScintClockRecAlg(name, **kwargs)
+    recoAlg.ClockReconstructionTool = tool
+    acc.addEventAlgo(recoAlg)
+
+    return acc
+
+# Return configured WaveformHit reconstruction algorithm
+# Specify data source (Veto, Trigger, Preshower, Calo, Test)
+def WaveformHitRecCfg(flags, name="WaveformRecAlg", source="", **kwargs):
+
+    acc = ComponentAccumulator()
+
+    tool = WaveformReconstructionTool(name=source+"WaveformRecTool", **kwargs)
+    kwargs.setdefault("WaveformContainerKey", source+"Waveforms")
+    kwargs.setdefault("WaveformHitContainerKey", source+"WaveformHits")
+    kwargs.setdefault("WaveformReconstructionTool", tool)
+              
+    recoAlg = CompFactory.ScintWaveformRecAlg(name, **kwargs)
+    recoAlg.WaveformReconstructionTool = tool
+    acc.addEventAlgo(recoAlg)
+
+    return acc
+
+def WaveformReconstructionOutputCfg(flags, **kwargs):
+    """ Return ComponentAccumulator with output for Waveform Reco"""
+    acc = ComponentAccumulator()
+    ItemList = [
+        "xAOD::WaveformHitContainer#*"
+        , "xAOD::WaveformHitAuxContainer#*"
+        , "xAOD::WaveformClock#*"
+        , "xAOD::WaveformClockAuxInfo#*"
+    ]
+    acc.merge(OutputStreamCfg(flags, "xAOD", ItemList))
+    # ostream = acc.getEventAlgo("OutputStreamRDO")
+    # ostream.TakeItemsFromInput = True # Don't know what this does
+    return acc

Scintillator/ScintRecAlgs/src/ScintClockRecAlg.cxx

+#include "ScintClockRecAlg.h"
+
+ScintClockRecAlg::ScintClockRecAlg(const std::string& name, 
+				   ISvcLocator* pSvcLocator)
+  : AthReentrantAlgorithm(name, pSvcLocator) { 
+}
+
+StatusCode 
+ScintClockRecAlg::initialize() {
+  ATH_MSG_INFO(name() << "::initalize()" );
+
+  // Initalize tools
+  ATH_CHECK( m_recoTool.retrieve() );
+
+  // Set key to read waveform from
+  ATH_CHECK( m_waveformContainerKey.initialize() );
+
+  // Set key to write container
+  ATH_CHECK( m_waveformClockKey.initialize() );
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode 
+ScintClockRecAlg::finalize() {
+  ATH_MSG_INFO(name() << "::finalize()");
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode 
+ScintClockRecAlg::execute(const EventContext& ctx) const {
+  ATH_MSG_INFO("Executing");
+
+  ATH_MSG_DEBUG("Run: " << ctx.eventID().run_number() 
+		<< " Event: " << ctx.eventID().event_number());
+
+  // Find the input waveform container
+  SG::ReadHandle<ScintWaveformContainer> waveformHandle(m_waveformContainerKey, ctx);
+
+  ATH_CHECK( waveformHandle.isValid() );
+  ATH_MSG_DEBUG("Found ReadHandle for Waveforms");
+
+  if (waveformHandle->size() == 0) {
+    ATH_MSG_INFO("Waveform container found with zero length!");
+    return StatusCode::SUCCESS;
+  }
+
+  // Check length (should be one)
+  if (waveformHandle->size() != 1) {
+    ATH_MSG_WARNING("Found waveform container " << m_waveformContainerKey 
+		    << " with size " << waveformHandle->size() << "!");
+  }
+
+  // Create the output clock container
+  SG::WriteHandle<xAOD::WaveformClock> clockHandle(m_waveformClockKey, ctx);
+  //ATH_CHECK( clockHandle.record( std::make_unique<WaveformClock>() ) );
+  ATH_CHECK( clockHandle.record( std::make_unique<xAOD::WaveformClock>(),
+				 std::make_unique<xAOD::WaveformClockAuxInfo>() ) );
+
+  auto clock = clockHandle.ptr();
+
+  // Reconstruct first element
+  for(const auto& wave : *waveformHandle) {
+
+    ATH_MSG_DEBUG("Reconstruct waveform for channel " << wave->channel());
+
+    // Make some sanity checks before going further
+    if (wave->adc_counts().size() == 0) {
+      ATH_MSG_ERROR( "Found waveform for channel " << wave->channel() 
+		     << " with size " << wave->adc_counts().size() << "!");
+      // Create dummy hit here, or just skip?
+      return StatusCode::FAILURE;
+    }
+    
+    if (wave->adc_counts().size() != wave->n_samples()) {
+      ATH_MSG_WARNING( "Found waveform for channel " << wave->channel() 
+		       << " with size " << wave->adc_counts().size() 
+		       << " not equal to number of samples " << wave->n_samples());
+      return StatusCode::SUCCESS;
+    }
+
+    // Reconstruct the hits
+    CHECK( m_recoTool->reconstruct(*wave, clock) );
+
+    // Only do one if there happen to be more
+    break;
+  }
+
+  ATH_MSG_DEBUG(waveformHandle.name() << " created");
+
+  return StatusCode::SUCCESS;
+}
+

Scintillator/ScintRecAlgs/src/ScintClockRecAlg.h

+#ifndef SCINTRECALGS_SCINTCLOCKRECALG_H
+#define SCINTRECALGS_SCINTCLOCKRECALG_H
+
+// Base class
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+
+// Input data
+#include "ScintRawEvent/ScintWaveformContainer.h"
+
+// Output data
+// #include "ScintRecEvent/WaveformClock.h"
+#include "xAODFaserWaveform/WaveformClock.h"
+#include "xAODFaserWaveform/WaveformClockAuxInfo.h"
+#include "ScintRecTools/IClockReconstructionTool.h"
+
+#include "StoreGate/ReadHandleKey.h"
+
+// Gaudi
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/ToolHandle.h"
+
+// STL
+#include <string>
+
+class ScintClockRecAlg : public AthReentrantAlgorithm {
+
+ public:
+  // Constructor
+  ScintClockRecAlg(const std::string& name, ISvcLocator* pSvcLocator);
+  virtual ~ScintClockRecAlg() = 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 */
+  //@{
+  ScintClockRecAlg() = delete;
+  ScintClockRecAlg(const ScintClockRecAlg&) = delete;
+  ScintClockRecAlg &operator=(const ScintClockRecAlg&) = delete;
+  //@}
+
+  /**
+   * @name Reconstruction tool
+   */
+  ToolHandle<IClockReconstructionTool> m_recoTool
+    {this, "ClockReconstructionTool", "ClockReconstructionTool"};
+
+  /**
+   * @name Input data using SG::ReadHandleKey
+   */
+  //@{
+  SG::ReadHandleKey<ScintWaveformContainer> m_waveformContainerKey
+    {this, "WaveformContainerKey", "ClockWaveforms"};
+  //@}
+
+  /**
+   * @name Output data using SG::WriteHandleKey
+   */
+  //@{
+  SG::WriteHandleKey<xAOD::WaveformClock> m_waveformClockKey
+    {this, "WaveformClockKey", "WaveformClock"};
+  //@}
+
+};
+
+#endif // SCINTRECALGS_SCINTCLOCKRECALG_H

Scintillator/ScintRecAlgs/src/ScintWaveformRecAlg.cxx

+#include "ScintWaveformRecAlg.h"
+
+ScintWaveformRecAlg::ScintWaveformRecAlg(const std::string& name, 
+					 ISvcLocator* pSvcLocator)
+  : AthReentrantAlgorithm(name, pSvcLocator) { 
+
+}
+
+StatusCode 
+ScintWaveformRecAlg::initialize() {
+  ATH_MSG_INFO(name() << "::initalize()" );
+
+  // Initalize tools
+  ATH_CHECK( m_recoTool.retrieve() );
+
+  // Set key to read waveform from
+  ATH_CHECK( m_waveformContainerKey.initialize() );
+
+  // Set key to read clock info
+  ATH_CHECK( m_clockKey.initialize() );
+
+  // Set key to write container
+  ATH_CHECK( m_waveformHitContainerKey.initialize() );
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode 
+ScintWaveformRecAlg::finalize() {
+  ATH_MSG_INFO(name() << "::finalize()");
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode 
+ScintWaveformRecAlg::execute(const EventContext& ctx) const {
+  ATH_MSG_INFO("Executing");
+
+  ATH_MSG_DEBUG("Run: " << ctx.eventID().run_number() 
+		<< " Event: " << ctx.eventID().event_number());
+
+  // Find the input waveform container
+  SG::ReadHandle<ScintWaveformContainer> waveformHandle(m_waveformContainerKey, ctx);
+
+  ATH_CHECK( waveformHandle.isValid() );
+  ATH_MSG_DEBUG("Found ReadHandle for ScintWaveformContainer " << m_waveformContainerKey);
+
+  if (waveformHandle->size() == 0) {
+    ATH_MSG_INFO("Waveform container found with zero length!");
+    return StatusCode::SUCCESS;
+  }
+
+  // Also find the clock information
+  SG::ReadHandle<xAOD::WaveformClock> clockHandle(m_clockKey, ctx);
+
+  ATH_CHECK( clockHandle.isValid() );
+  ATH_MSG_DEBUG("Found 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 each waveform
+  for( const auto& wave : *waveformHandle) {
+
+    ATH_MSG_DEBUG("Reconstruct waveform for channel " << wave->channel());
+
+    // Reconstruct the hits, may be more than one, so pass container
+    CHECK( m_recoTool->reconstruct(*wave, clockHandle.ptr(), 
+				   hitContainerHandle.ptr()) );
+
+  }
+
+  ATH_MSG_DEBUG("WaveformsHitContainer '" << hitContainerHandle.name() << "' filled with "<< hitContainerHandle->size() <<" items");
+
+  return StatusCode::SUCCESS;
+}
+

Scintillator/ScintRecAlgs/src/ScintWaveformRecAlg.h

+#ifndef SCINTRECALGS_SCINTWAVEFORMRECALG_H
+#define SCINTRECALGS_SCINTWAVEFORMRECALG_H
+
+// Base class
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+
+// Data classes
+#include "ScintRawEvent/ScintWaveformContainer.h"
+
+#include "xAODFaserWaveform/WaveformClock.h"
+#include "xAODFaserWaveform/WaveformClockAuxInfo.h"
+
+#include "xAODFaserWaveform/WaveformHit.h"
+#include "xAODFaserWaveform/WaveformHitContainer.h"
+#include "xAODFaserWaveform/WaveformHitAuxContainer.h"
+
+// Tool classes
+#include "ScintRecTools/IWaveformReconstructionTool.h"
+
+// Handles
+#include "StoreGate/ReadHandleKey.h"
+#include "StoreGate/WriteHandleKey.h"
+
+// Gaudi
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/ToolHandle.h"
+
+// STL
+#include <string>
+
+class ScintWaveformRecAlg : public AthReentrantAlgorithm {
+
+ public:
+  // Constructor
+  ScintWaveformRecAlg(const std::string& name, ISvcLocator* pSvcLocator);
+  virtual ~ScintWaveformRecAlg() = 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 */
+  //@{
+  ScintWaveformRecAlg() = delete;
+  ScintWaveformRecAlg(const ScintWaveformRecAlg&) = delete;
+  ScintWaveformRecAlg &operator=(const ScintWaveformRecAlg&) = delete;
+  //@}
+
+  /**
+   * @name Reconstruction tool
+   */
+  ToolHandle<IWaveformReconstructionTool> m_recoTool
+    {this, "WaveformReconstructionTool", "WaveformReconstructionTool"};
+
+  /**
+   * @name Input raw waveform data using SG::ReadHandleKey
+   */
+  //@{
+  SG::ReadHandleKey<ScintWaveformContainer> m_waveformContainerKey
+    {this, "WaveformContainerKey", ""};
+  //@}
+
+  /**
+   * @name Input WaveformClock data using SG::ReadHandleKey
+   */
+  //@{
+  SG::ReadHandleKey<xAOD::WaveformClock> m_clockKey
+    {this, "WaveformClockKey", "WaveformClock"};
+  //@}
+
+  /**
+   * @name Output data using SG::WriteHandleKey
+   */
+  //@{
+  SG::WriteHandleKey<xAOD::WaveformHitContainer> m_waveformHitContainerKey
+    {this, "WaveformHitContainerKey", ""};
+  //@}
+
+};
+
+#endif // SCINTRECALGS_SCINTWAVEFORMRECALG_H

Scintillator/ScintRecAlgs/src/components/ScintRecAlgs_entries.cxx

+#include "../ScintWaveformRecAlg.h"
+#include "../ScintClockRecAlg.h"
+
+DECLARE_COMPONENT( ScintWaveformRecAlg )
+DECLARE_COMPONENT( ScintClockRecAlg )

Scintillator/ScintRecTools/CMakeLists.txt

+################################################################################
+# Package: ScintRecTools
+################################################################################
+
+# Declare the package name:
+atlas_subdir( ScintRecTools )
+
+# External dependencies:
+find_package( ROOT )
+
+# Component(s) in the package:
+atlas_add_library( ScintRecToolsLib
+                   ScintRecTools/*.h src/*.cxx src/*.h
+                   PUBLIC_HEADERS ScintRecTools
+                   PRIVATE_INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
+                   LINK_LIBRARIES AthenaBaseComps AthenaKernel GeoPrimitives ScintRawEvent xAODFaserWaveform
+                   PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES}
+		   )
+
+atlas_add_component( ScintRecTools
+		     src/components/*.cxx 
+		     INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
+                     LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps GaudiKernel ScintRecToolsLib )
+
+

Scintillator/ScintRecTools/ScintRecTools/IClockReconstructionTool.h

+/*
+  Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+*/
+
+/**
+ * @file IClockReconstructionTool.h
+ * Header file for the IClockReconstructionTool class
+ * @author Eric Torrence, 2021
+ */
+
+
+#ifndef SCINTRECTOOLS_ICLOCKRECONSTRUCTIONTOOL_H
+#define SCINTRECTOOLS_ICLOCKRECONSTRUCTIONTOOL_H
+
+// Base class
+#include "GaudiKernel/IAlgTool.h"
+#include "GaudiKernel/ToolHandle.h"
+
+#include "xAODFaserWaveform/WaveformClock.h"
+
+class ScintWaveform;
+
+///Interface for Clock reco algorithms
+class IClockReconstructionTool : virtual public IAlgTool 
+{
+  public:
+
+  // InterfaceID
+  DeclareInterfaceID(IClockReconstructionTool, 1, 0);
+
+  virtual ~IClockReconstructionTool() = default;
+
+  // Reconstruct all peaks in a raw waveform
+  virtual StatusCode reconstruct(const ScintWaveform& wave, 
+				 xAOD::WaveformClock* clockdata) const = 0;
+
+};
+
+#endif // SCINTRECTOOLS_ICLOCKRECONSTRUCTIONTOOL_H

Scintillator/ScintRecTools/ScintRecTools/IWaveformReconstructionTool.h

+/*
+  Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+*/
+
+/**
+ * @file IWaveformReconstructionTool.h
+ * Header file for the IWaveformReconstructionTool class
+ * @author Eric Torrence, 2021
+ */
+
+
+#ifndef SCINTRECTOOLS_IWAVEFORMRECONSTRUCTIONTOOL_H
+#define SCINTRECTOOLS_IWAVEFORMRECONSTRUCTIONTOOL_H
+
+// Base class
+#include "GaudiKernel/IAlgTool.h"
+#include "GaudiKernel/ToolHandle.h"
+
+#include "xAODFaserWaveform/WaveformHitContainer.h"
+#include "xAODFaserWaveform/WaveformClock.h"
+
+class ScintWaveform;
+
+///Interface for Waveform reco algorithms
+class IWaveformReconstructionTool : virtual public IAlgTool 
+{
+  public:
+
+  // InterfaceID
+  DeclareInterfaceID(IWaveformReconstructionTool, 1, 0);
+
+  virtual ~IWaveformReconstructionTool() = default;
+
+  // Reconstruct all peaks in a raw waveform
+  virtual StatusCode reconstruct(const ScintWaveform& wave, 
+				 const xAOD::WaveformClock* clock, 
+				 xAOD::WaveformHitContainer* container) const = 0;
+
+};
+
+#endif // SCINTRECTOOLS_IWAVEFORMRECONSTRUCTIONTOOL_H

Scintillator/ScintRecTools/src/ClockReconstructionTool.cxx

+/*
+  Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+*/
+
+/**
+ * @file ClockReconstructionTool.cxx
+ * Implementation file for the ClockReconstructionTool.cxx
+ * @ author E. Torrence, 2021
+ **/
+
+#include "ClockReconstructionTool.h"
+
+#include "TVirtualFFT.h"
+
+#include <vector>
+
+// Constructor
+ClockReconstructionTool::ClockReconstructionTool(const std::string& type, const std::string& name, const IInterface* parent) :
+  base_class(type, name, parent)
+{
+}
+
+// Initialization
+StatusCode
+ClockReconstructionTool::initialize() {
+  ATH_MSG_INFO( name() << "::initalize()" );
+
+  return StatusCode::SUCCESS;
+}
+
+// Reconstruction step
+StatusCode
+ClockReconstructionTool::reconstruct(const ScintWaveform& raw_wave,
+				     xAOD::WaveformClock* clockdata) const {
+
+  ATH_MSG_DEBUG("Clock reconstruct called ");
+
+  // Check the output object
+  if (!clockdata) {
+    ATH_MSG_ERROR("WaveformClock passed to reconstruct() is null!");
+    return StatusCode::FAILURE;
+  }
+
+  // Set the trigger time
+  //clockdata->setTriggerTime(raw_wave.trigger_time_tag());
+  //ATH_MSG_DEBUG("Trigger time: " << raw_wave.trigger_time_tag());
+
+  // Digitized clock data, sampled at 500 MHz (2 ns)
+  auto counts = raw_wave.adc_counts();
+
+  // Need a double array for FFT
+  int N = counts.size();
+  std::vector<double> wave(N);
+  wave.assign(counts.begin(), counts.end());
+
+  ATH_MSG_DEBUG("Created double array with length " << wave.size() );
+  ATH_MSG_DEBUG("First 10 elements:");
+  for (unsigned int i=0; i< 10; i++)
+    ATH_MSG_DEBUG(" " << i << " " << wave[i]);
+
+  // delta_nu = 1/T where T is the total waveform length
+  // Multiplying (freqmult * i) gives the frequency of point i in the FFT
+  double freqmult = 1./(0.002 * N); // 2ns per point, so this is in MHz
+
+  TVirtualFFT *fftr2c = TVirtualFFT::FFT(1, &N, "R2C");
+  fftr2c->SetPoints(&wave[0]);
+  fftr2c->Transform();
+
+  // Get the coefficients
+  std::vector<double> re_full(N);
+  std::vector<double> im_full(N);
+  std::vector<double> magnitude(N/2); 
+  fftr2c->GetPointsComplex(&re_full[0],&im_full[0]);
+
+  // Normalize the magnitude (just using the positive complex frequencies)
+  unsigned int i=0;
+  magnitude[i] = sqrt(pow(re_full[i], 2) + pow(im_full[i], 2))/N;
+  for(i=1; i<magnitude.size(); i++) 
+    magnitude[i] = 2*sqrt(pow(re_full[i], 2) + pow(im_full[i], 2))/N;
+
+  // First, look at the DC offset
+  ATH_MSG_DEBUG("DC offset: " << magnitude[0]);
+
+  // Second, find max value (should be primary line at 40 MHz
+  unsigned int imax = max_element(magnitude.begin()+1, magnitude.end()) - magnitude.begin();
+
+  if (((int(imax)-1) < 0) || ((imax+1) >= magnitude.size())) {
+    ATH_MSG_WARNING("Found maximum frequency amplitude at postion " << imax << "!");
+  } else {
+    ATH_MSG_DEBUG("Magnitude array at peak:");
+    for(i = imax-1; i <= imax+1; i++)
+      ATH_MSG_DEBUG("Index: " << i << " Freq: " << i*freqmult << " Mag: " << magnitude[i]);
+  }
+
+  // Store results
+  clockdata->set_dc_offset(magnitude[0]);
+  clockdata->set_frequency(imax * freqmult);
+  clockdata->set_amplitude(magnitude[imax]);
+  clockdata->set_phase(atan2(im_full[imax], re_full[imax]));
+
+  ATH_MSG_DEBUG("Before correcting for finite resolution:");
+  ATH_MSG_DEBUG(*clockdata);
+
+  // Correct for the windowing to find the exact peak frequency
+  // Following: https://indico.cern.ch/event/132526/contributions/128902/attachments/99707/142376/Meeting1-06-11_FFT_corrections_for_tune_measurements.pdf
+
+  double dm;  // Fraction of integer frequency where real peak sits
+  if (magnitude[imax+1] >= magnitude[imax-1]) { // Past ipeak by dm
+    dm = magnitude[imax+1] / (magnitude[imax] + magnitude[imax + 1]);
+  } else { // Before ipeak by dm
+    dm = -magnitude[imax-1] / (magnitude[imax-1] + magnitude[imax]);
+  }
+
+  ATH_MSG_DEBUG("Found shift in frequency index: " << dm);
+
+  // Improved values
+
+  double phase = atan2(im_full[imax], re_full[imax]) - dm * M_PI;
+  // Fix any overflows
+  if (phase < M_PI) phase += (2*M_PI);
+  if (phase > M_PI) phase -= (2*M_PI);
+
+  clockdata->set_frequency( (imax+dm) * freqmult );
+  clockdata->set_phase (phase);
+  clockdata->set_amplitude( 2*M_PI*dm*magnitude[imax] / sin(M_PI * dm) );
+
+  ATH_MSG_DEBUG("After correcting for finite resolution:");
+  ATH_MSG_DEBUG(*clockdata);
+
+  delete fftr2c;
+
+  return StatusCode::SUCCESS;
+}

Scintillator/ScintRecTools/src/ClockReconstructionTool.h

+/*
+   Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+*/
+
+/** @file ClockReconstructionTool.h
+ *  Header file for ClockReconstructionTool.h
+ *
+ */
+#ifndef SCINTRECTOOLS_CLOCKRECONSTRUCTIONTOOL_H
+#define SCINTRECTOOLS_CLOCKRECONSTRUCTIONTOOL_H
+
+//Athena
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "ScintRecTools/IClockReconstructionTool.h"
+
+#include "ScintRawEvent/ScintWaveform.h"
+#include "xAODFaserWaveform/WaveformClock.h"
+
+//Gaudi
+#include "GaudiKernel/ToolHandle.h"
+
+//STL
+#include <string>
+#include <vector>
+
+class ClockReconstructionTool: public extends<AthAlgTool, IClockReconstructionTool> {
+ public:
+
+  /// Normal constructor for an AlgTool; 'properties' are also declared here
+  ClockReconstructionTool(const std::string& type, 
+			  const std::string& name, const IInterface* parent);
+
+  /// Retrieve the necessary services in initialize
+  StatusCode initialize();
+
+  /// Reconstruct hits from clock  
+  virtual StatusCode reconstruct(const ScintWaveform& wave,
+				 xAOD::WaveformClock* clockdata) const;
+
+ private:
+
+  //
+  // Baseline Estimation Parameters
+  //BooleanProperty m_useSimpleBaseline{this, "UseSimpleBaseline", false};
+  //IntegerProperty m_samplesForBaselineAverage{this, "SamplesForBaselineAverage", 40};
+  //FloatProperty m_baselineFitWindow{this, "BaselineFitWindow", 2.};
+
+};
+
+#endif // SCINTRECTOOLS_CLOCKRECONSTRUCTIONTOOL_H

Scintillator/ScintRecTools/src/WaveformBaselineData.cxx

+/*
+  Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+*/
+
+#include "WaveformBaselineData.h"
+
+WaveformBaselineData::WaveformBaselineData() {
+  clear();
+}
+
+WaveformBaselineData::WaveformBaselineData(bool isvalid) {
+  clear();
+  valid = isvalid;
+}
+
+void
+WaveformBaselineData::clear() {
+  channel = -1;
+  valid = false;
+  mean = 0.;
+  rms = 0.;
+  fit_status = -1;
+  peak = 0.;
+  chi2n = 0.;
+}
+
+std::ostream
+&operator<<( std::ostream& out, const WaveformBaselineData &base ) {
+  out << "WaveformBaselineData - Mean: " << base.mean << ", RMS: " << base.rms << std::endl;
+  return out;
+}
+
+
+
+
+

Scintillator/ScintRecTools/src/WaveformBaselineData.h

+/*
+  Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// WaveformBaselineData.h
+//   Header file for class WaveformBaselineData
+///////////////////////////////////////////////////////////////////
+// Class to handle the baseline information for a raw waveform
+///////////////////////////////////////////////////////////////////
+// Version 1.0   2/21/2021 Eric Torrence
+///////////////////////////////////////////////////////////////////
+#ifndef SCINTRECTOOLS_WAVEFORMBASELINEDATA_H
+#define SCINTRECTOOLS_WAVEFORMBASELINEDATA_H
+
+#include <iostream>
+
+/**
+ * @class WaveformBaselineData
+ * The baseline is the estimate of the waveform reading for no signal.
+ */
+  
+class WaveformBaselineData {
+      
+ public:
+      
+  /** Default constructor */
+  WaveformBaselineData();
+
+  // Set status directly
+  WaveformBaselineData(bool valid);
+
+  // move assignment defaulted
+  WaveformBaselineData & operator= (WaveformBaselineData &&) = default;
+  //assignment defaulted
+  WaveformBaselineData & operator= (const WaveformBaselineData &) = default;
+  //copy c'tor defaulted
+  WaveformBaselineData(const WaveformBaselineData &) = default;
+
+  /** Destructor */
+  virtual ~WaveformBaselineData() = default;
+
+  /** Clones */
+  virtual WaveformBaselineData* clone() const ;       
+
+  // Clear
+  void clear();
+            
+  // Data members
+  int channel;  // Channel number
+  bool valid;   // Overall valididty flag
+
+  double mean;  // Best value for mean baseline
+  double rms;   // Best value for baseline RMS
+  int fit_status; // Return value from fit
+  double peak;  // Amplitude of baseline fit
+  double chi2n; // Reduced chi2/Ndof from gauss fit
+
+ private:
+
+};
+
+std::ostream
+&operator<<(std::ostream &out, const WaveformBaselineData &wfm);
+
+///////////////////////////////////////////////////////////////////
+// Inline methods:
+///////////////////////////////////////////////////////////////////
+
+inline WaveformBaselineData* 
+WaveformBaselineData::clone() const { 
+  return new WaveformBaselineData(*this);  
+}
+
+#endif // SCINTRECTOOLS_WAVEFORMBASELINEDATA_H

Scintillator/ScintRecTools/src/WaveformFitResult.cxx

+/*
+  Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+*/
+
+#include "WaveformFitResult.h"
+
+WaveformFitResult::WaveformFitResult() {
+  clear();
+}
+
+WaveformFitResult::WaveformFitResult(bool isvalid) {
+  clear();
+  valid = isvalid;
+}
+
+void
+WaveformFitResult::clear() {
+  valid = false;
+  peak = 0.;
+  mean = 0.;
+  sigma = 0.;
+  fit_status = -1;
+  chi2n = 0.;
+}
+
+std::ostream
+&operator<<( std::ostream& out, const WaveformFitResult &fit ) {
+  out << "WaveformFitResult - Mean: " << fit.mean << ", Sigma: " << fit.sigma << ", Peak: " << fit.peak << std::endl;
+  return out;
+}
+
+
+
+
+

Scintillator/ScintRecTools/src/WaveformFitResult.h

+/*
+  Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// WaveformFitResult.h
+//   Header file for class WaveformFitResult
+///////////////////////////////////////////////////////////////////
+// Class to handle the fit output for a raw waveform
+///////////////////////////////////////////////////////////////////
+// Version 1.0   2/21/2021 Eric Torrence
+///////////////////////////////////////////////////////////////////
+#ifndef SCINTRECTOOLS_WAVEFORMFITRESULT_H
+#define SCINTRECTOOLS_WAVEFORMFITRESULT_H
+
+#include <iostream>
+
+/**
+ * @class WaveformFitResult
+ * The baseline is the estimate of the waveform reading for no signal.
+ */
+  
+class WaveformFitResult {
+      
+ public:
+      
+  /** Default constructor */
+  WaveformFitResult();
+
+  // Set status directly
+  WaveformFitResult(bool valid);
+
+  // move assignment defaulted
+  WaveformFitResult & operator= (WaveformFitResult &&) = default;
+  //assignment defaulted
+  WaveformFitResult & operator= (const WaveformFitResult &) = default;
+  //copy c'tor defaulted
+  WaveformFitResult(const WaveformFitResult &) = default;
+
+  /** Destructor */
+  virtual ~WaveformFitResult() = default;
+
+  // Clear
+  void clear();
+            
+  // Data members
+  bool valid;   // Overall valididty (from fit result == 0)
+
+  double peak;  // Best value for peak
+  double mean;  // Best value for mean
+  double sigma; // Best value for width
+  double alpha; // Crystal ball parameter
+  double nval;  // Crystal ball parameter
+  double chi2ndf; 
+  double integral;
+  double time;  // Time of leading edge of pulse
+
+  int fit_status; // Return value from fit
+
+  double chi2n; // Reduced chi2/Ndof from fit
+
+ private:
+
+};
+
+std::ostream
+&operator<<(std::ostream &out, const WaveformFitResult &wfm);
+
+///////////////////////////////////////////////////////////////////
+// Inline methods:
+///////////////////////////////////////////////////////////////////
+
+#endif // SCINTRECTOOLS_WAVEFORMFITRESULT_H

Scintillator/ScintRecTools/src/WaveformReconstructionTool.cxx

+/*
+  Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+*/
+
+/**
+ * @file WaveformReconstructionTool.cxx
+ * Implementation file for the WaveformReconstructionTool.cxx
+ * @ author E. Torrence, 2021
+ **/
+
+#include "WaveformReconstructionTool.h"
+
+#include "xAODFaserWaveform/WaveformHit.h"
+
+#include "TH1F.h"
+#include "TF1.h"
+#include "TFitResult.h"
+#include "TFitResultPtr.h"
+#include "TGraph.h"
+
+#include <vector>
+#include <tuple>
+
+// Constructor
+WaveformReconstructionTool::WaveformReconstructionTool(const std::string& type, const std::string& name, const IInterface* parent) :
+  base_class(type, name, parent)
+{
+}
+
+// Initialization
+StatusCode
+WaveformReconstructionTool::initialize() {
+  ATH_MSG_INFO( name() << "::initalize()" );
+
+  if (m_useSimpleBaseline.value()) {
+    ATH_MSG_INFO("Will use simple baseline estimation");
+  } else {
+    ATH_MSG_INFO("Will use fit to determine baseline");
+  }
+  return StatusCode::SUCCESS;
+}
+
+// Reconstruction step
+StatusCode
+WaveformReconstructionTool::reconstruct(const ScintWaveform& raw_wave,
+					const xAOD::WaveformClock* clock, 
+					xAOD::WaveformHitContainer* container) const {
+
+  ATH_MSG_DEBUG(" reconstruct called ");
+
+  // Check the container
+  if (!container) {
+    ATH_MSG_ERROR("WaveformHitCollection passed to reconstruct() is null!");
+    return StatusCode::FAILURE;
+  }
+
+  //
+  // We always want to create at least one hit, so create it here
+  xAOD::WaveformHit* hit = new xAOD::WaveformHit();
+  container->push_back(hit);
+  hit->set_channel(raw_wave.channel());
+
+  //
+  // Make some sanity checks on the waveform data
+  if (raw_wave.adc_counts().size() == 0) {
+    ATH_MSG_WARNING( "Found waveform for channel " << raw_wave.channel() 
+		     << " with size " << raw_wave.adc_counts().size() << "!");
+
+    hit->set_status_bit(xAOD::WaveformStatus::WAVEFORM_MISSING);
+    return StatusCode::SUCCESS;
+  }
+    
+  if (raw_wave.adc_counts().size() != raw_wave.n_samples()) {
+    ATH_MSG_WARNING( "Found waveform for channel " << raw_wave.channel() 
+		     << " with size " << raw_wave.adc_counts().size() 
+		     << " not equal to number of samples " << raw_wave.n_samples());
+
+    hit->set_status_bit(xAOD::WaveformStatus::WAVEFORM_INVALID);
+    return StatusCode::SUCCESS;
+  }
+
+  // Set channel
+  hit->set_channel(raw_wave.channel());
+
+  //
+  // Find baseline
+  WaveformBaselineData baseline;
+
+  if (m_useSimpleBaseline.value())
+    baseline = findSimpleBaseline(raw_wave);
+  else
+    baseline = findAdvancedBaseline(raw_wave);
+
+  if (!(baseline.valid)) {
+    ATH_MSG_WARNING("Failed to reconstruct baseline!");
+
+    hit->set_baseline_mean(0.);
+    hit->set_baseline_rms(-1.);
+    hit->set_status_bit(xAOD::WaveformStatus::BASELINE_FAILED);
+
+    return StatusCode::SUCCESS;
+  }
+
+  // Save baseline to hit collection object
+  hit->set_baseline_mean(baseline.mean);
+  hit->set_baseline_rms(baseline.rms);
+
+  //
+  // Create baseline-subtracted data array for both time and signal
+  // Time in ns from start of readout
+  unsigned int size = raw_wave.adc_counts().size();
+  std::vector<float> time(size);  
+  for (unsigned int i=0; i<size; i++)
+    time[i] = 2.*i;
+
+  // Baseline subtracted (and inverted) ADC waveform values
+  std::vector<float> wave(raw_wave.adc_counts().begin(), raw_wave.adc_counts().end());
+  for (auto& element : wave)
+    element = baseline.mean - element;
+
+  bool first = true;
+
+  // Now we iteratively find peaks and fit
+  while(true) {
+
+    //
+    // Find peak in array and return time and value arrays
+    // This range of data is also *removed* from original arrays
+    std::vector<float> wtime;
+    std::vector<float> wwave;
+
+    // All done if we don't have any peaks above threshold
+    // If we do find a significant peak, fill the window
+    if (! findPeak(baseline, time, wave, wtime, wwave) ) {
+      if (first) hit->set_status_bit(xAOD::WaveformStatus::THRESHOLD_FAILED);
+      break;
+    }
+
+    //
+    // Create new hit to fill
+    if (!first) {
+      hit = new xAOD::WaveformHit();
+      container->push_back(hit);
+      hit->set_status_bit(xAOD::WaveformStatus::SECONDARY);
+    }
+    first = false;
+
+    //
+    // Save windowed waveform to Hit object
+    hit->set_channel(raw_wave.channel());
+    hit->set_baseline_mean(baseline.mean);
+    hit->set_baseline_rms(baseline.rms);
+    hit->set_time_vector(wtime);
+    hit->set_wave_vector(wave);
+
+    //
+    // Find some raw values
+    WaveformFitResult raw = findRawHitValues(wtime, wwave);
+    hit->set_peak(raw.peak);
+    hit->set_mean(raw.mean);
+    hit->set_width(raw.sigma);
+    hit->set_integral(raw.integral);
+    hit->set_localtime(raw.mean);
+    hit->set_raw_peak(raw.peak);
+    hit->set_raw_integral(raw.integral);
+
+    //
+    // Perform Gaussian fit to waveform
+    WaveformFitResult gfit = fitGaussian(raw, wtime, wwave);
+    if (! gfit.valid) {
+      ATH_MSG_WARNING( " Gaussian waveform fit failed! ");
+      hit->set_status_bit(xAOD::WaveformStatus::GFIT_FAILED);
+    } 
+
+    // Fit results (or raw if it failed)
+    hit->set_peak(gfit.peak);
+    hit->set_mean(gfit.mean);
+    hit->set_width(gfit.sigma);
+    hit->set_integral(gfit.integral);
+    hit->set_localtime(gfit.time);
+
+    //
+    // Check for overflow
+    if (m_removeOverflow && findOverflow(baseline, wtime, wwave)) {
+      ATH_MSG_INFO("Found waveform overflow");
+      hit->set_status_bit(xAOD::WaveformStatus::WAVE_OVERFLOW);
+    }
+
+    //
+    // Perform CB fit
+    WaveformFitResult cbfit = fitCBall(gfit, wtime, wwave);
+    if (! cbfit.valid) {
+      ATH_MSG_WARNING("CrystalBall fit failed!");
+      // Still have gaussian parameters as an estimate
+      hit->set_status_bit(xAOD::WaveformStatus::CBFIT_FAILED);
+    } else {
+      hit->set_peak(cbfit.peak);
+      hit->set_mean(cbfit.mean);
+      hit->set_width(cbfit.sigma);
+      hit->set_integral(cbfit.integral);
+      hit->set_localtime(cbfit.time);
+
+      hit->set_alpha(cbfit.alpha);
+      hit->set_nval(cbfit.nval);
+    }
+
+    //
+    // Find time from clock
+    if (!clock) {
+      hit->set_status_bit(xAOD::WaveformStatus::CLOCK_INVALID);
+    } else {
+      hit->set_bcid_time(clock->time_from_clock(hit->localtime()));
+    }
+
+    if (! m_findMultipleHits) break;
+
+  } // End of loop over waveform data
+
+  ATH_MSG_DEBUG( "WaveformReconstructionTool finished for channel " 
+		 << raw_wave.channel() << " container size= " << container->size());
+
+  return StatusCode::SUCCESS;
+}
+
+bool
+WaveformReconstructionTool::findPeak(WaveformBaselineData& baseline, 
+   std::vector<float>& time, std::vector<float>& wave,
+   std::vector<float>& windowed_time, std::vector<float>& windowed_wave) const {
+
+  ATH_MSG_DEBUG("findPeak called");
+
+  // Find max value location in array
+  unsigned int imax = std::max_element(wave.begin(), wave.end()) - wave.begin();
+  float maxval = wave[imax];
+  ATH_MSG_DEBUG( "Found peak value " << maxval << " at position " << imax );
+
+  // Check if this is over threshold (in sigma)
+  if (maxval < m_peakThreshold*baseline.rms) {
+    ATH_MSG_DEBUG("Failed threshold");
+    return false;
+  }
+
+  // Make a window around this peak, values are in bins, so units of 2ns
+  // Ensure our window is within the vector range
+  int lo_edge = ((int(imax) + m_windowStart) >= 0 ? (imax + m_windowStart) : 0); 
+  int hi_edge = ((imax + m_windowStart + m_windowWidth) < wave.size() ? (imax + m_windowStart + m_windowWidth) : wave.size());
+  
+  ATH_MSG_DEBUG("Windowing waveform from " << lo_edge << " to " << hi_edge);
+  windowed_time = std::vector<float> (time.begin()+lo_edge, time.begin()+hi_edge);
+  windowed_wave = std::vector<float> (wave.begin()+lo_edge, wave.begin()+hi_edge);
+
+  // Remove these values from the original arrays so we can iterate
+  time.erase(time.begin()+lo_edge, time.begin()+hi_edge);
+  wave.erase(wave.begin()+lo_edge, wave.begin()+hi_edge);
+  
+  return true;
+}
+
+bool
+WaveformReconstructionTool::findOverflow(const WaveformBaselineData& base, 
+	      std::vector<float>& time, std::vector<float>& wave) const {
+
+  auto peakloc = std::max_element(wave.begin(), wave.end());
+
+  // If peak value is less than baseline, we have no overflow
+  if (*peakloc < int(base.mean)) return false;
+
+  ATH_MSG_DEBUG("Removing overflows from waveform with length " << wave.size());
+
+  // We have an overflow, remove all elements that are overflowing
+  unsigned int i = peakloc - wave.begin();
+  for (; i<wave.size(); i++) {
+    if (wave[i] < int(base.mean)) continue;
+
+    ATH_MSG_DEBUG("Removing position "<< i<< " with value " << wave[i] << " > " << int(base.mean));
+    // This is an overflow, remove elements
+    time.erase(time.begin() + i);
+    wave.erase(wave.begin() + i);
+    i--;  // Decrement so that loop itertaion points to next element
+  }
+
+  ATH_MSG_DEBUG("Removed overflows, length now " << wave.size());
+  return true;
+}
+
+WaveformBaselineData&
+WaveformReconstructionTool::findSimpleBaseline(const ScintWaveform& raw_wave) const {
+
+  ATH_MSG_DEBUG( "findSimpleBaseline called" );
+  //ATH_MSG_DEBUG( raw_wave );
+
+  // This must be static so we can return a reference
+  static WaveformBaselineData baseline;  
+  baseline.clear();
+
+  // Calling algorithm checks for proper data
+  // Here we just check algorithm-specific issues
+  if (int(raw_wave.n_samples()) < m_samplesForBaselineAverage.value()) {
+    ATH_MSG_WARNING( "Found waveform with " << raw_wave.n_samples() << " samples, not enough to find baseline!" );
+    return baseline;
+  }
+
+  const std::vector<unsigned int>& counts = raw_wave.adc_counts();
+
+  double sumx = 0.;
+  double sumx2 = 0.;
+
+  // Look at first n values in the waveform
+  unsigned int nvalues = m_samplesForBaselineAverage.value();
+
+  for (unsigned int i=0; i< nvalues; i++) {
+    sumx  += counts[i];
+    sumx2 += (counts[i] * counts[i]);
+  }
+  double mean = sumx / nvalues;
+  double mean2 = sumx2 / nvalues;
+
+  double rms = std::sqrt(mean2 - mean*mean);
+
+  baseline.valid = true;
+  baseline.channel = raw_wave.channel();
+  baseline.mean = mean;
+  baseline.rms = rms;
+  return baseline;
+}
+
+WaveformBaselineData&
+WaveformReconstructionTool::findAdvancedBaseline(const ScintWaveform& raw_wave) const {
+
+  ATH_MSG_DEBUG( "findAdvancedBaseline called" );
+
+  // This must be static so we can return a reference
+  static WaveformBaselineData baseline;  
+  baseline.clear();
+
+  // First histogram to find most likely value
+  // Turn these into configurables once this works
+  int nbins = m_baselineRangeBins;
+  double xlo = 0.;
+  double xhi = m_baselineRange;
+
+  TH1F h1("", "", nbins, xlo, xhi);
+
+  // Fill this histogram with the waveform
+  for (auto value : raw_wave.adc_counts()) {
+    //ATH_MSG_DEBUG( "Found value " << value );
+    h1.Fill(value);
+  }
+
+  // Find max bin
+  int maxbin = h1.GetMaximumBin();
+  double maxbinval = h1.GetXaxis()->GetBinCenter(maxbin);
+  ATH_MSG_DEBUG( "Found max bin at " << maxbinval << " counts");
+
+  // Fill second histogram with integer resolution around the peak
+  nbins = m_baselineFitRange;
+  xlo = int(maxbinval - nbins/2)-0.5;
+  xhi = xlo + nbins;
+  ATH_MSG_DEBUG( "Filling 2nd histogram from " << xlo << " to " << xhi);
+
+  TH1F h2("", "", nbins, xlo, xhi);
+  for (auto value : raw_wave.adc_counts()) {
+    h2.Fill(value);
+  }
+
+  // Start with full histogram range
+  double mean = h2.GetMean();
+  double rms = h2.GetRMS();
+  double peak = h2.GetBinContent(h2.GetMaximumBin());
+
+  ATH_MSG_DEBUG( "Initial Mean: " << mean << " RMS: " << rms << " Peak: " << peak );
+
+  // Restrict range to +/- 2 sigma of mean
+  double window = m_baselineFitWindow;  // Window range in sigma
+  h2.GetXaxis()->SetRangeUser(mean-window*rms, mean+window*rms);
+  mean = h2.GetMean();
+  rms = h2.GetRMS();
+  peak = h2.GetBinContent(h2.GetMaximumBin());
+
+  ATH_MSG_DEBUG( "2 Sigma Mean: " << mean << " RMS: " << rms << " Peak: " << peak);
+
+  // Finally, fit a Gaussian to the 2 sigma range
+  double fitlo = mean-window*rms;
+  double fithi = mean+window*rms;
+
+  // Define fit function and preset range
+  TF1 fitfunc("BaseFit", "gaus", fitlo, fithi);
+  fitfunc.SetParameters(peak, mean, rms);
+
+  // Fit Options:
+  // L - Likelihood fit
+  // Q - Quiet mode (V for verbose)
+  // N - Don't draw or store graphics
+  // S - Return fit results
+  // E - Better error estimation
+  // M - Use TMinuit IMPROVE to find a better minimum
+  // R - Use function range  
+  TFitResultPtr fit = h2.Fit(&fitfunc, "LQNSR", "");
+
+  int fitStatus = fit;
+  double chi2 = fit->Chi2();
+  unsigned int ndf = fit->Ndf();
+  if (ndf == 0) ndf = 1;
+
+  if (!fit->IsValid()) {
+    ATH_MSG_WARNING( " Baseline fit failed! ");
+  } else {
+    // Improve estimation with fit results
+    peak = fit->Parameter(0);
+    mean = fit->Parameter(1);
+    rms  = fit->Parameter(2);
+    ATH_MSG_DEBUG( "G Fit   Mean: " << mean << " RMS: " << rms << " Peak: " << peak << " Chi2/N: " << chi2/ndf);  
+  }
+
+  baseline.valid = true;
+  baseline.channel = raw_wave.channel();
+  baseline.mean = mean;
+  baseline.rms = rms;
+  baseline.fit_status = fitStatus; 
+  baseline.peak = peak;
+  baseline.chi2n = chi2/ndf;
+  return baseline;
+}
+
+WaveformFitResult&
+WaveformReconstructionTool::findRawHitValues(const std::vector<float> time, const std::vector<float> wave) const {
+
+  ATH_MSG_DEBUG("findRawHitValues called");
+
+  // This must be static so we can return a reference
+  static WaveformFitResult rfit;
+  rfit.clear();
+
+  // First, estimate starting values from input waveform
+  // Will use this to pre-set fit values, but also as a backup in case fit fails
+  double tot = 0.;
+  double sum = 0.;
+  double sum2 = 0.;
+  for (unsigned int i=0; i<time.size(); i++) {
+    tot += wave[i];
+    sum += time[i] * wave[i];
+    sum2 += time[i] * time[i] * wave[i];
+  }
+
+  // Pointer to peak
+  auto peakloc = std::max_element(wave.begin(), wave.end());
+
+  // Initial parameters from waveform
+  rfit.mean = sum/tot;
+  rfit.sigma = std::sqrt(sum2/tot - rfit.mean*rfit.mean);
+  rfit.peak = *peakloc;
+  rfit.integral = tot;
+  rfit.time = rfit.mean;
+
+  ATH_MSG_DEBUG( "Initial Mean: " << rfit.mean << " RMS: " << rfit.sigma << " Peak: " << rfit.peak);  
+
+  return rfit;
+}
+
+WaveformFitResult&
+WaveformReconstructionTool::fitGaussian(const WaveformFitResult& raw, const std::vector<float> time, const std::vector<float> wave) const {
+
+  ATH_MSG_DEBUG("fitGaussian called");
+
+  // This must be static so we can return a reference
+  static WaveformFitResult gfit;
+  gfit.clear();
+
+  // Start with raw values by default
+  gfit = raw;
+
+  // Start by creating a TGraph and fitting this with a Gaussian
+  // Must pass arrays to TGraph, use pointer to first vector element
+  TGraph tg(time.size(), &time[0], &wave[0]);
+
+  // Define fit function and preset range
+  TF1 gfunc("gfunc", "gaus");
+  gfunc.SetParameters(raw.peak, raw.mean, raw.sigma);
+  gfunc.SetParError(0, std::sqrt(raw.peak));
+  gfunc.SetParError(1, raw.sigma);
+  gfunc.SetParError(2, raw.sigma / 10.);
+  gfunc.SetParLimits(2, 0., 20.);  // Constrain width
+
+  TFitResultPtr gfitptr = tg.Fit(&gfunc, "QNS", "");
+
+  gfit.fit_status = gfitptr;
+  gfit.valid = (gfit.fit_status == 0);
+
+  if (!gfitptr->IsValid()) {
+    ATH_MSG_WARNING( " Gaussian waveform fit failed! ");
+  } else {
+    // Improve estimation with fit results
+    gfit.peak = gfitptr->Parameter(0);
+    gfit.mean = gfitptr->Parameter(1);
+    gfit.sigma = gfitptr->Parameter(2);
+    gfit.integral = gfunc.Integral(time.front(), time.back());
+
+    // Find time here
+    gfit.time = gfunc.GetX( gfit.peak * m_timingPeakFraction, time.front(), gfit.mean );
+
+    ATH_MSG_DEBUG("G Fit   Mean: " << gfit.mean << " RMS: " << gfit.sigma 
+		  << " Peak: " << gfit.peak << " Int: " << gfit.integral << " Time: " << gfit.time);
+  }
+
+  return gfit;
+}
+
+WaveformFitResult&
+WaveformReconstructionTool::fitCBall(const WaveformFitResult& gfit, 
+     const std::vector<float> time, const std::vector<float> wave) const {
+
+  ATH_MSG_DEBUG("fitCBall called");
+
+  // This must be static so we can return a reference
+  static WaveformFitResult cbfit;
+  cbfit.clear();
+
+  // Must pass arrays to TGraph, use pointer to first vector element
+  TGraph tg(time.size(), &time[0], &wave[0]);
+
+  // Values to preset CB fit
+  cbfit.peak = abs(gfit.peak);
+  cbfit.mean = gfit.mean;
+  cbfit.sigma = abs(gfit.sigma);
+  if (cbfit.sigma > 20.) cbfit.sigma = 5.;
+  cbfit.alpha = -0.25;  // Negative to get tail on high side
+  cbfit.nval = 3.;
+
+  // Define fit function and preset values
+  TF1 cbfunc("cbfunc", "crystalball");
+  cbfunc.SetParameter(0, cbfit.peak); // Peak height
+  cbfunc.SetParError(0, std::sqrt(cbfit.peak));
+  cbfunc.SetParameter(1, cbfit.mean); // Mean
+  cbfunc.SetParError(1, cbfit.sigma);
+  cbfunc.SetParameter(2, cbfit.sigma);  // Width
+  cbfunc.SetParError(2, cbfit.sigma/10.);
+  cbfunc.SetParLimits(2, 0., 20.);
+  cbfunc.SetParameter(3, cbfit.alpha); // Tail parameter (negative for high-side tail)
+  cbfunc.SetParError(3, -cbfit.alpha/10.);
+  cbfunc.SetParLimits(3, -10., 0.);
+  cbfunc.SetParameter(4, cbfit.nval);  // Tail power
+  cbfunc.SetParError(4, cbfit.nval/10.);
+  cbfunc.SetParLimits(4, 0., 1.E3);
+
+  TFitResultPtr cbfitptr = tg.Fit(&cbfunc, "QNS", "");
+
+  cbfit.fit_status = cbfitptr;
+  cbfit.valid = (cbfit.fit_status == 0);
+
+  if (!cbfitptr->IsValid()) {
+    ATH_MSG_WARNING( " Crystal Ball waveform fit failed! ");
+  } else {
+    // Improve estimation with fit results
+    cbfit.peak = cbfitptr->Parameter(0);
+    cbfit.mean = cbfitptr->Parameter(1);
+    cbfit.sigma = cbfitptr->Parameter(2);
+    cbfit.alpha = cbfitptr->Parameter(3);
+    cbfit.nval = cbfitptr->Parameter(4);
+
+    double chi2 = cbfitptr->Chi2();
+    unsigned int ndf = cbfitptr->Ndf();
+    if (ndf == 0) ndf = 1;
+    cbfit.chi2ndf = chi2/ndf;
+
+    cbfit.integral = cbfunc.Integral(time.front(), time.back());
+
+    // Find time here
+    cbfit.time = cbfunc.GetX( cbfit.peak * m_timingPeakFraction, time.front(), cbfit.mean );
+
+    ATH_MSG_DEBUG("CB Fit  Mean: " << cbfit.mean << " RMS: " << cbfit.sigma 
+		  << " Peak: " << cbfit.peak << " Int: " << cbfit.integral
+		  << " Time: " << cbfit.time 
+		  << " N: " << cbfit.nval << " Alpha: " << cbfit.alpha 
+		  << " Chi2/Ndf: " << cbfit.chi2ndf );
+
+  }
+
+  return cbfit;
+}

Scintillator/ScintRecTools/src/WaveformReconstructionTool.h

+/*
+   Copyright (C) 2021 CERN for the benefit of the FASER collaboration
+*/
+
+/** @file WaveformReconstructionTool.h
+ *  Header file for WaveformReconstructionTool.h
+ *
+ */
+#ifndef SCINTRECTOOLS_WAVEFORMRECONSTRUCTIONTOOL_H
+#define SCINTRECTOOLS_WAVEFORMRECONSTRUCTIONTOOL_H
+
+//Athena
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "ScintRecTools/IWaveformReconstructionTool.h"
+
+#include "ScintRawEvent/ScintWaveform.h"
+
+#include "WaveformBaselineData.h"
+#include "WaveformFitResult.h"
+
+//Gaudi
+#include "GaudiKernel/ToolHandle.h"
+
+//STL
+#include <string>
+#include <vector>
+
+class WaveformReconstructionTool: public extends<AthAlgTool, IWaveformReconstructionTool> {
+ public:
+
+  /// Normal constructor for an AlgTool; 'properties' are also declared here
+  WaveformReconstructionTool(const std::string& type, 
+			     const std::string& name, const IInterface* parent);
+
+  /// Retrieve the necessary services in initialize
+  StatusCode initialize();
+
+  /// Reconstruct hits from waveform
+  
+  virtual StatusCode reconstruct(const ScintWaveform& wave,
+				 const xAOD::WaveformClock* clock,
+				 xAOD::WaveformHitContainer* container) const;
+
+ private:
+
+  //
+  // Baseline Estimation Parameters
+  BooleanProperty m_useSimpleBaseline{this, "UseSimpleBaseline", false};
+  IntegerProperty m_samplesForBaselineAverage{this, "SamplesForBaselineAverage", 40};
+
+  //
+  // Parameters of initial histogram to find baseline max
+  // Range and bins to use, ratio should be an integer (bin width)
+  IntegerProperty m_baselineRange{this, "BaselineRange", 16000};
+  IntegerProperty m_baselineRangeBins{this, "BaselineRangeBins", 320};
+
+  // 
+  // Parameters for the Gaussian fit to the baseline peak (in counts)
+  // Range is total range to use to find truncated mean and width
+  IntegerProperty m_baselineFitRange{this, "BaselineFitRange", 200};
+  // Fit window is value (in sigma) of trucated width to use in final fit
+  FloatProperty m_baselineFitWindow{this, "BaselineFitWindow", 2.};
+
+  //
+  // Peak threshold (in sigma of baseline RMS) to find a hit
+  FloatProperty m_peakThreshold{this, "PeakThreshold", 10.};
+ 
+  //
+  // Window to define fitting range, in samples (2ns/sample)
+  IntegerProperty m_windowStart{this, "FitWindowStart", -10};
+  IntegerProperty m_windowWidth{this, "FitWindowWidth", 50};
+
+  //
+  // Remove overflow values from CB fit
+  BooleanProperty m_removeOverflow{this, "RemoveOverflow", true};
+
+  //
+  // Look for more than one hit in each channel
+  BooleanProperty m_findMultipleHits{this, "FindMultipleHits", false};
+
+  //
+  // Fraction of peak to set local hit time
+  FloatProperty m_timingPeakFraction{this, "TimingPeakFraction", 0.2};
+
+  // Baseline algorithms
+  WaveformBaselineData& findSimpleBaseline(const ScintWaveform& wave) const;
+  WaveformBaselineData& findAdvancedBaseline(const ScintWaveform& wave) const;
+
+  // Find peak in wave, return windowed region in windowed_time and windowed_wave
+  // Windowed region is removed from original vectors
+  // Returns true if peak found, false if not
+  bool findPeak(WaveformBaselineData& baseline,
+		std::vector<float>& time, std::vector<float>& wave,
+		std::vector<float>& windowed_time, std::vector<float>& windowed_wave) const;
+
+  // Get estimate from waveform data itself
+  WaveformFitResult& findRawHitValues(const std::vector<float> time, 
+				      const std::vector<float> wave) const;
+
+  // Fit windowed data to Gaussian (to get initial estimate of parameters
+  WaveformFitResult& fitGaussian(const WaveformFitResult& raw,
+				 const std::vector<float> time, 
+				 const std::vector<float> wave) const;
+
+  // Find overflows and remove points from arrays
+  bool findOverflow(const WaveformBaselineData& baseline, 
+		    std::vector<float>& time, std::vector<float>& wave) const;
+
+  // Fit windowed data to CrystalBall function
+  WaveformFitResult& fitCBall(const WaveformFitResult& gfit, 
+			      const std::vector<float> time, 
+			      const std::vector<float> wave) const;
+
+
+};
+
+#endif // SCINTRECTOOLS_WAVEFORMRECONSTRUCTIONTOOL_H