diff --git a/Calorimeter/CaloDigiAlgs/CMakeLists.txt b/Calorimeter/CaloDigiAlgs/CMakeLists.txt
index 42991aa03cc03af6e51ab161d36d1d61df1a7956..5a4de58a7469a2acb320de9c758ea1aef98ebcac 100644
--- a/Calorimeter/CaloDigiAlgs/CMakeLists.txt
+++ b/Calorimeter/CaloDigiAlgs/CMakeLists.txt
@@ -9,7 +9,7 @@ atlas_subdir( CaloDigiAlgs )
atlas_add_component( CaloDigiAlgs
src/*.cxx src/*.h
src/components/*.cxx
- LINK_LIBRARIES AthenaBaseComps Identifier StoreGateLib WaveRawEvent FaserCaloSimEvent WaveDigiToolsLib)
+ LINK_LIBRARIES AthenaBaseComps Identifier FaserCaloIdentifier StoreGateLib WaveRawEvent FaserCaloSimEvent WaveDigiToolsLib)
atlas_install_python_modules( python/*.py )
diff --git a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx
index 1abfcf9f11ecf71c135723f51d799828205ea2be..c54e1e37344d963f28a12f9c9f8b5d86781f93c4 100644
--- a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx
+++ b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx
@@ -2,13 +2,15 @@
#include "Identifier/Identifier.h"
-#include <vector>
+#include "FaserCaloSimEvent/CaloHitIdHelper.h"
+
#include <map>
#include <utility>
CaloWaveformDigiAlg::CaloWaveformDigiAlg(const std::string& name,
ISvcLocator* pSvcLocator)
- : AthReentrantAlgorithm(name, pSvcLocator) {
+ : AthReentrantAlgorithm(name, pSvcLocator)
+{
}
@@ -19,26 +21,25 @@ CaloWaveformDigiAlg::initialize() {
// Initalize tools
ATH_CHECK( m_digiTool.retrieve() );
-
// Set key to read waveform from
ATH_CHECK( m_caloHitContainerKey.initialize() );
// Set key to write container
ATH_CHECK( m_waveformContainerKey.initialize() );
- // Will eventually depend on the type of detector
- // TODO: Vary time at which centre it?
- // TODO: Change params compared to scint
- // m_kernel = new TF1("PDF", " ROOT::Math::crystalball_pdf(x, -0.9, 10, 4, 900)", 0, 1200);
-
+ // Set up helper
+ ATH_CHECK(detStore()->retrieve(m_ecalID, "EcalID"));
+
+ // Create CB time kernel and pre-evaluate for number of samples
m_kernel = new TF1("PDF", "[4] * ROOT::Math::crystalball_pdf(x, [0],[1],[2],[3])", 0, 1200);
- //m_kernel->SetParameters(-0.25,10,4,900);
m_kernel->SetParameter(0, m_CB_alpha);
m_kernel->SetParameter(1, m_CB_n);
m_kernel->SetParameter(2, m_CB_sigma);
m_kernel->SetParameter(3, m_CB_mean);
m_kernel->SetParameter(4, m_CB_norm);
+ // Pre-evaluate time kernel for each bin
+ m_timekernel = m_digiTool->evaluate_timekernel(m_kernel);
return StatusCode::SUCCESS;
}
@@ -55,12 +56,10 @@ CaloWaveformDigiAlg::finalize() {
StatusCode
CaloWaveformDigiAlg::execute(const EventContext& ctx) const {
ATH_MSG_DEBUG("Executing");
-
- ATH_MSG_DEBUG("Run: " << ctx.eventID().run_number()
- << " Event: " << ctx.eventID().event_number());
+ ATH_MSG_DEBUG("Run: " << ctx.eventID().run_number() << " Event: " << ctx.eventID().event_number());
// Find the input HIT collection
- SG::ReadHandle<CaloHitCollection> caloHitHandle(m_caloHitContainerKey, ctx);
+ SG::ReadHandle<CaloHitCollection> caloHitHandle(m_caloHitContainerKey, ctx);
ATH_CHECK( caloHitHandle.isValid() );
ATH_MSG_DEBUG("Found ReadHandle for CaloHitCollection " << m_caloHitContainerKey);
@@ -75,11 +74,50 @@ CaloWaveformDigiAlg::execute(const EventContext& ctx) const {
ATH_MSG_DEBUG("CaloHitCollection found with zero length!");
return StatusCode::SUCCESS;
}
-
- // Digitise the hits
- CHECK( m_digiTool->digitise<CaloHitCollection>(caloHitHandle.ptr(),
- waveformContainerHandle.ptr(), m_kernel,
- std::pair<float, float>(m_base_mean, m_base_rms)) );
+
+ // Create structure to store pulse for each channel
+ std::map<Identifier, std::vector<uint16_t>> waveforms = m_digiTool->create_waveform_map(m_ecalID);
+
+ for (const auto& tk : m_timekernel) {
+ std::map<unsigned int, float> counts;
+
+ // Convolve hit energy with evaluated kernel and sum for each hit id (i.e. channel)
+ for (const auto& hit : *caloHitHandle) {
+ counts[hit.identify()] += tk * hit.energyLoss();
+ }
+
+ // Subtract count from basleine and add result to correct waveform vector
+ for (const auto& c : counts) {
+
+ unsigned int baseline = m_digiTool->generate_baseline(m_base_mean, m_base_rms);
+ int value = baseline - c.second;
+
+ if (value < 0) {
+ ATH_MSG_WARNING("Found pulse " << c.second << " larger than baseline " << c.first);
+ value = 0; // Protect against scaling signal above baseline
+ }
+
+ // Convert hit id to Identifier and store
+ Identifier id = CaloHitIdHelper::GetHelper()->getIdentifier(c.first);
+ waveforms[id].push_back(value);
+ }
+ }
+
+
+ //m_chrono->chronoStop("Digit");
+ //m_chrono->chronoStart("Write");
+
+ // Loop over wavefrom vectors to make and store waveform
+ unsigned int nsamples = m_digiTool->nsamples();
+ for (const auto& w : waveforms) {
+ RawWaveform* wfm = new RawWaveform();
+ wfm->setWaveform(0, w.second);
+ wfm->setIdentifier(w.first);
+ wfm->setSamples(nsamples);
+ waveformContainerHandle->push_back(wfm);
+ }
+
+ //m_chrono->chronoStop("Write");
ATH_MSG_DEBUG("WaveformsHitContainer " << waveformContainerHandle.name() << "' filled with "<< waveformContainerHandle->size() <<" items");
diff --git a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h
index f99d13b5565712045bcef7507fd04bca6bf0d07a..8eeb0a44eba0dde0ede044036b57406cfcf25efc 100644
--- a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h
+++ b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h
@@ -19,11 +19,15 @@
#include "GaudiKernel/ServiceHandle.h"
#include "GaudiKernel/ToolHandle.h"
+// Helpers
+#include "FaserCaloIdentifier/EcalID.h"
+
// ROOT
#include "TF1.h"
// STL
#include <string>
+#include <vector>
class CaloWaveformDigiAlg : public AthReentrantAlgorithm {
@@ -48,6 +52,8 @@ class CaloWaveformDigiAlg : public AthReentrantAlgorithm {
CaloWaveformDigiAlg &operator=(const CaloWaveformDigiAlg&) = delete;
//@}
+ /** @name Steerable pameters for crystal ball and baseline **/
+ //@{
Gaudi::Property<double> m_CB_alpha {this, "CB_alpha", 0, "Alpha of the crystal ball function"};
Gaudi::Property<double> m_CB_n {this, "CB_n", 0, "n of the crystal ball function"};
Gaudi::Property<double> m_CB_mean {this, "CB_mean", 0, "Mean of the crystal ball function"};
@@ -56,10 +62,16 @@ class CaloWaveformDigiAlg : public AthReentrantAlgorithm {
Gaudi::Property<double> m_base_mean {this, "base_mean", 0, "Mean of the baseline"};
Gaudi::Property<double> m_base_rms {this, "base_rms", 0, "RMS of the baseline"};
+ //@}
- /// Kernel PDF
- TF1* m_kernel;
+ /** Kernel PDF and evaluated values **/
+ //@{
+ TF1* m_kernel;
+ std::vector<float> m_timekernel;
+ //@}
+ /// Detector ID helper
+ const EcalID* m_ecalID{nullptr};
/**
* @name Digitisation tool
@@ -89,4 +101,6 @@ class CaloWaveformDigiAlg : public AthReentrantAlgorithm {
};
+
+
#endif // CALODIGIALGS_CALODIGIALG_H
diff --git a/Calorimeter/FaserCaloSimEvent/CMakeLists.txt b/Calorimeter/FaserCaloSimEvent/CMakeLists.txt
index 9777d4178e6e465f0f57de0c82173b0620fbe46a..2b034a3189207f186984b163c0c5a034bee7f956 100644
--- a/Calorimeter/FaserCaloSimEvent/CMakeLists.txt
+++ b/Calorimeter/FaserCaloSimEvent/CMakeLists.txt
@@ -18,11 +18,11 @@ atlas_add_library( FaserCaloSimEvent
PRIVATE_INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${GEANT4_INCLUDE_DIRS}
DEFINITIONS ${CLHEP_DEFINITIONS}
LINK_LIBRARIES ${CLHEP_LIBRARIES} AthAllocators AthenaKernel CxxUtils GeneratorObjects HitManagement StoreGateLib SGtests
- PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES} FaserCaloIdentifier )
+ PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES} FaserCaloIdentifier Identifier)
atlas_add_dictionary( FaserCaloSimEventDict
FaserCaloSimEvent/CaloSimEventDict.h
FaserCaloSimEvent/selection.xml
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS} ${GEANT4_INCLUDE_DIRS}
- LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} AthAllocators CxxUtils GeneratorObjects HitManagement StoreGateLib SGtests FaserCaloIdentifier FaserCaloSimEvent )
+ LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} AthAllocators CxxUtils GeneratorObjects HitManagement StoreGateLib SGtests FaserCaloIdentifier FaserCaloSimEvent Identifier)
diff --git a/Calorimeter/FaserCaloSimEvent/FaserCaloSimEvent/CaloHitIdHelper.h b/Calorimeter/FaserCaloSimEvent/FaserCaloSimEvent/CaloHitIdHelper.h
index 5500649bc0c789aa5119503ee357839ee7e11b9e..e94581470c0fa76dddbec99c44aa93ba61484cdb 100644
--- a/Calorimeter/FaserCaloSimEvent/FaserCaloSimEvent/CaloHitIdHelper.h
+++ b/Calorimeter/FaserCaloSimEvent/FaserCaloSimEvent/CaloHitIdHelper.h
@@ -23,6 +23,11 @@
// This class is singleton and static method and variable are used.
#include "CxxUtils/checker_macros.h"
+
+#include "Identifier/Identifier.h"
+
+#include "FaserCaloIdentifier/EcalID.h"
+
ATLAS_NO_CHECK_FILE_THREAD_SAFETY;
class CaloHitIdHelper : HitIdHelper {
@@ -43,6 +48,8 @@ class CaloHitIdHelper : HitIdHelper {
// Left or Right
int getModule(const int& hid) const;
+ Identifier getIdentifier(const int& hid) const;
+
//
// Info packing:
int buildHitId(const int, const int) const;
@@ -54,6 +61,9 @@ class CaloHitIdHelper : HitIdHelper {
//
// Initialize the helper, only called by the constructor
void Initialize();
+
+ /// Detector ID helper
+ const EcalID* m_ecalID{nullptr};
};
#endif // CALOSIMEVENT_CALOHITIDHELPER
diff --git a/Calorimeter/FaserCaloSimEvent/src/CaloHitIdHelper.cxx b/Calorimeter/FaserCaloSimEvent/src/CaloHitIdHelper.cxx
index eec88553dae53bc6e415b385494b32a7f83645ad..0e98fbd6aa3c1be718b049ab91e74ae92b682ed4 100644
--- a/Calorimeter/FaserCaloSimEvent/src/CaloHitIdHelper.cxx
+++ b/Calorimeter/FaserCaloSimEvent/src/CaloHitIdHelper.cxx
@@ -6,8 +6,6 @@
#include "FaserCaloSimEvent/CaloHitIdHelper.h"
#include "StoreGate/StoreGateSvc.h"
-#include "StoreGate/StoreGateSvc.h"
-#include "FaserCaloIdentifier/EcalID.h"
#include "G4Types.hh"
#ifdef G4MULTITHREADED
@@ -42,10 +40,9 @@ void CaloHitIdHelper::Initialize() {
// determine whether hits were created with an SLHC dictionary
// in which case eta module field is expanded.
// Need to lock this thread-unsafe retrieval
- const EcalID* pix;
ServiceHandle<StoreGateSvc> detStore ("DetectorStore", "CaloHitIdHelper");
if (detStore.retrieve().isSuccess()) {
- if (detStore->retrieve(pix, "EcalID").isFailure()) { pix = 0; }
+ if (detStore->retrieve(m_ecalID, "EcalID").isFailure()) { m_ecalID = 0; }
}
InitializeField("Row", 0, 2);
@@ -64,6 +61,14 @@ int CaloHitIdHelper::getModule(const int& hid) const
return this->GetFieldValue("Module", hid);
}
+// identifier
+Identifier CaloHitIdHelper::getIdentifier(const int& hid) const
+{
+ return m_ecalID->pmt_id(getRow(hid), getModule(hid), 0);
+
+}
+
+
//
// Info packing:
int CaloHitIdHelper::buildHitId( const int row,
diff --git a/Scintillator/ScintDigiAlgs/CMakeLists.txt b/Scintillator/ScintDigiAlgs/CMakeLists.txt
index 43695c5f90d96053a8b44fa3ac41c0f4c159ad99..3ebe3d691decc51518eb6f64ecae03efd76aaf73 100644
--- a/Scintillator/ScintDigiAlgs/CMakeLists.txt
+++ b/Scintillator/ScintDigiAlgs/CMakeLists.txt
@@ -9,7 +9,7 @@ atlas_subdir( ScintDigiAlgs )
atlas_add_component( ScintDigiAlgs
src/*.cxx src/*.h
src/components/*.cxx
- LINK_LIBRARIES AthenaBaseComps Identifier StoreGateLib WaveRawEvent ScintSimEvent WaveDigiToolsLib)
+ LINK_LIBRARIES AthenaBaseComps Identifier ScintIdentifier StoreGateLib WaveRawEvent ScintSimEvent WaveDigiToolsLib)
atlas_install_python_modules( python/*.py )
diff --git a/Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py b/Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py
index 7f2ba8530575eeb6a6c9b51779c482f84bdfc648..292a19c0890672bb3a1530b080d936c8716a2fff 100644
--- a/Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py
+++ b/Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py
@@ -15,6 +15,13 @@ dict_CB_param["Timing"]=dict(CB_alpha=-0.32, CB_n=65, CB_mean=846, CB_sigma=5.3,
dict_CB_param["Veto"]=dict(CB_alpha=-0.38, CB_n=25, CB_mean=815, CB_sigma=7.7, CB_norm = 1000) # copy from Trigger; Veto was not in TestBeam, but in sim "Veto" is the TestBeam Trigger component
dict_CB_param["Preshower"]=dict(CB_alpha=-0.32, CB_n=65, CB_mean=846, CB_sigma=5.3, CB_norm = 500)
+dict_baseline_params = {
+ "Trigger" : {"mean" : 15000, "rms" : 3},
+ "Timing" : {"mean" : 15000, "rms" : 3},
+ "Veto" : {"mean" : 15000, "rms" : 3},
+ "Preshower" : {"mean" : 15000, "rms" : 3},
+ }
+
# One stop shopping for normal FASER data
def ScintWaveformDigitizationCfg(flags):
""" Return all algorithms and tools for Waveform digitization """
@@ -48,8 +55,8 @@ def ScintWaveformDigiCfg(flags, name="ScintWaveformDigiAlg", source="", **kwargs
digiAlg.CB_sigma = dict_CB_param[source]["CB_sigma"]
digiAlg.CB_norm = dict_CB_param[source]["CB_norm"]
- digiAlg.base_mean = 15000
- digiAlg.base_rms = 3
+ digiAlg.base_mean = dict_baseline_params[source]["mean"]
+ digiAlg.base_rms = dict_baseline_params[source]["rms"]
kwargs.setdefault("WaveformDigitisationTool", tool)
diff --git a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
index 91cffd93552d222b90c200cf50fa7842e2b7177a..d5235a92ccdacbc2631e4dafc66a3f9dea48dc9d 100644
--- a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
+++ b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
@@ -1,11 +1,8 @@
-
#include "ScintWaveformDigiAlg.h"
-#include "Identifier/Identifier.h"
+#include "ScintSimEvent/ScintHitIdHelper.h"
-#include <vector>
#include <map>
-#include <utility>
ScintWaveformDigiAlg::ScintWaveformDigiAlg(const std::string& name,
@@ -21,25 +18,28 @@ ScintWaveformDigiAlg::initialize() {
// Initalize tools
ATH_CHECK( m_digiTool.retrieve() );
-
// Set key to read waveform from
ATH_CHECK( m_scintHitContainerKey.initialize() );
// Set key to write container
ATH_CHECK( m_waveformContainerKey.initialize() );
- // Will eventually depend on the type of detector
- // TODO: Vary time at which centre it?
- // TODO: Better parameters
-
+ // Set up helpers
+ ATH_CHECK(detStore()->retrieve(m_vetoID, "VetoID"));
+ ATH_CHECK(detStore()->retrieve(m_triggerID, "TriggerID"));
+ ATH_CHECK(detStore()->retrieve(m_preshowerID, "PreshowerID"));
-
+ // Create CB time kernel and pre-evaluate for number of samples
m_kernel = new TF1("PDF", "[4] * ROOT::Math::crystalball_pdf(x, [0],[1],[2],[3])", 0, 1200);
m_kernel->SetParameter(0, m_CB_alpha);
m_kernel->SetParameter(1, m_CB_n);
m_kernel->SetParameter(2, m_CB_sigma);
m_kernel->SetParameter(3, m_CB_mean);
m_kernel->SetParameter(4, m_CB_norm);
+
+ // Pre-evaluate time kernel for each bin
+ m_timekernel = m_digiTool->evaluate_timekernel(m_kernel);
+
return StatusCode::SUCCESS;
}
@@ -55,11 +55,9 @@ ScintWaveformDigiAlg::finalize() {
StatusCode
ScintWaveformDigiAlg::execute(const EventContext& ctx) const {
ATH_MSG_DEBUG("Executing");
+ ATH_MSG_DEBUG("Run: " << ctx.eventID().run_number() << " Event: " << ctx.eventID().event_number());
- ATH_MSG_DEBUG("Run: " << ctx.eventID().run_number()
- << " Event: " << ctx.eventID().event_number());
-
- // Find the input HIT collection
+ // Find the input HITS collection
SG::ReadHandle<ScintHitCollection> scintHitHandle(m_scintHitContainerKey, ctx);
ATH_CHECK( scintHitHandle.isValid() );
@@ -70,16 +68,59 @@ ScintWaveformDigiAlg::execute(const EventContext& ctx) const {
ATH_CHECK( waveformContainerHandle.record( std::make_unique<RawWaveformContainer>()) );
ATH_MSG_DEBUG("WaveformsContainer '" << waveformContainerHandle.name() << "' initialized");
-
+
if (scintHitHandle->size() == 0) {
ATH_MSG_DEBUG("ScintHitCollection found with zero length!");
return StatusCode::SUCCESS;
}
+
+ // Create structure to store pulse for each channel
+ std::map<Identifier, std::vector<uint16_t>> waveforms;
+
+ auto first = *scintHitHandle->begin();
+ if (first.isVeto()) {
+ waveforms = m_digiTool->create_waveform_map(m_vetoID);
+ } else if (first.isTrigger()) {
+ waveforms = m_digiTool->create_waveform_map(m_triggerID);
+ } else if (first.isPreshower()) {
+ waveforms = m_digiTool->create_waveform_map(m_preshowerID);
+ }
+
+ // Loop over time samples
+ for (const auto& tk : m_timekernel) {
+ std::map<unsigned int, float> counts;
+
+ // Convolve hit energy with evaluated kernel and sum for each hit id (i.e. channel)
+ for (const auto& hit : *scintHitHandle) {
+ counts[hit.identify()] += tk * hit.energyLoss();
+ }
+
+ // Subtract count from basleine and add result to correct waveform vector
+ for (const auto& c : counts) {
+
+ unsigned int baseline = m_digiTool->generate_baseline(m_base_mean, m_base_rms);
+ int value = baseline - c.second;
+
+ if (value < 0) {
+ ATH_MSG_WARNING("Found pulse " << c.second << " larger than baseline " << c.first);
+ value = 0; // Protect against scaling signal above baseline
+ }
+
+ // Convert hit id to Identifier and store
+ Identifier id = ScintHitIdHelper::GetHelper()->getIdentifier(c.first);
+ waveforms[id].push_back(value);
+ }
+ }
- // Digitise the hits
- CHECK( m_digiTool->digitise<ScintHitCollection>(scintHitHandle.ptr(),
- waveformContainerHandle.ptr(), m_kernel,
- std::pair<float, float>(m_base_mean, m_base_rms)) );
+ // Loop over wavefrom vectors to make and store raw waveform
+ unsigned int nsamples = m_digiTool->nsamples();
+ for (const auto& w : waveforms) {
+ RawWaveform* wfm = new RawWaveform();
+ wfm->setWaveform(0, w.second);
+ wfm->setIdentifier(w.first);
+ wfm->setSamples(nsamples);
+ waveformContainerHandle->push_back(wfm);
+ }
ATH_MSG_DEBUG("WaveformsHitContainer " << waveformContainerHandle.name() << "' filled with "<< waveformContainerHandle->size() <<" items");
diff --git a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h
index 9bf7843a1867eeffa4587a39875913ea0685d773..de03370698620986377c695852af07c7d8c79b18 100644
--- a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h
+++ b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h
@@ -19,12 +19,21 @@
#include "GaudiKernel/ServiceHandle.h"
#include "GaudiKernel/ToolHandle.h"
+//Helpers
+#include "ScintIdentifier/VetoID.h"
+#include "ScintIdentifier/TriggerID.h"
+#include "ScintIdentifier/PreshowerID.h"
+#include "ScintSimEvent/ScintHitIdHelper.h"
+#include "Identifier/Identifier.h"
+
+
// ROOT
#include "TF1.h"
// STL
#include <string>
+#include <vector>
class ScintWaveformDigiAlg : public AthReentrantAlgorithm {
@@ -50,6 +59,10 @@ class ScintWaveformDigiAlg : public AthReentrantAlgorithm {
ScintWaveformDigiAlg &operator=(const ScintWaveformDigiAlg&) = delete;
//@}
+ ///
+
+ /** @name Steerable pameters for crystal ball and baseline **/
+ //@{
Gaudi::Property<double> m_CB_alpha {this, "CB_alpha", 0, "Alpha of the crystal ball function"};
Gaudi::Property<double> m_CB_n {this, "CB_n", 0, "n of the crystal ball function"};
Gaudi::Property<double> m_CB_mean {this, "CB_mean", 0, "Mean of the crystal ball function"};
@@ -58,10 +71,20 @@ class ScintWaveformDigiAlg : public AthReentrantAlgorithm {
Gaudi::Property<double> m_base_mean {this, "base_mean", 0, "Mean of the baseline"};
Gaudi::Property<double> m_base_rms {this, "base_rms", 0, "RMS of the baseline"};
+ //@}
- /// Kernel PDF
- TF1* m_kernel;
+ /** Kernel PDF and evaluated values **/
+ //@{
+ TF1* m_kernel;
+ std::vector<float> m_timekernel;
+ //@}
+
+
+ /// Detector ID helpers
+ const VetoID* m_vetoID{nullptr};
+ const TriggerID* m_triggerID{nullptr};
+ const PreshowerID* m_preshowerID{nullptr};
/**
* @name Digitisation tool
@@ -91,4 +114,5 @@ class ScintWaveformDigiAlg : public AthReentrantAlgorithm {
};
+
#endif // SCINTDIGIALGS_SCINTDIGIALG_H
diff --git a/Scintillator/ScintSimEvent/CMakeLists.txt b/Scintillator/ScintSimEvent/CMakeLists.txt
index 27b63871df2aa95ea230c6be47454573a7c14d16..07224a272694548a6c39740e081059d7652a9d91 100644
--- a/Scintillator/ScintSimEvent/CMakeLists.txt
+++ b/Scintillator/ScintSimEvent/CMakeLists.txt
@@ -18,11 +18,11 @@ atlas_add_library( ScintSimEvent
PRIVATE_INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${GEANT4_INCLUDE_DIRS}
DEFINITIONS ${CLHEP_DEFINITIONS}
LINK_LIBRARIES ${CLHEP_LIBRARIES} AthAllocators AthenaKernel CxxUtils GeneratorObjects HitManagement StoreGateLib SGtests
- PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES} ScintIdentifier )
+ PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES} ScintIdentifier Identifier)
atlas_add_dictionary( ScintSimEventDict
ScintSimEvent/ScintSimEventDict.h
ScintSimEvent/selection.xml
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS} ${GEANT4_INCLUDE_DIRS}
- LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} AthAllocators CxxUtils GeneratorObjects HitManagement StoreGateLib SGtests ScintIdentifier ScintSimEvent )
+ LINK_LIBRARIES ${ROOT_LIBRARIES} ${CLHEP_LIBRARIES} AthAllocators CxxUtils GeneratorObjects HitManagement StoreGateLib SGtests ScintIdentifier ScintSimEvent Identifier)
diff --git a/Scintillator/ScintSimEvent/ScintSimEvent/ScintHitIdHelper.h b/Scintillator/ScintSimEvent/ScintSimEvent/ScintHitIdHelper.h
index 2874bac894b15684252de5a531081f338750c10c..233acb37e9a01ff47d8e3f93924e473eebe37f62 100644
--- a/Scintillator/ScintSimEvent/ScintSimEvent/ScintHitIdHelper.h
+++ b/Scintillator/ScintSimEvent/ScintSimEvent/ScintHitIdHelper.h
@@ -23,6 +23,15 @@
// This class is singleton and static method and variable are used.
#include "CxxUtils/checker_macros.h"
+
+//Helpers
+#include "ScintIdentifier/VetoID.h"
+#include "ScintIdentifier/TriggerID.h"
+#include "ScintIdentifier/PreshowerID.h"
+
+#include "Identifier/Identifier.h"
+
+
ATLAS_NO_CHECK_FILE_THREAD_SAFETY;
class ScintHitIdHelper : HitIdHelper {
@@ -43,6 +52,9 @@ class ScintHitIdHelper : HitIdHelper {
// Layer/Disk
int getPlate(const int& hid) const;
+ // identifier
+ Identifier getIdentifier(const int& hid) const;
+
//
// Info packing:
int buildHitId(const int, const int, const int) const;
@@ -54,6 +66,12 @@ class ScintHitIdHelper : HitIdHelper {
//
// Initialize the helper, only called by the constructor
void Initialize();
+
+ /// Detector ID helpers
+ const VetoID* m_vetoID{nullptr};
+ const TriggerID* m_triggerID{nullptr};
+ const PreshowerID* m_preshowerID{nullptr};
+
};
#endif // SCINTSIMEVENT_SCINTHITIDHELPER
diff --git a/Scintillator/ScintSimEvent/src/ScintHitIdHelper.cxx b/Scintillator/ScintSimEvent/src/ScintHitIdHelper.cxx
index cacdb5d720bfdc7759188481346899d48069d856..5f115d47109690e545267e6cbc561fc5b8e7afdf 100644
--- a/Scintillator/ScintSimEvent/src/ScintHitIdHelper.cxx
+++ b/Scintillator/ScintSimEvent/src/ScintHitIdHelper.cxx
@@ -6,7 +6,6 @@
#include "ScintSimEvent/ScintHitIdHelper.h"
#include "StoreGate/StoreGateSvc.h"
-#include "ScintIdentifier/VetoID.h"
#include "G4Types.hh"
#ifdef G4MULTITHREADED
@@ -44,7 +43,9 @@ void ScintHitIdHelper::Initialize() {
const VetoID* pix;
ServiceHandle<StoreGateSvc> detStore ("DetectorStore", "ScitHitIdHelper");
if (detStore.retrieve().isSuccess()) {
- if (detStore->retrieve(pix, "VetoID").isFailure()) { pix = 0; }
+ if (detStore->retrieve(m_vetoID, "VetoID").isFailure()) { m_vetoID = 0; }
+ if (detStore->retrieve(m_triggerID, "TriggerID").isFailure()) { m_triggerID = 0; }
+ if (detStore->retrieve(m_preshowerID, "PreshowerID").isFailure()) { m_preshowerID = 0; }
}
InitializeField("VetoTriggerPreshower", 0, 2);
@@ -87,6 +88,20 @@ int ScintHitIdHelper::getPlate(const int& hid) const
return this->GetFieldValue("Plate", hid);
}
+
+// identifier
+Identifier ScintHitIdHelper::getIdentifier(const int& hid) const
+{
+ if (isVeto(hid)) {
+ return m_vetoID->pmt_id(getStation(hid), getPlate(hid), 0);
+ } else if (isTrigger(hid)) {
+ return m_triggerID->pmt_id(getStation(hid), getPlate(hid), 0);
+ } else if (isPreshower(hid)) {
+ return m_preshowerID->pmt_id(getStation(hid), getPlate(hid), 0);
+ }
+ return Identifier();
+}
+
//
// Info packing:
int ScintHitIdHelper::buildHitId(const int veto_trigger_preshower,
diff --git a/Waveform/WaveDigiTools/CMakeLists.txt b/Waveform/WaveDigiTools/CMakeLists.txt
index 692fdb69bc14451ba5a6a660d011bad5b14b66e5..d7e9fd857b273c764dad780aceebd0d3f58acc3a 100644
--- a/Waveform/WaveDigiTools/CMakeLists.txt
+++ b/Waveform/WaveDigiTools/CMakeLists.txt
@@ -13,13 +13,12 @@ atlas_add_library( WaveDigiToolsLib
WaveDigiTools/*.h src/*.cxx src/*.h
PUBLIC_HEADERS WaveDigiTools
PRIVATE_INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
- LINK_LIBRARIES AthenaBaseComps AthenaKernel GeoPrimitives WaveRawEvent
+ LINK_LIBRARIES AthenaBaseComps AthenaKernel GeoPrimitives WaveRawEvent Identifier
PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES}
)
atlas_add_component( WaveDigiTools
src/components/*.cxx
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
- LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps GaudiKernel WaveDigiToolsLib )
-
+ LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps GaudiKernel WaveDigiToolsLib)
diff --git a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h
index 7241e8167d36facc42748141b58f5b640b4d33ee..c30351902a1b79e9a58ad33d5f58a5d04ed13faa 100644
--- a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h
+++ b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h
@@ -23,10 +23,15 @@
#include "WaveRawEvent/RawWaveformContainer.h"
#include "WaveRawEvent/RawWaveform.h"
+#include "Identifier/Identifier.h"
+
#include "TF1.h"
#include "TRandom3.h"
#include <utility>
+#include <map>
+#include <vector>
+
///Interface for waveform digitisation tools
class IWaveformDigitisationTool : virtual public IAlgTool
@@ -42,19 +47,25 @@ public:
virtual ~IWaveformDigitisationTool() = default;
- // Digitise HITS to Raw waveform
- template<class CONT>
- StatusCode digitise(const CONT* hitCollection,
- RawWaveformContainer* waveContainer,
- TF1* kernel, std::pair<float, float> base
- ) const;
+ /// Evaluate time kernel over time samples
+ virtual std::vector<float> evaluate_timekernel(TF1* kernel) const = 0;
+
+ /// Generate random baseline
+ virtual unsigned int generate_baseline(int mean, int rms) const = 0;
+
+ /// Create structure to store pulse for each channel
+ template <class T>
+ std::map<Identifier, std::vector<uint16_t>> create_waveform_map(const T* idHelper) const;
+
+ /// Number of time samples
+ unsigned int nsamples() const { return m_nsamples; }
private:
ServiceHandle<IMessageSvc> m_msgSvc;
protected:
TRandom3* m_random;
-
+ unsigned int m_nsamples;
};
#include "WaveDigiTools/IWaveformDigitisationTool.icc"
diff --git a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.icc b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.icc
index b07e1836ac5cedfefa87a8fcbe8c910de886693c..41b8c2650319a448df63d11e098fc8d0784dc056 100644
--- a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.icc
+++ b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.icc
@@ -1,62 +1,17 @@
-#include <vector>
-#include <map>
+#include "Identifier/Identifier.h"
+#include "Identifier/ExpandedIdentifier.h"
-template<class CONT>
-StatusCode IWaveformDigitisationTool::digitise(const CONT* hitCollection,
- RawWaveformContainer* container, TF1* kernel,
- std::pair<float, float> base) const {
+template <class ID>
+std::map<Identifier, std::vector<uint16_t>> IWaveformDigitisationTool::create_waveform_map(const ID* idHelper) const {
+ std::map<Identifier, std::vector<uint16_t>> waveforms;
- // Check the container
- if (!container) {
- MsgStream log(&(*m_msgSvc), name());
- log << MSG::ERROR << "HitCollection passed to digitise() is null!" << endmsg;
- return StatusCode::FAILURE;
+ for (auto itr = idHelper->pmt_begin(); itr != idHelper->pmt_end(); ++itr) {
+ const ExpandedIdentifier& extId = *itr;
+ Identifier id = idHelper->pmt_id(extId);
+ waveforms[id] = std::vector<uint16_t>();
+ waveforms[id].reserve(m_nsamples);
}
- unsigned int size = 600; // TODO: how know the correct number of time samples?
- std::vector<float> time(size);
- for (unsigned int i=0; i<size; i++) time[i] = 2.*i;
-
- std::map<unsigned int, std::vector<uint16_t>> waveforms;
- //unsigned int baseline = 800;
-
- // TODO: varying time to centre of edge of bin (odd = centre, even = edge)
-
- // Loop over time samples
- for (const auto& t : time) {
- std::map<unsigned int, float> counts;
-
- unsigned int baseline = m_random->Gaus(base.first, base.second);
-
- // Convolve hit energy with kernel and sum for each ID (i.e. channel)
- for (const auto& hit : *hitCollection) {
- counts[hit.identify()] += kernel->Eval(t) * hit.energyLoss();
- //std::cout << "HIT " << hit.identify() << " @ " << t << ": " << kernel->Eval(t) << " " << hit.energyLoss() << " -> " << counts[hit.identify()] << std::endl;
- }
-
- // Add count to correct waveform vec
- for (const auto& c : counts) {
- int value = baseline - c.second;
- if (value < 0) {
- MsgStream log(&(*m_msgSvc), name());
- log << MSG::WARNING << "Found pulse " << c.second << " larger than baseline " << c.first << endmsg;
- value = 0; // Protect against scaling signal above baseline
- }
- waveforms[c.first].push_back(value);
- //std::cout << "ADC " << c.first << " @ " << t << ": " << baseline - c.second << std::endl;
- }
- }
-
- // Loop over wavefrom vecs to make and store waveform
- for (const auto& w : waveforms) {
- RawWaveform* wfm = new RawWaveform();
- wfm->setWaveform(0, w.second);
- wfm->setIdentifier(Identifier(w.first));
- wfm->setSamples(size);
- container->push_back(wfm);
- }
-
-
- return StatusCode::SUCCESS;
+ return waveforms;
}
diff --git a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx
index 1c71fe193105f089bec81a5821d4a5558cdb39ca..c62d7f8a753490907bd1e8a513cad5341a2a46be 100644
--- a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx
+++ b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx
@@ -20,8 +20,28 @@ WaveformDigitisationTool::WaveformDigitisationTool(const std::string& type, cons
StatusCode
WaveformDigitisationTool::initialize() {
ATH_MSG_INFO( name() << "::initalize()" );
+
+ m_nsamples = 600;
m_random = new TRandom3();
+
return StatusCode::SUCCESS;
}
+std::vector<float>
+WaveformDigitisationTool::evaluate_timekernel(TF1* kernel) const {
+
+ std::vector<float> timekernel;
+ timekernel.reserve(m_nsamples);
+
+ for (unsigned int i=0; i<m_nsamples; i++) {
+ timekernel.push_back(kernel->Eval(2.*i));
+ }
+
+ return timekernel;
+}
+
+unsigned int
+WaveformDigitisationTool::generate_baseline(int mean, int rms) const {
+ return m_random->Gaus(mean, rms);
+}
diff --git a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h
index 8a5ba71f3dd124fcdd2c6b4b8124ee96591512da..e2dd5169152845824927baeeae7ce8fc36ab46f8 100644
--- a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h
+++ b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h
@@ -28,6 +28,13 @@ class WaveformDigitisationTool: public extends<AthAlgTool, IWaveformDigitisation
/// Retrieve the necessary services in initialize
StatusCode initialize();
+ /// Evaluate time kernel over samples
+ std::vector<float> evaluate_timekernel(TF1* kernel) const;
+
+ /// Generate random baseline
+ unsigned int generate_baseline(int mean, int rms) const;
+
+
private:
// None
diff --git a/Waveform/WaveRecAlgs/src/RawWaveformRecAlg.cxx b/Waveform/WaveRecAlgs/src/RawWaveformRecAlg.cxx
index 97f122e50f51c17f2ac3f0e0ff3174cf2f337347..3248939a0879178789e7c770fb584fc0b85cce07 100644
--- a/Waveform/WaveRecAlgs/src/RawWaveformRecAlg.cxx
+++ b/Waveform/WaveRecAlgs/src/RawWaveformRecAlg.cxx
@@ -52,6 +52,13 @@ RawWaveformRecAlg::execute(const EventContext& ctx) const {
ATH_CHECK( waveformHandle.isValid() );
ATH_MSG_DEBUG("Found ReadHandle for RawWaveformContainer " << m_waveformContainerKey);
+ // Find the output waveform container
+ SG::WriteHandle<xAOD::WaveformHitContainer> hitContainerHandle(m_waveformHitContainerKey, ctx);
+ ATH_CHECK( hitContainerHandle.record( std::make_unique<xAOD::WaveformHitContainer>(),
+ std::make_unique<xAOD::WaveformHitAuxContainer>() ) );
+
+ ATH_MSG_DEBUG("WaveformsHitContainer '" << hitContainerHandle.name() << "' initialized");
+
if (waveformHandle->size() == 0) {
ATH_MSG_DEBUG("Waveform container found with zero length!");
return StatusCode::SUCCESS;
@@ -69,12 +76,6 @@ RawWaveformRecAlg::execute(const EventContext& ctx) const {
ATH_MSG_WARNING("Didn't find ReadHandle for WaveformClock!");
}
- // Find the output waveform container
- SG::WriteHandle<xAOD::WaveformHitContainer> hitContainerHandle(m_waveformHitContainerKey, ctx);
- ATH_CHECK( hitContainerHandle.record( std::make_unique<xAOD::WaveformHitContainer>(),
- std::make_unique<xAOD::WaveformHitAuxContainer>() ) );
-
- ATH_MSG_DEBUG("WaveformsHitContainer '" << hitContainerHandle.name() << "' initialized");
// Reconstruct all waveforms
CHECK( m_recoTool->reconstructAll(*waveformHandle, clockptr, hitContainerHandle.ptr()) );