diff --git a/JugDigi/CMakeLists.txt b/JugDigi/CMakeLists.txt index 2f33672..bc70eac 100644 --- a/JugDigi/CMakeLists.txt +++ b/JugDigi/CMakeLists.txt @@ -14,10 +14,12 @@ gaudi_add_module(JugDigiPlugins JugBase JugAlgo EICrecon::algorithms_calorimetry_library + EICrecon::algorithms_digi_library ROOT::Core ROOT::RIO ROOT::Tree EDM4HEP::edm4hep EDM4EIC::edm4eic k4FWCore::k4FWCore + GENCONF_OPTIONS --load-library=JugAlgoPreload.so ) target_include_directories(JugDigiPlugins PUBLIC diff --git a/JugDigi/src/components/PhotoMultiplierDigi.cpp b/JugDigi/src/components/PhotoMultiplierDigi.cpp index 82160a0..d8a41cf 100644 --- a/JugDigi/src/components/PhotoMultiplierDigi.cpp +++ b/JugDigi/src/components/PhotoMultiplierDigi.cpp @@ -1,221 +1,8 @@ // SPDX-License-Identifier: LGPL-3.0-or-later -// Copyright (C) 2022 Chao Peng +// Copyright (C) 2024 Wouter Deconinck -/* General PhotoMultiplier Digitization - * - * Apply the given quantum efficiency for photon detection - * Converts the number of detected photons to signal amplitude - * - * Author: Chao Peng (ANL) - * Date: 10/02/2020 - */ - -#include -#include -#include -#include - -#include "GaudiAlg/GaudiAlgorithm.h" -#include "GaudiAlg/Transformer.h" -#include "GaudiAlg/GaudiTool.h" -#include "GaudiKernel/RndmGenerators.h" -#include "GaudiKernel/PhysicalConstants.h" - -#include - -// Event Model related classes -#include "edm4eic/RawTrackerHitCollection.h" -#include "edm4hep/MCParticleCollection.h" -#include "edm4hep/SimTrackerHitCollection.h" - - -using namespace Gaudi::Units; - -namespace Jug::Digi { - -/** PhotoMultiplierDigi. - * - * \ingroup digi - */ -class PhotoMultiplierDigi : public GaudiAlgorithm -{ -public: - DataHandle - m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this}; - DataHandle - m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this}; - Gaudi::Property>> - u_quantumEfficiency{this, "quantumEfficiency", {{2.6*eV, 0.3}, {7.0*eV, 0.3}}}; - Gaudi::Property m_hitTimeWindow{this, "hitTimeWindow", 20.0*ns}; - Gaudi::Property m_timeStep{this, "timeStep", 0.0625*ns}; - Gaudi::Property m_speMean{this, "speMean", 80.0}; - Gaudi::Property m_speError{this, "speError", 16.0}; - Gaudi::Property m_pedMean{this, "pedMean", 200.0}; - Gaudi::Property m_pedError{this, "pedError", 3.0}; - Rndm::Numbers m_rngUni, m_rngNorm; - - // constructor - PhotoMultiplierDigi(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; - } - - auto randSvc = svc("RndmGenSvc", true); - auto sc1 = m_rngUni.initialize(randSvc, Rndm::Flat(0., 1.)); - auto sc2 = m_rngNorm.initialize(randSvc, Rndm::Gauss(0., 1.)); - if (!sc1.isSuccess() || !sc2.isSuccess()) { - error() << "Cannot initialize random generator!" << endmsg; - return StatusCode::FAILURE; - } - - qe_init(); - - return StatusCode::SUCCESS; - } - - StatusCode execute() override - { - // input collection - const auto &sim = *m_inputHitCollection.get(); - // Create output collections - auto &raw = *m_outputHitCollection.createAndPut(); - - struct HitData { int npe; double signal; double time; }; - std::unordered_map> hit_groups; - // collect the photon hit in the same cell - // calculate signal - for(const auto& ahit : sim) { - // quantum efficiency - if (!qe_pass(ahit.getEDep(), m_rngUni())) { - continue; - } - // cell id, time, signal amplitude - uint64_t id = ahit.getCellID(); - double time = ahit.getMCParticle().getTime(); - double amp = m_speMean + m_rngNorm()*m_speError; - - // group hits - auto it = hit_groups.find(id); - if (it != hit_groups.end()) { - size_t i = 0; - for (auto git = it->second.begin(); git != it->second.end(); ++git, ++i) { - if (std::abs(time - git->time) <= (m_hitTimeWindow/ns)) { - git->npe += 1; - git->signal += amp; - break; - } - } - // no hits group found - if (i >= it->second.size()) { - it->second.emplace_back(HitData{1, amp + m_pedMean + m_pedError*m_rngNorm(), time}); - } - } else { - hit_groups[id] = {HitData{1, amp + m_pedMean + m_pedError*m_rngNorm(), time}}; - } - } - - // build hit - for (auto &it : hit_groups) { - for (auto &data : it.second) { - raw.create( - it.first, - static_cast(data.signal), - static_cast(data.time/(m_timeStep/ns)) - ); - } - } - - return StatusCode::SUCCESS; - } - -private: - void qe_init() - { - auto &qeff = u_quantumEfficiency.value(); - - // sort quantum efficiency data first - std::sort(qeff.begin(), qeff.end(), - [] (const std::pair &v1, const std::pair &v2) { - return v1.first < v2.first; - }); - - // sanity checks - if (qeff.empty()) { - qeff = {{2.6*eV, 0.3}, {7.0*eV, 0.3}}; - warning() << "Invalid quantum efficiency data provided, using default values: " << qeff << endmsg; - } - if (qeff.front().first > 3.0*eV) { - warning() << "Quantum efficiency data start from " << qeff.front().first/eV - << " eV, maybe you are using wrong units?" << endmsg; - } - if (qeff.back().first < 6.0*eV) { - warning() << "Quantum efficiency data end at " << qeff.back().first/eV - << " eV, maybe you are using wrong units?" << endmsg; - } - } - - // helper function for linear interpolation - // Comp return is defined as: equal, 0; greater, > 0; less, < 0 - template - RndmIter interval_search(RndmIter beg, RndmIter end, const T &val, Compare comp) const - { - // special cases - auto dist = std::distance(beg, end); - if ((dist < 2) || (comp(*beg, val) > 0) || (comp(*std::prev(end), val) < 0)) { - return end; - } - auto mid = std::next(beg, dist / 2); - - while (mid != end) { - if (comp(*mid, val) == 0) { - return mid; - } else if (comp(*mid, val) > 0) { - end = mid; - } else { - beg = std::next(mid); - } - mid = std::next(beg, std::distance(beg, end)/2); - } - - if (mid == end || comp(*mid, val) > 0) { - return std::prev(mid); - } - return mid; - } - - bool qe_pass(double ev, double rand) const - { - const auto &qeff = u_quantumEfficiency.value(); - auto it = interval_search(qeff.begin(), qeff.end(), ev, - [] (const std::pair &vals, double val) { - return vals.first - val; - }); - - if (it == qeff.end()) { - // info() << ev/eV << " eV is out of QE data range, assuming 0% efficiency" << endmsg; - return false; - } - - double prob = it->second; - auto itn = std::next(it); - if (itn != qeff.end() && (itn->first - it->first != 0)) { - prob = (it->second*(itn->first - ev) + itn->second*(ev - it->first)) / (itn->first - it->first); - } - - // info() << ev/eV << " eV, QE: " << prob*100. << "%" << endmsg; - return rand <= prob; - } -}; +#include +#include // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) -DECLARE_COMPONENT(PhotoMultiplierDigi) - -} // namespace Jug::Digi +JUGALGO_DEFINE_ALGORITHM(PhotoMultiplierDigi, eicrecon::PhotoMultiplierHitDigi, Jug::Digi) diff --git a/JugDigi/src/components/SiliconTrackerDigi.cpp b/JugDigi/src/components/SiliconTrackerDigi.cpp index c54dbfe..9537392 100644 --- a/JugDigi/src/components/SiliconTrackerDigi.cpp +++ b/JugDigi/src/components/SiliconTrackerDigi.cpp @@ -1,104 +1,8 @@ // SPDX-License-Identifier: LGPL-3.0-or-later -// Copyright (C) 2022 Whitney Armstrong, Wouter Deconinck, Sylvester Joosten +// Copyright (C) 2024 Wouter Deconinck -#include -#include +#include +#include -#include "Gaudi/Property.h" -#include "GaudiAlg/GaudiAlgorithm.h" -#include "GaudiAlg/GaudiTool.h" -#include "GaudiAlg/Transformer.h" -#include "GaudiKernel/PhysicalConstants.h" -#include "GaudiKernel/RndmGenerators.h" - -#include - -// Event Model related classes -// edm4hep's tracker hit is the input collectiopn -#include "edm4hep/MCParticle.h" -#include "edm4hep/SimTrackerHitCollection.h" -// edm4eic's RawTrackerHit is the output -#include "edm4eic/RawTrackerHitCollection.h" - -namespace Jug::Digi { - -/** Silicon detector digitization. - * - * \ingroup digi - */ -class SiliconTrackerDigi : public GaudiAlgorithm { -private: - Gaudi::Property m_timeResolution{this, "timeResolution", 10}; // todo : add units - Gaudi::Property m_threshold{this, "threshold", 0. * Gaudi::Units::keV}; - Rndm::Numbers m_gaussDist; - DataHandle m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, - this}; - DataHandle m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, - this}; - -public: - // ill-formed: using GaudiAlgorithm::GaudiAlgorithm; - SiliconTrackerDigi(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; - } - IRndmGenSvc* randSvc = svc("RndmGenSvc", true); - StatusCode sc = m_gaussDist.initialize(randSvc, Rndm::Gauss(0.0, m_timeResolution.value())); - if (!sc.isSuccess()) { - return StatusCode::FAILURE; - } - return StatusCode::SUCCESS; - } - StatusCode execute() override { - // input collection - const auto* const simhits = m_inputHitCollection.get(); - // Create output collections - auto* rawhits = m_outputHitCollection.createAndPut(); - // edm4eic::RawTrackerHitCollection* rawHitCollection = new edm4eic::RawTrackerHitCollection(); - std::map cell_hit_map; - for (const auto& ahit : *simhits) { - if (msgLevel(MSG::DEBUG)) { - debug() << "--------------------" << ahit.getCellID() << endmsg; - debug() << "Hit in cellID = " << ahit.getCellID() << endmsg; - debug() << " position = (" << ahit.getPosition().x << "," << ahit.getPosition().y << "," - << ahit.getPosition().z << ")" << endmsg; - debug() << " xy_radius = " << std::hypot(ahit.getPosition().x, ahit.getPosition().y) << endmsg; - debug() << " momentum = (" << ahit.getMomentum().x << "," << ahit.getMomentum().y << "," - << ahit.getMomentum().z << ")" << endmsg; - } - if (ahit.getEDep() * Gaudi::Units::keV < m_threshold) { - if (msgLevel(MSG::DEBUG)) { - debug() << " edep = " << ahit.getEDep() << " (below threshold of " << m_threshold / Gaudi::Units::keV - << " keV)" << endmsg; - } - continue; - } else { - if (msgLevel(MSG::DEBUG)) { - debug() << " edep = " << ahit.getEDep() << endmsg; - } - } - if (cell_hit_map.count(ahit.getCellID()) == 0) { - cell_hit_map[ahit.getCellID()] = rawhits->size(); - rawhits->create( - ahit.getCellID(), - ahit.getMCParticle().getTime() * 1e6 + m_gaussDist() * 1e3, // ns->fs - std::llround(ahit.getEDep() * 1e6) - ); - } else { - auto hit = (*rawhits)[cell_hit_map[ahit.getCellID()]]; - hit.setTimeStamp(ahit.getMCParticle().getTime() * 1e6 + m_gaussDist() * 1e3); - auto ch = hit.getCharge(); - hit.setCharge(ch + std::llround(ahit.getEDep() * 1e6)); - } - } - return StatusCode::SUCCESS; - } -}; // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) -DECLARE_COMPONENT(SiliconTrackerDigi) - -} // namespace Jug::Digi +JUGALGO_DEFINE_ALGORITHM(SiliconTrackerDigi, eicrecon::SiliconTrackerDigi, Jug::Digi)