From 5b1e1e91d15d7488222f733d3b7e39d37e47d3c3 Mon Sep 17 00:00:00 2001
From: Eric Torrence <eric.torrence@cern.ch>
Date: Wed, 6 Dec 2023 11:48:11 +0100
Subject: [PATCH] Add tof to waveform times
---
Calorimeter/CaloDigiAlgs/CMakeLists.txt | 2 +-
.../CaloDigiAlgs/python/CaloDigiAlgsConfig.py | 3 +
.../CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx | 223 +++++++++++++++---
.../CaloDigiAlgs/src/CaloWaveformDigiAlg.h | 10 +
.../Digitization/scripts/faser_digi.py | 6 +-
.../Digitization/scripts/faser_digi_merge.py | 6 +-
Scintillator/ScintDigiAlgs/CMakeLists.txt | 2 +-
.../python/ScintDigiAlgsConfig.py | 14 +-
.../src/ScintWaveformDigiAlg.cxx | 142 ++++++++---
.../ScintDigiAlgs/src/ScintWaveformDigiAlg.h | 6 +-
.../WaveDigiTools/IWaveformDigitisationTool.h | 6 +-
.../IWaveformDigitisationTool.icc | 3 +-
.../src/WaveformDigitisationTool.cxx | 18 +-
.../src/WaveformDigitisationTool.h | 17 ++
14 files changed, 367 insertions(+), 91 deletions(-)
diff --git a/Calorimeter/CaloDigiAlgs/CMakeLists.txt b/Calorimeter/CaloDigiAlgs/CMakeLists.txt
index 89d20eae8..baa67edad 100644
--- a/Calorimeter/CaloDigiAlgs/CMakeLists.txt
+++ b/Calorimeter/CaloDigiAlgs/CMakeLists.txt
@@ -11,7 +11,7 @@ atlas_add_component( CaloDigiAlgs
src/components/*.cxx
LINK_LIBRARIES AthenaBaseComps Identifier FaserCaloIdentifier
WaveformConditionsToolsLib StoreGateLib WaveRawEvent
- FaserCaloSimEvent WaveDigiToolsLib)
+ CaloReadoutGeometry FaserCaloSimEvent WaveDigiToolsLib)
atlas_install_python_modules( python/*.py )
diff --git a/Calorimeter/CaloDigiAlgs/python/CaloDigiAlgsConfig.py b/Calorimeter/CaloDigiAlgs/python/CaloDigiAlgsConfig.py
index 85b56fe61..9f5b1355d 100644
--- a/Calorimeter/CaloDigiAlgs/python/CaloDigiAlgsConfig.py
+++ b/Calorimeter/CaloDigiAlgs/python/CaloDigiAlgsConfig.py
@@ -44,6 +44,9 @@ def CaloWaveformDigiCfg(flags, name="CaloWaveformDigiAlg", **kwargs):
acc = ComponentAccumulator()
+ # Need this later
+ #adv_time = kwargs.pop('advancedTiming')
+
tool = CompFactory.WaveformDigitisationTool(name="CaloWaveformDigtisationTool")
kwargs.setdefault("WaveformDigitisationTool", tool)
diff --git a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx
index c3f5313c0..7e9b6f485 100644
--- a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx
+++ b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx
@@ -1,9 +1,12 @@
#include "CaloWaveformDigiAlg.h"
-#include "Identifier/Identifier.h"
-
#include "FaserCaloSimEvent/CaloHitIdHelper.h"
+#include "CaloReadoutGeometry/EcalDetectorManager.h"
+
+#include "CaloReadoutGeometry/CaloDetectorElement.h"
+#include "GaudiKernel/PhysicalConstants.h"
+
#include <map>
#include <utility>
#include <cmath>
@@ -95,49 +98,178 @@ CaloWaveformDigiAlg::execute(const EventContext& ctx) const {
// Also save the baseline parameters
m_base_mean = m_digiCondTool->base_mean(ctx);
m_base_rms = m_digiCondTool->base_rms(ctx);
+
+ // And print out our timing simulation
+ if (m_advancedTiming) {
+ ATH_MSG_INFO("Using improved digitization timing");
+ } else {
+ ATH_MSG_INFO("Using simple digitization timing");
+ }
+ }
+
+ // Create TOF values
+ if (!m_time_of_flight.size()) {
+
+ // Detector manager to find global location
+ const CaloDD::EcalDetectorManager* detman;
+ ATH_CHECK(detStore()->retrieve(detman, "Ecal"));
+ m_time_of_flight = create_tof_map(m_ecalID, detman);
+
}
// Create structure to store pulse for each channel
- std::map<Identifier, std::vector<uint16_t>> waveforms = m_digiTool->create_waveform_map(m_ecalID);
+ std::map<Identifier, std::vector<uint16_t>> waveforms;
+
+ // This gives us empty vectors
+ waveforms = m_digiTool->create_waveform_map(m_ecalID);
+
+ // Identifiers we have seen hits from
+ std::set<Identifier> id_set;
// Sum energy for each channel (i.e. identifier)
- std::map<unsigned int, float> esum;
+ std::map<Identifier, float> esum;
+
+ // Histograms of E vs. t, make sure this is empty
+ m_digiTool->clear_evst_hist(m_evst_hist);
+
for (const auto& hit : *caloHitHandle) {
- esum[hit.identify()] += hit.energyLoss();
- }
- // Loop over time samples
- for (const auto& tk : m_timekernel) {
- std::map<unsigned int, float> counts;
+ float tof;
- // 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();
- //}
+ Identifier id = CaloHitIdHelper::GetHelper()->getIdentifier(hit.identify());
+ //id = hit.identify();
+ id_set.insert(id);
+ esum[id] += hit.energyLoss();
+ tof = m_time_of_flight[id];
- // Convolve summed energy with evaluated kernel for each hit id (i.e. channel)
- for (const auto& e : esum) {
- counts[e.first] = tk * e.second;
- }
+ m_digiTool->fill_evst_hist(m_evst_hist, id, hit.meanTime()-tof, hit.energyLoss());
- // Subtract count from basleine and add result to correct waveform vector
- for (const auto& c : counts) {
+ // Print out hit times
+ ATH_MSG_DEBUG("Scint Hit: " << id << " E: " << hit.energyLoss()
+ << " t: " << hit.meanTime()
+ << " x: " << hit.localStartPosition().x()
+ << " y: " << hit.localStartPosition().y()
+ << " z: " << hit.localStartPosition().z()
+ << " tof: " << tof
+ );
- float baseline = m_digiTool->generate_baseline(m_base_mean, m_base_rms);
- int value = std::round(baseline - c.second);
+ } // Done with loop over hits
- if (value < 0) {
- ATH_MSG_WARNING("Found pulse " << c.second << " larger than baseline " << c.first);
- value = 0; // Protect against scaling signal above baseline
+ // Choose how we construct the waveform
+ if (m_advancedTiming) {
+
+ // Loop over seen identifiers
+ for (const auto& id : id_set) {
+
+ // Convolve histogram with time kernel
+
+ // Use this vector to sum our values
+ std::vector<double> dwave;
+
+ // Loop over bins and fill with baseline values
+ int i = m_digiTool->digitizer_samples();
+ while(i--) { // Use while to avoid unused variable warning
+
+ // Push back baseline value so we can index this
+ float baseline = m_digiTool->generate_baseline(m_base_mean, m_base_rms);
+ // vector starts empty, push back each baseline element here
+ dwave.push_back(baseline);
+ } // End of loop over bins
+
+ // Get the histogram
+ auto it = m_evst_hist.find(id);
+ if (it == m_evst_hist.end()) {
+ ATH_MSG_WARNING("Didn't find EvsT hist for id " << id <<"!");
+ continue;
}
- // Convert hit id to Identifier and store
- Identifier id = CaloHitIdHelper::GetHelper()->getIdentifier(c.first);
- waveforms[id].push_back(value);
+ // Histogram of E vs. t filled above
+ auto h = it->second;
+
+ // Now, do it again and convolute any signals with the time kernal
+ for(unsigned int ibin=1; int(ibin) <= h->GetNbinsX(); ibin++) {
+
+ float value = h->GetBinContent(ibin);
+
+ // If bin is empty, move on
+ if (value <= 0.) continue;
+
+ // Location of histogram bin in digitizer array
+ // Must divide by oversampling
+ unsigned int jbin = (ibin-1)/m_digiTool->over_samples();
+
+ // Histogram doesn't start at t=0, must subtract pre_samples
+ int ik = jbin - m_digiTool->pre_samples();
+
+ // Convolute the histogram value with our time kernel
+ for (const auto& tk : m_timekernel) {
+
+ // Don't write off the end
+ if (ik < 0) {
+ ik++;
+ continue;
+ }
+
+ if (ik >= int(dwave.size())) break;
+
+ dwave[ik] -= (value * tk);
+
+ ik++; // Iterate
+ }
+ } // End of loop over histogram bins
+
+ // Finally protect against underflows and 'digitize' the values
+ for(unsigned int i=0; i<dwave.size(); i++) {
+
+ // waveforms is initialized empty, must push back each element here
+ if (dwave[i] < 0.) {
+ waveforms[id].push_back(0);
+ } else {
+ waveforms[id].push_back(std::round(dwave[i]));
+ }
+
+ // Debug weird overflow values
+ // if (waveforms[id][i] > m_base_mean + 5*m_base_rms) {
+ // ATH_MSG_WARNING("Found waveform id " << id << " bin " << i << " value " << waveforms[id][i] << "!");
+ // }
+
+ }
+
+ // Check what we have done
+ ATH_MSG_DEBUG("Waveform for ID " << id);
+ for(unsigned int i=400; i<450; i++)
+ ATH_MSG_DEBUG(2*i << "ns -> " << waveforms[id][i]);
+
+ } // End of loop over ids
+
+ } else {
+
+ // Loop over time samples
+ for (const auto& tk : m_timekernel) {
+ std::map<Identifier, float> counts;
+
+ // Convolve summed energy with evaluated kernel for each hit id (i.e. channel)
+ for (const auto& e : esum) {
+ counts[e.first] = tk * e.second;
+ }
+
+ // Subtract count from basleine and add result to correct waveform vector
+ for (const auto& c : counts) {
+
+ float baseline = m_digiTool->generate_baseline(m_base_mean, m_base_rms);
+ int value = std::round(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
+ waveforms[c.first].push_back(value);
+ }
}
}
-
// This is a bit of a hack to make sure all waveforms have
// at least baseline entries. Should really add this to the
// logic above
@@ -147,7 +279,7 @@ CaloWaveformDigiAlg::execute(const EventContext& ctx) const {
// Waveform was empty, fill with baseline
int channel = m_mappingTool->getChannelMapping(w.first);
ATH_MSG_DEBUG("Writing baseline into empty waveform in channel "<< channel);
- int i = m_digiTool->nsamples();
+ int i = m_digiTool->digitizer_samples();
while(i--) { // Use while to avoid unused variable warning with for
int baseline = m_digiTool->generate_baseline(m_base_mean, m_base_rms);
waveforms[w.first].push_back(baseline);
@@ -158,13 +290,13 @@ CaloWaveformDigiAlg::execute(const EventContext& ctx) const {
//m_chrono->chronoStart("Write");
// Loop over wavefrom vectors to make and store waveform
- unsigned int nsamples = m_digiTool->nsamples();
+ unsigned int digitizer_samples = m_digiTool->digitizer_samples();
for (const auto& w : waveforms) {
RawWaveform* wfm = new RawWaveform();
wfm->setWaveform(0, w.second);
wfm->setIdentifier(w.first);
wfm->setChannel(m_mappingTool->getChannelMapping(w.first));
- wfm->setSamples(nsamples);
+ wfm->setSamples(digitizer_samples);
waveformContainerHandle->push_back(wfm);
}
@@ -172,5 +304,34 @@ CaloWaveformDigiAlg::execute(const EventContext& ctx) const {
ATH_MSG_DEBUG("WaveformsHitContainer " << waveformContainerHandle.name() << "' filled with "<< waveformContainerHandle->size() <<" items");
+ // Debugging printout
+ for (auto it : m_evst_hist) {
+ ATH_MSG_DEBUG("Nonzero bins Evst hist id = " << it.first);
+ for (unsigned int ibin=0; int(ibin) <= it.second->GetNbinsX(); ibin++) {
+ if (it.second->GetBinContent(ibin) == 0) continue;
+ ATH_MSG_DEBUG(" bin: " << ibin << " t:" << it.second->GetBinCenter(ibin) << " val: " << it.second->GetBinContent(ibin));
+ }
+ }
+
+ // Cleanup any histograms we made
+ m_digiTool->clear_evst_hist(m_evst_hist);
+
return StatusCode::SUCCESS;
}
+
+template <class ID, class DM>
+std::map<Identifier, float>
+CaloWaveformDigiAlg::create_tof_map(const ID* idHelper, const DM* detman) const {
+
+ std::map<Identifier, float> time_of_flight;
+ for (auto itr = idHelper->pmt_begin(); itr != idHelper->pmt_end(); ++itr) {
+ const ExpandedIdentifier& extId = *itr;
+ Identifier pmtid = idHelper->pmt_id(extId);
+
+ CaloDD::CaloDetectorElement* element = detman->getDetectorElement(pmtid);
+ time_of_flight[pmtid] = element->center().z() / CLHEP::c_light;
+
+ ATH_MSG_DEBUG("Found TOF for " << pmtid << " = " << time_of_flight[pmtid]);
+ }
+ return time_of_flight;
+}
diff --git a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h
index 7ee25a271..c8d381ec2 100644
--- a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h
+++ b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h
@@ -23,6 +23,7 @@
// Helpers
#include "FaserCaloIdentifier/EcalID.h"
+#include "Identifier/Identifier.h"
// ROOT
#include "TF1.h"
@@ -45,6 +46,10 @@ class CaloWaveformDigiAlg: public AthReentrantAlgorithm {
virtual StatusCode finalize() override;
//@}
+ //
+ // Simulate detailed timing of waveforms
+ BooleanProperty m_advancedTiming{this, "AdvancedTiming", true};
+
private:
/** @name Disallow default instantiation, copy, assignment */
@@ -63,6 +68,8 @@ class CaloWaveformDigiAlg: public AthReentrantAlgorithm {
//@{
mutable TF1* m_kernel;
mutable std::vector<float> m_timekernel;
+ mutable EvstHistMap m_evst_hist;
+ mutable std::map<Identifier, float> m_time_of_flight;
//@}
/// Detector ID helper
@@ -105,6 +112,9 @@ class CaloWaveformDigiAlg: public AthReentrantAlgorithm {
{this, "WaveformContainerKey", ""};
//@}
+ template <class ID, class DM>
+ std::map<Identifier, float> create_tof_map(const ID* idHelper, const DM* detman) const;
+
};
diff --git a/Control/CalypsoExample/Digitization/scripts/faser_digi.py b/Control/CalypsoExample/Digitization/scripts/faser_digi.py
index 0da1361ce..81a097617 100755
--- a/Control/CalypsoExample/Digitization/scripts/faser_digi.py
+++ b/Control/CalypsoExample/Digitization/scripts/faser_digi.py
@@ -27,6 +27,8 @@ parser.add_argument("-t", "--tag", default="",
help="Specify digi tag (to append to output filename)")
parser.add_argument("--subtractTime", type=float, default=-999.,
help="Subtract time parameter for SCT RDOs")
+parser.add_argument("--simpleTiming", action="store_true",
+ help="Use simple waveform time digitization")
parser.add_argument("--digiTag", default="",
help="Specify tag for waveform digi folder")
parser.add_argument("--short", default="",
@@ -169,10 +171,10 @@ pualg.PileUpTools['FaserSCT_DigitizationTool'].SurfaceChargesGenerator.SubtractT
# Pass something to set folder tag
from CaloDigiAlgs.CaloDigiAlgsConfig import CaloWaveformDigitizationCfg
-acc.merge(CaloWaveformDigitizationCfg(ConfigFlags, digiTag=args.digiTag))
+acc.merge(CaloWaveformDigitizationCfg(ConfigFlags, digiTag=args.digiTag, AdvancedTiming=(not args.simpleTiming)))
from ScintDigiAlgs.ScintDigiAlgsConfig import ScintWaveformDigitizationCfg
-acc.merge(ScintWaveformDigitizationCfg(ConfigFlags, digiTag=args.digiTag))
+acc.merge(ScintWaveformDigitizationCfg(ConfigFlags, digiTag=args.digiTag, AdvancedTiming=(not args.simpleTiming)))
# Configure verbosity
if args.verbose:
diff --git a/Control/CalypsoExample/Digitization/scripts/faser_digi_merge.py b/Control/CalypsoExample/Digitization/scripts/faser_digi_merge.py
index b945ffbbd..9d931d3d2 100755
--- a/Control/CalypsoExample/Digitization/scripts/faser_digi_merge.py
+++ b/Control/CalypsoExample/Digitization/scripts/faser_digi_merge.py
@@ -35,6 +35,8 @@ parser.add_argument("-t", "--tag", default="",
help="Specify digi tag (to append to output filename)")
parser.add_argument("--subtractTime", type=float, default=-999.,
help="Subtract time parameter for SCT RDOs")
+parser.add_argument("--simpleTiming", action="store_true",
+ help="Use simple waveform time digitization")
parser.add_argument("--digiTag", default="",
help="Specify tag for waveform digi folder")
parser.add_argument("--short", default="",
@@ -239,10 +241,10 @@ pualg.PileUpTools['FaserSCT_DigitizationTool'].SurfaceChargesGenerator.SubtractT
# Pass something to set folder tag
from CaloDigiAlgs.CaloDigiAlgsConfig import CaloWaveformDigitizationCfg
-acc.merge(CaloWaveformDigitizationCfg(ConfigFlags, digiTag=args.digiTag))
+acc.merge(CaloWaveformDigitizationCfg(ConfigFlags, digiTag=args.digiTag, AdvancedTiming=(not args.simpleTiming)))
from ScintDigiAlgs.ScintDigiAlgsConfig import ScintWaveformDigitizationCfg
-acc.merge(ScintWaveformDigitizationCfg(ConfigFlags, digiTag=args.digiTag))
+acc.merge(ScintWaveformDigitizationCfg(ConfigFlags, digiTag=args.digiTag, AdvancedTiming=(not args.simpleTiming)))
# Configure verbosity
if args.verbose:
diff --git a/Scintillator/ScintDigiAlgs/CMakeLists.txt b/Scintillator/ScintDigiAlgs/CMakeLists.txt
index 5c1874e25..31950e073 100644
--- a/Scintillator/ScintDigiAlgs/CMakeLists.txt
+++ b/Scintillator/ScintDigiAlgs/CMakeLists.txt
@@ -11,7 +11,7 @@ atlas_add_component( ScintDigiAlgs
src/components/*.cxx
LINK_LIBRARIES AthenaBaseComps Identifier ScintIdentifier
WaveformConditionsToolsLib StoreGateLib WaveRawEvent
- ScintSimEvent WaveDigiToolsLib)
+ ScintReadoutGeometry ScintSimEvent WaveDigiToolsLib)
atlas_install_python_modules( python/*.py )
diff --git a/Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py b/Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py
index 538f8d2a3..5ed331105 100644
--- a/Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py
+++ b/Scintillator/ScintDigiAlgs/python/ScintDigiAlgsConfig.py
@@ -20,13 +20,13 @@ def ScintWaveformDigitizationCfg(flags, **kwargs):
kwargs.pop("digiTag")
if "TB" in flags.GeoModel.FaserVersion:
- acc.merge(ScintWaveformDigiCfg(flags, "VetoWaveformDigiAlg", "Veto"))
- acc.merge(ScintWaveformDigiCfg(flags, "PreshowerWaveformDigiAlg", "Preshower"))
+ acc.merge(ScintWaveformDigiCfg(flags, "VetoWaveformDigiAlg", "Veto", **kwargs))
+ acc.merge(ScintWaveformDigiCfg(flags, "PreshowerWaveformDigiAlg", "Preshower", **kwargs))
else:
- acc.merge(ScintWaveformDigiCfg(flags, "TriggerWaveformDigiAlg", "Trigger"))
- acc.merge(ScintWaveformDigiCfg(flags, "VetoWaveformDigiAlg", "Veto"))
- acc.merge(ScintWaveformDigiCfg(flags, "VetoNuWaveformDigiAlg", "VetoNu"))
- acc.merge(ScintWaveformDigiCfg(flags, "PreshowerWaveformDigiAlg", "Preshower"))
+ acc.merge(ScintWaveformDigiCfg(flags, "TriggerWaveformDigiAlg", "Trigger", **kwargs))
+ acc.merge(ScintWaveformDigiCfg(flags, "VetoWaveformDigiAlg", "Veto", **kwargs))
+ acc.merge(ScintWaveformDigiCfg(flags, "VetoNuWaveformDigiAlg", "VetoNu", **kwargs))
+ acc.merge(ScintWaveformDigiCfg(flags, "PreshowerWaveformDigiAlg", "Preshower", **kwargs))
acc.merge(ScintWaveformDigitizationOutputCfg(flags))
acc.merge(WaveformCableMappingCfg(flags))
@@ -39,7 +39,7 @@ def ScintWaveformDigiCfg(flags, name="ScintWaveformDigiAlg", source="", **kwargs
acc = ComponentAccumulator()
- tool = CompFactory.WaveformDigitisationTool(name=source+"WaveformDigtisationTool", **kwargs)
+ tool = CompFactory.WaveformDigitisationTool(name=source+"WaveformDigtisationTool")
kwargs.setdefault("WaveformDigitisationTool", tool)
kwargs.setdefault("ScintHitContainerKey", source+"Hits")
diff --git a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
index 33bd45d7d..c007f61ba 100644
--- a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
+++ b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
@@ -1,6 +1,13 @@
#include "ScintWaveformDigiAlg.h"
-#include "ScintSimEvent/ScintHitIdHelper.h"
+#include "ScintReadoutGeometry/VetoDetectorManager.h"
+#include "ScintReadoutGeometry/VetoNuDetectorManager.h"
+#include "ScintReadoutGeometry/TriggerDetectorManager.h"
+#include "ScintReadoutGeometry/PreshowerDetectorManager.h"
+
+#include "ScintReadoutGeometry/ScintDetectorElement.h"
+
+#include "GaudiKernel/PhysicalConstants.h"
#include <map>
#include <cmath>
@@ -103,14 +110,13 @@ ScintWaveformDigiAlg::execute(const EventContext& ctx) const {
ATH_MSG_INFO("Using simple digitization timing");
}
}
-
+
// Create structure to store pulse for each channel
std::map<Identifier, std::vector<uint16_t>> waveforms;
-
- // This gives us empty vectors
+
auto first = *scintHitHandle->begin();
if (first.isVeto()) {
- waveforms = m_digiTool->create_waveform_map(m_vetoID);
+ waveforms = m_digiTool->create_waveform_map(m_vetoID);
} else if (first.isVetoNu()) {
waveforms = m_digiTool->create_waveform_map(m_vetoNuID);
} else if (first.isTrigger()) {
@@ -119,7 +125,35 @@ ScintWaveformDigiAlg::execute(const EventContext& ctx) const {
waveforms = m_digiTool->create_waveform_map(m_preshowerID);
}
+ // Create TOF values
+ if (!m_time_of_flight.size()) {
+
+ // Detector manager to find global location
+ if (first.isVeto()) {
+ const ScintDD::VetoDetectorManager* detman;
+ ATH_CHECK(detStore()->retrieve(detman, "Veto"));
+ m_time_of_flight = create_tof_map(m_vetoID, detman);
+
+ } else if (first.isVetoNu()) {
+ const ScintDD::VetoNuDetectorManager* detman;
+ ATH_CHECK(detStore()->retrieve(detman, "VetoNu"));
+ m_time_of_flight = create_tof_map(m_vetoNuID, detman);
+
+ } else if (first.isTrigger()) {
+ const ScintDD::TriggerDetectorManager* detman;
+ ATH_CHECK(detStore()->retrieve(detman, "Trigger"));
+ m_time_of_flight = create_tof_map(m_triggerID, detman);
+
+ } else if (first.isPreshower()) {
+ const ScintDD::PreshowerDetectorManager* detman;
+ ATH_CHECK(detStore()->retrieve(detman, "Preshower"));
+ m_time_of_flight = create_tof_map(m_preshowerID, detman);
+
+ }
+ }
+
// Identifiers we have seen hits from
+ Identifier id;
std::set<Identifier> id_set;
// Sum energy for each channel (i.e. identifier)
@@ -128,12 +162,11 @@ ScintWaveformDigiAlg::execute(const EventContext& ctx) const {
// Histograms of E vs. t, make sure this is empty
m_digiTool->clear_evst_hist(m_evst_hist);
- Identifier id;
-
// Loop over all hits and make histograms of E vs. t for each channel
for (const auto& hit : *scintHitHandle) {
float dtime; // Time offset due to scint. geometry
+ float tof;
// Trigger plane, horizontal readout
if (first.isTrigger()) {
@@ -148,13 +181,16 @@ ScintWaveformDigiAlg::execute(const EventContext& ctx) const {
id = m_triggerID->pmt_id(plate_id, 0); // ID for PMT 0
id_set.insert(id);
esum[id] += hit.energyLoss();
- m_digiTool->fill_evst_hist(m_evst_hist, id, hit.meanTime()+dtime, hit.energyLoss());
- id_set.insert(id);
+ tof = m_time_of_flight[id];
+
+ m_digiTool->fill_evst_hist(m_evst_hist, id, hit.meanTime()+dtime-tof, hit.energyLoss());
id = m_triggerID->pmt_id(plate_id, 1); // ID for PMT 1
id_set.insert(id);
esum[id] += hit.energyLoss();
- m_digiTool->fill_evst_hist(m_evst_hist, id, hit.meanTime()-dtime, hit.energyLoss());
+ tof = m_time_of_flight[id];
+
+ m_digiTool->fill_evst_hist(m_evst_hist, id, hit.meanTime()-dtime-tof, hit.energyLoss());
} else {
// All others have only 1 PMT
@@ -162,19 +198,22 @@ ScintWaveformDigiAlg::execute(const EventContext& ctx) const {
id = hit.getIdentifier();
id_set.insert(id);
esum[id] += hit.energyLoss();
+ tof = m_time_of_flight[id];
// +y means shorter times
dtime = -m_digiCondTool->time_slope(ctx) * (hit.localStartPosition().y() + hit.localEndPosition().y()) / 2.;
- m_digiTool->fill_evst_hist(m_evst_hist, id, hit.meanTime()+dtime, hit.energyLoss());
+ m_digiTool->fill_evst_hist(m_evst_hist, id, hit.meanTime()+dtime-tof, hit.energyLoss());
}
+
// Print out hit times
ATH_MSG_DEBUG("Scint Hit: " << id << " E: " << hit.energyLoss()
- << " t: " << hit.meanTime()
- << " x: " << hit.localStartPosition().x()
- << " y: " << hit.localStartPosition().y()
- << " dt: " << dtime
- );
+ << " t: " << hit.meanTime()
+ << " x: " << hit.localStartPosition().x()
+ << " y: " << hit.localStartPosition().y()
+ << " dt: " << dtime
+ << " tof: " << tof
+ );
} // Done with loop over hits
@@ -186,43 +225,57 @@ ScintWaveformDigiAlg::execute(const EventContext& ctx) const {
// Convolve histogram with time kernel
- // Get the histogram
- auto it = m_evst_hist.find(id);
- if (it == m_evst_hist.end()) {
- ATH_MSG_WARNING("Didn't find EvsT hist for id " << id <<"!");
- continue;
- }
- auto h = it->second;
-
// Use this vector to sum our values
std::vector<double> dwave;
// Loop over bins and fill with baseline values
- for(unsigned int ibin=1; ibin <= h->GetNbinsX(); ibin++) {
+ int i = m_digiTool->digitizer_samples();
+ while(i--) { // Use while to avoid unused variable warning with for
// Push back baseline value so we can index this
float baseline = m_digiTool->generate_baseline(m_base_mean, m_base_rms);
- // Starts empty, push back each element here
+ // vector starts empty, push back each baseline element here
dwave.push_back(baseline);
-
} // End of loop over bins
+ // Get the histogram
+ auto it = m_evst_hist.find(id);
+ if (it == m_evst_hist.end()) {
+ ATH_MSG_WARNING("Didn't find EvsT hist for id " << id <<"!");
+ continue;
+ }
+
+ // Histogram of E vs. t filled above
+ auto h = it->second;
+
// Now, do it again and convolute any signals with the time kernal
- // Because root histos start with bin=1, subtract 1 from ibin
- for(unsigned int ibin=1; ibin <= h->GetNbinsX(); ibin++) {
+ for(unsigned int ibin=1; int(ibin) <= h->GetNbinsX(); ibin++) {
float value = h->GetBinContent(ibin);
// If bin is empty, move on
if (value <= 0.) continue;
+ // Location of histogram bin in digitizer array
+ // Must divide by oversampling
+ unsigned int jbin = (ibin-1)/m_digiTool->over_samples();
+
+ // Histogram doesn't start at t=0, must subtract pre_samples
+ int ik = jbin - m_digiTool->pre_samples();
+
// Convolute the histogram value with our time kernel
- unsigned int ik = 0;
for (const auto& tk : m_timekernel) {
+
// Don't write off the end
- if ((ibin + ik - 1) >= dwave.size()) break;
+ if (ik < 0) {
+ ik++;
+ continue;
+ }
+
+ if (ik >= int(dwave.size())) break;
- dwave[ibin+ik-1] -= (value * tk);
+ dwave[ik] -= (value * tk);
+
ik++; // Iterate
}
} // End of loop over histogram bins
@@ -291,7 +344,7 @@ ScintWaveformDigiAlg::execute(const EventContext& ctx) const {
// Waveform was empty, fill with baseline
int channel = m_mappingTool->getChannelMapping(w.first);
ATH_MSG_DEBUG("Writing baseline into empty waveform in channel "<< channel);
- int i = m_digiTool->nsamples();
+ int i = m_digiTool->digitizer_samples();
while(i--) { // Use while to avoid unused variable warning with for
float baseline = m_digiTool->generate_baseline(m_base_mean, m_base_rms);
waveforms[w.first].push_back(std::round(baseline));
@@ -299,7 +352,7 @@ ScintWaveformDigiAlg::execute(const EventContext& ctx) const {
}
// Loop over waveform vectors to make and store raw waveform
- unsigned int nsamples = m_digiTool->nsamples();
+ unsigned int nsamples = m_digiTool->digitizer_samples();
for (const auto& w : waveforms) {
RawWaveform* wfm = new RawWaveform();
wfm->setWaveform(0, w.second);
@@ -321,9 +374,9 @@ ScintWaveformDigiAlg::execute(const EventContext& ctx) const {
// Debugging printout
for (auto it : m_evst_hist) {
ATH_MSG_DEBUG("Nonzero bins Evst hist id = " << it.first);
- for (unsigned int ibin=0; ibin <= it.second->GetNbinsX(); ibin++) {
+ for (unsigned int ibin=0; int(ibin) <= it.second->GetNbinsX(); ibin++) {
if (it.second->GetBinContent(ibin) == 0) continue;
- ATH_MSG_DEBUG(" bin: " << ibin << " val: " << it.second->GetBinContent(ibin));
+ ATH_MSG_DEBUG(" bin: " << ibin << " t:" << it.second->GetBinCenter(ibin) << " val: " << it.second->GetBinContent(ibin));
}
}
@@ -334,3 +387,20 @@ ScintWaveformDigiAlg::execute(const EventContext& ctx) const {
}
+template <class ID, class DM>
+std::map<Identifier, float>
+ScintWaveformDigiAlg::create_tof_map(const ID* idHelper, const DM* detman) const {
+
+ std::map<Identifier, float> time_of_flight;
+ for (auto itr = idHelper->pmt_begin(); itr != idHelper->pmt_end(); ++itr) {
+ const ExpandedIdentifier& extId = *itr;
+ Identifier pmtid = idHelper->pmt_id(extId);
+
+ ScintDD::ScintDetectorElement* element = detman->getDetectorElement(pmtid);
+ time_of_flight[pmtid] = element->center().z() / CLHEP::c_light;
+
+ ATH_MSG_DEBUG("Found TOF for " << pmtid << " = " << time_of_flight[pmtid]);
+ }
+ return time_of_flight;
+}
+
diff --git a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h
index 7f76328a3..f439e982f 100644
--- a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h
+++ b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h
@@ -26,10 +26,8 @@
#include "ScintIdentifier/VetoNuID.h"
#include "ScintIdentifier/TriggerID.h"
#include "ScintIdentifier/PreshowerID.h"
-#include "ScintSimEvent/ScintHitIdHelper.h"
#include "Identifier/Identifier.h"
-
// ROOT
#include "TF1.h"
@@ -75,6 +73,7 @@ class ScintWaveformDigiAlg : public AthReentrantAlgorithm {
mutable TF1* m_kernel;
mutable std::vector<float> m_timekernel;
mutable EvstHistMap m_evst_hist;
+ mutable std::map<Identifier, float> m_time_of_flight;
//mutable std::map<Identifier, TH1D*> m_evst_hist;
//@}
@@ -121,6 +120,9 @@ class ScintWaveformDigiAlg : public AthReentrantAlgorithm {
{this, "WaveformContainerKey", ""};
//@}
+ template <class ID, class DM>
+ std::map<Identifier, float> create_tof_map(const ID* idHelper, const DM* detman) const;
+
};
diff --git a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h
index 2009b495d..814b25cd6 100644
--- a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h
+++ b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h
@@ -60,7 +60,10 @@ public:
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; }
+ virtual unsigned int digitizer_samples() const = 0;
+ virtual float digitizer_period() const = 0; // in ns
+ virtual unsigned int over_samples() const = 0;
+ virtual unsigned int pre_samples() const = 0;
/// Fill histograms
virtual void fill_evst_hist(EvstHistMap& map, const Identifier& id, float time, float energy) const = 0;
@@ -73,7 +76,6 @@ private:
protected:
TRandom3* m_random;
- unsigned int m_nsamples;
};
diff --git a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.icc b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.icc
index fd50ecfcf..140e4c5eb 100644
--- a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.icc
+++ b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.icc
@@ -10,10 +10,9 @@ std::map<Identifier, std::vector<uint16_t>> IWaveformDigitisationTool::create_w
const ExpandedIdentifier& extId = *itr;
Identifier id = idHelper->pmt_id(extId);
waveforms[id] = std::vector<uint16_t>();
- waveforms[id].reserve(m_nsamples);
+ waveforms[id].reserve(digitizer_samples());
}
return waveforms;
}
-
diff --git a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx
index 1b9ad0f44..a8f6207eb 100644
--- a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx
+++ b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx
@@ -21,7 +21,10 @@ StatusCode
WaveformDigitisationTool::initialize() {
ATH_MSG_INFO( name() << "::initalize()" );
- m_nsamples = 600;
+ ATH_MSG_INFO( m_digitizerPeriod);
+ ATH_MSG_INFO( m_digitizerSamples);
+ ATH_MSG_INFO( m_overSamples);
+
m_random = new TRandom3(0);
return StatusCode::SUCCESS;
@@ -32,10 +35,10 @@ std::vector<float>
WaveformDigitisationTool::evaluate_timekernel(TF1* kernel) const {
std::vector<float> timekernel;
- timekernel.reserve(m_nsamples);
+ timekernel.reserve(m_digitizerSamples);
- for (unsigned int i=0; i<m_nsamples; i++) {
- timekernel.push_back(kernel->Eval(2.*i));
+ for (unsigned int i=0; int(i) < m_digitizerSamples; i++) {
+ timekernel.push_back(kernel->Eval(m_digitizerPeriod*i));
}
return timekernel;
@@ -48,14 +51,19 @@ WaveformDigitisationTool::generate_baseline(float mean, float rms) const {
/// Fill EvsT histograms
+/// Time should be in ns
void
WaveformDigitisationTool::fill_evst_hist(EvstHistMap& map, const Identifier& id, float time, float energy) const {
if (map.find(id) == map.end()) {
ATH_MSG_DEBUG("Creating EvsT hist for channel" << id);
- map[id] = new TH1D("", "", m_nsamples, 0., 2.*m_nsamples);
+ // Start histogram before zero to pick up any early arriving hits
+ map[id] = new TH1D("", "", m_overSamples * (m_digitizerSamples+m_preSamples), -m_digitizerPeriod*m_preSamples, m_digitizerPeriod*m_digitizerSamples);
}
+ if (time < (-m_digitizerPeriod*m_preSamples))
+ ATH_MSG_WARNING("ID " << id << " hit found with time " << time << " E " << energy << "!");
+
map[id]->Fill(time, energy);
}
diff --git a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h
index ba73b3c8d..b36765493 100644
--- a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h
+++ b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h
@@ -39,6 +39,23 @@ class WaveformDigitisationTool: public extends<AthAlgTool, IWaveformDigitisation
void clear_evst_hist(EvstHistMap& map) const;
+ /// Waveform digitization parameters
+ FloatProperty m_digitizerPeriod { this, "DigitizerPeriod", 2.0, "Duration of one sample (in ns)" };
+ IntegerProperty m_digitizerSamples { this, "DigitizerSamples", 600, "Number of digitizer samples" };
+
+ /// This is used to set the time scale over which the hit energy
+ /// is histogramed before convoluting with the time kernel
+ IntegerProperty m_overSamples { this, "OverSamples", 8, "Oversampling of hit times compared to digitizer period" };
+
+ /// Digitizer samples to start before t0 (time = preSamples * digiPeriod)
+ IntegerProperty m_preSamples { this, "PreSamples", 10, "Presamples for energy histogram" };
+
+ /// Access functions for parameters
+ unsigned int digitizer_samples() const { return m_digitizerSamples; }
+ float digitizer_period() const { return m_digitizerPeriod; }
+ unsigned int over_samples() const { return m_overSamples; }
+ unsigned int pre_samples() const {return m_preSamples; }
+
private:
/// None
};
--
GitLab