From 0922c11483a0d5a6607005fa77ca8cfd1823ff22 Mon Sep 17 00:00:00 2001
From: Eric Torrence <eric.torrence@cern.ch>
Date: Tue, 5 Dec 2023 05:34:48 +0100
Subject: [PATCH] Add hit timing to waveform dititization
---
.../src/ScintWaveformDigiAlg.cxx | 213 +++++++++++++-----
.../ScintDigiAlgs/src/ScintWaveformDigiAlg.h | 8 +-
.../WaveDigiTools/IWaveformDigitisationTool.h | 9 +
.../IWaveformDigitisationTool.icc | 2 +
.../src/WaveformDigitisationTool.cxx | 27 +++
.../src/WaveformDigitisationTool.h | 7 +-
.../IWaveformDigiConditionsTool.h | 5 +
.../src/WaveformDigiConditionsTool.cxx | 11 +
.../src/WaveformDigiConditionsTool.h | 3 +
9 files changed, 229 insertions(+), 56 deletions(-)
diff --git a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
index 1869860e8..33bd45d7d 100644
--- a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
+++ b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
@@ -95,11 +95,19 @@ ScintWaveformDigiAlg::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 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);
@@ -109,79 +117,168 @@ ScintWaveformDigiAlg::execute(const EventContext& ctx) const {
waveforms = m_digiTool->create_waveform_map(m_triggerID);
} else if (first.isPreshower()) {
waveforms = m_digiTool->create_waveform_map(m_preshowerID);
- }
+ }
+
+ // Identifiers we have seen hits from
+ std::set<Identifier> id_set;
// Sum energy for each channel (i.e. identifier)
std::map<Identifier, float> esum;
+ // 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
+
+ // Trigger plane, horizontal readout
if (first.isTrigger()) {
Identifier plate_id = m_triggerID->plate_id(hit.getIdentifier());
+
+ // Time delay based on horizontal distance
+ // Don't know if I have the sign right here
+ // Assume +x is further from PMT 0
+ dtime = m_digiCondTool->time_slope(ctx) * (hit.localStartPosition().x()+hit.localEndPosition().x()) / 2.;
+
// Need to do something better here, but for now just fill both
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);
+
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());
+
} else {
// All others have only 1 PMT
// Use detector id (not hit id) throughout
id = hit.getIdentifier();
+ id_set.insert(id);
esum[id] += hit.energyLoss();
- }
- }
- // Loop over time samples
- for (const auto& tk : m_timekernel) {
- std::map<Identifier, float> counts;
-
-// // Convolve hit energy with evaluated kernel and sum for each hit id (i.e. channel)
-// Identifier id;
-// for (const auto& hit : *scintHitHandle) {
-//
-// // Special handling for trigger scintillator
-// if (first.isTrigger()) {
-// // Hits are created per plate
-// // Need to split between 2 PMTs
-//
-// // Identifier from hit is plate ID
-// Identifier plate_id = m_triggerID->plate_id(hit.getIdentifier());
-//
-// // Need to do something better here, but for now just fill both
-// id = m_triggerID->pmt_id(plate_id, 0); // ID for PMT 0
-// counts[id] += tk * hit.energyLoss();
-// id = m_triggerID->pmt_id(plate_id, 1); // ID for PMT 1
-// counts[id] += tk * hit.energyLoss();
-//
-// } else {
-// // All others have only 1 PMT
-// // Use detector id (not hit id) throughout
-// id = hit.getIdentifier();
-// counts[id] += tk * hit.energyLoss();
-// }
-// }
-
- // Convolve summed energy with evaluated kernel for each hit id (i.e. channel)
- for (const auto& e : esum) {
- // Convert hit id to Identifier and store
- id = e.first;
- counts[id] = tk * e.second;
+ // +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());
}
- // 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()
+ << " dt: " << dtime
+ );
- float baseline = m_digiTool->generate_baseline(m_base_mean, m_base_rms);
- int value = std::round(baseline - c.second);
+ } // Done with loop over hits
+
+ // Choose how we construct the waveform
+ if (m_advancedTiming) {
+
+ // Loop over seen identifiers
+ for (const auto& id : id_set) {
+
+ // 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++) {
+
+ // 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
+ dwave.push_back(baseline);
+
+ } // End of loop over bins
+
+ // 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++) {
+
+ float value = h->GetBinContent(ibin);
+
+ // If bin is empty, move on
+ if (value <= 0.) continue;
+
+ // 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;
+
+ dwave[ibin+ik-1] -= (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] << "!");
+ // }
- 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
- id = c.first;
- waveforms[id].push_back(value);
+ // 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) {
+ // Convert hit id to Identifier and store
+ id = e.first;
+ counts[id] = 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
+ id = c.first;
+ waveforms[id].push_back(value);
+ }
}
}
@@ -196,12 +293,12 @@ ScintWaveformDigiAlg::execute(const EventContext& ctx) const {
ATH_MSG_DEBUG("Writing baseline into empty waveform in channel "<< channel);
int i = m_digiTool->nsamples();
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);
+ float baseline = m_digiTool->generate_baseline(m_base_mean, m_base_rms);
+ waveforms[w.first].push_back(std::round(baseline));
}
}
- // Loop over wavefrom vectors to make and store raw waveform
+ // Loop over waveform vectors to make and store raw waveform
unsigned int nsamples = m_digiTool->nsamples();
for (const auto& w : waveforms) {
RawWaveform* wfm = new RawWaveform();
@@ -221,5 +318,19 @@ ScintWaveformDigiAlg::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; ibin <= it.second->GetNbinsX(); ibin++) {
+ if (it.second->GetBinContent(ibin) == 0) continue;
+ ATH_MSG_DEBUG(" bin: " << ibin << " val: " << it.second->GetBinContent(ibin));
+ }
+ }
+
+ // Cleanup any histograms we made
+ m_digiTool->clear_evst_hist(m_evst_hist);
+
return StatusCode::SUCCESS;
}
+
+
diff --git a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h
index 99c374f55..7f76328a3 100644
--- a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h
+++ b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h
@@ -33,7 +33,6 @@
// ROOT
#include "TF1.h"
-
// STL
#include <string>
#include <vector>
@@ -52,6 +51,9 @@ class ScintWaveformDigiAlg : public AthReentrantAlgorithm {
virtual StatusCode finalize() override;
//@}
+ //
+ // Simulate detailed timing of waveforms
+ BooleanProperty m_advancedTiming{this, "AdvancedTiming", true};
private:
@@ -68,13 +70,13 @@ class ScintWaveformDigiAlg : public AthReentrantAlgorithm {
mutable float m_base_mean;
mutable float m_base_rms;
-
/** Kernel PDF and evaluated values **/
//@{
mutable TF1* m_kernel;
mutable std::vector<float> m_timekernel;
+ mutable EvstHistMap m_evst_hist;
+ //mutable std::map<Identifier, TH1D*> m_evst_hist;
//@}
-
/// Detector ID helpers
const VetoID* m_vetoID{nullptr};
diff --git a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h
index 37cb7af18..2009b495d 100644
--- a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h
+++ b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h
@@ -27,11 +27,13 @@
#include "TF1.h"
#include "TRandom3.h"
+#include "TH1D.h"
#include <utility>
#include <map>
#include <vector>
+typedef std::map<Identifier, TH1D*> EvstHistMap;
///Interface for waveform digitisation tools
class IWaveformDigitisationTool : virtual public IAlgTool
@@ -60,12 +62,19 @@ public:
/// Number of time samples
unsigned int nsamples() const { return m_nsamples; }
+ /// Fill histograms
+ virtual void fill_evst_hist(EvstHistMap& map, const Identifier& id, float time, float energy) const = 0;
+
+ // Remove old histograms
+ virtual void clear_evst_hist(EvstHistMap& map) const = 0;
+
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 41b8c2650..fd50ecfcf 100644
--- a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.icc
+++ b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.icc
@@ -15,3 +15,5 @@ std::map<Identifier, std::vector<uint16_t>> IWaveformDigitisationTool::create_w
return waveforms;
}
+
+
diff --git a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx
index 3a48027d5..1b9ad0f44 100644
--- a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx
+++ b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx
@@ -45,3 +45,30 @@ float
WaveformDigitisationTool::generate_baseline(float mean, float rms) const {
return m_random->Gaus(mean, rms);
}
+
+
+/// Fill EvsT histograms
+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);
+ }
+
+ map[id]->Fill(time, energy);
+}
+
+/// Delete histograms and clear map
+void
+WaveformDigitisationTool::clear_evst_hist(EvstHistMap& map) const {
+
+ // Delete histogram objects
+ for (auto& it : map)
+ delete it.second;
+
+ // Also remove map elements (now null pointers)
+ map.clear();
+}
+
+
diff --git a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h
index b6554c337..ba73b3c8d 100644
--- a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h
+++ b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h
@@ -34,10 +34,13 @@ class WaveformDigitisationTool: public extends<AthAlgTool, IWaveformDigitisation
/// Generate random baseline
float generate_baseline(float mean, float rms) const;
+ /// Fill EvsT histograms
+ void fill_evst_hist(EvstHistMap& map, const Identifier& id, float time, float energy) const;
- private:
- // None
+ void clear_evst_hist(EvstHistMap& map) const;
+ private:
+ /// None
};
#endif // WAVEDIGITOOLS_WAVEFORMDIGITISATIONTOOL_H
diff --git a/Waveform/WaveformConditions/WaveformConditionsTools/WaveformConditionsTools/IWaveformDigiConditionsTool.h b/Waveform/WaveformConditions/WaveformConditionsTools/WaveformConditionsTools/IWaveformDigiConditionsTool.h
index 1bfaf0dd9..3f531beca 100644
--- a/Waveform/WaveformConditions/WaveformConditionsTools/WaveformConditionsTools/IWaveformDigiConditionsTool.h
+++ b/Waveform/WaveformConditions/WaveformConditionsTools/WaveformConditionsTools/IWaveformDigiConditionsTool.h
@@ -57,6 +57,11 @@ class IWaveformDigiConditionsTool: virtual public IAlgTool {
virtual float cb_n(void) const = 0;
virtual float cb_n(const EventContext& ctx) const = 0;
+
+ // Time dependence parameters
+ virtual float time_slope(void) const = 0;
+ virtual float time_slope(const EventContext& ctx) const = 0;
+
};
//----------------------------------------------------------------------
diff --git a/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformDigiConditionsTool.cxx b/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformDigiConditionsTool.cxx
index 4e2fedd97..d16f72242 100644
--- a/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformDigiConditionsTool.cxx
+++ b/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformDigiConditionsTool.cxx
@@ -93,6 +93,12 @@ WaveformDigiConditionsTool::cb_n(const EventContext& ctx) const {
return get_value(ctx, "cb_n");
}
+float
+WaveformDigiConditionsTool::time_slope(const EventContext& ctx) const {
+ return 1/230.; // 23 cm/ns
+ // return get_value(ctx, "time_slope");
+}
+
//----------------------------------------------------------------------
float
WaveformDigiConditionsTool::base_mean(void) const {
@@ -136,6 +142,11 @@ WaveformDigiConditionsTool::cb_n(void) const {
return cb_n(ctx);
}
+float
+WaveformDigiConditionsTool::time_slope(void) const {
+ const EventContext& ctx(Gaudi::Hive::currentContext());
+ return time_slope(ctx);
+}
//----------------------------------------------------------------------
float
diff --git a/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformDigiConditionsTool.h b/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformDigiConditionsTool.h
index e5e823208..eab1a65d7 100644
--- a/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformDigiConditionsTool.h
+++ b/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformDigiConditionsTool.h
@@ -67,6 +67,9 @@ class WaveformDigiConditionsTool: public extends<AthAlgTool, IWaveformDigiCondit
virtual float cb_n(void) const override;
virtual float cb_n(const EventContext& ctx) const override;
+ // Time dependence parameters
+ virtual float time_slope(void) const override;
+ virtual float time_slope(const EventContext& ctx) const override;
private:
--
GitLab