diff --git a/JugAlgo/src/components/AlgoServiceSvc.cpp b/JugAlgo/src/components/AlgoServiceSvc.cpp index 4af79335ad335f3c5018b87e2576becb6d6e3a15..0e9c2bf929d7c6e3d015a092b263b1d2c21972fe 100644 --- a/JugAlgo/src/components/AlgoServiceSvc.cpp +++ b/JugAlgo/src/components/AlgoServiceSvc.cpp @@ -9,6 +9,7 @@ #include #include +#include #include // too many services? :P @@ -83,6 +84,17 @@ StatusCode AlgoServiceSvc::initialize() { info() << "Setting up algorithms::GeoSvc" << endmsg; auto* geo = static_cast(svc); geo->init(m_geoSvc->detector()); + } else if (name == algorithms::RandomSvc::kName) { + // Setup random service + //m_rndmGenSvc = service("RndmGenSvc"); + //if (!m_rndmGenSvc) { + // error() << "Unable to locate Random Service. " + // << "Make sure you have RndmGenSvc in the right order in the configuration." << endmsg; + // return StatusCode::FAILURE; + //} + info() << "Setting up algorithms::RandomSvc" << endmsg; + auto* random = static_cast(svc); + random->init(); // FIXME: algorithms::core::random service is not using framework } else { fatal() << "Unknown service encountered, please implement the necessary framework hooks" << endmsg; diff --git a/JugDigi/CMakeLists.txt b/JugDigi/CMakeLists.txt index 41d082b7ac5d9cb04e803f478fbccf499ce3cdf4..a177c3973eb3c5178422af7b4dd6349eaa8bb49e 100644 --- a/JugDigi/CMakeLists.txt +++ b/JugDigi/CMakeLists.txt @@ -11,7 +11,8 @@ gaudi_add_module(JugDigiPlugins ${JugDigiPlugins_sources} LINK Gaudi::GaudiAlgLib Gaudi::GaudiKernel - JugBase + JugBase JugAlgo + algocore algocalorimetry ROOT::Core ROOT::RIO ROOT::Tree EDM4HEP::edm4hep EDM4EIC::edm4eic diff --git a/JugDigi/src/components/CalorimeterHitDigi.cpp b/JugDigi/src/components/CalorimeterHitDigi.cpp index a108d7656b022ab4fad08521a3c12895a00a2e29..c93bbf2b4aab9ed6c8f8a349b747c03fe10037b9 100644 --- a/JugDigi/src/components/CalorimeterHitDigi.cpp +++ b/JugDigi/src/components/CalorimeterHitDigi.cpp @@ -1,264 +1,65 @@ // SPDX-License-Identifier: LGPL-3.0-or-later // Copyright (C) 2022 Chao Peng, Wouter Deconinck, Sylvester Joosten -// 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 +#include +#include -#include -#include -#include - -#include "GaudiAlg/GaudiTool.h" -#include "GaudiAlg/Transformer.h" -#include "GaudiKernel/PhysicalConstants.h" #include "Gaudi/Property.h" -#include "GaudiKernel/RndmGenerators.h" - -#include "DDRec/CellIDPositionConverter.h" -#include "DDSegmentation/BitFieldCoder.h" - -#include "JugBase/IGeoSvc.h" -#include "JugBase/DataHandle.h" - -#include "fmt/format.h" -#include "fmt/ranges.h" - -// Event Model related classes -#include "edm4hep/SimCalorimeterHitCollection.h" -#include "edm4eic/RawCalorimeterHitCollection.h" -#include "edm4eic/RawCalorimeterHitData.h" - - -using namespace Gaudi::Units; namespace Jug::Digi { - /** Generic calorimeter hit digitiziation. - * - * \ingroup digi - * \ingroup calorimetry - */ - class CalorimeterHitDigi : public GaudiAlgorithm { - public: - // additional smearing resolutions - Gaudi::Property> u_eRes{this, "energyResolutions", {}}; // a/sqrt(E/GeV) + b + c/(E/GeV) - Gaudi::Property m_tRes{this, "timeResolution", 0.0 * ns}; - // single hit energy deposition threshold - Gaudi::Property m_threshold{this, "threshold", 1. * keV}; - - // digitization settings - Gaudi::Property m_capADC{this, "capacityADC", 8096}; - Gaudi::Property m_dyRangeADC{this, "dynamicRangeADC", 100 * MeV}; - Gaudi::Property m_pedMeanADC{this, "pedestalMean", 400}; - Gaudi::Property m_pedSigmaADC{this, "pedestalSigma", 3.2}; - Gaudi::Property m_resolutionTDC{this, "resolutionTDC", 10 * ps}; - - Gaudi::Property m_corrMeanScale{this, "scaleResponse", 1.0}; - // These are better variable names for the "energyResolutions" array which is a bit - // magic @FIXME - //Gaudi::Property m_corrSigmaCoeffE{this, "responseCorrectionSigmaCoeffE", 0.0}; - //Gaudi::Property m_corrSigmaCoeffSqrtE{this, "responseCorrectionSigmaCoeffSqrtE", 0.0}; - - // signal sums - // @TODO: implement signal sums with timing - // field names to generate id mask, the hits will be grouped by masking the field - Gaudi::Property> u_fields{this, "signalSumFields", {}}; - // ref field ids are used for the merged hits, 0 is used if nothing provided - Gaudi::Property> u_refs{this, "fieldRefNumbers", {}}; - Gaudi::Property m_geoSvcName{this, "geoServiceName", "GeoSvc"}; - Gaudi::Property m_readout{this, "readoutClass", ""}; - - // unitless counterparts of inputs - double dyRangeADC{0}, stepTDC{0}, tRes{0}, eRes[3] = {0., 0., 0.}; - Rndm::Numbers m_normDist; - SmartIF m_geoSvc; - uint64_t id_mask{0}, ref_mask{0}; - - DataHandle m_inputHitCollection{ - "inputHitCollection", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputHitCollection{ - "outputHitCollection", Gaudi::DataHandle::Writer, this}; - - // ill-formed: using GaudiAlgorithm::GaudiAlgorithm; - CalorimeterHitDigi(const std::string& name, ISvcLocator* svcLoc) - : GaudiAlgorithm(name, svcLoc) { - declareProperty("inputHitCollection", m_inputHitCollection, ""); - declareProperty("outputHitCollection", m_outputHitCollection, ""); - } - - StatusCode initialize() override - { - if (GaudiAlgorithm::initialize().isFailure()) { - return StatusCode::FAILURE; - } - // random number generator from service - auto randSvc = svc("RndmGenSvc", true); - auto sc = m_normDist.initialize(randSvc, Rndm::Gauss(0.0, 1.0)); - if (!sc.isSuccess()) { - return StatusCode::FAILURE; - } - // set energy resolution numbers - 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) - dyRangeADC = m_dyRangeADC.value() / GeV; - tRes = m_tRes.value() / ns; - stepTDC = ns / m_resolutionTDC.value(); - - // need signal sum - if (!u_fields.value().empty()) { - m_geoSvc = service(m_geoSvcName); - // sanity checks - if (!m_geoSvc) { - error() << "Unable to locate Geometry Service. " - << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." - << endmsg; - return StatusCode::FAILURE; - } - if (m_readout.value().empty()) { - error() << "readoutClass is not provided, it is needed to know the fields in readout ids" - << endmsg; - return StatusCode::FAILURE; - } - - // get decoders - try { - auto id_desc = m_geoSvc->detector()->readout(m_readout).idSpec(); - 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); - } - ref_mask = id_desc.encode(ref_fields); - // debug() << fmt::format("Referece id mask for the fields {:#064b}", ref_mask) << endmsg; - } catch (...) { - error() << "Failed to load ID decoder for " << m_readout << endmsg; - return StatusCode::FAILURE; - } - id_mask = ~id_mask; - info() << fmt::format("ID mask in {:s}: {:#064b}", m_readout, id_mask) << endmsg; - return StatusCode::SUCCESS; - } - - return StatusCode::SUCCESS; - } - - StatusCode execute() override - { - if (!u_fields.value().empty()) { - signal_sum_digi(); - } else { - single_hits_digi(); - } - return StatusCode::SUCCESS; - } - - private: - void single_hits_digi() { - // input collections - const auto* const simhits = m_inputHitCollection.get(); - // Create output collections - auto* rawhits = m_outputHitCollection.createAndPut(); - for (const auto& ahit : *simhits) { - // Note: juggler internal unit of energy is GeV - const double eDep = ahit.getEnergy(); - - // apply additional calorimeter noise to corrected energy deposit - const double eResRel = (eDep > m_threshold) - ? m_normDist() * std::sqrt( - std::pow(eRes[0] / std::sqrt(eDep), 2) + - std::pow(eRes[1], 2) + - std::pow(eRes[2] / (eDep), 2) - ) - : 0; - - const double ped = m_pedMeanADC + m_normDist() * m_pedSigmaADC; - const long long adc = std::llround(ped + eDep * (m_corrMeanScale + eResRel) / dyRangeADC * m_capADC); - - double time = std::numeric_limits::max(); - for (const auto& c : ahit.getContributions()) { - if (c.getTime() <= time) { - time = c.getTime(); - } - } - const long long tdc = std::llround((time + m_normDist() * tRes) * stepTDC); - - edm4eic::RawCalorimeterHit rawhit( - ahit.getCellID(), - (adc > m_capADC.value() ? m_capADC.value() : adc), - tdc - ); - rawhits->push_back(rawhit); - } - } - - void signal_sum_digi() { - const auto* const simhits = m_inputHitCollection.get(); - auto* rawhits = m_outputHitCollection.createAndPut(); - - // find the hits that belong to the same group (for merging) - std::unordered_map> merge_map; - for (const 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(); - // 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(); - for (const auto& c : hits[i].getContributions()) { - if (c.getTime() <= time) { - time = c.getTime(); - } - } - } - } - - // safety check - const double eResRel = (edep > m_threshold) - ? m_normDist() * eRes[0] / std::sqrt(edep) + - m_normDist() * eRes[1] + - m_normDist() * eRes[2] / edep - : 0; - - double ped = m_pedMeanADC + m_normDist() * m_pedSigmaADC; - unsigned long long adc = std::llround(ped + edep * (1. + eResRel) / dyRangeADC * m_capADC); - unsigned long long tdc = std::llround((time + m_normDist() * tRes) * stepTDC); - - edm4eic::RawCalorimeterHit rawhit( - id, - (adc > m_capADC.value() ? m_capADC.value() : adc), - tdc - ); - rawhits->push_back(rawhit); - } - } - }; - // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) - DECLARE_COMPONENT(CalorimeterHitDigi) +namespace { + using AlgoBase = Jug::Algo::Algorithm; +} + +class CalorimeterHitDigi : public AlgoBase { + +public: + CalorimeterHitDigi(const std::string& name, ISvcLocator* svcLoc) : AlgoBase(name, svcLoc) {} + + virtual StatusCode configure() { + setAlgoProp("energyResolutions", u_eRes.value()); + setAlgoProp("timeResolution", m_tRes.value()); + setAlgoProp("threshold", m_threshold.value()); + setAlgoProp("capacityADC", m_capADC.value()); + setAlgoProp("dynamicRangeADC", m_dyRangeADC.value()); + setAlgoProp("pedMeanADC", m_pedMeanADC.value()); + setAlgoProp("pedSigmaADC", m_pedSigmaADC.value()); + setAlgoProp("resolutionTDC", m_resolutionTDC.value()); + setAlgoProp("scaleResponse", m_corrMeanScale.value()); + setAlgoProp("fieldRefNumbers", u_refs.value()); + setAlgoProp("geoServiceName", m_geoSvcName.value()); + setAlgoProp("readoutClass", m_readout.value()); + return StatusCode::SUCCESS; + } + +private: + // additional smearing resolutions + Gaudi::Property> u_eRes{this, "energyResolutions", {}}; // a/sqrt(E/GeV) + b + c/(E/GeV) + Gaudi::Property m_tRes{this, "timeResolution", 0.0 * dd4hep::ns}; + // single hit energy deposition threshold + Gaudi::Property m_threshold{this, "threshold", 1. * dd4hep::keV}; + + // digitization settings + Gaudi::Property m_capADC{this, "capacityADC", 8096}; + Gaudi::Property m_dyRangeADC{this, "dynamicRangeADC", 100 * dd4hep::MeV}; + Gaudi::Property m_pedMeanADC{this, "pedestalMean", 400}; + Gaudi::Property m_pedSigmaADC{this, "pedestalSigma", 3.2}; + Gaudi::Property m_resolutionTDC{this, "resolutionTDC", 0.010 * dd4hep::ns}; + Gaudi::Property m_corrMeanScale{this, "scaleResponse", 1.0}; + + // signal sums + // @TODO: implement signal sums with timing + // field names to generate id mask, the hits will be grouped by masking the field + Gaudi::Property> u_fields{this, "signalSumFields", {}}; + // ref field ids are used for the merged hits, 0 is used if nothing provided + Gaudi::Property> u_refs{this, "fieldRefNumbers", {}}; + Gaudi::Property m_geoSvcName{this, "geoServiceName", "GeoSvc"}; + Gaudi::Property m_readout{this, "readoutClass", ""}; + +}; + +// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) +DECLARE_COMPONENT(CalorimeterHitDigi) } // namespace Jug::Digi diff --git a/external/algorithms/calorimetry/include/algorithms/calorimetry/CalorimeterHitDigi.cpp b/external/algorithms/calorimetry/include/algorithms/calorimetry/CalorimeterHitDigi.cpp deleted file mode 100644 index a108d7656b022ab4fad08521a3c12895a00a2e29..0000000000000000000000000000000000000000 --- a/external/algorithms/calorimetry/include/algorithms/calorimetry/CalorimeterHitDigi.cpp +++ /dev/null @@ -1,264 +0,0 @@ -// SPDX-License-Identifier: LGPL-3.0-or-later -// Copyright (C) 2022 Chao Peng, Wouter Deconinck, Sylvester Joosten - -// 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 - -#include -#include -#include - -#include "GaudiAlg/GaudiTool.h" -#include "GaudiAlg/Transformer.h" -#include "GaudiKernel/PhysicalConstants.h" -#include "Gaudi/Property.h" -#include "GaudiKernel/RndmGenerators.h" - -#include "DDRec/CellIDPositionConverter.h" -#include "DDSegmentation/BitFieldCoder.h" - -#include "JugBase/IGeoSvc.h" -#include "JugBase/DataHandle.h" - -#include "fmt/format.h" -#include "fmt/ranges.h" - -// Event Model related classes -#include "edm4hep/SimCalorimeterHitCollection.h" -#include "edm4eic/RawCalorimeterHitCollection.h" -#include "edm4eic/RawCalorimeterHitData.h" - - -using namespace Gaudi::Units; - -namespace Jug::Digi { - - /** Generic calorimeter hit digitiziation. - * - * \ingroup digi - * \ingroup calorimetry - */ - class CalorimeterHitDigi : public GaudiAlgorithm { - public: - // additional smearing resolutions - Gaudi::Property> u_eRes{this, "energyResolutions", {}}; // a/sqrt(E/GeV) + b + c/(E/GeV) - Gaudi::Property m_tRes{this, "timeResolution", 0.0 * ns}; - // single hit energy deposition threshold - Gaudi::Property m_threshold{this, "threshold", 1. * keV}; - - // digitization settings - Gaudi::Property m_capADC{this, "capacityADC", 8096}; - Gaudi::Property m_dyRangeADC{this, "dynamicRangeADC", 100 * MeV}; - Gaudi::Property m_pedMeanADC{this, "pedestalMean", 400}; - Gaudi::Property m_pedSigmaADC{this, "pedestalSigma", 3.2}; - Gaudi::Property m_resolutionTDC{this, "resolutionTDC", 10 * ps}; - - Gaudi::Property m_corrMeanScale{this, "scaleResponse", 1.0}; - // These are better variable names for the "energyResolutions" array which is a bit - // magic @FIXME - //Gaudi::Property m_corrSigmaCoeffE{this, "responseCorrectionSigmaCoeffE", 0.0}; - //Gaudi::Property m_corrSigmaCoeffSqrtE{this, "responseCorrectionSigmaCoeffSqrtE", 0.0}; - - // signal sums - // @TODO: implement signal sums with timing - // field names to generate id mask, the hits will be grouped by masking the field - Gaudi::Property> u_fields{this, "signalSumFields", {}}; - // ref field ids are used for the merged hits, 0 is used if nothing provided - Gaudi::Property> u_refs{this, "fieldRefNumbers", {}}; - Gaudi::Property m_geoSvcName{this, "geoServiceName", "GeoSvc"}; - Gaudi::Property m_readout{this, "readoutClass", ""}; - - // unitless counterparts of inputs - double dyRangeADC{0}, stepTDC{0}, tRes{0}, eRes[3] = {0., 0., 0.}; - Rndm::Numbers m_normDist; - SmartIF m_geoSvc; - uint64_t id_mask{0}, ref_mask{0}; - - DataHandle m_inputHitCollection{ - "inputHitCollection", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputHitCollection{ - "outputHitCollection", Gaudi::DataHandle::Writer, this}; - - // ill-formed: using GaudiAlgorithm::GaudiAlgorithm; - CalorimeterHitDigi(const std::string& name, ISvcLocator* svcLoc) - : GaudiAlgorithm(name, svcLoc) { - declareProperty("inputHitCollection", m_inputHitCollection, ""); - declareProperty("outputHitCollection", m_outputHitCollection, ""); - } - - StatusCode initialize() override - { - if (GaudiAlgorithm::initialize().isFailure()) { - return StatusCode::FAILURE; - } - // random number generator from service - auto randSvc = svc("RndmGenSvc", true); - auto sc = m_normDist.initialize(randSvc, Rndm::Gauss(0.0, 1.0)); - if (!sc.isSuccess()) { - return StatusCode::FAILURE; - } - // set energy resolution numbers - 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) - dyRangeADC = m_dyRangeADC.value() / GeV; - tRes = m_tRes.value() / ns; - stepTDC = ns / m_resolutionTDC.value(); - - // need signal sum - if (!u_fields.value().empty()) { - m_geoSvc = service(m_geoSvcName); - // sanity checks - if (!m_geoSvc) { - error() << "Unable to locate Geometry Service. " - << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." - << endmsg; - return StatusCode::FAILURE; - } - if (m_readout.value().empty()) { - error() << "readoutClass is not provided, it is needed to know the fields in readout ids" - << endmsg; - return StatusCode::FAILURE; - } - - // get decoders - try { - auto id_desc = m_geoSvc->detector()->readout(m_readout).idSpec(); - 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); - } - ref_mask = id_desc.encode(ref_fields); - // debug() << fmt::format("Referece id mask for the fields {:#064b}", ref_mask) << endmsg; - } catch (...) { - error() << "Failed to load ID decoder for " << m_readout << endmsg; - return StatusCode::FAILURE; - } - id_mask = ~id_mask; - info() << fmt::format("ID mask in {:s}: {:#064b}", m_readout, id_mask) << endmsg; - return StatusCode::SUCCESS; - } - - return StatusCode::SUCCESS; - } - - StatusCode execute() override - { - if (!u_fields.value().empty()) { - signal_sum_digi(); - } else { - single_hits_digi(); - } - return StatusCode::SUCCESS; - } - - private: - void single_hits_digi() { - // input collections - const auto* const simhits = m_inputHitCollection.get(); - // Create output collections - auto* rawhits = m_outputHitCollection.createAndPut(); - for (const auto& ahit : *simhits) { - // Note: juggler internal unit of energy is GeV - const double eDep = ahit.getEnergy(); - - // apply additional calorimeter noise to corrected energy deposit - const double eResRel = (eDep > m_threshold) - ? m_normDist() * std::sqrt( - std::pow(eRes[0] / std::sqrt(eDep), 2) + - std::pow(eRes[1], 2) + - std::pow(eRes[2] / (eDep), 2) - ) - : 0; - - const double ped = m_pedMeanADC + m_normDist() * m_pedSigmaADC; - const long long adc = std::llround(ped + eDep * (m_corrMeanScale + eResRel) / dyRangeADC * m_capADC); - - double time = std::numeric_limits::max(); - for (const auto& c : ahit.getContributions()) { - if (c.getTime() <= time) { - time = c.getTime(); - } - } - const long long tdc = std::llround((time + m_normDist() * tRes) * stepTDC); - - edm4eic::RawCalorimeterHit rawhit( - ahit.getCellID(), - (adc > m_capADC.value() ? m_capADC.value() : adc), - tdc - ); - rawhits->push_back(rawhit); - } - } - - void signal_sum_digi() { - const auto* const simhits = m_inputHitCollection.get(); - auto* rawhits = m_outputHitCollection.createAndPut(); - - // find the hits that belong to the same group (for merging) - std::unordered_map> merge_map; - for (const 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(); - // 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(); - for (const auto& c : hits[i].getContributions()) { - if (c.getTime() <= time) { - time = c.getTime(); - } - } - } - } - - // safety check - const double eResRel = (edep > m_threshold) - ? m_normDist() * eRes[0] / std::sqrt(edep) + - m_normDist() * eRes[1] + - m_normDist() * eRes[2] / edep - : 0; - - double ped = m_pedMeanADC + m_normDist() * m_pedSigmaADC; - unsigned long long adc = std::llround(ped + edep * (1. + eResRel) / dyRangeADC * m_capADC); - unsigned long long tdc = std::llround((time + m_normDist() * tRes) * stepTDC); - - edm4eic::RawCalorimeterHit rawhit( - id, - (adc > m_capADC.value() ? m_capADC.value() : adc), - tdc - ); - rawhits->push_back(rawhit); - } - } - }; - // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) - DECLARE_COMPONENT(CalorimeterHitDigi) - -} // namespace Jug::Digi diff --git a/external/algorithms/calorimetry/include/algorithms/calorimetry/CalorimeterHitDigi.h b/external/algorithms/calorimetry/include/algorithms/calorimetry/CalorimeterHitDigi.h new file mode 100644 index 0000000000000000000000000000000000000000..c86431f9f78693371bdf1df225d8f077ab391b59 --- /dev/null +++ b/external/algorithms/calorimetry/include/algorithms/calorimetry/CalorimeterHitDigi.h @@ -0,0 +1,90 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Chao Peng, Wouter Deconinck, Sylvester Joosten + +// 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 + +#include +#include +#include +#include + +#include "DDRec/CellIDPositionConverter.h" +#include "DDSegmentation/BitFieldCoder.h" + +#include "fmt/format.h" +#include "fmt/ranges.h" + +// Event Model related classes +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "edm4eic/RawCalorimeterHitCollection.h" + +namespace algorithms::calorimetry { + + using CalorimeterHitDigiAlgorithm = Algorithm< + Input, + Output>; + + /** Generic calorimeter hit digitiziation. + * + * \ingroup digi + * \ingroup calorimetry + */ + class CalorimeterHitDigi : public CalorimeterHitDigiAlgorithm { + public: + using Input = CalorimeterHitDigiAlgorithm::Input; + using Output = CalorimeterHitDigiAlgorithm::Output; + + CalorimeterHitDigi(std::string_view name) + : CalorimeterHitDigiAlgorithm{name, + {"inputHitCollection"}, + {"outputHitCollection"}} {}; + + void init(); + void process(const Input&, const Output&); + + private: + void single_hits_digi(const Input&, const Output&); + void signal_sum_digi(const Input&, const Output&); + + // additional smearing resolutions + Property> u_eRes{this, "energyResolutions", {}}; // a/sqrt(E/GeV) + b + c/(E/GeV) + Property m_tRes{this, "timeResolution", 0.0 * dd4hep::ns}; + // single hit energy deposition threshold + Property m_threshold{this, "threshold", 1. * dd4hep::keV}; + + // digitization settings + Property m_capADC{this, "capacityADC", 8096}; + Property m_dyRangeADC{this, "dynamicRangeADC", 100 * dd4hep::MeV}; + Property m_pedMeanADC{this, "pedestalMean", 400}; + Property m_pedSigmaADC{this, "pedestalSigma", 3.2}; + Property m_resolutionTDC{this, "resolutionTDC", 0.010 * dd4hep::ns}; + + Property m_corrMeanScale{this, "scaleResponse", 1.0}; + // These are better variable names for the "energyResolutions" array which is a bit + // magic @FIXME + //Property m_corrSigmaCoeffE{this, "responseCorrectionSigmaCoeffE", 0.0}; + //Property m_corrSigmaCoeffSqrtE{this, "responseCorrectionSigmaCoeffSqrtE", 0.0}; + + // signal sums + // @TODO: implement signal sums with timing + // field names to generate id mask, the hits will be grouped by masking the field + Property> u_fields{this, "signalSumFields", {}}; + // ref field ids are used for the merged hits, 0 is used if nothing provided + Property> u_refs{this, "fieldRefNumbers", {}}; + Property m_readout{this, "readoutClass", ""}; + + double eRes[3] = {0., 0., 0.}; + uint64_t id_mask{0}, ref_mask{0}; + + const GeoSvc& m_geoSvc = GeoSvc::instance(); + const RandomSvc& m_randomSvc = RandomSvc::instance(); + }; + +} // namespace algoriths::calorimetry diff --git a/external/algorithms/calorimetry/src/CalorimeterHitDigi.cpp b/external/algorithms/calorimetry/src/CalorimeterHitDigi.cpp index a108d7656b022ab4fad08521a3c12895a00a2e29..deae3f9590260382acac888bc9c1aec3826fcaf0 100644 --- a/external/algorithms/calorimetry/src/CalorimeterHitDigi.cpp +++ b/external/algorithms/calorimetry/src/CalorimeterHitDigi.cpp @@ -10,22 +10,12 @@ // Author: Chao Peng // Date: 06/02/2021 +#include + #include #include #include -#include "GaudiAlg/GaudiTool.h" -#include "GaudiAlg/Transformer.h" -#include "GaudiKernel/PhysicalConstants.h" -#include "Gaudi/Property.h" -#include "GaudiKernel/RndmGenerators.h" - -#include "DDRec/CellIDPositionConverter.h" -#include "DDSegmentation/BitFieldCoder.h" - -#include "JugBase/IGeoSvc.h" -#include "JugBase/DataHandle.h" - #include "fmt/format.h" #include "fmt/ranges.h" @@ -34,231 +24,152 @@ #include "edm4eic/RawCalorimeterHitCollection.h" #include "edm4eic/RawCalorimeterHitData.h" - -using namespace Gaudi::Units; - -namespace Jug::Digi { - - /** Generic calorimeter hit digitiziation. - * - * \ingroup digi - * \ingroup calorimetry - */ - class CalorimeterHitDigi : public GaudiAlgorithm { - public: - // additional smearing resolutions - Gaudi::Property> u_eRes{this, "energyResolutions", {}}; // a/sqrt(E/GeV) + b + c/(E/GeV) - Gaudi::Property m_tRes{this, "timeResolution", 0.0 * ns}; - // single hit energy deposition threshold - Gaudi::Property m_threshold{this, "threshold", 1. * keV}; - - // digitization settings - Gaudi::Property m_capADC{this, "capacityADC", 8096}; - Gaudi::Property m_dyRangeADC{this, "dynamicRangeADC", 100 * MeV}; - Gaudi::Property m_pedMeanADC{this, "pedestalMean", 400}; - Gaudi::Property m_pedSigmaADC{this, "pedestalSigma", 3.2}; - Gaudi::Property m_resolutionTDC{this, "resolutionTDC", 10 * ps}; - - Gaudi::Property m_corrMeanScale{this, "scaleResponse", 1.0}; - // These are better variable names for the "energyResolutions" array which is a bit - // magic @FIXME - //Gaudi::Property m_corrSigmaCoeffE{this, "responseCorrectionSigmaCoeffE", 0.0}; - //Gaudi::Property m_corrSigmaCoeffSqrtE{this, "responseCorrectionSigmaCoeffSqrtE", 0.0}; - - // signal sums - // @TODO: implement signal sums with timing - // field names to generate id mask, the hits will be grouped by masking the field - Gaudi::Property> u_fields{this, "signalSumFields", {}}; - // ref field ids are used for the merged hits, 0 is used if nothing provided - Gaudi::Property> u_refs{this, "fieldRefNumbers", {}}; - Gaudi::Property m_geoSvcName{this, "geoServiceName", "GeoSvc"}; - Gaudi::Property m_readout{this, "readoutClass", ""}; - - // unitless counterparts of inputs - double dyRangeADC{0}, stepTDC{0}, tRes{0}, eRes[3] = {0., 0., 0.}; - Rndm::Numbers m_normDist; - SmartIF m_geoSvc; - uint64_t id_mask{0}, ref_mask{0}; - - DataHandle m_inputHitCollection{ - "inputHitCollection", Gaudi::DataHandle::Reader, this}; - DataHandle m_outputHitCollection{ - "outputHitCollection", Gaudi::DataHandle::Writer, this}; - - // ill-formed: using GaudiAlgorithm::GaudiAlgorithm; - CalorimeterHitDigi(const std::string& name, ISvcLocator* svcLoc) - : GaudiAlgorithm(name, svcLoc) { - declareProperty("inputHitCollection", m_inputHitCollection, ""); - declareProperty("outputHitCollection", m_outputHitCollection, ""); +namespace algorithms::calorimetry { + +void CalorimeterHitDigi::init() { + // set energy resolution numbers + for (size_t i = 0; i < u_eRes.size() && i < 3; ++i) { + eRes[i] = u_eRes[i]; + } + + // need signal sum + if (!u_fields.value().empty()) { + // sanity checks + if (m_readout.value().empty()) { + error() << "readoutClass is not provided, it is needed to know the fields in readout ids" + << endmsg; + return; } - StatusCode initialize() override - { - if (GaudiAlgorithm::initialize().isFailure()) { - return StatusCode::FAILURE; - } - // random number generator from service - auto randSvc = svc("RndmGenSvc", true); - auto sc = m_normDist.initialize(randSvc, Rndm::Gauss(0.0, 1.0)); - if (!sc.isSuccess()) { - return StatusCode::FAILURE; - } - // set energy resolution numbers - 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) - dyRangeADC = m_dyRangeADC.value() / GeV; - tRes = m_tRes.value() / ns; - stepTDC = ns / m_resolutionTDC.value(); - - // need signal sum - if (!u_fields.value().empty()) { - m_geoSvc = service(m_geoSvcName); - // sanity checks - if (!m_geoSvc) { - error() << "Unable to locate Geometry Service. " - << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." - << endmsg; - return StatusCode::FAILURE; - } - if (m_readout.value().empty()) { - error() << "readoutClass is not provided, it is needed to know the fields in readout ids" - << endmsg; - return StatusCode::FAILURE; - } - - // get decoders - try { - auto id_desc = m_geoSvc->detector()->readout(m_readout).idSpec(); - 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); - } - ref_mask = id_desc.encode(ref_fields); - // debug() << fmt::format("Referece id mask for the fields {:#064b}", ref_mask) << endmsg; - } catch (...) { - error() << "Failed to load ID decoder for " << m_readout << endmsg; - return StatusCode::FAILURE; - } - id_mask = ~id_mask; - info() << fmt::format("ID mask in {:s}: {:#064b}", m_readout, id_mask) << endmsg; - return StatusCode::SUCCESS; + // get decoders + try { + auto id_desc = m_geoSvc.detector()->readout(m_readout).idSpec(); + 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); } - - return StatusCode::SUCCESS; + ref_mask = id_desc.encode(ref_fields); + // debug() << fmt::format("Referece id mask for the fields {:#064b}", ref_mask) << endmsg; + } catch (...) { + error() << "Failed to load ID decoder for " << m_readout << endmsg; + return; } - - StatusCode execute() override - { - if (!u_fields.value().empty()) { - signal_sum_digi(); - } else { - single_hits_digi(); + id_mask = ~id_mask; + info() << fmt::format("ID mask in {:s}: {:#064b}", m_readout, id_mask) << endmsg; + } +} + +void CalorimeterHitDigi::process( + const CalorimeterHitDigi::Input& input, + const CalorimeterHitDigi::Output& output +) { + if (!u_fields.value().empty()) { + signal_sum_digi(input, output); + } else { + single_hits_digi(input, output); + } +} + +void CalorimeterHitDigi::single_hits_digi( + const CalorimeterHitDigi::Input& input, + const CalorimeterHitDigi::Output& output +) { + const auto [simhits] = input; + auto [rawhits] = output; + + for (const auto& ahit : *simhits) { + // Note: juggler internal unit of energy is GeV + const double eDep = ahit.getEnergy(); + + // apply additional calorimeter noise to corrected energy deposit + const double eResRel = (eDep > m_threshold) + ? m_randomSvc.normal() * std::sqrt( + std::pow(eRes[0] / std::sqrt(eDep), 2) + + std::pow(eRes[1], 2) + + std::pow(eRes[2] / (eDep), 2) + ) + : 0; + + const double ped = m_pedMeanADC + m_randomSvc.normal() * m_pedSigmaADC; + const long long adc = std::llround(ped + eDep * (m_corrMeanScale + eResRel) / m_dyRangeADC * m_capADC); + + double time = std::numeric_limits::max(); + for (const auto& c : ahit.getContributions()) { + if (c.getTime() <= time) { + time = c.getTime(); } - return StatusCode::SUCCESS; } - - private: - void single_hits_digi() { - // input collections - const auto* const simhits = m_inputHitCollection.get(); - // Create output collections - auto* rawhits = m_outputHitCollection.createAndPut(); - for (const auto& ahit : *simhits) { - // Note: juggler internal unit of energy is GeV - const double eDep = ahit.getEnergy(); - - // apply additional calorimeter noise to corrected energy deposit - const double eResRel = (eDep > m_threshold) - ? m_normDist() * std::sqrt( - std::pow(eRes[0] / std::sqrt(eDep), 2) + - std::pow(eRes[1], 2) + - std::pow(eRes[2] / (eDep), 2) - ) - : 0; - - const double ped = m_pedMeanADC + m_normDist() * m_pedSigmaADC; - const long long adc = std::llround(ped + eDep * (m_corrMeanScale + eResRel) / dyRangeADC * m_capADC); - - double time = std::numeric_limits::max(); - for (const auto& c : ahit.getContributions()) { + const long long tdc = std::llround((time + m_randomSvc.normal() * m_tRes) / m_resolutionTDC); + + edm4eic::RawCalorimeterHit rawhit( + ahit.getCellID(), + (adc > m_capADC.value() ? m_capADC.value() : adc), + tdc + ); + rawhits->push_back(rawhit); + } +} + +void CalorimeterHitDigi::signal_sum_digi( + const CalorimeterHitDigi::Input& input, + const CalorimeterHitDigi::Output& output +) { + const auto [simhits] = input; + auto [rawhits] = output; + + // find the hits that belong to the same group (for merging) + std::unordered_map> merge_map; + for (const 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(); + // 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(); + for (const auto& c : hits[i].getContributions()) { if (c.getTime() <= time) { time = c.getTime(); } } - const long long tdc = std::llround((time + m_normDist() * tRes) * stepTDC); - - edm4eic::RawCalorimeterHit rawhit( - ahit.getCellID(), - (adc > m_capADC.value() ? m_capADC.value() : adc), - tdc - ); - rawhits->push_back(rawhit); - } - } - - void signal_sum_digi() { - const auto* const simhits = m_inputHitCollection.get(); - auto* rawhits = m_outputHitCollection.createAndPut(); - - // find the hits that belong to the same group (for merging) - std::unordered_map> merge_map; - for (const 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(); - // 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(); - for (const auto& c : hits[i].getContributions()) { - if (c.getTime() <= time) { - time = c.getTime(); - } - } - } - } - - // safety check - const double eResRel = (edep > m_threshold) - ? m_normDist() * eRes[0] / std::sqrt(edep) + - m_normDist() * eRes[1] + - m_normDist() * eRes[2] / edep - : 0; - - double ped = m_pedMeanADC + m_normDist() * m_pedSigmaADC; - unsigned long long adc = std::llround(ped + edep * (1. + eResRel) / dyRangeADC * m_capADC); - unsigned long long tdc = std::llround((time + m_normDist() * tRes) * stepTDC); - - edm4eic::RawCalorimeterHit rawhit( - id, - (adc > m_capADC.value() ? m_capADC.value() : adc), - tdc - ); - rawhits->push_back(rawhit); } } - }; - // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) - DECLARE_COMPONENT(CalorimeterHitDigi) -} // namespace Jug::Digi + // safety check + const double eResRel = (edep > m_threshold) + ? m_randomSvc.normal() * eRes[0] / std::sqrt(edep) + + m_randomSvc.normal() * eRes[1] + + m_randomSvc.normal() * eRes[2] / edep + : 0; + + double ped = m_pedMeanADC + m_randomSvc.normal() * m_pedSigmaADC; + unsigned long long adc = std::llround(ped + edep * (1. + eResRel) / m_dyRangeADC * m_capADC); + unsigned long long tdc = std::llround((time + m_randomSvc.normal() * m_tRes) / m_resolutionTDC); + + edm4eic::RawCalorimeterHit rawhit( + id, + (adc > m_capADC.value() ? m_capADC.value() : adc), + tdc + ); + rawhits->push_back(rawhit); + } +} + +} // namespace algorithms::calorimetry diff --git a/external/algorithms/core/include/algorithms/random.h b/external/algorithms/core/include/algorithms/random.h new file mode 100644 index 0000000000000000000000000000000000000000..09ed60417b3f51211ce84633c7a348ec3df02d48 --- /dev/null +++ b/external/algorithms/core/include/algorithms/random.h @@ -0,0 +1,27 @@ +#pragma once + +#include +#include +#include + +#include +#include + +namespace algorithms { + +class RandomSvc : public LoggedService { +public: + // Initialize the random service + void init(); + + double uniform(double a = 0.0, double b = 1.0) const; + double normal(double mu = 0.0, double sigma = 1.0) const; + +private: + // Standard mersenne_twister_engine seeded with rd() + mutable std::mt19937 rng{std::random_device{}()}; + + ALGORITHMS_DEFINE_LOGGED_SERVICE(RandomSvc) +}; + +} // namespace algorithms diff --git a/external/algorithms/core/src/random.cpp b/external/algorithms/core/src/random.cpp new file mode 100644 index 0000000000000000000000000000000000000000..97e337815dbfc67e389304d9f62a6a889ec08cae --- /dev/null +++ b/external/algorithms/core/src/random.cpp @@ -0,0 +1,41 @@ +#include + +#include +#include +#include +#include + +namespace algorithms { + +void RandomSvc::init() { } + +double RandomSvc::uniform(double a, double b) const { + // initialize the random uniform number generator (runif) in a range 0 to 1 + static std::uniform_real_distribution<> runif(a, b); + return runif(rng); +} + +double RandomSvc::normal(double mu, double sigma) const { + // https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform#Implementation + constexpr double epsilon = std::numeric_limits::epsilon(); + constexpr double two_pi = 2.0 * M_PI; + + static std::uniform_real_distribution<> runif(0.0, 1.0); + + double u1, u2; + do + { + u1 = runif(rng); + } + while (u1 <= epsilon); + u2 = runif(rng); + + //compute z0 and z1 + auto mag = sigma * sqrt(-2.0 * log(u1)); + auto z0 = mag * cos(two_pi * u2) + mu; + //auto z1 = mag * sin(two_pi * u2) + mu; + + return z0; +} + +} // namespace algorithms