src/components/SiliconTrackerDigi.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Whitney Armstrong, Wouter Deconinck, Sylvester Joosten
0003
0004 #include <algorithm>
0005 #include <cmath>
0006
0007 #include "Gaudi/Property.h"
0008 #include "Gaudi/Algorithm.h"
0009 #include "GaudiKernel/PhysicalConstants.h"
0010 #include "GaudiKernel/RndmGenerators.h"
0011
0012 #include <k4FWCore/DataHandle.h>
0013
0014 // Event Model related classes
0015 // edm4hep's tracker hit is the input collectiopn
0016 #include "edm4hep/EDM4hepVersion.h"
0017 #include "edm4hep/MCParticleCollection.h"
0018 #include "edm4hep/SimTrackerHitCollection.h"
0019 // edm4eic's RawTrackerHit is the output
0020 #include "edm4eic/RawTrackerHitCollection.h"
0021
0022 namespace Jug::Digi {
0023
0024 /** Silicon detector digitization.
0025  *
0026  * \ingroup digi
0027  */
0028 class SiliconTrackerDigi : public Gaudi::Algorithm {
0029 private:
0030   Gaudi::Property<double> m_timeResolution{this, "timeResolution", 10}; // todo : add units
0031   Gaudi::Property<double> m_threshold{this, "threshold", 0. * Gaudi::Units::keV};
0032   Rndm::Numbers m_gaussDist;
0033   mutable DataHandle<edm4hep::SimTrackerHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader,
0034                                                                     this};
0035   mutable DataHandle<edm4eic::RawTrackerHitCollection> m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer,
0036                                                                   this};
0037
0038 public:
0039   //  ill-formed: using Gaudi::Algorithm::GaudiAlgorithm;
0040   SiliconTrackerDigi(const std::string& name, ISvcLocator* svcLoc) : Gaudi::Algorithm(name, svcLoc) {
0041     declareProperty("inputHitCollection", m_inputHitCollection, "");
0042     declareProperty("outputHitCollection", m_outputHitCollection, "");
0043   }
0044   StatusCode initialize() override {
0045     if (Gaudi::Algorithm::initialize().isFailure()) {
0046       return StatusCode::FAILURE;
0047     }
0048     IRndmGenSvc* randSvc = Gaudi::svcLocator()->service<IRndmGenSvc>("RndmGenSvc", true);
0049     StatusCode sc        = m_gaussDist.initialize(randSvc, Rndm::Gauss(0.0, m_timeResolution.value()));
0050     if (!sc.isSuccess()) {
0051       return StatusCode::FAILURE;
0052     }
0053     return StatusCode::SUCCESS;
0054   }
0055   StatusCode execute(const EventContext&) const override {
0056     // input collection
0057     const auto* const simhits = m_inputHitCollection.get();
0058     // Create output collections
0059     auto* rawhits = m_outputHitCollection.createAndPut();
0060     // edm4eic::RawTrackerHitCollection* rawHitCollection = new edm4eic::RawTrackerHitCollection();
0061     std::map<long long, int> cell_hit_map;
0062     for (const auto& ahit : *simhits) {
0063       if (msgLevel(MSG::DEBUG)) {
0064         debug() << "--------------------" << ahit.getCellID() << endmsg;
0065         debug() << "Hit in cellID = " << ahit.getCellID() << endmsg;
0066         debug() << "     position = (" << ahit.getPosition().x << "," << ahit.getPosition().y << ","
0067                 << ahit.getPosition().z << ")" << endmsg;
0068         debug() << "    xy_radius = " << std::hypot(ahit.getPosition().x, ahit.getPosition().y) << endmsg;
0069         debug() << "     momentum = (" << ahit.getMomentum().x << "," << ahit.getMomentum().y << ","
0070                 << ahit.getMomentum().z << ")" << endmsg;
0071       }
0072       if (ahit.getEDep() * Gaudi::Units::keV < m_threshold) {
0073         if (msgLevel(MSG::DEBUG)) {
0074           debug() << "         edep = " << ahit.getEDep() << " (below threshold of " << m_threshold / Gaudi::Units::keV
0075                   << " keV)" << endmsg;
0076         }
0077         continue;
0078       } else {
0079         if (msgLevel(MSG::DEBUG)) {
0080           debug() << "         edep = " << ahit.getEDep() << endmsg;
0081         }
0082       }
0083       if (cell_hit_map.count(ahit.getCellID()) == 0) {
0084         cell_hit_map[ahit.getCellID()] = rawhits->size();
0085         rawhits->create(
0086           ahit.getCellID(),
0087 #if EDM4HEP_BUILD_VERSION >= EDM4HEP_VERSION(0, 99, 0)
0088           ahit.getParticle().getTime() * 1e6 + m_gaussDist() * 1e3, // ns->fs
0089 #else
0090           ahit.getMCParticle().getTime() * 1e6 + m_gaussDist() * 1e3, // ns->fs
0091 #endif
0092           std::llround(ahit.getEDep() * 1e6)
0093         );
0094       } else {
0095         auto hit = (*rawhits)[cell_hit_map[ahit.getCellID()]];
0096 #if EDM4HEP_BUILD_VERSION >= EDM4HEP_VERSION(0, 99, 0)
0097         hit.setTimeStamp(ahit.getParticle().getTime() * 1e6 + m_gaussDist() * 1e3);
0098 #else
0099         hit.setTimeStamp(ahit.getMCParticle().getTime() * 1e6 + m_gaussDist() * 1e3);
0100 #endif
0101         auto ch = hit.getCharge();
0102         hit.setCharge(ch + std::llround(ahit.getEDep() * 1e6));
0103       }
0104     }
0105     return StatusCode::SUCCESS;
0106   }
0107 };
0108 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0109 DECLARE_COMPONENT(SiliconTrackerDigi)
0110
0111 } // namespace Jug::Digi