diff --git a/Scintillator/ScintEventCnv/ScintByteStream/src/ScintWaveformDecoderTool.cxx b/Scintillator/ScintEventCnv/ScintByteStream/src/ScintWaveformDecoderTool.cxx
index de415ef295d11ca6ca9fec2724edfa231f852560..55c9bd960477cb6fbdf7871098a8a5cae9f9c3a8 100644
--- a/Scintillator/ScintEventCnv/ScintByteStream/src/ScintWaveformDecoderTool.cxx
+++ b/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;
}
diff --git a/Scintillator/ScintRawEvent/src/ScintWaveform.cxx b/Scintillator/ScintRawEvent/src/ScintWaveform.cxx
index b525f2e0b51f7ea6e41ee283633ab0c55e3c085e..bc9f5eb912e4bd7cd75fefda9b4a829699f92578 100644
--- a/Scintillator/ScintRawEvent/src/ScintWaveform.cxx
+++ b/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;
}
diff --git a/Scintillator/ScintRecAlgs/CMakeLists.txt b/Scintillator/ScintRecAlgs/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..7ad1776170385a4f2203030dbac8fc5f7f5f43ab
--- /dev/null
+++ b/Scintillator/ScintRecAlgs/CMakeLists.txt
@@ -0,0 +1,15 @@
+################################################################################
+# 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 )
+
diff --git a/Scintillator/ScintRecAlgs/python/ScintRecAlgsConfig.py b/Scintillator/ScintRecAlgs/python/ScintRecAlgsConfig.py
new file mode 100644
index 0000000000000000000000000000000000000000..6319e9a01771a042f162c11f51d4f5b8d3bba181
--- /dev/null
+++ b/Scintillator/ScintRecAlgs/python/ScintRecAlgsConfig.py
@@ -0,0 +1,70 @@
+""" 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
diff --git a/Scintillator/ScintRecAlgs/src/ScintClockRecAlg.cxx b/Scintillator/ScintRecAlgs/src/ScintClockRecAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..45ceaa30ba91d814072a02d312755d03d3bd1b85
--- /dev/null
+++ b/Scintillator/ScintRecAlgs/src/ScintClockRecAlg.cxx
@@ -0,0 +1,94 @@
+#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;
+}
+
diff --git a/Scintillator/ScintRecAlgs/src/ScintClockRecAlg.h b/Scintillator/ScintRecAlgs/src/ScintClockRecAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..3fde8295da9b5ca4cdecf68bceae205df2a7ef8a
--- /dev/null
+++ b/Scintillator/ScintRecAlgs/src/ScintClockRecAlg.h
@@ -0,0 +1,72 @@
+#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
diff --git a/Scintillator/ScintRecAlgs/src/ScintWaveformRecAlg.cxx b/Scintillator/ScintRecAlgs/src/ScintWaveformRecAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..e3d42247e30d9b295e037338f93f092bcee32281
--- /dev/null
+++ b/Scintillator/ScintRecAlgs/src/ScintWaveformRecAlg.cxx
@@ -0,0 +1,81 @@
+#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;
+}
+
diff --git a/Scintillator/ScintRecAlgs/src/ScintWaveformRecAlg.h b/Scintillator/ScintRecAlgs/src/ScintWaveformRecAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..7d43e97523bdba7ddc8749dddc40108dafa29f35
--- /dev/null
+++ b/Scintillator/ScintRecAlgs/src/ScintWaveformRecAlg.h
@@ -0,0 +1,86 @@
+#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
diff --git a/Scintillator/ScintRecAlgs/src/components/ScintRecAlgs_entries.cxx b/Scintillator/ScintRecAlgs/src/components/ScintRecAlgs_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..9634fc8f9a4adbd9d6950b1727fc5f19f72d7ec0
--- /dev/null
+++ b/Scintillator/ScintRecAlgs/src/components/ScintRecAlgs_entries.cxx
@@ -0,0 +1,5 @@
+#include "../ScintWaveformRecAlg.h"
+#include "../ScintClockRecAlg.h"
+
+DECLARE_COMPONENT( ScintWaveformRecAlg )
+DECLARE_COMPONENT( ScintClockRecAlg )
diff --git a/Scintillator/ScintRecTools/CMakeLists.txt b/Scintillator/ScintRecTools/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a128a73845fdfdfea13d2c43e6085d2f5e99b2d1
--- /dev/null
+++ b/Scintillator/ScintRecTools/CMakeLists.txt
@@ -0,0 +1,25 @@
+################################################################################
+# 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 )
+
+
diff --git a/Scintillator/ScintRecTools/ScintRecTools/IClockReconstructionTool.h b/Scintillator/ScintRecTools/ScintRecTools/IClockReconstructionTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..46712a75ca41eceabf1331fd2ae1412d9f7f429b
--- /dev/null
+++ b/Scintillator/ScintRecTools/ScintRecTools/IClockReconstructionTool.h
@@ -0,0 +1,39 @@
+/*
+ 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
diff --git a/Scintillator/ScintRecTools/ScintRecTools/IWaveformReconstructionTool.h b/Scintillator/ScintRecTools/ScintRecTools/IWaveformReconstructionTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..00218de1cf251e239be6cd5e066e02721154be61
--- /dev/null
+++ b/Scintillator/ScintRecTools/ScintRecTools/IWaveformReconstructionTool.h
@@ -0,0 +1,41 @@
+/*
+ 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
diff --git a/Scintillator/ScintRecTools/src/ClockReconstructionTool.cxx b/Scintillator/ScintRecTools/src/ClockReconstructionTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..92442bed407a6b62d653d803d42ce7c18a7b6ccf
--- /dev/null
+++ b/Scintillator/ScintRecTools/src/ClockReconstructionTool.cxx
@@ -0,0 +1,133 @@
+/*
+ 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;
+}
diff --git a/Scintillator/ScintRecTools/src/ClockReconstructionTool.h b/Scintillator/ScintRecTools/src/ClockReconstructionTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..c7f76d6b6934ecd54339e0e762d1245b52d4a0dd
--- /dev/null
+++ b/Scintillator/ScintRecTools/src/ClockReconstructionTool.h
@@ -0,0 +1,50 @@
+/*
+ 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
diff --git a/Scintillator/ScintRecTools/src/WaveformBaselineData.cxx b/Scintillator/ScintRecTools/src/WaveformBaselineData.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..f2e92625cff312592eeeadbde13715cc04a69775
--- /dev/null
+++ b/Scintillator/ScintRecTools/src/WaveformBaselineData.cxx
@@ -0,0 +1,36 @@
+/*
+ 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;
+}
+
+
+
+
+
diff --git a/Scintillator/ScintRecTools/src/WaveformBaselineData.h b/Scintillator/ScintRecTools/src/WaveformBaselineData.h
new file mode 100644
index 0000000000000000000000000000000000000000..13144a4f92da45982a32a7b25053d075a2a8bb2b
--- /dev/null
+++ b/Scintillator/ScintRecTools/src/WaveformBaselineData.h
@@ -0,0 +1,75 @@
+/*
+ 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
diff --git a/Scintillator/ScintRecTools/src/WaveformFitResult.cxx b/Scintillator/ScintRecTools/src/WaveformFitResult.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..45dcdae23e6760a28617e5507893e44a36be6a38
--- /dev/null
+++ b/Scintillator/ScintRecTools/src/WaveformFitResult.cxx
@@ -0,0 +1,35 @@
+/*
+ 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;
+}
+
+
+
+
+
diff --git a/Scintillator/ScintRecTools/src/WaveformFitResult.h b/Scintillator/ScintRecTools/src/WaveformFitResult.h
new file mode 100644
index 0000000000000000000000000000000000000000..92ecf0fc6e69d1b481ffdbc84b0959be2049fd6a
--- /dev/null
+++ b/Scintillator/ScintRecTools/src/WaveformFitResult.h
@@ -0,0 +1,73 @@
+/*
+ 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
diff --git a/Scintillator/ScintRecTools/src/WaveformReconstructionTool.cxx b/Scintillator/ScintRecTools/src/WaveformReconstructionTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..037cde51247e2c70b7f5db4728d336e28233d8ff
--- /dev/null
+++ b/Scintillator/ScintRecTools/src/WaveformReconstructionTool.cxx
@@ -0,0 +1,579 @@
+/*
+ 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;
+}
diff --git a/Scintillator/ScintRecTools/src/WaveformReconstructionTool.h b/Scintillator/ScintRecTools/src/WaveformReconstructionTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..45a1e9c79d1229e337e6223090959d6b955f7b03
--- /dev/null
+++ b/Scintillator/ScintRecTools/src/WaveformReconstructionTool.h
@@ -0,0 +1,117 @@
+/*
+ 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
diff --git a/Scintillator/ScintRecTools/src/components/ScintRecTools_entries.cxx b/Scintillator/ScintRecTools/src/components/ScintRecTools_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..a96a42c9b5e4a0007558650314ef6e1144281877
--- /dev/null
+++ b/Scintillator/ScintRecTools/src/components/ScintRecTools_entries.cxx
@@ -0,0 +1,5 @@
+#include "../ClockReconstructionTool.h"
+#include "../WaveformReconstructionTool.h"
+
+DECLARE_COMPONENT( ClockReconstructionTool )
+DECLARE_COMPONENT( WaveformReconstructionTool )
diff --git a/xAOD/xAODFaserWaveform/CMakeLists.txt b/xAOD/xAODFaserWaveform/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..1185251139a20ea833a672371173ad00c2987ccd
--- /dev/null
+++ b/xAOD/xAODFaserWaveform/CMakeLists.txt
@@ -0,0 +1,26 @@
+# Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+
+# Declare the package name.
+atlas_subdir( xAODFaserWaveform )
+
+# External dependencies.
+find_package( xAODUtilities )
+
+# Component(s) in the package.
+atlas_add_library( xAODFaserWaveform
+ xAODFaserWaveform/*.h xAODFaserWaveform/versions/*.h Root/*.cxx
+ PUBLIC_HEADERS xAODFaserWaveform
+ LINK_LIBRARIES xAODCore )
+
+atlas_add_xaod_smart_pointer_dicts(
+ INPUT xAODFaserWaveform/selection.xml
+ OUTPUT _selectionFile
+ OBJECTS "xAOD::WaveformClock_v1"
+ CONTAINERS "xAOD::WaveformHitContainer_v1")
+
+atlas_add_dictionary( xAODFaserWaveformDict
+ xAODFaserWaveform/xAODFaserWaveformDict.h
+ ${_selectionFile}
+ LINK_LIBRARIES xAODCore xAODFaserWaveform
+ EXTRA_FILES Root/dict/*.cxx )
+
diff --git a/xAOD/xAODFaserWaveform/Root/WaveformClockAuxInfo_v1.cxx b/xAOD/xAODFaserWaveform/Root/WaveformClockAuxInfo_v1.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..61d5dc8e8d736c44ceba480964551bed85ccf981
--- /dev/null
+++ b/xAOD/xAODFaserWaveform/Root/WaveformClockAuxInfo_v1.cxx
@@ -0,0 +1,18 @@
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+// Local include(s):
+#include "xAODFaserWaveform/versions/WaveformClockAuxInfo_v1.h"
+
+namespace xAOD {
+
+ WaveformClockAuxInfo_v1::WaveformClockAuxInfo_v1()
+ : AuxInfoBase(),
+ frequency(0), phase(0), dc_offset(0), amplitude(0) {
+ AUX_VARIABLE( frequency );
+ AUX_VARIABLE( phase );
+ AUX_VARIABLE( dc_offset );
+ AUX_VARIABLE( amplitude );
+ }
+} // namespace xAOD
diff --git a/xAOD/xAODFaserWaveform/Root/WaveformClock_v1.cxx b/xAOD/xAODFaserWaveform/Root/WaveformClock_v1.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..04498ccc8ecfc6e0a2c62d1c404318988229b891
--- /dev/null
+++ b/xAOD/xAODFaserWaveform/Root/WaveformClock_v1.cxx
@@ -0,0 +1,47 @@
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+// EDM include(s):
+#include "xAODCore/AuxStoreAccessorMacros.h"
+
+// Local include(s):
+#include "xAODFaserWaveform/versions/WaveformClock_v1.h"
+
+namespace xAOD {
+
+ WaveformClock_v1::WaveformClock_v1() : SG::AuxElement() {
+ }
+
+ AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformClock_v1, double, frequency, set_frequency )
+
+ AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformClock_v1, double, phase, set_phase )
+
+ AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformClock_v1, double, dc_offset, set_dc_offset )
+
+ AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformClock_v1, double, amplitude, set_amplitude )
+
+ float WaveformClock_v1::time_from_clock(float time) const {
+ // Figure out which integer cycle we are on, must add 1/4 of a cycle
+ // because FFT finds phase of cosine, which would be middle of
+ // clock HI region. frequency is in MHz, so convert here to GHz
+ int ncycle = int(time*frequency() / 1.E3 + phase() / (2*M_PI) + 0.25);
+ float dt = time - (ncycle - phase() / (2*M_PI) - 0.25) * 1.E3 / frequency();
+ return dt;
+ }
+
+} // namespace xAOD
+
+namespace xAOD {
+
+ std::ostream& operator<<(std::ostream& s, const xAOD::WaveformClock_v1& clk) {
+ s << "xAODWaveformClock: frequency=" << clk.frequency()
+ << " phase=" << clk.phase()
+ << " amplitude=" << clk.amplitude()
+ << " offset=" << clk.dc_offset()
+ << std::endl;
+
+ return s;
+ }
+
+} // namespace xAOD
diff --git a/xAOD/xAODFaserWaveform/Root/WaveformHitAuxContainer_v1.cxx b/xAOD/xAODFaserWaveform/Root/WaveformHitAuxContainer_v1.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..76aed56247567f13e39fbcbf7471950dc93b3c23
--- /dev/null
+++ b/xAOD/xAODFaserWaveform/Root/WaveformHitAuxContainer_v1.cxx
@@ -0,0 +1,35 @@
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+// Local include(s):
+#include "xAODFaserWaveform/versions/WaveformHitAuxContainer_v1.h"
+
+namespace xAOD {
+
+ WaveformHitAuxContainer_v1::WaveformHitAuxContainer_v1()
+ : AuxContainerBase() {
+
+ AUX_VARIABLE(channel);
+ AUX_VARIABLE(localtime);
+ AUX_VARIABLE(peak);
+ AUX_VARIABLE(width);
+ AUX_VARIABLE(integral);
+ AUX_VARIABLE(bcid_time);
+ AUX_VARIABLE(raw_peak);
+ AUX_VARIABLE(raw_integral);
+
+ AUX_VARIABLE(baseline_mean);
+ AUX_VARIABLE(baseline_rms);
+
+ AUX_VARIABLE(status);
+ AUX_VARIABLE(mean);
+ AUX_VARIABLE(alpha);
+ AUX_VARIABLE(nval);
+ AUX_VARIABLE(time_vector);
+ AUX_VARIABLE(wave_vector);
+
+ }
+
+} // namespace xAOD
+
diff --git a/xAOD/xAODFaserWaveform/Root/WaveformHit_v1.cxx b/xAOD/xAODFaserWaveform/Root/WaveformHit_v1.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..4d4ce278978e68024637c117bc1e9456c2e22730
--- /dev/null
+++ b/xAOD/xAODFaserWaveform/Root/WaveformHit_v1.cxx
@@ -0,0 +1,61 @@
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+// EDM include(s):
+#include "xAODCore/AuxStoreAccessorMacros.h"
+
+// Local include(s):
+#include "xAODFaserWaveform/versions/WaveformHit_v1.h"
+
+namespace xAOD {
+
+ WaveformHit_v1::WaveformHit_v1() : SG::AuxElement() {
+ }
+
+ AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, unsigned int, channel, set_channel )
+
+ AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, float, localtime, set_localtime )
+
+ AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, float, peak, set_peak )
+
+ AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, float, width, set_width )
+
+ AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, float, integral, set_integral )
+
+ AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, float, bcid_time, set_bcid_time )
+
+ AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, float, raw_peak, set_raw_peak )
+
+ AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, float, raw_integral, set_raw_integral )
+
+ AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, float, baseline_mean, set_baseline_mean )
+
+ AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, float, baseline_rms, set_baseline_rms )
+
+ AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, unsigned int, status, set_status )
+
+ AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, float, mean, set_mean )
+
+ AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, float, alpha, set_alpha )
+
+ AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, float, nval, set_nval )
+
+ AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, std::vector<float>, time_vector, set_time_vector )
+
+ AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( WaveformHit_v1, std::vector<float>, wave_vector, set_wave_vector )
+
+} // namespace xAOD
+
+namespace xAOD {
+
+ std::ostream& operator<<(std::ostream& s, const xAOD::WaveformHit_v1& hit) {
+ s << "xAODWaveformHit: channel=" << hit.channel()
+ << " local time=" << hit.localtime()
+ << " peak=" << hit.peak()
+ << std::endl;
+
+ return s;
+ }
+
+} // namespace xAOD
diff --git a/xAOD/xAODFaserWaveform/Root/dict/ContainerProxies.cxx b/xAOD/xAODFaserWaveform/Root/dict/ContainerProxies.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..5ef691ac9c9cd4e7091a9ed1d6640090d69bc016
--- /dev/null
+++ b/xAOD/xAODFaserWaveform/Root/dict/ContainerProxies.cxx
@@ -0,0 +1,8 @@
+// EDM include(s):
+#include "xAODCore/AddDVProxy.h"
+
+// Local include(s):
+#include "xAODFaserWaveform/WaveformHitContainer.h"
+
+// Set up the collection proxies:
+ADD_NS_DV_PROXY( xAOD, WaveformHitContainer );
diff --git a/xAOD/xAODFaserWaveform/Root/xAODFaserWaveformCLIDs.cxx b/xAOD/xAODFaserWaveform/Root/xAODFaserWaveformCLIDs.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..acd57bce0d941b16c815e9805b7da8e76bf2ef27
--- /dev/null
+++ b/xAOD/xAODFaserWaveform/Root/xAODFaserWaveformCLIDs.cxx
@@ -0,0 +1,10 @@
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+//simple includes to force the CLASS_DEF etc to be encountered during compile
+
+#include "xAODFaserWaveform/WaveformHitContainer.h"
+#include "xAODFaserWaveform/WaveformHitAuxContainer.h"
+#include "xAODFaserWaveform/WaveformClock.h"
+#include "xAODFaserWaveform/WaveformClockAuxInfo.h"
diff --git a/xAOD/xAODFaserWaveform/xAODFaserWaveform/WaveformClock.h b/xAOD/xAODFaserWaveform/xAODFaserWaveform/WaveformClock.h
new file mode 100644
index 0000000000000000000000000000000000000000..951d21a7a706e01c11b642490e4d71f7d019a8af
--- /dev/null
+++ b/xAOD/xAODFaserWaveform/xAODFaserWaveform/WaveformClock.h
@@ -0,0 +1,22 @@
+// Dear emacs, this is -*- c++ -*-
+
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+#ifndef XAODFASERWAVEFORM_WAVEFORMCLOCK_H
+#define XAODFASERWAVEFORM_WAVEFORMCLOCK_H
+
+// Local include(s):
+#include "xAODFaserWaveform/versions/WaveformClock_v1.h"
+
+namespace xAOD {
+ /// Declare the latest version of the class
+ typedef WaveformClock_v1 WaveformClock;
+}
+
+// Set up a CLID for the container:
+#include "xAODCore/CLASS_DEF.h"
+CLASS_DEF( xAOD::WaveformClock, 58376762, 1 )
+
+#endif // XAODFASERWAVEFORM_WAVEFORMCLOCK_H
diff --git a/xAOD/xAODFaserWaveform/xAODFaserWaveform/WaveformClockAuxInfo.h b/xAOD/xAODFaserWaveform/xAODFaserWaveform/WaveformClockAuxInfo.h
new file mode 100644
index 0000000000000000000000000000000000000000..5cdc5294f041b7370ca745e09dfe55223b47b3f9
--- /dev/null
+++ b/xAOD/xAODFaserWaveform/xAODFaserWaveform/WaveformClockAuxInfo.h
@@ -0,0 +1,22 @@
+// Dear emacs, this is -*- c++ -*-
+
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+#ifndef XAODFASERWAVEFORM_WAVEFORMCLOCKAUXINFO_H
+#define XAODFASERWAVEFORM_WAVEFORMCLOCKAUXINFO_H
+
+// Local include(s):
+#include "xAODFaserWaveform/versions/WaveformClockAuxInfo_v1.h"
+
+namespace xAOD {
+ /// Declare the latest version of the class
+ typedef WaveformClockAuxInfo_v1 WaveformClockAuxInfo;
+}
+
+// Set up a CLID for the container:
+#include "xAODCore/CLASS_DEF.h"
+CLASS_DEF( xAOD::WaveformClockAuxInfo, 150155999, 1 )
+
+#endif // XAODFASERWAVEFORM_WAVEFORMCLOCKAUXINFO_H
diff --git a/xAOD/xAODFaserWaveform/xAODFaserWaveform/WaveformHit.h b/xAOD/xAODFaserWaveform/xAODFaserWaveform/WaveformHit.h
new file mode 100644
index 0000000000000000000000000000000000000000..6f01f4fc47ca8e213dbd8634ea7444b5db0270b1
--- /dev/null
+++ b/xAOD/xAODFaserWaveform/xAODFaserWaveform/WaveformHit.h
@@ -0,0 +1,22 @@
+// Dear emacs, this is -*- c++ -*-
+
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+#ifndef XAODFASERWAVEFORM_WAVEFORMHIT_H
+#define XAODFASERWAVEFORM_WAVEFORMHIT_H
+
+// Local include(s):
+#include "xAODFaserWaveform/versions/WaveformHit_v1.h"
+
+namespace xAOD {
+ /// Declare the latest version of the class
+ typedef WaveformHit_v1 WaveformHit;
+}
+
+// Set up a CLID for the container:
+#include "xAODCore/CLASS_DEF.h"
+CLASS_DEF( xAOD::WaveformHit, 131600577, 1 )
+
+#endif // XAODFASERWAVEFORM_WAVEFORMHIT_H
diff --git a/xAOD/xAODFaserWaveform/xAODFaserWaveform/WaveformHitAuxContainer.h b/xAOD/xAODFaserWaveform/xAODFaserWaveform/WaveformHitAuxContainer.h
new file mode 100644
index 0000000000000000000000000000000000000000..014da4c02d3c80f2f2068f5994579deb33ecf056
--- /dev/null
+++ b/xAOD/xAODFaserWaveform/xAODFaserWaveform/WaveformHitAuxContainer.h
@@ -0,0 +1,22 @@
+// Dear emacs, this is -*- c++ -*-
+
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+#ifndef XAODFASERWAVEFORM_WAVEFORMHITAUXCONTAINER_H
+#define XAODFASERWAVEFORM_WAVEFORMHITAUXCONTAINER_H
+
+// Local include(s):
+#include "xAODFaserWaveform/versions/WaveformHitAuxContainer_v1.h"
+
+namespace xAOD {
+ /// Declare the latest version of the class
+ typedef WaveformHitAuxContainer_v1 WaveformHitAuxContainer;
+}
+
+// Set up a CLID for the container:
+#include "xAODCore/CLASS_DEF.h"
+CLASS_DEF( xAOD::WaveformHitAuxContainer, 1262669778, 1 )
+
+#endif // XAODFASERWAVEFORM_WAVEFORMHITAUXCONTAINER_H
diff --git a/xAOD/xAODFaserWaveform/xAODFaserWaveform/WaveformHitContainer.h b/xAOD/xAODFaserWaveform/xAODFaserWaveform/WaveformHitContainer.h
new file mode 100644
index 0000000000000000000000000000000000000000..41de79d21d8fa88389710a9dbef214653bc6edf1
--- /dev/null
+++ b/xAOD/xAODFaserWaveform/xAODFaserWaveform/WaveformHitContainer.h
@@ -0,0 +1,22 @@
+// Dear emacs, this is -*- c++ -*-
+
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+#ifndef XAODFASERWAVEFORM_WAVEFORMHITCONTAINER_H
+#define XAODFASERWAVEFORM_WAVEFORMHITCONTAINER_H
+
+// Local include(s):
+#include "xAODFaserWaveform/versions/WaveformHitContainer_v1.h"
+
+namespace xAOD {
+ /// Declare the latest version of the class
+ typedef WaveformHitContainer_v1 WaveformHitContainer;
+}
+
+// Set up a CLID for the container:
+#include "xAODCore/CLASS_DEF.h"
+CLASS_DEF( xAOD::WaveformHitContainer, 1156844391, 1 )
+
+#endif // XAODFASERWAVEFORM_WAVEFORMHITCONTAINER_H
diff --git a/xAOD/xAODFaserWaveform/xAODFaserWaveform/selection.xml b/xAOD/xAODFaserWaveform/xAODFaserWaveform/selection.xml
new file mode 100644
index 0000000000000000000000000000000000000000..44146cd55fad3798ba00f365c671949e36e38423
--- /dev/null
+++ b/xAOD/xAODFaserWaveform/xAODFaserWaveform/selection.xml
@@ -0,0 +1,24 @@
+<!-- Copyright (C) 2020 CERN for the benefit of the FASER collaboration -->
+<lcgdict>
+
+ <class name="xAOD::WaveformHit_v1" />
+ <typedef name="xAOD::WaveformHit" />
+
+ <class name="xAOD::WaveformHitContainer_v1"
+ id="fb49e082-7130-462f-9500-de9a2aef8a24" />
+ <typedef name="xAOD::WaveformHitContainer" />
+
+ <class name="xAOD::WaveformHitAuxContainer_v1"
+ id="091007bd-f4ab-4862-b905-1a0ed96c962f" />
+ <typedef name="xAOD::WaveformHitAuxContainer" />
+
+ <class name="xAOD::WaveformClock_v1"
+ id="461aa5ca-6c72-46ee-b5d5-5bd6df4e2848" />
+ <typedef name="xAOD::WaveformClock" />
+
+ <class name="xAOD::WaveformClockAuxInfo_v1"
+ id="fe0600cc-f2f3-4be5-9723-257646dbb8f7" />
+ <typedef name="xAOD::WaveformClockAuxInfo" />
+
+
+</lcgdict>
diff --git a/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformClockAuxInfo_v1.h b/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformClockAuxInfo_v1.h
new file mode 100644
index 0000000000000000000000000000000000000000..24741e2476ab952f316a8ceeaf7611361bdd291d
--- /dev/null
+++ b/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformClockAuxInfo_v1.h
@@ -0,0 +1,45 @@
+// Dear emacs, this is -*- c++ -*-
+
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+#ifndef XAODFASERWAVEFORM_VERSIONS_WAVEFORMCLOCKAUXINFO_V1_H
+#define XAODFASERWAVEFORM_VERSIONS_WAVEFORMCLOCKAUXINFO_V1_H
+
+// System include(s):
+extern "C" {
+# include "stdint.h"
+}
+
+// xAOD include(s):
+#include "xAODCore/AuxInfoBase.h"
+
+namespace xAOD {
+
+ /// Class holding the data handled by WaveformClock_v1
+
+ class WaveformClockAuxInfo_v1 : public AuxInfoBase {
+
+ public:
+ /// Default constructor
+ WaveformClockAuxInfo_v1();
+
+ private:
+ /// @name Basic variables
+ ///@ {
+ double frequency;
+ double phase;
+ double dc_offset;
+ double amplitude;
+ ///@}
+
+ }; // class WaveformClockAuxInfo_v1
+
+} // namespace xAOD
+
+// Set up a CLID and StoreGate inheritance for the class:
+#include "xAODCore/BaseInfo.h"
+SG_BASE( xAOD::WaveformClockAuxInfo_v1, xAOD::AuxInfoBase );
+
+#endif // XAODFASERWAVEFORM_VERSIONS_WAVEFORMCLOCKAUXINFO_V1_H
diff --git a/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformClock_v1.h b/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformClock_v1.h
new file mode 100644
index 0000000000000000000000000000000000000000..437b3ecabf1f71c9be0c353d27bebec1d0916d11
--- /dev/null
+++ b/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformClock_v1.h
@@ -0,0 +1,61 @@
+// Dear emacs, this is -*- c++ -*-
+
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+#ifndef XAODFASERWAVEFORM_VERSIONS_WAVEFORMCLOCK_V1_H
+#define XAODFASERWAVEFORM_VERSIONS_WAVEFORMCLOCK_V1_H
+
+// System include(s):
+extern "C" {
+# include "stdint.h"
+}
+
+// Core EDM include(s):
+#include "AthContainers/AuxElement.h"
+
+namespace xAOD {
+
+ // Cllss describing pulses in the waveform digitizer
+ class WaveformClock_v1 : public SG::AuxElement {
+
+ public:
+ /// Defaullt constructor
+ WaveformClock_v1();
+
+ /// @name Access WaveformClock elements
+ /// @{
+
+ /// Clock Frequency (in MHz)
+ double frequency() const;
+ void set_frequency(double value);
+
+ /// Clock Phase (in Radians)
+ double phase() const;
+ void set_phase(double value);
+
+ /// DC Clock offset
+ double dc_offset() const;
+ void set_dc_offset(double value);
+
+ /// Amplitude of primary freq. component
+ double amplitude() const;
+ void set_amplitude(double value);
+
+ /// Distance of time (in ns) from previous rising clock edge
+ float time_from_clock(float time) const;
+
+ /// @}
+
+ }; // class WaveformClock_v1
+
+ std::ostream& operator<<(std::ostream& s, const xAOD::WaveformClock_v1& clk);
+
+} // namespace xAOD
+
+// Declare the inheritance of the type:
+#include "xAODCore/BaseInfo.h"
+SG_BASE( xAOD::WaveformClock_v1, SG::AuxElement );
+
+#endif // XAODFASERWAVEFORM_VERSIONS_WAVEFORMCLOCK_V1_H
diff --git a/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHitAuxContainer_v1.h b/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHitAuxContainer_v1.h
new file mode 100644
index 0000000000000000000000000000000000000000..0ea61b39a567a61882a390aa35277363771e970f
--- /dev/null
+++ b/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHitAuxContainer_v1.h
@@ -0,0 +1,63 @@
+// Dear emacs, this is -*- c++ -*-
+
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+#ifndef XAODFASERWAVEFORM_VERSIONS_WAVEFORMHITAUXCONTAINER_V1_H
+#define XAODFASERWAVEFORM_VERSIONS_WAVEFORMHITAUXCONTAINER_V1_H
+
+// STL include(s):
+#include <vector>
+
+// EDM include(s):
+#include "xAODCore/AuxContainerBase.h"
+
+namespace xAOD {
+
+ /// Auxiliary container for WaveformHit containers
+
+ class WaveformHitAuxContainer_v1 : public AuxContainerBase {
+
+ public:
+ /// Default constructor
+ WaveformHitAuxContainer_v1();
+ /// Destructor
+ ~WaveformHitAuxContainer_v1() {}
+
+ private:
+ /// @name Basic variables
+ ///@ {
+ std::vector<unsigned int> channel;
+ std::vector<float> localtime;
+ std::vector<float> peak;
+ std::vector<float> width;
+ std::vector<float> integral;
+ std::vector<float> bcid_time;
+
+ std::vector<float> raw_peak;
+ std::vector<float> raw_integral;
+
+ std::vector<float> baseline_mean;
+ std::vector<float> baseline_rms;
+
+ std::vector<unsigned int> status;
+
+ std::vector<float> mean;
+ std::vector<float> alpha;
+ std::vector<float> nval;
+
+ std::vector<std::vector<float>> time_vector;
+ std::vector<std::vector<float>> wave_vector;
+
+ ///@}
+
+ }; // class WaveformHitAuxContainer_v1
+
+} // namespace xAOD
+
+// Set up a CLID and StoreGate inheritance for the class:
+#include "xAODCore/BaseInfo.h"
+SG_BASE( xAOD::WaveformHitAuxContainer_v1, xAOD::AuxContainerBase );
+
+#endif // XAODFASERWAVEFORM_VERSIONS_WAVEFORMHITAUXCONTAINER_V1_H
diff --git a/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHitContainer_v1.h b/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHitContainer_v1.h
new file mode 100644
index 0000000000000000000000000000000000000000..180cf9719cacaded61292aa67d401d38472016cd
--- /dev/null
+++ b/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHitContainer_v1.h
@@ -0,0 +1,26 @@
+// Dear emacs, this is -*- c++ -*-
+
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+#ifndef XAODFASERWAVEFORM_VERSIONS_WAVEFORMHITCONTAINER_V1_H
+#define XAODFASERWAVEFORM_VERSIONS_WAVEFORMHITCONTAINER_V1_H
+
+// System include(s):
+extern "C" {
+# include "stdint.h"
+}
+
+// EDM include(s):
+#include "AthContainers/DataVector.h"
+
+// Local includes:
+#include "xAODFaserWaveform/versions/WaveformHit_v1.h"
+
+namespace xAOD {
+ // Define the container as a simple DataVector
+ typedef DataVector<WaveformHit_v1> WaveformHitContainer_v1;
+}
+
+#endif // XAODFASERWAVEFORM_VERSIONS_WAVEFORMHITCONTAINER_V1_H
diff --git a/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHit_v1.h b/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHit_v1.h
new file mode 100644
index 0000000000000000000000000000000000000000..976c51512ae734902db9a5b3e30a309c5610723e
--- /dev/null
+++ b/xAOD/xAODFaserWaveform/xAODFaserWaveform/versions/WaveformHit_v1.h
@@ -0,0 +1,135 @@
+// Dear emacs, this is -*- c++ -*-
+
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+#ifndef XAODFASERWAVEFORM_VERSIONS_WAVEFORMHIT_V1_H
+#define XAODFASERWAVEFORM_VERSIONS_WAVEFORMHIT_V1_H
+
+// System include(s):
+extern "C" {
+# include "stdint.h"
+}
+
+#include <vector>
+
+// Core EDM include(s):
+#include "AthContainers/AuxElement.h"
+
+namespace xAOD {
+
+ // Meaning of status bits
+ enum WaveformStatus {
+ THRESHOLD_FAILED=0, // No hit found over threshold
+ SECONDARY, // Second or more hit found in raw waveform
+ WAVE_OVERFLOW, // Overflows removed from fit
+ BASELINE_FAILED, // Baseline determination failed
+ GFIT_FAILED, // Gaussian fit failed
+ CBFIT_FAILED, // CrystalBall fit failed
+ CLOCK_INVALID,
+ WAVEFORM_MISSING, // Input waveform data missing
+ WAVEFORM_INVALID // Input waveform data length mismatch
+ };
+
+ // Cllss describing pulses in the waveform digitizer
+ class WaveformHit_v1 : public SG::AuxElement {
+
+ public:
+ /// Defaullt constructor
+ WaveformHit_v1();
+
+ /// @name Access WaveformHit elements
+ /// @{
+
+ /// Waveform channel
+ unsigned int channel() const;
+ void set_channel(unsigned int value);
+
+ /// Best results
+ float localtime() const;
+ void set_localtime(float value);
+
+ float peak() const;
+ void set_peak(float value);
+
+ float width() const;
+ void set_width(float value);
+
+ float integral() const;
+ void set_integral(float value);
+
+ /// Time from previous clock edge
+ float bcid_time() const;
+ void set_bcid_time(float value);
+
+ /// Raw values from waveform
+ float raw_peak() const;
+ void set_raw_peak(float value);
+
+ float raw_integral() const;
+ void set_raw_integral(float value);
+
+ /// Bsaeline mean
+ float baseline_mean() const;
+ void set_baseline_mean(float value);
+
+ /// Baseline rms
+ float baseline_rms() const;
+ void set_baseline_rms(float value);
+
+ /// Status word
+ unsigned int status() const;
+ void set_status(unsigned int value);
+
+ /// Other fit results
+ float mean() const;
+ void set_mean(float value);
+
+ /// Crystal Ball fit parameters
+ float alpha() const;
+ void set_alpha(float value);
+
+ float nval() const;
+ void set_nval(float value);
+
+ /// Raw time and waveform data
+ std::vector<float> time_vector() const;
+ void set_time_vector(std::vector<float> value);
+
+ std::vector<float> wave_vector() const;
+ void set_wave_vector(std::vector<float> value);
+
+ /// Status bit access functions
+ void set_status_bit(WaveformStatus bit) {
+ this->set_status(this->status() | (1<<bit));
+ }
+ bool status_bit(WaveformStatus bit) const {
+ return (this->status() | (1<<bit));
+ }
+
+ bool threshold() const {
+ return !(this->status_bit(xAOD::WaveformStatus::THRESHOLD_FAILED));
+ }
+
+ bool overflow() const {
+ return (this->status_bit(xAOD::WaveformStatus::WAVE_OVERFLOW));
+ }
+
+ bool secondary() const {
+ return this->status_bit(xAOD::WaveformStatus::SECONDARY);
+ }
+
+ /// @}
+
+ }; // class WaveformHit_v1
+
+ std::ostream& operator<<(std::ostream& s, const xAOD::WaveformHit_v1& hit);
+}
+
+// Declare the inheritance of the type:
+#include "xAODCore/BaseInfo.h"
+SG_BASE( xAOD::WaveformHit_v1, SG::AuxElement );
+
+
+#endif // XAODFASERWAVEFORM_VERSIONS_WAVEFORMHIT_V1_H
diff --git a/xAOD/xAODFaserWaveform/xAODFaserWaveform/xAODFaserWaveformDict.h b/xAOD/xAODFaserWaveform/xAODFaserWaveform/xAODFaserWaveformDict.h
new file mode 100644
index 0000000000000000000000000000000000000000..5bb2eceb5665288f05ed3481de0d201d825989ff
--- /dev/null
+++ b/xAOD/xAODFaserWaveform/xAODFaserWaveform/xAODFaserWaveformDict.h
@@ -0,0 +1,33 @@
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+#ifndef XAODFASERWAVEFORM_XAODFASERWAVEFORMDICT_H
+#define XAODFASERWAVEFORM_XAODFASERWAVEFORMDICT_H
+
+// Local includes
+#include "xAODFaserWaveform/WaveformHit.h"
+#include "xAODFaserWaveform/WaveformHitContainer.h"
+#include "xAODFaserWaveform/WaveformHitAuxContainer.h"
+
+#include "xAODFaserWaveform/WaveformClock.h"
+#include "xAODFaserWaveform/WaveformClockAuxInfo.h"
+
+#include "xAODFaserWaveform/versions/WaveformHit_v1.h"
+#include "xAODFaserWaveform/versions/WaveformHitContainer_v1.h"
+#include "xAODFaserWaveform/versions/WaveformHitAuxContainer_v1.h"
+
+#include "xAODFaserWaveform/versions/WaveformClock_v1.h"
+#include "xAODFaserWaveform/versions/WaveformClockAuxInfo_v1.h"
+
+// EDM include(s).
+#include "xAODCore/tools/DictHelpers.h"
+
+namespace {
+ struct GCCXML_DUMMY_INSTANTIATION_XAODFASERWAVEFORM {
+ XAOD_INSTANTIATE_NS_CONTAINER_TYPES( xAOD, WaveformHitContainer_v1 );
+ XAOD_INSTANTIATE_NS_OBJECT_TYPES( xAOD, WaveformClock_v1 );
+ };
+}
+
+#endif // XAODFASERWAVEFORM_XAODFASERWAVEFORMDICT_H
diff --git a/xAOD/xAODFaserWaveformAthenaPool/CMakeLists.txt b/xAOD/xAODFaserWaveformAthenaPool/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..99f9ee00fdbfd45e353401e1c6480e1ec0a0d370
--- /dev/null
+++ b/xAOD/xAODFaserWaveformAthenaPool/CMakeLists.txt
@@ -0,0 +1,14 @@
+# Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+
+# Declare the package name.
+atlas_subdir( xAODFaserWaveformAthenaPool )
+
+# Component(s) in the package:
+atlas_add_poolcnv_library( xAODFaserWaveformAthenaPoolPoolCnv
+ src/*.h src/*.cxx
+ FILES xAODFaserWaveform/WaveformClock.h xAODFaserWaveform/WaveformClockAuxInfo.h xAODFaserWaveform/WaveformHitContainer.h xAODFaserWaveform/WaveformHitAuxContainer.h
+ TYPES_WITH_NAMESPACE xAOD::WaveformClock xAOD::WaveformClockAuxInfo xAOD::WaveformHitContainer xAOD::WaveformHitAuxContainer
+ CNV_PFX xAOD
+ LINK_LIBRARIES AthenaPoolCnvSvcLib AthenaPoolUtilities xAODFaserWaveform )
+
+
diff --git a/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformClockAuxInfoCnv.cxx b/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformClockAuxInfoCnv.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..b5024b1ae3058eccec0eaa0f4ab4e349cecf465f
--- /dev/null
+++ b/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformClockAuxInfoCnv.cxx
@@ -0,0 +1,6 @@
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+// Dummy source file so that cmake will know this is a custom converter.
+// xAODWaveformClockAuxInfoCnv.cxx
diff --git a/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformClockAuxInfoCnv.h b/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformClockAuxInfoCnv.h
new file mode 100644
index 0000000000000000000000000000000000000000..5ac8b4c0db48ff15f3647f71f37f072682e55fa8
--- /dev/null
+++ b/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformClockAuxInfoCnv.h
@@ -0,0 +1,15 @@
+// Dear emacs, this is -*- c++ -*-
+
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+#ifndef XAODFASERWAVEFORMATHENAPOOL_XAODWAVEFORMCLOCKAUXINFOCNV_H
+#define XAODFASERWAVEFORMATHENAPOOL_XAODWAVEFORMCLOCKAUXINFOCNV_H
+
+#include "xAODFaserWaveform/WaveformClockAuxInfo.h"
+#include "AthenaPoolCnvSvc/T_AthenaPoolAuxContainerCnv.h"
+
+typedef T_AthenaPoolAuxContainerCnv<xAOD::WaveformClockAuxInfo> xAODWaveformClockAuxInfoCnv;
+
+#endif // XAODFASERWAVEFORMATHENAPOOL_XAODFASERWAVEFORMCLOCKAUXINFOCNV_H
diff --git a/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformClockCnv.cxx b/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformClockCnv.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..a175dccc75b115fad4d53637e900d97d67552435
--- /dev/null
+++ b/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformClockCnv.cxx
@@ -0,0 +1,6 @@
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+// Dummy source file so that cmake will know this is a custom converter.
+// xAODWaveformClockCnv.cxx
diff --git a/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformClockCnv.h b/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformClockCnv.h
new file mode 100644
index 0000000000000000000000000000000000000000..2f591d164011f6ec3c4e83b5a89d55b0fc725c68
--- /dev/null
+++ b/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformClockCnv.h
@@ -0,0 +1,15 @@
+// Dear emacs, this is -*- c++ -*-
+
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+#ifndef XAODFASERWAVEFORMATHENAPOOL_XAODWAVEFORMCLOCKCNV_H
+#define XAODFASERWAVEFORMATHENAPOOL_XAODWAVEFORMCLOCKCNV_H
+
+#include "xAODFaserWaveform/WaveformClock.h"
+#include "AthenaPoolCnvSvc/T_AthenaPoolxAODCnv.h"
+
+typedef T_AthenaPoolxAODCnv<xAOD::WaveformClock> xAODWaveformClockCnv;
+
+#endif // XAODFASERWAVEFORMATHENAPOOL_XAODWAVEFORMCLOCKCNV_H
diff --git a/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformHitAuxContainerCnv.cxx b/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformHitAuxContainerCnv.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..31d86ebb89db520de8029ae21581a7ffca95ddd0
--- /dev/null
+++ b/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformHitAuxContainerCnv.cxx
@@ -0,0 +1,6 @@
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+// Dummy source file so that cmake will know this is a custom converter.
+// xAODWaveformHitAuxContainerCnv.cxx
diff --git a/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformHitAuxContainerCnv.h b/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformHitAuxContainerCnv.h
new file mode 100644
index 0000000000000000000000000000000000000000..bffe321a489859d89c98a13d29e24ea4a5d7c338
--- /dev/null
+++ b/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformHitAuxContainerCnv.h
@@ -0,0 +1,15 @@
+// Dear emacs, this is -*- c++ -*-
+
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+#ifndef XAODFASERWAVEFORMATHENAPOOL_XAODWAVEFORMHITAUXCONTAINERCNV_H
+#define XAODFASERWAVEFORMATHENAPOOL_XAODWAVEFORMHITAUXCONTAINERCNV_H
+
+#include "xAODFaserWaveform/WaveformHitAuxContainer.h"
+#include "AthenaPoolCnvSvc/T_AthenaPoolAuxContainerCnv.h"
+
+typedef T_AthenaPoolAuxContainerCnv<xAOD::WaveformHitAuxContainer> xAODWaveformHitAuxContainerCnv;
+
+#endif // XAODFASERWAVEFORMATHENAPOOL_XAODFASERWAVEFORMHITAUXCONTAINERCNV_H
diff --git a/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformHitContainerCnv.cxx b/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformHitContainerCnv.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..2cd6731a4992cdc4b8ec2be6fcc7eeded4690727
--- /dev/null
+++ b/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformHitContainerCnv.cxx
@@ -0,0 +1,6 @@
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+// Dummy source file so that cmake will know this is a custom converter.
+// xAODWaveformHitContainerCnv.cxx
diff --git a/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformHitContainerCnv.h b/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformHitContainerCnv.h
new file mode 100644
index 0000000000000000000000000000000000000000..770bb1625a210f23efb850bda8223a06e6496b5e
--- /dev/null
+++ b/xAOD/xAODFaserWaveformAthenaPool/src/xAODWaveformHitContainerCnv.h
@@ -0,0 +1,15 @@
+// Dear emacs, this is -*- c++ -*-
+
+/*
+ Copyright (C) 2020 CERN for the benefit of the FASER collaboration
+*/
+
+#ifndef XAODFASERWAVEFORMATHENAPOOL_XAODWAVEFORMHITCONTAINERCNV_H
+#define XAODFASERWAVEFORMATHENAPOOL_XAODWAVEFORMHITCONTAINERCNV_H
+
+#include "xAODFaserWaveform/WaveformHitContainer.h"
+#include "AthenaPoolCnvSvc/T_AthenaPoolxAODCnv.h"
+
+typedef T_AthenaPoolxAODCnv<xAOD::WaveformHitContainer> xAODWaveformHitContainerCnv;
+
+#endif // XAODFASERWAVEFORMATHENAPOOL_XAODWAVEFORMHITCONTAINERCNV_H