diff --git a/Calorimeter/CaloDigiAlgs/CMakeLists.txt b/Calorimeter/CaloDigiAlgs/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..42991aa03cc03af6e51ab161d36d1d61df1a7956 --- /dev/null +++ b/Calorimeter/CaloDigiAlgs/CMakeLists.txt @@ -0,0 +1,15 @@ +################################################################################ +# Package: CaloDigiAlgs +################################################################################ + +# Declare the package name: +atlas_subdir( CaloDigiAlgs ) + +# Component(s) in the package: +atlas_add_component( CaloDigiAlgs + src/*.cxx src/*.h + src/components/*.cxx + LINK_LIBRARIES AthenaBaseComps Identifier StoreGateLib WaveRawEvent FaserCaloSimEvent WaveDigiToolsLib) + +atlas_install_python_modules( python/*.py ) + diff --git a/Calorimeter/CaloDigiAlgs/python/CaloDigiAlgsConfig.py b/Calorimeter/CaloDigiAlgs/python/CaloDigiAlgsConfig.py new file mode 100644 index 0000000000000000000000000000000000000000..6a0b9e35717e92eabed64b111b811ddf6283c85d --- /dev/null +++ b/Calorimeter/CaloDigiAlgs/python/CaloDigiAlgsConfig.py @@ -0,0 +1,52 @@ +# Copyright (C) 2020-2021 CERN for the benefit of the FASER collaboration + +from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator +from AthenaConfiguration.ComponentFactory import CompFactory + +from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg + + +# One stop shopping for normal FASER data +def CaloWaveformDigitizationCfg(flags): + """ Return all algorithms and tools for Waveform digitization """ + acc = ComponentAccumulator() + + if not flags.Input.isMC: + return + + acc.merge(CaloWaveformDigiCfg(flags, "CaloWaveformDigiAlg")) + return acc + +# Return configured digitization algorithm from SIM hits +def CaloWaveformDigiCfg(flags, name="CaloWaveformDigiAlg", **kwargs): + + acc = ComponentAccumulator() + + tool = CompFactory.WaveformDigitisationTool(name="CaloWaveformDigtisationTool", **kwargs) + + kwargs.setdefault("CaloHitContainerKey", "EcalHits") + kwargs.setdefault("WaveformContainerKey", "CaloWaveforms") + + digiAlg = CompFactory.CaloWaveformDigiAlg(name, **kwargs) + kwargs.setdefault("WaveformDigitisationTool", tool) + + digiAlg.CB_alpha = -0.9 + digiAlg.CB_n = 10 + digiAlg.CB_sigma = 4 + digiAlg.CB_mean = 820 + + + acc.addEventAlgo(digiAlg) + + return acc + +def CaloWaveformDigitizationOutputCfg(flags, **kwargs): + """ Return ComponentAccumulator with output for Waveform Digi""" + acc = ComponentAccumulator() + ItemList = [ + "RawWaveformContainer#*" + ] + acc.merge(OutputStreamCfg(flags, "RDO", ItemList)) + # ostream = acc.getEventAlgo("OutputStreamRDO") + # ostream.TakeItemsFromInput = True # Don't know what this does + return acc diff --git a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx new file mode 100644 index 0000000000000000000000000000000000000000..718dcb4f398aa2ad110510c622dd159e81969a54 --- /dev/null +++ b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx @@ -0,0 +1,84 @@ +#include "CaloWaveformDigiAlg.h" + +#include "Identifier/Identifier.h" + +#include <vector> +#include <map> + +CaloWaveformDigiAlg::CaloWaveformDigiAlg(const std::string& name, + ISvcLocator* pSvcLocator) + : AthReentrantAlgorithm(name, pSvcLocator) { + +} + +StatusCode +CaloWaveformDigiAlg::initialize() { + ATH_MSG_INFO(name() << "::initalize()" ); + + // Initalize tools + ATH_CHECK( m_digiTool.retrieve() ); + + + // Set key to read waveform from + ATH_CHECK( m_caloHitContainerKey.initialize() ); + + // Set key to write container + ATH_CHECK( m_waveformContainerKey.initialize() ); + + // Will eventually depend on the type of detector + // TODO: Vary time at which centre it? + // TODO: Change params compared to scint + // m_kernel = new TF1("PDF", " ROOT::Math::crystalball_pdf(x, -0.9, 10, 4, 900)", 0, 1200); + + m_kernel = new TF1("PDF", "ROOT::Math::crystalball_pdf(x, [0],[1],[2],[3])", 0, 1200); + //m_kernel->SetParameters(-0.25,10,4,900); + m_kernel->SetParameter(0, m_CB_alpha); + m_kernel->SetParameter(1, m_CB_n); + m_kernel->SetParameter(2, m_CB_sigma); + m_kernel->SetParameter(3, m_CB_mean); + + + return StatusCode::SUCCESS; +} + +StatusCode +CaloWaveformDigiAlg::finalize() { + ATH_MSG_INFO(name() << "::finalize()"); + + delete m_kernel; + + return StatusCode::SUCCESS; +} + +StatusCode +CaloWaveformDigiAlg::execute(const EventContext& ctx) const { + ATH_MSG_DEBUG("Executing"); + + ATH_MSG_DEBUG("Run: " << ctx.eventID().run_number() + << " Event: " << ctx.eventID().event_number()); + + // Find the input HIT collection + SG::ReadHandle<CaloHitCollection> caloHitHandle(m_caloHitContainerKey, ctx); + + ATH_CHECK( caloHitHandle.isValid() ); + ATH_MSG_DEBUG("Found ReadHandle for CaloHitCollection " << m_caloHitContainerKey); + + if (caloHitHandle->size() == 0) { + ATH_MSG_DEBUG("CaloHitCollection found with zero length!"); + return StatusCode::SUCCESS; + } + + // Find the output waveform container + SG::WriteHandle<RawWaveformContainer> waveformContainerHandle(m_waveformContainerKey, ctx); + ATH_CHECK( waveformContainerHandle.record( std::make_unique<RawWaveformContainer>()) ); + + ATH_MSG_DEBUG("WaveformsContainer '" << waveformContainerHandle.name() << "' initialized"); + + // Digitise the hits + CHECK( m_digiTool->digitise<CaloHitCollection>(caloHitHandle.ptr(), + waveformContainerHandle.ptr(), m_kernel) ); + + ATH_MSG_DEBUG("WaveformsHitContainer " << waveformContainerHandle.name() << "' filled with "<< waveformContainerHandle->size() <<" items"); + + return StatusCode::SUCCESS; +} diff --git a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h new file mode 100644 index 0000000000000000000000000000000000000000..9de68257e540890d29143df6a054578334f12503 --- /dev/null +++ b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h @@ -0,0 +1,92 @@ +#ifndef CALODIGIALGS_CALOWAVEFORMDIGIALG_H +#define CALODIGIALGS_CALOWAVEFORMDIGIALG_H + +// Base class +#include "AthenaBaseComps/AthReentrantAlgorithm.h" + +// Data classes +#include "WaveRawEvent/RawWaveformContainer.h" +#include "FaserCaloSimEvent/CaloHitCollection.h" + +// Tool classes +#include "WaveDigiTools/IWaveformDigitisationTool.h" + +// Handles +#include "StoreGate/ReadHandleKey.h" +#include "StoreGate/WriteHandleKey.h" + +// Gaudi +#include "GaudiKernel/ServiceHandle.h" +#include "GaudiKernel/ToolHandle.h" + +// ROOT +#include "TF1.h" + +// STL +#include <string> + +class CaloWaveformDigiAlg : public AthReentrantAlgorithm { + + public: + // Constructor + CaloWaveformDigiAlg(const std::string& name, ISvcLocator* pSvcLocator); + virtual ~CaloWaveformDigiAlg() = default; + + /** @name Usual algorithm methods */ + //@{ + virtual StatusCode initialize() override; + virtual StatusCode execute(const EventContext& ctx) const override; + virtual StatusCode finalize() override; + //@} + + private: + + /** @name Disallow default instantiation, copy, assignment */ + //@{ + CaloWaveformDigiAlg() = delete; + CaloWaveformDigiAlg(const CaloWaveformDigiAlg&) = delete; + CaloWaveformDigiAlg &operator=(const CaloWaveformDigiAlg&) = delete; + //@} + + Gaudi::Property<double> m_CB_alpha {this, "CB_alpha", 0, "Alpha of the crystal ball function"}; + Gaudi::Property<double> m_CB_n {this, "CB_n", 0, "n of the crystal ball function"}; + Gaudi::Property<double> m_CB_mean {this, "CB_mean", 0, "Mean of the crystal ball function"}; + Gaudi::Property<double> m_CB_sigma {this, "CB_sigma", 0, "Sigma of the crystal ball function"}; + + + + + + /// Kernel PDF + TF1* m_kernel; + + + /** + * @name Digitisation tool + */ + ToolHandle<IWaveformDigitisationTool> m_digiTool + {this, "WaveformDigitisationTool", "WaveformDigitisationTool"}; + + + /** + * @name Input HITS using SG::ReadHandleKey + */ + //@{ + + SG::ReadHandleKey<CaloHitCollection> m_caloHitContainerKey + {this, "CaloHitContainerKey", ""}; + + //@} + + + /** + * @name Output data using SG::WriteHandleKey + */ + //@{ + SG::WriteHandleKey<RawWaveformContainer> m_waveformContainerKey + {this, "WaveformContainerKey", ""}; + //@} + +}; + +#endif // CALODIGIALGS_CALODIGIALG_H diff --git a/Calorimeter/CaloDigiAlgs/src/components/CaloDigiAlgs_entries.cxx b/Calorimeter/CaloDigiAlgs/src/components/CaloDigiAlgs_entries.cxx new file mode 100644 index 0000000000000000000000000000000000000000..832cfbb096ffdb13430ba92dc2bfc7b1125bd53c --- /dev/null +++ b/Calorimeter/CaloDigiAlgs/src/components/CaloDigiAlgs_entries.cxx @@ -0,0 +1,2 @@ +#include "../CaloWaveformDigiAlg.h" +DECLARE_COMPONENT( CaloWaveformDigiAlg ) diff --git a/Scintillator/ScintDigiAlgs/CMakeLists.txt b/Scintillator/ScintDigiAlgs/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..43695c5f90d96053a8b44fa3ac41c0f4c159ad99 --- /dev/null +++ b/Scintillator/ScintDigiAlgs/CMakeLists.txt @@ -0,0 +1,15 @@ +################################################################################ +# Package: ScintDigiAlgs +################################################################################ + +# Declare the package name: +atlas_subdir( ScintDigiAlgs ) + +# Component(s) in the package: +atlas_add_component( ScintDigiAlgs + src/*.cxx src/*.h + src/components/*.cxx + LINK_LIBRARIES AthenaBaseComps Identifier StoreGateLib WaveRawEvent ScintSimEvent WaveDigiToolsLib) + +atlas_install_python_modules( python/*.py ) + diff --git a/Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py b/Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py new file mode 100644 index 0000000000000000000000000000000000000000..7c322060ab68f8580cc6b517743edf46209dee31 --- /dev/null +++ b/Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py @@ -0,0 +1,64 @@ +# Copyright (C) 2020-2021 CERN for the benefit of the FASER collaboration + +from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator +from AthenaConfiguration.ComponentFactory import CompFactory + +from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg + +# Crystallball parameter dictionary used in simulated digitized wave reconstruction. +# Crystalball function Parameters estimated from Deion's slides uploaded at +# https://indico.cern.ch/event/1099652/contributions/4626975/attachments/2352595/4013927/Faser-Physics-run3933-plots.pdf (20/01/2022) +# Parameters are per scintillator source, but not per channel. +dict_CB_param = {} +dict_CB_param["Trigger"]=dict(CB_alpha=-0.38, CB_n=25, CB_mean=815, CB_sigma=7.7) +dict_CB_param["Timing"]=dict(CB_alpha=-0.32, CB_n=65, CB_mean=846, CB_sigma=5.3) # copy from Preshower; Timing was not in TestBeam +dict_CB_param["Veto"]=dict(CB_alpha=-0.38, CB_n=25, CB_mean=815, CB_sigma=7.7) # copy from Trigger; Veto was not in TestBeam, but in sim "Veto" is the TestBeam Trigger component +dict_CB_param["Preshower"]=dict(CB_alpha=-0.32, CB_n=65, CB_mean=846, CB_sigma=5.3) + +# One stop shopping for normal FASER data +def ScintWaveformDigitizationCfg(flags): + """ Return all algorithms and tools for Waveform digitization """ + acc = ComponentAccumulator() + + if not flags.Input.isMC: + return + + if "TB" not in flags.GeoModel.FaserVersion: + acc.merge(ScintWaveformDigiCfg(flags, "TimingWaveformDigiAlg", "Trigger")) + acc.merge(ScintWaveformDigiCfg(flags, "VetoWaveformDigiAlg", "Veto")) + acc.merge(ScintWaveformDigiCfg(flags, "PreshowerWaveformDigiAlg", "Preshower")) + return acc + +# Return configured digitization algorithm from SIM hits +# Specify data source (Veto, Trigger, Preshower) +def ScintWaveformDigiCfg(flags, name="ScintWaveformDigiAlg", source="", **kwargs): + + acc = ComponentAccumulator() + + tool = CompFactory.WaveformDigitisationTool(name=source+"WaveformDigtisationTool", **kwargs) + + kwargs.setdefault("ScintHitContainerKey", source+"Hits") + kwargs.setdefault("WaveformContainerKey", source+"Waveforms") + + digiAlg = CompFactory.ScintWaveformDigiAlg(name, **kwargs) + digiAlg.CB_alpha = dict_CB_param[source]["CB_alpha"] + digiAlg.CB_n = dict_CB_param[source]["CB_n"] + digiAlg.CB_mean = dict_CB_param[source]["CB_mean"] + digiAlg.CB_sigma = dict_CB_param[source]["CB_sigma"] + + kwargs.setdefault("WaveformDigitisationTool", tool) + + acc.addEventAlgo(digiAlg) + + return acc + +def ScintWaveformDigitizationOutputCfg(flags, **kwargs): + """ Return ComponentAccumulator with output for Waveform Digi""" + acc = ComponentAccumulator() + ItemList = [ + "RawWaveformContainer#*" + ] + acc.merge(OutputStreamCfg(flags, "RDO", ItemList)) + # ostream = acc.getEventAlgo("OutputStreamRDO") + # ostream.TakeItemsFromInput = True # Don't know what this does + return acc diff --git a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx new file mode 100644 index 0000000000000000000000000000000000000000..03c1873d11d4ecd83d2d8ad22eef07a76b7338a2 --- /dev/null +++ b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx @@ -0,0 +1,85 @@ + +#include "ScintWaveformDigiAlg.h" + +#include "Identifier/Identifier.h" + +#include <vector> +#include <map> + + +ScintWaveformDigiAlg::ScintWaveformDigiAlg(const std::string& name, + ISvcLocator* pSvcLocator) + : AthReentrantAlgorithm(name, pSvcLocator) { + +} + +StatusCode +ScintWaveformDigiAlg::initialize() { + ATH_MSG_INFO(name() << "::initalize()" ); + + // Initalize tools + ATH_CHECK( m_digiTool.retrieve() ); + + + // Set key to read waveform from + ATH_CHECK( m_scintHitContainerKey.initialize() ); + + // Set key to write container + ATH_CHECK( m_waveformContainerKey.initialize() ); + + // Will eventually depend on the type of detector + // TODO: Vary time at which centre it? + // TODO: Better parameters + + + + m_kernel = new TF1("PDF", "ROOT::Math::crystalball_pdf(x, [0],[1],[2],[3])", 0, 1200); + m_kernel->SetParameter(0, m_CB_alpha); + m_kernel->SetParameter(1, m_CB_n); + m_kernel->SetParameter(2, m_CB_sigma); + m_kernel->SetParameter(3, m_CB_mean); + + return StatusCode::SUCCESS; +} + +StatusCode +ScintWaveformDigiAlg::finalize() { + ATH_MSG_INFO(name() << "::finalize()"); + + delete m_kernel; + + return StatusCode::SUCCESS; +} + +StatusCode +ScintWaveformDigiAlg::execute(const EventContext& ctx) const { + ATH_MSG_DEBUG("Executing"); + + ATH_MSG_DEBUG("Run: " << ctx.eventID().run_number() + << " Event: " << ctx.eventID().event_number()); + + // Find the input HIT collection + SG::ReadHandle<ScintHitCollection> scintHitHandle(m_scintHitContainerKey, ctx); + + ATH_CHECK( scintHitHandle.isValid() ); + ATH_MSG_DEBUG("Found ReadHandle for ScintHitCollection " << m_scintHitContainerKey); + + if (scintHitHandle->size() == 0) { + ATH_MSG_DEBUG("ScintHitCollection found with zero length!"); + return StatusCode::SUCCESS; + } + + // Find the output waveform container + SG::WriteHandle<RawWaveformContainer> waveformContainerHandle(m_waveformContainerKey, ctx); + ATH_CHECK( waveformContainerHandle.record( std::make_unique<RawWaveformContainer>()) ); + + ATH_MSG_DEBUG("WaveformsContainer '" << waveformContainerHandle.name() << "' initialized"); + + // Digitise the hits + CHECK( m_digiTool->digitise<ScintHitCollection>(scintHitHandle.ptr(), + waveformContainerHandle.ptr(), m_kernel) ); + + ATH_MSG_DEBUG("WaveformsHitContainer " << waveformContainerHandle.name() << "' filled with "<< waveformContainerHandle->size() <<" items"); + + return StatusCode::SUCCESS; +} diff --git a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h new file mode 100644 index 0000000000000000000000000000000000000000..af60e02d8169194ccc2459f0b7a69727b2544eba --- /dev/null +++ b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h @@ -0,0 +1,91 @@ +#ifndef SCINTDIGIALGS_SCINTWAVEFORMDIGIALG_H +#define SCINTDIGIALGS_SCINTWAVEFORMDIGIALG_H + +// Base class +#include "AthenaBaseComps/AthReentrantAlgorithm.h" + +// Data classes +#include "WaveRawEvent/RawWaveformContainer.h" +#include "ScintSimEvent/ScintHitCollection.h" + +// Tool classes +#include "WaveDigiTools/IWaveformDigitisationTool.h" + +// Handles +#include "StoreGate/ReadHandleKey.h" +#include "StoreGate/WriteHandleKey.h" + +// Gaudi +#include "GaudiKernel/ServiceHandle.h" +#include "GaudiKernel/ToolHandle.h" + +// ROOT +#include "TF1.h" + + +// STL +#include <string> + +class ScintWaveformDigiAlg : public AthReentrantAlgorithm { + + public: + // Constructor + ScintWaveformDigiAlg(const std::string& name, ISvcLocator* pSvcLocator); + virtual ~ScintWaveformDigiAlg() = default; + + /** @name Usual algorithm methods */ + //@{ + virtual StatusCode initialize() override; + virtual StatusCode execute(const EventContext& ctx) const override; + virtual StatusCode finalize() override; + //@} + + + private: + + /** @name Disallow default instantiation, copy, assignment */ + //@{ + ScintWaveformDigiAlg() = delete; + ScintWaveformDigiAlg(const ScintWaveformDigiAlg&) = delete; + ScintWaveformDigiAlg &operator=(const ScintWaveformDigiAlg&) = delete; + //@} + + Gaudi::Property<double> m_CB_alpha {this, "CB_alpha", 0, "Alpha of the crystal ball function"}; + Gaudi::Property<double> m_CB_n {this, "CB_n", 0, "n of the crystal ball function"}; + Gaudi::Property<double> m_CB_mean {this, "CB_mean", 0, "Mean of the crystal ball function"}; + Gaudi::Property<double> m_CB_sigma {this, "CB_sigma", 0, "Sigma of the crystal ball function"}; + + + /// Kernel PDF + TF1* m_kernel; + + + /** + * @name Digitisation tool + */ + ToolHandle<IWaveformDigitisationTool> m_digiTool + {this, "WaveformDigitisationTool", "WaveformDigitisationTool"}; + + + /** + * @name Input HITS using SG::ReadHandleKey + */ + //@{ + + SG::ReadHandleKey<ScintHitCollection> m_scintHitContainerKey + {this, "ScintHitContainerKey", ""}; + + //@} + + + /** + * @name Output data using SG::WriteHandleKey + */ + //@{ + SG::WriteHandleKey<RawWaveformContainer> m_waveformContainerKey + {this, "WaveformContainerKey", ""}; + //@} + +}; + +#endif // SCINTDIGIALGS_SCINTDIGIALG_H diff --git a/Scintillator/ScintDigiAlgs/src/components/ScintDigiAlgs_entries.cxx b/Scintillator/ScintDigiAlgs/src/components/ScintDigiAlgs_entries.cxx new file mode 100644 index 0000000000000000000000000000000000000000..3a8deff1f06b692e1cf0928472fb12426274f059 --- /dev/null +++ b/Scintillator/ScintDigiAlgs/src/components/ScintDigiAlgs_entries.cxx @@ -0,0 +1,2 @@ +#include "../ScintWaveformDigiAlg.h" +DECLARE_COMPONENT( ScintWaveformDigiAlg ) diff --git a/Waveform/WaveDigiTools/CMakeLists.txt b/Waveform/WaveDigiTools/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..692fdb69bc14451ba5a6a660d011bad5b14b66e5 --- /dev/null +++ b/Waveform/WaveDigiTools/CMakeLists.txt @@ -0,0 +1,25 @@ +################################################################################ +# Package: WaveDigiTools +################################################################################ + +# Declare the package name: +atlas_subdir( WaveDigiTools ) + +# External dependencies: +find_package( ROOT ) + +# Component(s) in the package: +atlas_add_library( WaveDigiToolsLib + WaveDigiTools/*.h src/*.cxx src/*.h + PUBLIC_HEADERS WaveDigiTools + PRIVATE_INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} + LINK_LIBRARIES AthenaBaseComps AthenaKernel GeoPrimitives WaveRawEvent + PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES} + ) + +atlas_add_component( WaveDigiTools + src/components/*.cxx + INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} + LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps GaudiKernel WaveDigiToolsLib ) + + diff --git a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h new file mode 100644 index 0000000000000000000000000000000000000000..3e85ff839a04f6d02c794fa74567841450c9cfeb --- /dev/null +++ b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h @@ -0,0 +1,55 @@ +/* + Copyright (C) 2021 CERN for the benefit of the FASER collaboration +*/ + +/** + * @file IWaveformDigitisationTool.h + * Header file for the IWaveformDigitisationTool class + * @author Carl Gwilliam, 2021 + */ + + +#ifndef WAVEDIGITOOLS_IWAVEFORMDIGITISATIONTOOL_H +#define WAVEDIGITOOLS_IWAVEFORMDIGITISATIONTOOL_H + +// Base class +#include "GaudiKernel/IAlgTool.h" +#include "GaudiKernel/ToolHandle.h" + +#include "GaudiKernel/ServiceHandle.h" +#include "GaudiKernel/IMessageSvc.h" +#include "GaudiKernel/MsgStream.h" + +#include "WaveRawEvent/RawWaveformContainer.h" +#include "WaveRawEvent/RawWaveform.h" + +#include "TF1.h" + +///Interface for waveform digitisation tools +class IWaveformDigitisationTool : virtual public IAlgTool +{ +public: + + // InterfaceID + DeclareInterfaceID(IWaveformDigitisationTool, 1, 0); + + IWaveformDigitisationTool(): + m_msgSvc ( "MessageSvc", "ITrkEventCnvTool" ) + {} + + virtual ~IWaveformDigitisationTool() = default; + + // Digitise HITS to Raw waveform + template<class CONT> + StatusCode digitise(const CONT* hitCollection, + RawWaveformContainer* waveContainer, TF1* kernel) const; + +private: + ServiceHandle<IMessageSvc> m_msgSvc; + +}; + +#include "WaveDigiTools/IWaveformDigitisationTool.icc" + + +#endif //WAVEDIGITOOLS_IWAVEFORMDIGITISATIONTOOL_H diff --git a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.icc b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.icc new file mode 100644 index 0000000000000000000000000000000000000000..57d4839bda5f2d30286ec412e3f01c92ce353b11 --- /dev/null +++ b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.icc @@ -0,0 +1,51 @@ +#include <vector> +#include <map> + +template<class CONT> +StatusCode IWaveformDigitisationTool::digitise(const CONT* hitCollection, + RawWaveformContainer* container, TF1* kernel) const { + + + // Check the container + if (!container) { + MsgStream log(&(*m_msgSvc), name()); + log << MSG::ERROR << "HitCollection passed to digitise() is null!" << endmsg; + return StatusCode::FAILURE; + } + + unsigned int size = 600; // TODO: how know the correct number of time samples? + std::vector<float> time(size); + for (unsigned int i=0; i<size; i++) time[i] = 2.*i; + + std::map<unsigned int, std::vector<uint16_t>> waveforms; + unsigned int baseline = 8000; // TODO: vary this + add noise + + // Loop over time samples + for (const auto& t : time) { + std::map<unsigned int, float> counts; + + // Convolve hit energy with kernel and sum for each ID (i.e. channel) + for (const auto& hit : *hitCollection) { + counts[hit.identify()] += kernel->Eval(t) * hit.energyLoss(); + //std::cout << "HIT " << hit.identify() << " @ " << t << ": " << kernel->Eval(t) << " " << hit.energyLoss() << " -> " << counts[hit.identify()] << std::endl; + } + + // Add count to correct waveform vec + for (const auto& c : counts) { + waveforms[c.first].push_back(baseline - c.second); + //std::cout << "ADC " << c.first << " @ " << t << ": " << baseline - c.second << std::endl; + } + } + + // Loop over wavefrom vecs to make and store waveform + for (const auto& w : waveforms) { + RawWaveform* wfm = new RawWaveform(); + wfm->setWaveform(0, w.second); + wfm->setIdentifier(Identifier(w.first)); + wfm->setSamples(size); + container->push_back(wfm); + } + + + return StatusCode::SUCCESS; +} diff --git a/Waveform/WaveDigiTools/share/WaveformDigiAndRecoExample_jobOptions.py b/Waveform/WaveDigiTools/share/WaveformDigiAndRecoExample_jobOptions.py new file mode 100644 index 0000000000000000000000000000000000000000..779031828c9c5736f9403e403514e19285029499 --- /dev/null +++ b/Waveform/WaveDigiTools/share/WaveformDigiAndRecoExample_jobOptions.py @@ -0,0 +1,87 @@ +#!/usr/bin/env python + +import sys + +if __name__ == "__main__": + + fileroot = "output" + if len(sys.argv) > 1: + filename = sys.argv[1] + + doRDO = False + if len(sys.argv) > 2: + filename = sys.argv[2] + + from AthenaCommon.Logging import log, logging + from AthenaCommon.Constants import DEBUG, VERBOSE, INFO + from CalypsoConfiguration.AllConfigFlags import ConfigFlags + from AthenaCommon.Configurable import Configurable + + Configurable.configurableRun3Behavior = True + + log.setLevel(VERBOSE) + + ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" # Always needed; must match FaserVersion + ConfigFlags.IOVDb.DatabaseInstance = "OFLP200" # Use MC conditions for now + ConfigFlags.Input.ProjectName = "mc21" # Needed to bypass autoconfig + ConfigFlags.Input.isMC = True # Needed to bypass autoconfig + ConfigFlags.GeoModel.FaserVersion = "FASER-TB00" # FASER geometry + ConfigFlags.Common.isOnline = False + ConfigFlags.GeoModel.Align.Dynamic = False + + ConfigFlags.Input.Files = [ + "/eos/project-f/faser-commissioning/Simulation/Test/TB.Elec.200GeV.SIM.root" + ] + + if doRDO: + ConfigFlags.Output.RDOFileName = f"{fileroot}.RDO.root" + else: + ConfigFlags.addFlag("Output.xAODFileName", f"{fileroot}.xAOD.root") + ConfigFlags.Output.ESDFileName = f"{fileroot}.ESD.root" + + ConfigFlags.lock() + + from CalypsoConfiguration.MainServicesConfig import MainServicesCfg + acc = MainServicesCfg(ConfigFlags) + + from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg + from AthenaPoolCnvSvc.PoolWriteConfig import PoolWriteCfg + + acc.merge(PoolReadCfg(ConfigFlags)) + acc.merge(PoolWriteCfg(ConfigFlags)) + + if doRDO: + from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg + itemList = [ + "RawWaveformContainer#*" + ] + acc.merge(OutputStreamCfg(ConfigFlags, "RDO", itemList,disableEventTag=True)) + + + else: + from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg + itemList = [ + "xAOD::EventInfo#*", + "xAOD::WaveformHitContainer#*", + "xAOD::WaveformHitAuxContainer#*", + ] + + acc.merge(OutputStreamCfg(ConfigFlags, "xAOD", itemList, disableEventTag=True)) + + from ScintDigiAlgs.ScintDigiAlgsConfig import ScintWaveformDigitizationCfg + acc.merge(ScintWaveformDigitizationCfg(ConfigFlags)) + + from CaloDigiAlgs.CaloDigiAlgsConfig import CaloWaveformDigitizationCfg + acc.merge(CaloWaveformDigitizationCfg(ConfigFlags)) + + if not doRDO: + from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionCfg + acc.merge(WaveformReconstructionCfg(ConfigFlags)) + + #acc.foreach_component("*").OutputLevel = VERBOSE + + # Execute and finish + sc = acc.run(maxEvents=100) + + # Success should be 0 + sys.exit(not sc.isSuccess()) diff --git a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx new file mode 100644 index 0000000000000000000000000000000000000000..e4776da3bf1b140fa1c614642d939f46557b0ba6 --- /dev/null +++ b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx @@ -0,0 +1,26 @@ +/* + Copyright (C) 2021 CERN for the benefit of the FASER collaboration +*/ + +/** + * @file WaveformDigitisationTool.cxx + * Implementation file for the WaveformDigitisationTool class + * @ author C. Gwilliam, 2021 + **/ + +#include "WaveformDigitisationTool.h" + +// Constructor +WaveformDigitisationTool::WaveformDigitisationTool(const std::string& type, const std::string& name, const IInterface* parent) : + base_class(type, name, parent) +{ +} + +// Initialization +StatusCode +WaveformDigitisationTool::initialize() { + ATH_MSG_INFO( name() << "::initalize()" ); + return StatusCode::SUCCESS; +} + + diff --git a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h new file mode 100644 index 0000000000000000000000000000000000000000..8a5ba71f3dd124fcdd2c6b4b8124ee96591512da --- /dev/null +++ b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h @@ -0,0 +1,36 @@ +/* + Copyright (C) 2021 CERN for the benefit of the FASER collaboration +*/ + +/** @file WaveformDigitisationTool.h + * Header file for WaveformDigitisationTool.h + * + */ +#ifndef WAVEDIGITOOLS_WAVEFORMDIGITISATIONTOOL_H +#define WAVEDIGITOOLS_WAVEFORMDIGITISATIONTOOL_H + +//Athena +#include "AthenaBaseComps/AthAlgTool.h" +#include "WaveDigiTools/IWaveformDigitisationTool.h" + +//Gaudi +#include "GaudiKernel/ToolHandle.h" + +//STL + +class WaveformDigitisationTool: public extends<AthAlgTool, IWaveformDigitisationTool> { + public: + + /// Normal constructor for an AlgTool; 'properties' are also declared here + WaveformDigitisationTool(const std::string& type, + const std::string& name, const IInterface* parent); + + /// Retrieve the necessary services in initialize + StatusCode initialize(); + + private: + // None + +}; + +#endif // WAVEDIGITOOLS_WAVEFORMDIGITISATIONTOOL_H diff --git a/Waveform/WaveDigiTools/src/components/WaveDigiTools_entries.cxx b/Waveform/WaveDigiTools/src/components/WaveDigiTools_entries.cxx new file mode 100644 index 0000000000000000000000000000000000000000..169476d905bb3e3a56de0d15aac7b7d13c3bdeb1 --- /dev/null +++ b/Waveform/WaveDigiTools/src/components/WaveDigiTools_entries.cxx @@ -0,0 +1,3 @@ +#include "../WaveformDigitisationTool.h" + +DECLARE_COMPONENT( WaveformDigitisationTool ) diff --git a/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py b/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py index 39a99ed0c39d01e08c5f555e4e3efd98e0bc70a1..64764cd1c30ca7f039a28799dd16cc22bb50720c 100644 --- a/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py +++ b/Waveform/WaveRecAlgs/python/WaveRecAlgsConfig.py @@ -11,24 +11,27 @@ WaveformReconstructionTool = CompFactory.WaveformReconstructionTool ClockReconstructionTool = CompFactory.ClockReconstructionTool # One stop shopping for normal FASER data -def WaveformReconstructionCfg(flags, naive = True): +def WaveformReconstructionCfg(flags, naive = False): """ Return all algorithms and tools for Waveform reconstruction """ acc = ComponentAccumulator() + if not flags.Input.isMC: + acc.merge(WaveformClockRecCfg(flags, "ClockRecAlg")) + if flags.Input.isMC and naive: - if not "TB" in flags.GeoModel.FaserVersion: + if "TB" not in flags.GeoModel.FaserVersion: acc.merge(PseudoScintHitToWaveformRecCfg(flags, "PseudoTimingHitWaveformRecAlg", "Trigger")) acc.merge(PseudoScintHitToWaveformRecCfg(flags, "PseudoVetoHitToWaveformRecAlg", "Veto")) acc.merge(PseudoScintHitToWaveformRecCfg(flags, "PseudoPresehowerHitWaveformRecAlg", "Preshower")) acc.merge(PseudoCaloHitToWaveformRecCfg(flags, "PseudoCaloHitWaveformRecAlg")) return acc - acc.merge(WaveformClockRecCfg(flags, "ClockRecAlg")) - + if "TB" not in flags.GeoModel.FaserVersion: + acc.merge(WaveformHitRecCfg(flags, "TimingWaveformRecAlg", "Trigger")) acc.merge(WaveformHitRecCfg(flags, "VetoWaveformRecAlg", "Veto")) - acc.merge(WaveformHitRecCfg(flags, "TimingWaveformRecAlg", "Trigger")) acc.merge(WaveformHitRecCfg(flags, "PreshowerWaveformRecAlg", "Preshower")) acc.merge(WaveformHitRecCfg(flags, "CaloWaveformRecAlg", "Calo")) + return acc # Return configured WaveformClock reconstruction algorithm @@ -52,7 +55,13 @@ def WaveformHitRecCfg(flags, name="WaveformRecAlg", source="", **kwargs): acc = ComponentAccumulator() + if flags.Input.isMC: + kwargs.setdefault("PeakThreshold", 5) + tool = WaveformReconstructionTool(name=source+"WaveformRecTool", **kwargs) + + if "PeakThreshold" in kwargs: kwargs.pop("PeakThreshold") + kwargs.setdefault("WaveformContainerKey", source+"Waveforms") kwargs.setdefault("WaveformHitContainerKey", source+"WaveformHits") kwargs.setdefault("WaveformReconstructionTool", tool)