From 386c281a6ea6bfcb4e2d420282109447c0184046 Mon Sep 17 00:00:00 2001 From: Chao Peng Date: Wed, 15 Feb 2023 23:57:06 -0600 Subject: [PATCH 1/3] remove ScFiDigi --- .../calorimetry/CalorimeterScFiDigi.cc | 192 ------------------ .../calorimetry/CalorimeterScFiDigi.h | 91 --------- ...rimeterHit_factory_EcalBarrelScFiRawHits.h | 6 +- 3 files changed, 3 insertions(+), 286 deletions(-) delete mode 100644 src/algorithms/calorimetry/CalorimeterScFiDigi.cc delete mode 100644 src/algorithms/calorimetry/CalorimeterScFiDigi.h diff --git a/src/algorithms/calorimetry/CalorimeterScFiDigi.cc b/src/algorithms/calorimetry/CalorimeterScFiDigi.cc deleted file mode 100644 index 42e6800aa3..0000000000 --- a/src/algorithms/calorimetry/CalorimeterScFiDigi.cc +++ /dev/null @@ -1,192 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -// Copyright (C) 2022 Chao Peng, Wouter Deconinck, Sylvester Joosten, Barak Schmookler, David Lawrence - -// A digitization digitization algorithm specifically for CalorimeterHit from Scintillating Fibers -// 1. Fiber signals are merged if they are within the same light guide -// 2. Digitize the energy with dynamic ADC range and add pedestal (mean +- sigma) -// 3. Time conversion with smearing resolution (absolute value) -// TODO: Sum with timing; waveform implementation -// -// Author: Chao Peng -// Date: 02/12/2023 - - -#include "CalorimeterScFiDigi.h" - -#include -#include -#include -#include -using namespace dd4hep; - -// -// TODO: -// - Array type configuration parameters are not yet supported in JANA (needs to be added) -// - Random number service needs to be resolved (on global scale) - - - -//------------------------ -// AlgorithmInit -//------------------------ -void CalorimeterScFiDigi::AlgorithmInit(std::shared_ptr& logger) { - - // Assume all configuration parameter data members have been filled in already. - - // Gaudi implments a random number generator service. It is not clear to me how this - // can work. There are multiple race conditions that occur in parallel event processing: - // 1. The exact same events processed by a given thread in one invocation will not - // neccessarily be the combination of events any thread sees in a subsequest - // invocation. Thus, you can't rely on thread_local storage. - // 2. Its possible for the factory execution order to be modified by the presence of - // a processor (e.g. monitoring plugin). This is not as serious since changing the - // command line should cause one not to expect reproducibility. Still, one may - // expect the inclusion of an "observer" plugin not to have such side affects. - // - // More information will be needed. In the meantime, we implement a local random number - // generator. Ideally, this would be seeded with the run number+event number, but for - // now, just use default values defined in header file. - - // set energy resolution numbers - m_log=logger; - for (size_t i = 0; i < u_eRes.size() && i < 3; ++i) { - eRes[i] = u_eRes[i]; - } - - // using juggler internal units (GeV, mm, radian, ns) - tRes = m_tRes / dd4hep::ns; - stepTDC = dd4hep::ns / m_resolutionTDC; - - // need signal sum - if (!u_fields.empty()) { - - // sanity checks - if (!m_geoSvc) { - m_log->error("Unable to locate Geometry Service."); - throw std::runtime_error("Unable to locate Geometry Service."); - } - if (m_readout.empty()) { - m_log->error("readoutClass is not provided, it is needed to know the fields in readout ids."); - throw std::runtime_error("readoutClass is not provided."); - } - - // get decoders - try { - auto id_desc = m_geoSvc->detector()->readout(m_readout).idSpec(); - id_dec = id_desc.decoder(); - id_mask = 0; - std::vector> ref_fields; - for (size_t i = 0; i < u_fields.size(); ++i) { - id_mask |= id_desc.field(u_fields[i])->mask(); - // use the provided id number to find ref cell, or use 0 - int ref = i < u_refs.size() ? u_refs[i] : 0; - ref_fields.emplace_back(u_fields[i], ref); - } - // need to also merge z - if (!m_zsegment.empty()) { - z_idx = id_dec->index(m_zsegment); - // treat z segment as a field to merge - id_mask |= id_desc.field(m_zsegment)->mask(); - // but new z value will be determined in merging, 0 is a placeholder here - ref_fields.emplace_back(m_zsegment, 0); - } - ref_mask = id_desc.encode(ref_fields); - // debug() << fmt::format("Referece id mask for the fields {:#064b}", ref_mask) << endmsg; - // update z field index from the decoder - } catch (...) { - m_log->warn("Failed to load ID decoder for {}", m_readout); - japp->Quit(); - return; - } - id_mask = ~id_mask; - //LOG_INFO(default_cout_logger) << fmt::format("ID mask in {:s}: {:#064b}", m_readout, id_mask) << LOG_END; - m_log->info("ID mask in {:s}: {:#064b}", m_readout, id_mask); - } -} - - - -//------------------------ -// AlgorithmChangeRun -//------------------------ -void CalorimeterScFiDigi::AlgorithmChangeRun() { - /// This is automatically run before Process, when a new run number is seen - /// Usually we update our calibration constants by asking a JService - /// to give us the latest data for this run number -} - -//------------------------ -// AlgorithmProcess -//------------------------ -void CalorimeterScFiDigi::AlgorithmProcess() { - - // Delete any output objects left from last event. - // (Should already have been done for us, but just to be bullet-proof.) - for( auto h : rawhits ) delete h; - rawhits.clear(); - - light_guide_digi(); -} - -//------------------------ -// light_guide_digi -//------------------------ -void CalorimeterScFiDigi::light_guide_digi( void ){ - - // find the hits that belong to the same group (for merging) - std::unordered_map> merge_map; - for (auto ahit : simhits) { - int64_t hid = (ahit->getCellID() & id_mask) | ref_mask; - auto it = merge_map.find(hid); - - if (it == merge_map.end()) { - merge_map[hid] = {ahit}; - } else { - it->second.push_back(ahit); - } - } - - // signal sum - for (auto &[id, hits] : merge_map) { - double edep = hits[0]->getEnergy(); - double time = hits[0]->getContributions(0).getTime(); - double max_edep = hits[0]->getEnergy(); - int64_t mid = hits[0]->getCellID(); - - // sum energy, take time from the most energetic hit - // TODO, implement a timing window to sum or split the hits group - for (size_t i = 1; i < hits.size(); ++i) { - int64_t ztmp = id_dec->get(hits[i]->getCellID(), z_idx); - edep += hits[i]->getEnergy(); - if (hits[i]->getEnergy() > max_edep) { - max_edep = hits[i]->getEnergy(); - mid = hits[i]->getCellID(); - for (const auto& c : hits[i]->getContributions()) { - if (c.getTime() <= time) { - time = c.getTime(); - } - } - } - } - - // safety check - const double eResRel = (edep > m_threshold) - ? m_normDist(generator) * eRes[0] / std::sqrt(edep) + - m_normDist(generator) * eRes[1] + - m_normDist(generator) * eRes[2] / edep - : 0; - // digitize - double ped = m_pedMeanADC + m_normDist(generator) * m_pedSigmaADC; - unsigned long long adc = std::llround(ped + edep * (m_corrMeanScale + eResRel) / m_dyRangeADC * m_capADC); - unsigned long long tdc = std::llround((time + m_normDist(generator) * tRes) * stepTDC); - - // use the cellid from the most energetic hit in this group - auto rawhit = new edm4hep::RawCalorimeterHit( - mid, - (adc > m_capADC ? m_capADC : adc), - tdc - ); - rawhits.push_back(rawhit); - } - m_log->debug("size before digi: {:d}, after: {:d}", simhits.size(), rawhits.size()); -} diff --git a/src/algorithms/calorimetry/CalorimeterScFiDigi.h b/src/algorithms/calorimetry/CalorimeterScFiDigi.h deleted file mode 100644 index b952e66dc4..0000000000 --- a/src/algorithms/calorimetry/CalorimeterScFiDigi.h +++ /dev/null @@ -1,91 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -// Copyright (C) 2022 Chao Peng, Wouter Deconinck, Sylvester Joosten, Barak Schmookler, David Lawrence - -// A general digitization for CalorimeterHit from simulation -// 1. Smear energy deposit with a/sqrt(E/GeV) + b + c/E or a/sqrt(E/GeV) (relative value) -// 2. Digitize the energy with dynamic ADC range and add pedestal (mean +- sigma) -// 3. Time conversion with smearing resolution (absolute value) -// 4. Signal is summed if the SumFields are provided -// -// Author: Chao Peng -// Date: 06/02/2021 - - -#pragma once - -#include - -#include - -#include -#include -#include - -class CalorimeterScFiDigi { - - // Insert any member variables here - -public: - CalorimeterScFiDigi() = default; - ~CalorimeterScFiDigi(){for( auto h : rawhits ) delete h;} // better to use smart pointer? - virtual void AlgorithmInit(std::shared_ptr& logger); - virtual void AlgorithmChangeRun() ; - virtual void AlgorithmProcess() ; - - //-------- Configuration Parameters ------------ - //instantiate new spdlog logger - std::shared_ptr m_log; - - // Name of input data type (collection) - std::string m_input_tag; - - // additional smearing resolutions - std::vector u_eRes; - double m_tRes; - - // single hit energy deposition threshold - double m_threshold=1.0*dd4hep::keV; // {this, "threshold", 1. * keV}; - - // digitization settings - unsigned int m_capADC; - double m_dyRangeADC; - unsigned int m_pedMeanADC; - double m_pedSigmaADC; - double m_resolutionTDC; - double m_corrMeanScale; - - // signal sums - std::vector u_fields; - std::vector u_refs; - std::string m_geoSvcName; - std::string m_readout; - std::string m_zsegment; - - // This may be used to declare the data members as JANA configuration parameters. - // This should compile OK even without JANA so long as you don't try using it. - // To use it, do something like the following: - // - // mycalohitdigi->SetJANAConfigParameters( japp, "BEMC"); - // - // The above will register config. parameters like: "BEMC:tag". - // The configuration parameter members of this class should be set to thier - // defaults *before* calling this. - //----------------------------------------------- - - // unitless counterparts of inputs - double dyRangeADC{0}, stepTDC{0}, tRes{0}, eRes[3] = {0., 0., 0.}; - //Rndm::Numbers m_normDist; - std::shared_ptr m_geoSvc; - dd4hep::BitFieldCoder* id_dec = nullptr; - uint64_t id_mask{0}, ref_mask{0}, z_idx{0}; - - // inputs/outputs - std::vector simhits; - std::vector rawhits; - -private: - std::default_random_engine generator; // TODO: need something more appropriate here - std::normal_distribution m_normDist; // defaults to mean=0, sigma=1 - - void light_guide_digi(); -}; diff --git a/src/detectors/BEMC/RawCalorimeterHit_factory_EcalBarrelScFiRawHits.h b/src/detectors/BEMC/RawCalorimeterHit_factory_EcalBarrelScFiRawHits.h index bec8344e3a..22bdee5ef9 100644 --- a/src/detectors/BEMC/RawCalorimeterHit_factory_EcalBarrelScFiRawHits.h +++ b/src/detectors/BEMC/RawCalorimeterHit_factory_EcalBarrelScFiRawHits.h @@ -9,7 +9,7 @@ #include #include #include -#include +#include #include #include #include @@ -18,7 +18,7 @@ -class RawCalorimeterHit_factory_EcalBarrelScFiRawHits : public JFactoryT, CalorimeterScFiDigi { +class RawCalorimeterHit_factory_EcalBarrelScFiRawHits : public JFactoryT, CalorimeterHitDigi { public: @@ -34,7 +34,7 @@ class RawCalorimeterHit_factory_EcalBarrelScFiRawHits : public JFactoryT Date: Thu, 16 Feb 2023 00:17:11 -0600 Subject: [PATCH 2/3] remove m_zsegment --- .../BEMC/RawCalorimeterHit_factory_EcalBarrelScFiRawHits.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/detectors/BEMC/RawCalorimeterHit_factory_EcalBarrelScFiRawHits.h b/src/detectors/BEMC/RawCalorimeterHit_factory_EcalBarrelScFiRawHits.h index 22bdee5ef9..604a93bf2b 100644 --- a/src/detectors/BEMC/RawCalorimeterHit_factory_EcalBarrelScFiRawHits.h +++ b/src/detectors/BEMC/RawCalorimeterHit_factory_EcalBarrelScFiRawHits.h @@ -48,7 +48,6 @@ class RawCalorimeterHit_factory_EcalBarrelScFiRawHits : public JFactoryTGetService(); // TODO: implement named geometry service? // This is another option for exposing the data members as JANA configuration parameters. @@ -66,7 +65,6 @@ class RawCalorimeterHit_factory_EcalBarrelScFiRawHits : public JFactoryTSetDefaultParameter("BEMC:EcalBarrelScFiRawHits:fieldRefNumbers", u_refs); app->SetDefaultParameter("BEMC:EcalBarrelScFiRawHits:geoServiceName", m_geoSvcName); app->SetDefaultParameter("BEMC:EcalBarrelScFiRawHits:readoutClass", m_readout); - app->SetDefaultParameter("BEMC:EcalBarrelScFiRawHits:zSegment", m_zsegment); // Call Init for generic algorithm AlgorithmInit(m_log); From 7c6179e28dd25a1fb80f1d849de145512c2e4922 Mon Sep 17 00:00:00 2001 From: Chao Peng Date: Thu, 16 Feb 2023 11:39:55 -0600 Subject: [PATCH 3/3] add some flexibility to take wrong field names --- .../calorimetry/CalorimeterHitDigi.cc | 19 +++++++++++++------ .../calorimetry/CalorimeterHitDigi.h | 3 ++- ...rimeterHit_factory_EcalBarrelScFiRawHits.h | 4 ++-- 3 files changed, 17 insertions(+), 9 deletions(-) diff --git a/src/algorithms/calorimetry/CalorimeterHitDigi.cc b/src/algorithms/calorimetry/CalorimeterHitDigi.cc index a8fb5c7ff7..18739dc080 100644 --- a/src/algorithms/calorimetry/CalorimeterHitDigi.cc +++ b/src/algorithms/calorimetry/CalorimeterHitDigi.cc @@ -64,9 +64,9 @@ void CalorimeterHitDigi::AlgorithmInit(std::shared_ptr& logger) tRes = m_tRes / dd4hep::ns; stepTDC = dd4hep::ns / m_resolutionTDC; - // need signal sum + // all these are for signal sum at digitization level + merge_hits = false; if (!u_fields.empty()) { - // sanity checks if (!m_geoSvc) { m_log->error("Unable to locate Geometry Service."); @@ -91,13 +91,17 @@ void CalorimeterHitDigi::AlgorithmInit(std::shared_ptr& logger) ref_mask = id_desc.encode(ref_fields); // debug() << fmt::format("Referece id mask for the fields {:#064b}", ref_mask) << endmsg; } catch (...) { - m_log->warn("Failed to load ID decoder for {}", m_readout); - japp->Quit(); + // a workaround to avoid breaking the whole analysis if a field is not in some configurations + // TODO: it should be a fatal error to not cause unexpected analysis results + m_log->warn("Failed to load ID decoder for {}, hits will not be merged.", m_readout); + // japp->Quit(); return; } id_mask = ~id_mask; //LOG_INFO(default_cout_logger) << fmt::format("ID mask in {:s}: {:#064b}", m_readout, id_mask) << LOG_END; m_log->info("ID mask in {:s}: {:#064b}", m_readout, id_mask); + // all checks passed + merge_hits = true; } } @@ -122,7 +126,7 @@ void CalorimeterHitDigi::AlgorithmProcess() { for( auto h : rawhits ) delete h; rawhits.clear(); - if (!u_fields.empty()) { + if (merge_hits) { signal_sum_digi(); } else { single_hits_digi(); @@ -191,15 +195,18 @@ void CalorimeterHitDigi::signal_sum_digi( void ){ } // signal sum + // NOTE: we take the cellID of the most energetic hit in this group so it is a real cellID from an MC hit for (auto &[id, hits] : merge_map) { double edep = hits[0]->getEnergy(); double time = hits[0]->getContributions(0).getTime(); double max_edep = hits[0]->getEnergy(); + auto mid = hits[0]->getCellID(); // sum energy, take time from the most energetic hit for (size_t i = 1; i < hits.size(); ++i) { edep += hits[i]->getEnergy(); if (hits[i]->getEnergy() > max_edep) { max_edep = hits[i]->getEnergy(); + mid = hits[i]->getCellID(); for (const auto& c : hits[i]->getContributions()) { if (c.getTime() <= time) { time = c.getTime(); @@ -225,7 +232,7 @@ void CalorimeterHitDigi::signal_sum_digi( void ){ unsigned long long tdc = std::llround((time + m_normDist(generator) * tRes) * stepTDC); auto rawhit = new edm4hep::RawCalorimeterHit( - id, + mid, (adc > m_capADC ? m_capADC : adc), tdc ); diff --git a/src/algorithms/calorimetry/CalorimeterHitDigi.h b/src/algorithms/calorimetry/CalorimeterHitDigi.h index 0525c4c6b6..19555e7d7f 100644 --- a/src/algorithms/calorimetry/CalorimeterHitDigi.h +++ b/src/algorithms/calorimetry/CalorimeterHitDigi.h @@ -73,7 +73,8 @@ class CalorimeterHitDigi { // unitless counterparts of inputs double dyRangeADC{0}, stepTDC{0}, tRes{0}, eRes[3] = {0., 0., 0.}; - //Rndm::Numbers m_normDist; + // variables for merging at digitization step + bool merge_hits = false; std::shared_ptr m_geoSvc; uint64_t id_mask{0}, ref_mask{0}; diff --git a/src/detectors/BEMC/RawCalorimeterHit_factory_EcalBarrelScFiRawHits.h b/src/detectors/BEMC/RawCalorimeterHit_factory_EcalBarrelScFiRawHits.h index 604a93bf2b..89753ad7da 100644 --- a/src/detectors/BEMC/RawCalorimeterHit_factory_EcalBarrelScFiRawHits.h +++ b/src/detectors/BEMC/RawCalorimeterHit_factory_EcalBarrelScFiRawHits.h @@ -46,8 +46,8 @@ class RawCalorimeterHit_factory_EcalBarrelScFiRawHits : public JFactoryTGetService(); // TODO: implement named geometry service? // This is another option for exposing the data members as JANA configuration parameters.