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