diff --git a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx
index 7e9b6f4852eb2db02b458acf6af579fd1ada2b84..4d122f2e8ec6844496c4bff025779fca9045da1c 100644
--- a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx
+++ b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.cxx
@@ -33,11 +33,12 @@ CaloWaveformDigiAlg::initialize() {
// Set key to write container
ATH_CHECK( m_waveformContainerKey.initialize() );
- // Set up helper
- ATH_CHECK(detStore()->retrieve(m_ecalID, "EcalID"));
-
- // Need to set this on first event
- m_kernel = 0;
+ // And print out our timing simulation
+ if (m_advancedTiming) {
+ ATH_MSG_INFO("Using improved digitization timing");
+ } else {
+ ATH_MSG_INFO("Using simple digitization timing");
+ }
return StatusCode::SUCCESS;
}
@@ -46,9 +47,6 @@ StatusCode
CaloWaveformDigiAlg::finalize() {
ATH_MSG_INFO(name() << "::finalize()");
- if (m_kernel)
- delete m_kernel;
-
return StatusCode::SUCCESS;
}
@@ -74,221 +72,29 @@ CaloWaveformDigiAlg::execute(const EventContext& ctx) const {
return StatusCode::SUCCESS;
}
- if (!m_kernel) {
-
- ATH_MSG_INFO(name() << ": initialize waveform digitization kernel");
- ATH_MSG_INFO(" Norm: " << m_digiCondTool->cb_norm(ctx));
- ATH_MSG_INFO(" Mean: " << m_digiCondTool->cb_mean(ctx));
- ATH_MSG_INFO(" Sigma: " << m_digiCondTool->cb_sigma(ctx));
- ATH_MSG_INFO(" Alpha: " << m_digiCondTool->cb_alpha(ctx));
- ATH_MSG_INFO(" N: " << m_digiCondTool->cb_n(ctx));
-
- // Set up waveform shape
- m_kernel = new TF1("PDF", "[4] * ROOT::Math::crystalball_pdf(x, [0],[1],[2],[3])", 0, 1200);
-
- m_kernel->SetParameter(0, m_digiCondTool->cb_alpha(ctx));
- m_kernel->SetParameter(1, m_digiCondTool->cb_n(ctx));
- m_kernel->SetParameter(2, m_digiCondTool->cb_sigma(ctx));
- m_kernel->SetParameter(3, m_digiCondTool->cb_mean(ctx));
- m_kernel->SetParameter(4, m_digiCondTool->cb_norm(ctx));
-
- // Pre-evaluate time kernel for each bin
- m_timekernel = m_digiTool->evaluate_timekernel(m_kernel);
-
- // 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;
+ waveforms.clear();
- // This gives us empty vectors
- waveforms = m_digiTool->create_waveform_map(m_ecalID);
+ if (m_first) {
+ m_first = false;
- // 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);
-
- for (const auto& hit : *caloHitHandle) {
-
- float tof;
-
- Identifier id = CaloHitIdHelper::GetHelper()->getIdentifier(hit.identify());
- //id = hit.identify();
- id_set.insert(id);
- esum[id] += hit.energyLoss();
- tof = m_time_of_flight[id];
-
- m_digiTool->fill_evst_hist(m_evst_hist, id, hit.meanTime()-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()
- << " z: " << hit.localStartPosition().z()
- << " tof: " << tof
- );
-
- } // Done with loop over hits
+ // Write kernel parameters
+ ATH_MSG_INFO(name() << ": initialize waveform digitization kernel");
+ ATH_MSG_INFO(" Norm: " << m_digiCondTool->cb_norm());
+ ATH_MSG_INFO(" Mean: " << m_digiCondTool->cb_mean());
+ ATH_MSG_INFO(" Sigma: " << m_digiCondTool->cb_sigma());
+ ATH_MSG_INFO(" Alpha: " << m_digiCondTool->cb_alpha());
+ ATH_MSG_INFO(" N: " << m_digiCondTool->cb_n());
+ }
- // Choose how we construct the waveform
+ // Try our new tool instead
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;
- }
-
- // 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
-
+ waveforms = m_digiTool->generate_calo_timing_waveforms(ctx, m_digiCondTool, caloHitHandle.get());
} 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
- for (const auto& w : waveforms) {
- if (w.second.size() > 0) continue;
-
- // 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->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);
- }
+ waveforms = m_digiTool->generate_calo_waveforms(ctx, m_digiCondTool, caloHitHandle.get());
}
- //m_chrono->chronoStop("Digit");
- //m_chrono->chronoStart("Write");
-
// Loop over wavefrom vectors to make and store waveform
unsigned int digitizer_samples = m_digiTool->digitizer_samples();
for (const auto& w : waveforms) {
@@ -304,34 +110,9 @@ 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);
+ // Cleanup
+ waveforms.clear();
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 c8d381ec2c76c26b6913d9a5a13e8192ee597261..f874f174f992bc35bf9fa3db021bfdd5f30ef130 100644
--- a/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h
+++ b/Calorimeter/CaloDigiAlgs/src/CaloWaveformDigiAlg.h
@@ -21,13 +21,6 @@
#include "GaudiKernel/ServiceHandle.h"
#include "GaudiKernel/ToolHandle.h"
-// Helpers
-#include "FaserCaloIdentifier/EcalID.h"
-#include "Identifier/Identifier.h"
-
-// ROOT
-#include "TF1.h"
-
// STL
#include <string>
#include <vector>
@@ -59,51 +52,38 @@ class CaloWaveformDigiAlg: public AthReentrantAlgorithm {
CaloWaveformDigiAlg &operator=(const CaloWaveformDigiAlg&) = delete;
//@}
- // We use these a lot, so read them once from the conditions tool
- mutable float m_base_mean;
- mutable float m_base_rms;
-
- /** Kernel PDF and evaluated values **/
- // Mark these mutable, as they must be changed on first event...
- //@{
- 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
- const EcalID* m_ecalID{nullptr};
-
/**
* @name Digitisation tool
*/
+ //@{
ToolHandle<IWaveformDigitisationTool> m_digiTool
{this, "WaveformDigitisationTool", "WaveformDigitisationTool"};
+ //@}
/**
* @name Mapping tool
*/
+ //@{
ToolHandle<IWaveformCableMappingTool> m_mappingTool
{this, "WaveformCableMappingTool", "WaveformCableMappingTool"};
+ //@}
/**
* @name Digitization parameters tool
*/
+ //@{
ToolHandle<IWaveformDigiConditionsTool> m_digiCondTool
{this, "DigiConditionsTool", "CaloDigiConditionsTool"};
+ //@}
/**
* @name Input HITS using SG::ReadHandleKey
*/
//@{
-
SG::ReadHandleKey<CaloHitCollection> m_caloHitContainerKey
{this, "CaloHitContainerKey", ""};
-
//@}
-
/**
* @name Output data using SG::WriteHandleKey
*/
@@ -112,9 +92,7 @@ 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;
-
+ mutable bool m_first {true};
};
diff --git a/Calorimeter/CaloG4/EcalG4_SD/src/EcalSensorSDTool.cxx b/Calorimeter/CaloG4/EcalG4_SD/src/EcalSensorSDTool.cxx
index d59abfdae394dfe48196df0bddcc3aad1922b294..b9d9d3efae0ba53151d9edeaaa794102ea9315e6 100644
--- a/Calorimeter/CaloG4/EcalG4_SD/src/EcalSensorSDTool.cxx
+++ b/Calorimeter/CaloG4/EcalG4_SD/src/EcalSensorSDTool.cxx
@@ -21,7 +21,7 @@ EcalSensorSDTool::EcalSensorSDTool(const std::string& type, const std::string& n
, m_birk_c1 ( 0.013 * CLHEP::g/CLHEP::MeV/CLHEP::cm2 ) // 1st coef. of Birk's
, m_birk_c2 ( 9.6E-6 * CLHEP::g*CLHEP::g/CLHEP::MeV/CLHEP::MeV/CLHEP::cm2/CLHEP::cm2 ) // 2nd coef. of Birk's law
, m_birk_c1correction ( 0.57142857) //correction of c1 for 2e charged part.
- , m_a_local_outer_ecal ( 0.01 )
+ , m_a_local_outer_ecal ( 0.05 ) // sinusoidal non uniformity amplitude
, m_a_global_outer_ecal ( 0.03 ) // global non uniformity amplitude
, m_a_reflection_height ( 0.09 ) // reflection on the edges - height
, m_a_reflection_width ( 6. * CLHEP::mm ) // reflection on the edges - width
diff --git a/Scintillator/ScintDigiAlgs/CMakeLists.txt b/Scintillator/ScintDigiAlgs/CMakeLists.txt
index 31950e073b68352e35450d769e04521d13660731..d2de59c66529f8e9fb68e68a48259e764ecd96be 100644
--- a/Scintillator/ScintDigiAlgs/CMakeLists.txt
+++ b/Scintillator/ScintDigiAlgs/CMakeLists.txt
@@ -9,9 +9,9 @@ atlas_subdir( ScintDigiAlgs )
atlas_add_component( ScintDigiAlgs
src/*.cxx src/*.h
src/components/*.cxx
- LINK_LIBRARIES AthenaBaseComps Identifier ScintIdentifier
- WaveformConditionsToolsLib StoreGateLib WaveRawEvent
- ScintReadoutGeometry ScintSimEvent WaveDigiToolsLib)
+ LINK_LIBRARIES AthenaBaseComps Identifier StoreGateLib
+ ScintIdentifier ScintReadoutGeometry ScintSimEvent
+ WaveRawEvent WaveformConditionsToolsLib WaveDigiToolsLib)
atlas_install_python_modules( python/*.py )
diff --git a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
index c007f61ba17f48dddf1a0d48e8cec604f2b08be9..00d498b31ff73084f008048f6a643c2b3fabc699 100644
--- a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
+++ b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.cxx
@@ -34,15 +34,13 @@ ScintWaveformDigiAlg::initialize() {
// Set key to write container
ATH_CHECK( m_waveformContainerKey.initialize() );
- // Set up helpers
- ATH_CHECK(detStore()->retrieve(m_vetoID, "VetoID"));
- ATH_CHECK(detStore()->retrieve(m_vetoNuID, "VetoNuID"));
- ATH_CHECK(detStore()->retrieve(m_triggerID, "TriggerID"));
- ATH_CHECK(detStore()->retrieve(m_preshowerID, "PreshowerID"));
-
- // Need to set this on first event
- m_kernel = 0;
-
+ // And print out our timing simulation
+ if (m_advancedTiming) {
+ ATH_MSG_INFO("Using improved digitization timing");
+ } else {
+ ATH_MSG_INFO("Using simple digitization timing");
+ }
+
return StatusCode::SUCCESS;
}
@@ -50,9 +48,6 @@ StatusCode
ScintWaveformDigiAlg::finalize() {
ATH_MSG_INFO(name() << "::finalize()");
- if (m_kernel)
- delete m_kernel;
-
return StatusCode::SUCCESS;
}
@@ -78,277 +73,27 @@ ScintWaveformDigiAlg::execute(const EventContext& ctx) const {
return StatusCode::SUCCESS;
}
- if (!m_kernel) {
-
- ATH_MSG_INFO(name() << ": initialize waveform digitization kernel");
- ATH_MSG_INFO(" Norm: " << m_digiCondTool->cb_norm(ctx));
- ATH_MSG_INFO(" Mean: " << m_digiCondTool->cb_mean(ctx));
- ATH_MSG_INFO(" Sigma: " << m_digiCondTool->cb_sigma(ctx));
- ATH_MSG_INFO(" Alpha: " << m_digiCondTool->cb_alpha(ctx));
- ATH_MSG_INFO(" N: " << m_digiCondTool->cb_n(ctx));
-
- // Set up waveform shape
- m_kernel = new TF1("PDF", "[4] * ROOT::Math::crystalball_pdf(x, [0],[1],[2],[3])", 0, 1200);
-
- m_kernel->SetParameter(0, m_digiCondTool->cb_alpha(ctx));
- m_kernel->SetParameter(1, m_digiCondTool->cb_n(ctx));
- m_kernel->SetParameter(2, m_digiCondTool->cb_sigma(ctx));
- m_kernel->SetParameter(3, m_digiCondTool->cb_mean(ctx));
- m_kernel->SetParameter(4, m_digiCondTool->cb_norm(ctx));
-
- // Pre-evaluate time kernel for each bin
- m_timekernel = m_digiTool->evaluate_timekernel(m_kernel);
-
- // 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;
-
- auto first = *scintHitHandle->begin();
- if (first.isVeto()) {
- waveforms = m_digiTool->create_waveform_map(m_vetoID);
- } else if (first.isVetoNu()) {
- waveforms = m_digiTool->create_waveform_map(m_vetoNuID);
- } else if (first.isTrigger()) {
- waveforms = m_digiTool->create_waveform_map(m_triggerID);
- } else if (first.isPreshower()) {
- 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);
+ waveforms.clear();
- } else if (first.isTrigger()) {
- const ScintDD::TriggerDetectorManager* detman;
- ATH_CHECK(detStore()->retrieve(detman, "Trigger"));
- m_time_of_flight = create_tof_map(m_triggerID, detman);
+ if (m_first) {
+ m_first = false;
- } else if (first.isPreshower()) {
- const ScintDD::PreshowerDetectorManager* detman;
- ATH_CHECK(detStore()->retrieve(detman, "Preshower"));
- m_time_of_flight = create_tof_map(m_preshowerID, detman);
-
- }
+ // Write kernel parameters
+ ATH_MSG_INFO(name() << ": initialize waveform digitization kernel");
+ ATH_MSG_INFO(" Norm: " << m_digiCondTool->cb_norm());
+ ATH_MSG_INFO(" Mean: " << m_digiCondTool->cb_mean());
+ ATH_MSG_INFO(" Sigma: " << m_digiCondTool->cb_sigma());
+ ATH_MSG_INFO(" Alpha: " << m_digiCondTool->cb_alpha());
+ ATH_MSG_INFO(" N: " << m_digiCondTool->cb_n());
}
-
- // Identifiers we have seen hits from
- Identifier id;
- 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);
-
- // 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()) {
- 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();
- 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();
- 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
- // Use detector id (not hit id) throughout
- 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-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
- << " tof: " << tof
- );
-
- } // Done with loop over hits
-
- // Choose how we construct the waveform
+
+ // Try our new tool instead
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 with for
-
- // 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;
- }
-
- // 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
-
+ waveforms = m_digiTool->generate_scint_timing_waveforms(ctx, m_digiCondTool, scintHitHandle.get());
} 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);
- }
- }
- }
-
- // 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
- for (const auto& w : waveforms) {
- if (w.second.size() > 0) continue;
-
- // 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->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));
- }
+ waveforms = m_digiTool->generate_scint_waveforms(ctx, m_digiCondTool, scintHitHandle.get());
}
// Loop over waveform vectors to make and store raw waveform
@@ -371,36 +116,10 @@ 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; 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);
+ // Cleanup
+ waveforms.clear();
return StatusCode::SUCCESS;
}
-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 f439e982f0315859bb7a8a5eb7bdc76d32771988..feaba91ae4de7e35dc750cb6d91bb46b519a1d3e 100644
--- a/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h
+++ b/Scintillator/ScintDigiAlgs/src/ScintWaveformDigiAlg.h
@@ -21,16 +21,6 @@
#include "GaudiKernel/ServiceHandle.h"
#include "GaudiKernel/ToolHandle.h"
-//Helpers
-#include "ScintIdentifier/VetoID.h"
-#include "ScintIdentifier/VetoNuID.h"
-#include "ScintIdentifier/TriggerID.h"
-#include "ScintIdentifier/PreshowerID.h"
-#include "Identifier/Identifier.h"
-
-// ROOT
-#include "TF1.h"
-
// STL
#include <string>
#include <vector>
@@ -62,56 +52,38 @@ class ScintWaveformDigiAlg : public AthReentrantAlgorithm {
ScintWaveformDigiAlg &operator=(const ScintWaveformDigiAlg&) = delete;
//@}
- ///
-
- // We use these a lot, so read them once from the conditions tool
- 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, float> m_time_of_flight;
- //mutable std::map<Identifier, TH1D*> m_evst_hist;
- //@}
-
- /// Detector ID helpers
- const VetoID* m_vetoID{nullptr};
- const VetoNuID* m_vetoNuID{nullptr};
- const TriggerID* m_triggerID{nullptr};
- const PreshowerID* m_preshowerID{nullptr};
-
/**
* @name Digitisation tool
*/
+ //@{
ToolHandle<IWaveformDigitisationTool> m_digiTool
{this, "WaveformDigitisationTool", "WaveformDigitisationTool"};
+ //@}
/**
* @name Mapping tool
*/
+ //@{
ToolHandle<IWaveformCableMappingTool> m_mappingTool
{this, "WaveformCableMappingTool", "WaveformCableMappingTool"};
+ //@}
/**
* @name Digitization parameters tool
*/
+ //@{
ToolHandle<IWaveformDigiConditionsTool> m_digiCondTool
{this, "DigiConditionsTool", ""};
+ //@}
/**
* @name Input HITS using SG::ReadHandleKey
*/
//@{
-
SG::ReadHandleKey<ScintHitCollection> m_scintHitContainerKey
{this, "ScintHitContainerKey", ""};
-
//@}
-
/**
* @name Output data using SG::WriteHandleKey
*/
@@ -120,9 +92,7 @@ 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;
-
+ mutable bool m_first {true};
};
diff --git a/Waveform/WaveDigiTools/CMakeLists.txt b/Waveform/WaveDigiTools/CMakeLists.txt
index d7e9fd857b273c764dad780aceebd0d3f58acc3a..64329d8fcc00e60316bdf97f76d7f1b2619bb2cd 100644
--- a/Waveform/WaveDigiTools/CMakeLists.txt
+++ b/Waveform/WaveDigiTools/CMakeLists.txt
@@ -13,12 +13,16 @@ 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 Identifier
+ LINK_LIBRARIES AthenaBaseComps AthenaKernel GeoPrimitives
+ WaveRawEvent Identifier ScintSimEvent ScintReadoutGeometry
+ FaserCaloSimEvent CaloReadoutGeometry
+ WaveformConditionsToolsLib
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 814b25cd64855c1b0f054bf8ff92c12821663e27..765f7f0e85d18c69354eb910a9324da653d352aa 100644
--- a/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h
+++ b/Waveform/WaveDigiTools/WaveDigiTools/IWaveformDigitisationTool.h
@@ -22,6 +22,8 @@
#include "WaveRawEvent/RawWaveformContainer.h"
#include "WaveRawEvent/RawWaveform.h"
+#include "ScintSimEvent/ScintHitCollection.h"
+#include "FaserCaloSimEvent/CaloHitCollection.h"
#include "Identifier/Identifier.h"
@@ -33,6 +35,9 @@
#include <map>
#include <vector>
+// Forward declaration
+class IWaveformDigiConditionsTool;
+
typedef std::map<Identifier, TH1D*> EvstHistMap;
///Interface for waveform digitisation tools
@@ -49,15 +54,40 @@ public:
virtual ~IWaveformDigitisationTool() = default;
- /// Evaluate time kernel over time samples
- virtual std::vector<float> evaluate_timekernel(TF1* kernel) const = 0;
+ /// Create waveforms - old method
+ virtual std::map<Identifier, std::vector<uint16_t>>& generate_scint_waveforms(
+ const EventContext& ctx,
+ const ToolHandle<IWaveformDigiConditionsTool>& digiCondTool,
+ const ScintHitCollection* hitCollection) const = 0;
+
+ /// Create waveforms using hit timing
+ virtual std::map<Identifier, std::vector<uint16_t>>& generate_scint_timing_waveforms(
+ const EventContext& ctx,
+ const ToolHandle<IWaveformDigiConditionsTool>& digiCondTool,
+ const ScintHitCollection* hitCollection) const = 0;
+
+ /// Create waveforms - old method
+ virtual std::map<Identifier, std::vector<uint16_t>>& generate_calo_waveforms(
+ const EventContext& ctx,
+ const ToolHandle<IWaveformDigiConditionsTool>& digiCondTool,
+ const CaloHitCollection* hitCollection) const = 0;
+
+ /// Create waveforms using hit timing
+ virtual std::map<Identifier, std::vector<uint16_t>>& generate_calo_timing_waveforms(
+ const EventContext& ctx,
+ const ToolHandle<IWaveformDigiConditionsTool>& digiCondTool,
+ const CaloHitCollection* hitCollection) const = 0;
+
+
+ /* /// Evaluate time kernel over time samples */
+ /* virtual std::vector<float> evaluate_timekernel(TF1* kernel) const = 0; */
- /// Generate random baseline
- virtual float generate_baseline(float mean, float rms) const = 0;
+ /* /// Generate random baseline */
+ /* virtual float generate_baseline(float mean, float rms) const = 0; */
- /// Create structure to store pulse for each channel
+ /* /// 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;
+ std::map<Identifier, std::vector<uint16_t>> create_waveform_map(const T* idHelper) const;
/// Number of time samples
virtual unsigned int digitizer_samples() const = 0;
@@ -66,10 +96,10 @@ public:
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;
+ //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;
+ //virtual void clear_evst_hist(EvstHistMap& map) const = 0;
private:
ServiceHandle<IMessageSvc> m_msgSvc;
diff --git a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx
index a8f6207eb77ae2b05cae6adde3f73f314dcbb235..9e9500d07f0f81e6e9d99484ace9dcd30af8eed8 100644
--- a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx
+++ b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.cxx
@@ -10,6 +10,9 @@
#include "WaveformDigitisationTool.h"
+#include "ScintReadoutGeometry/ScintDetectorElement.h"
+#include "CaloReadoutGeometry/CaloDetectorElement.h"
+
// Constructor
WaveformDigitisationTool::WaveformDigitisationTool(const std::string& type, const std::string& name, const IInterface* parent) :
base_class(type, name, parent)
@@ -21,16 +24,30 @@ StatusCode
WaveformDigitisationTool::initialize() {
ATH_MSG_INFO( name() << "::initalize()" );
- ATH_MSG_INFO( m_digitizerPeriod);
- ATH_MSG_INFO( m_digitizerSamples);
- ATH_MSG_INFO( m_overSamples);
+ // Setup ID helpers
+ ATH_CHECK(detStore()->retrieve(m_vetoID, "VetoID"));
+ ATH_CHECK(detStore()->retrieve(m_vetoNuID, "VetoNuID"));
+ ATH_CHECK(detStore()->retrieve(m_triggerID, "TriggerID"));
+ ATH_CHECK(detStore()->retrieve(m_preshowerID, "PreshowerID"));
+ ATH_CHECK(detStore()->retrieve(m_ecalID, "EcalID"));
+
+ // Setup detector managers
+ ATH_CHECK(detStore()->retrieve(m_vetoDetMan, "Veto"));
+ ATH_CHECK(detStore()->retrieve(m_vetoNuDetMan, "VetoNu"));
+ ATH_CHECK(detStore()->retrieve(m_triggerDetMan, "Trigger"));
+ ATH_CHECK(detStore()->retrieve(m_preshowerDetMan, "Preshower"));
+ ATH_CHECK(detStore()->retrieve(m_caloDetMan, "Ecal"));
+
+ // Show our parameters (for debugging for now)
+ ATH_MSG_INFO(m_digitizerPeriod);
+ ATH_MSG_INFO(m_digitizerSamples);
+ ATH_MSG_INFO(m_overSamples);
m_random = new TRandom3(0);
return StatusCode::SUCCESS;
}
-
std::vector<float>
WaveformDigitisationTool::evaluate_timekernel(TF1* kernel) const {
@@ -79,4 +96,528 @@ WaveformDigitisationTool::clear_evst_hist(EvstHistMap& map) const {
map.clear();
}
+//template <class HitCollection>
+//std::map<Identifier, std::vector<uint16_t>>&
+//WaveformDigitisationTool::generateWaveforms(
+// const EventContext& ctx,
+// const ToolHandle<IWaveformDigiConditionsTool>& digiCondTool,
+// const HitCollection* hitCollection) const {
+
+//std::map<Identifier, std::vector<uint16_t>>&
+//WaveformDigitisationTool::test_context(const EventContext& ctx) const {
+// static std::map<Identifier, std::vector<uint16_t>> waveforms;
+// return waveforms;
+//}
+
+TF1*
+WaveformDigitisationTool::create_kernel(const EventContext& ctx,
+ const ToolHandle<IWaveformDigiConditionsTool>& digiCondTool) const {
+
+ TF1* kernel = new TF1("PDF", "[4] * ROOT::Math::crystalball_pdf(x, [0],[1],[2],[3])", 0, 1200);
+
+ kernel->SetParameter(0, digiCondTool->cb_alpha(ctx));
+ kernel->SetParameter(1, digiCondTool->cb_n(ctx));
+ kernel->SetParameter(2, digiCondTool->cb_sigma(ctx));
+ kernel->SetParameter(3, digiCondTool->cb_mean(ctx));
+ kernel->SetParameter(4, digiCondTool->cb_norm(ctx));
+
+ return kernel;
+}
+
+std::map<Identifier, std::vector<uint16_t>>&
+WaveformDigitisationTool::generate_scint_waveforms(
+ const EventContext& ctx,
+ const ToolHandle<IWaveformDigiConditionsTool>& digiCondTool,
+ const ScintHitCollection* hitCollection) const {
+
+ // Create digitization kernel
+ TF1* kernel = create_kernel(ctx, digiCondTool);
+
+ // Helpers, to be filled below as needed
+
+ // Create structure to store waveform for each channel
+ // We return a reference to this, so it must be static
+ static std::map<Identifier, std::vector<uint16_t>> waveforms;
+
+ // Make sure this is empty
+ waveforms.clear();
+
+ auto first = hitCollection->begin();
+
+ if (first->isVeto()) {
+ waveforms = create_waveform_map(m_vetoID);
+ } else if (first->isVetoNu()) {
+ waveforms = create_waveform_map(m_vetoNuID);
+ } else if (first->isTrigger()) {
+ waveforms = create_waveform_map(m_triggerID);
+ } else if (first->isPreshower()) {
+ waveforms = create_waveform_map(m_preshowerID);
+ }
+
+ // Make energy sums
+ Identifier id;
+ std::map<Identifier, float> esum;
+
+ for (const auto& hit : *hitCollection) {
+
+ if (hit.isTrigger()) {
+ Identifier plate_id = m_triggerID->plate_id(hit.getIdentifier());
+
+ id = m_triggerID->pmt_id(plate_id, 0); // ID for PMT 0
+ esum[id] += hit.energyLoss();
+
+ id = m_triggerID->pmt_id(plate_id, 1); // ID for PMT 1
+ esum[id] += hit.energyLoss();
+ } else {
+ id = hit.getIdentifier();
+ esum[id] += hit.energyLoss();
+ }
+ } // Done with loop over hits
+
+ // Now make waveforms
+ int value;
+ float signal, baseline;
+
+ for (const auto& w : waveforms) {
+
+ auto id = w.first;
+
+ // Now fill each entry with the sum of the kernel
+ // and some random background
+ for (unsigned int i=0; int(i) < m_digitizerSamples; i++) {
+ signal = kernel->Eval(m_digitizerPeriod * i) * esum[id];
+ baseline = generate_baseline(digiCondTool->base_mean(ctx),
+ digiCondTool->base_rms(ctx));
+
+ // Waveform is some random baseline plus a negative signal
+ value = std::round(baseline - signal);
+
+ // Don't let value go below zero
+ if (value < 0)
+ waveforms[id].push_back(0);
+ else
+ waveforms[id].push_back(value);
+
+ } // End of loop over bins
+
+ } // End of loop over channels
+
+ // Clean up after ourselves
+ delete kernel;
+
+ return waveforms;
+}
+
+std::map<Identifier, std::vector<uint16_t>>&
+WaveformDigitisationTool::generate_scint_timing_waveforms(
+ const EventContext& ctx,
+ const ToolHandle<IWaveformDigiConditionsTool>& digiCondTool,
+ const ScintHitCollection* hitCollection) const {
+
+ // Create structure to store waveform for each channel
+ // We return a reference to this, so it must be static
+ static std::map<Identifier, std::vector<uint16_t>> waveforms;
+
+ // Make sure this is empty
+ waveforms.clear();
+
+ // Create TOF values
+ std::map<Identifier, float> time_of_flight;
+
+ auto first = hitCollection->begin();
+
+ if (first->isVeto()) {
+ waveforms = create_waveform_map(m_vetoID);
+ time_of_flight = create_tof_map(m_vetoID, m_vetoDetMan);
+ } else if (first->isVetoNu()) {
+ waveforms = create_waveform_map(m_vetoNuID);
+ time_of_flight = create_tof_map(m_vetoNuID, m_vetoNuDetMan);
+ } else if (first->isTrigger()) {
+ waveforms = create_waveform_map(m_triggerID);
+ time_of_flight = create_tof_map(m_triggerID, m_triggerDetMan);
+ } else if (first->isPreshower()) {
+ waveforms = create_waveform_map(m_preshowerID);
+ time_of_flight = create_tof_map(m_preshowerID, m_preshowerDetMan);
+ }
+
+ // Make energy sums
+ Identifier id;
+
+ EvstHistMap evst_hist;
+
+ // Make sure we clear this
+ clear_evst_hist(evst_hist);
+
+ for (const auto& hit : *hitCollection) {
+
+ float dtime; // Time offset due to light propogation time in scintillator
+ float tof; // Time offset due to detector z position
+
+ // Trigger plane, horizontal readout
+ if (hit.isTrigger()) {
+ Identifier plate_id = m_triggerID->plate_id(hit.getIdentifier());
+
+ // Time delay based on horizontal distance
+ // Checked sign with Deion
+ dtime = digiCondTool->time_slope(ctx) * (hit.localStartPosition().x()+hit.localEndPosition().x()) / 2.;
+ id = m_triggerID->pmt_id(plate_id, 0); // ID for PMT 0
+ tof = time_of_flight[id];
+
+ // Fill our histogram for this channel
+ fill_evst_hist(evst_hist, id, hit.meanTime()-dtime-tof, hit.energyLoss());
+
+ id = m_triggerID->pmt_id(plate_id, 1); // ID for PMT 1
+ tof = time_of_flight[id];
+ fill_evst_hist(evst_hist, id, hit.meanTime()+dtime-tof, hit.energyLoss());
+
+ // Could also change amplitude based on distance to PMT, but not yet
+
+ } else {
+
+ // All others have only 1 PMT, use vertical position for dtime
+
+ id = hit.getIdentifier();
+ tof = time_of_flight[id];
+
+ // +y means shorter times
+ dtime = -digiCondTool->time_slope(ctx) * (hit.localStartPosition().y() + hit.localEndPosition().y()) / 2.;
+ fill_evst_hist(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
+ << " tof: " << tof
+ );
+
+ } // Done with loop over hits
+
+ // Debugging printout
+ for (auto it : 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));
+ }
+ }
+
+ // Now make waveforms
+ convolve_waveforms(ctx, digiCondTool, evst_hist, waveforms);
+
+ // Check what we have done
+ for (const auto& w : waveforms) {
+ ATH_MSG_DEBUG("Waveform for ID " << w.first);
+ for(unsigned int i=400; i<450; i++)
+ ATH_MSG_DEBUG(digitizer_period()*i << "ns -> " << w.second[i]);
+ }
+
+ // Cleanup
+ clear_evst_hist(evst_hist);
+
+ return waveforms;
+}
+
+std::map<Identifier, std::vector<uint16_t>>&
+WaveformDigitisationTool::generate_calo_waveforms(
+ const EventContext& ctx,
+ const ToolHandle<IWaveformDigiConditionsTool>& digiCondTool,
+ const CaloHitCollection* hitCollection) const {
+
+ // Create digitization kernel
+ TF1* kernel = create_kernel(ctx, digiCondTool);
+
+ // Helpers, to be filled below as needed
+
+ // Create structure to store waveform for each channel
+ // We return a reference to this, so it must be static
+ static std::map<Identifier, std::vector<uint16_t>> waveforms;
+
+ // Make sure this is empty
+ waveforms.clear();
+
+ waveforms = create_waveform_map(m_ecalID);
+
+ // Make energy sums
+ Identifier id;
+ std::map<Identifier, float> esum;
+
+ for (const auto& hit : *hitCollection) {
+ id = hit.getIdentifier();
+ esum[id] += hit.energyLoss();
+ } // Done with loop over hits
+
+ // Now make waveforms
+ int value;
+ float signal, baseline;
+
+ for (const auto& w : waveforms) {
+
+ auto id = w.first;
+
+ // Now fill each entry with the sum of the kernel
+ // and some random background
+ for (unsigned int i=0; int(i) < m_digitizerSamples; i++) {
+ signal = kernel->Eval(m_digitizerPeriod * i) * esum[id];
+ baseline = generate_baseline(digiCondTool->base_mean(ctx),
+ digiCondTool->base_rms(ctx));
+
+ // Waveform is some random baseline plus a negative signal
+ value = std::round(baseline - signal);
+
+ // Don't let value go below zero
+ if (value < 0)
+ waveforms[id].push_back(0);
+ else
+ waveforms[id].push_back(value);
+
+ } // End of loop over bins
+
+ } // End of loop over channels
+
+ // Clean up after ourselves
+ delete kernel;
+
+ return waveforms;
+}
+
+std::map<Identifier, std::vector<uint16_t>>&
+WaveformDigitisationTool::generate_calo_timing_waveforms(
+ const EventContext& ctx,
+ const ToolHandle<IWaveformDigiConditionsTool>& digiCondTool,
+ const CaloHitCollection* hitCollection) const {
+
+ // Create structure to store waveform for each channel
+ // We return a reference to this, so it must be static
+ static std::map<Identifier, std::vector<uint16_t>> waveforms;
+
+ // Make sure this is empty
+ waveforms.clear();
+
+ // Create TOF values
+ std::map<Identifier, float> time_of_flight;
+
+ waveforms = create_waveform_map(m_ecalID);
+ time_of_flight = create_calo_tof_map(m_ecalID, m_caloDetMan);
+
+ // Make energy sums
+ Identifier id;
+
+ EvstHistMap evst_hist;
+
+ // Make sure we clear this
+ clear_evst_hist(evst_hist);
+
+ for (const auto& hit : *hitCollection) {
+
+ float tof; // Time offset due to detector z position
+
+ id = hit.getIdentifier();
+ tof = time_of_flight[id];
+
+ fill_evst_hist(evst_hist, id, hit.meanTime()-tof, hit.energyLoss());
+
+ // Print out hit times
+ ATH_MSG_DEBUG("Calo Hit: " << id << " E: " << hit.energyLoss()
+ << " t: " << hit.meanTime()
+ << " x: " << hit.localStartPosition().x()
+ << " y: " << hit.localStartPosition().y()
+ << " z: " << hit.localStartPosition().z()
+ << " tof: " << tof
+ );
+
+ } // Done with loop over hits
+
+ // Debugging printout
+ for (auto it : 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));
+ }
+ }
+
+ // Now make waveforms
+ convolve_waveforms(ctx, digiCondTool, evst_hist, waveforms);
+
+ // Check what we have done
+ for (const auto& w : waveforms) {
+ ATH_MSG_DEBUG("Waveform for ID " << w.first);
+ for(unsigned int i=400; i<450; i++)
+ ATH_MSG_DEBUG(digitizer_period()*i << "ns -> " << w.second[i]);
+ }
+
+ // Cleanup
+ clear_evst_hist(evst_hist);
+
+ return waveforms;
+}
+
+
+template <class ID, class DM>
+std::map<Identifier, float>
+WaveformDigitisationTool::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;
+}
+
+std::map<Identifier, float>
+WaveformDigitisationTool::create_calo_tof_map(const EcalID* idHelper, const CaloDD::EcalDetectorManager* 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;
+}
+
+void
+WaveformDigitisationTool::convolve_waveforms(
+ const EventContext& ctx,
+ const ToolHandle<IWaveformDigiConditionsTool>& digiCondTool,
+ const EvstHistMap& evst_hist,
+ std::map<Identifier, std::vector<uint16_t>>& waveforms
+ ) const {
+
+ // Create digitization kernel
+ TF1* kernel = create_kernel(ctx, digiCondTool);
+
+ // Convolve Evst histogram vs. kernel for each waveform
+ for (const auto& w : waveforms) {
+
+ auto id = w.first;
+
+ // Convolve Evst histogram with time kernel
+
+ // Get the histogram
+ auto it = evst_hist.find(id);
+ if (it == evst_hist.end()) {
+ ATH_MSG_DEBUG("Didn't find EvsT hist for id " << id <<"!");
+
+ // Fill waveform with background only as Evst is empty
+ for (unsigned int i=0; i < digitizer_samples(); i++) {
+
+ float baseline = generate_baseline(digiCondTool->base_mean(ctx),
+ digiCondTool->base_rms(ctx));
+
+ // Waveform is some random baseline
+ int ivalue = std::round(baseline);
+
+ // Don't let value go below zero
+ if (ivalue < 0)
+ waveforms[id].push_back(0);
+ else
+ waveforms[id].push_back(ivalue);
+
+ } // End of loop over digitizer samples
+
+ // Done with this waveform
+ continue;
+ }
+
+ auto h = it->second;
+
+ // Use this vector to sum our values
+ // We will oversample by the same amount as the Evst histogram
+ std::vector<double> dwave(digitizer_samples() * over_samples());
+ // ATH_MSG_DEBUG("Creating dwave with length " << dwave.size());
+
+ // Convolute evst signals with the time kernel
+ for(unsigned int ih=1; int(ih) < h->GetNbinsX(); ih++) {
+
+ float value = h->GetBinContent(ih);
+
+ // If bin is empty, move on
+ if (value <= 0.) continue;
+
+ // Bin to write into in dwave
+ // Histogram doesn't start at t=0, so subtract pre_samples
+ int iw = ih - 1 - (pre_samples() * over_samples());
+
+ // Convolute the histogram with the time kernel
+ for (unsigned int ik=0; ik < dwave.size(); ik++) {
+ // Figure out destination bin first, skip if out of range
+ if ((iw + int(ik)) < 0) continue; // add int() to fix comp warnings
+ if ((iw + ik) >= dwave.size()) continue;
+
+ float ktime = ik * digitizer_period()/over_samples();
+ float kval = kernel->Eval(ktime) * value;
+ dwave[iw+ik] += kval;
+
+ // if (kval > 0.01)
+ // ATH_MSG_DEBUG(" Bin " << ih << " kernel bin " << ik << " time " << ktime << " value " << kval);
+ }
+
+ } // Done with convolution
+
+ // Debug printout
+ // ATH_MSG_DEBUG("dwave array for " << id);
+ // unsigned int iw=0;
+ // for (auto v : dwave) {
+ // iw++;
+ // if (v <= 0.1) continue;
+ // ATH_MSG_DEBUG(" bin "<< iw << " val " << v);
+ // }
+
+ // Need to down-sample dwave to fit waveform spacing
+ // Take average (or sum / over_sampling) to get amplitude correct
+ std::vector<double> dswave(digitizer_samples());
+ // ATH_MSG_DEBUG("Creating dswave with length " << dswave.size());
+
+ for (unsigned int i=0; i < dwave.size(); i++) {
+ dswave[i/over_samples()] += (dwave[i] / over_samples());
+ }
+
+ // Debug printout
+ // ATH_MSG_DEBUG("dswave array for " << id);
+ // iw=0;
+ // for (auto v : dswave) {
+ // iw++;
+ // if (v <= 0.1) continue;
+ // ATH_MSG_DEBUG(" bin "<< iw << " val " << v);
+ // }
+
+ // Now fill each entry of the waveform with the down-sampled signal
+ // and some random background
+ for (unsigned int i=0; i < digitizer_samples(); i++) {
+
+ float signal = dswave[i];
+ float baseline = generate_baseline(digiCondTool->base_mean(ctx),
+ digiCondTool->base_rms(ctx));
+
+ // Waveform is some random baseline plus a negative signal
+ int ivalue = std::round(baseline - signal);
+
+ // Don't let value go below zero
+ if (ivalue < 0)
+ waveforms[id].push_back(0);
+ else
+ waveforms[id].push_back(ivalue);
+
+ } // End of loop over dswave
+
+ } // End of loop over channels
+
+ // Clean up after ourselves
+ delete kernel;
+
+}
diff --git a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h
index b367654938984214f1d0de2d9f8d9b8a313eceb2..d71fa9ed6ef248fbb204a2200db5bb0daf4c1da4 100644
--- a/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h
+++ b/Waveform/WaveDigiTools/src/WaveformDigitisationTool.h
@@ -12,9 +12,31 @@
//Athena
#include "AthenaBaseComps/AthAlgTool.h"
#include "WaveDigiTools/IWaveformDigitisationTool.h"
+#include "WaveformConditionsTools/IWaveformDigiConditionsTool.h"
//Gaudi
#include "GaudiKernel/ToolHandle.h"
+#include "StoreGate/ReadHandleKey.h"
+
+// Data classes
+#include "ScintSimEvent/ScintHitCollection.h"
+
+// Helpers
+#include "ScintIdentifier/VetoID.h"
+#include "ScintIdentifier/VetoNuID.h"
+#include "ScintIdentifier/TriggerID.h"
+#include "ScintIdentifier/PreshowerID.h"
+#include "FaserCaloIdentifier/EcalID.h"
+
+// Detector Managers
+#include "ScintReadoutGeometry/VetoDetectorManager.h"
+#include "ScintReadoutGeometry/VetoNuDetectorManager.h"
+#include "ScintReadoutGeometry/TriggerDetectorManager.h"
+#include "ScintReadoutGeometry/PreshowerDetectorManager.h"
+#include "CaloReadoutGeometry/EcalDetectorManager.h"
+
+// ROOT
+#include "TF1.h"
//STL
@@ -45,7 +67,7 @@ class WaveformDigitisationTool: public extends<AthAlgTool, IWaveformDigitisation
/// 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" };
+ IntegerProperty m_overSamples { this, "OverSamples", 20, "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" };
@@ -56,8 +78,72 @@ class WaveformDigitisationTool: public extends<AthAlgTool, IWaveformDigitisation
unsigned int over_samples() const { return m_overSamples; }
unsigned int pre_samples() const {return m_preSamples; }
+ ///
+ /// New functions to handle entire digitization process
+
+ //template <class HitCollection>
+ //std::map<Identifier, std::vector<uint16_t>>& generateWaveforms(
+ // const EventContext& ctx,
+ // const ToolHandle<IWaveformDigiConditionsTool>&,
+ // const HitCollection*) const;
+
+ /// Simple (old-style) digitization
+ std::map<Identifier, std::vector<uint16_t>>& generate_scint_waveforms(
+ const EventContext& ctx,
+ const ToolHandle<IWaveformDigiConditionsTool>&,
+ const ScintHitCollection*) const;
+
+ /// Less simple (new-style) digitization with timing
+ std::map<Identifier, std::vector<uint16_t>>& generate_scint_timing_waveforms(
+ const EventContext& ctx,
+ const ToolHandle<IWaveformDigiConditionsTool>&,
+ const ScintHitCollection*) const;
+
+ /// Simple (old-style) digitization
+ std::map<Identifier, std::vector<uint16_t>>& generate_calo_waveforms(
+ const EventContext& ctx,
+ const ToolHandle<IWaveformDigiConditionsTool>&,
+ const CaloHitCollection*) const;
+
+ /// Less simple (new-style) digitization with timing
+ std::map<Identifier, std::vector<uint16_t>>& generate_calo_timing_waveforms(
+ const EventContext& ctx,
+ const ToolHandle<IWaveformDigiConditionsTool>&,
+ const CaloHitCollection*) const;
+
+ private:
+
+ template <class ID, class DM>
+ std::map<Identifier, float>
+ create_tof_map(const ID* idHelper, const DM* detman) const;
+
+ std::map<Identifier, float>
+ create_calo_tof_map(const EcalID* idHelper, const CaloDD::EcalDetectorManager* detman) const;
+
+ TF1* create_kernel(const EventContext& ctx,
+ const ToolHandle<IWaveformDigiConditionsTool>& digiCondTool) const;
+
+ void convolve_waveforms(
+ const EventContext& ctx,
+ const ToolHandle<IWaveformDigiConditionsTool>& digiCondTool,
+ const EvstHistMap& evst_hist,
+ std::map<Identifier, std::vector<uint16_t>>& waveforms) const;
+
private:
- /// None
+
+ /// ID helpers
+ const VetoID* m_vetoID{nullptr};
+ const VetoNuID* m_vetoNuID{nullptr};
+ const TriggerID* m_triggerID{nullptr};
+ const PreshowerID* m_preshowerID{nullptr};
+ const EcalID* m_ecalID{nullptr};
+
+ /// Detector manager helpers
+ const ScintDD::VetoDetectorManager* m_vetoDetMan{nullptr};
+ const ScintDD::VetoNuDetectorManager* m_vetoNuDetMan{nullptr};
+ const ScintDD::TriggerDetectorManager* m_triggerDetMan{nullptr};
+ const ScintDD::PreshowerDetectorManager* m_preshowerDetMan{nullptr};
+ const CaloDD::EcalDetectorManager* m_caloDetMan{nullptr};
};
#endif // WAVEDIGITOOLS_WAVEFORMDIGITISATIONTOOL_H
diff --git a/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformDigiConditionsTool.cxx b/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformDigiConditionsTool.cxx
index d16f72242af7e25143be762341bd6b89741d79c2..0eacb0a966c2aff1e2e82129ffe0c4f8c153eb7a 100644
--- a/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformDigiConditionsTool.cxx
+++ b/Waveform/WaveformConditions/WaveformConditionsTools/src/WaveformDigiConditionsTool.cxx
@@ -167,7 +167,7 @@ WaveformDigiConditionsTool::get_value(const EventContext& ctx, std::string arg)
if (payload.exists(arg) and not payload[arg].isNull()) {
val = payload[arg].data<float>();
- ATH_MSG_DEBUG("Found digi COOL channel " << m_cool_channel << " " << arg << " as " << val);
+ // ATH_MSG_DEBUG("Found digi COOL channel " << m_cool_channel << " " << arg << " as " << val);
} else {
ATH_MSG_WARNING("No valid " << arg << " value found for digi COOL channel "<< m_cool_channel<<"!");
}