File indexing completed on 2025-02-24 09:23:26
0001
0002
0003
0004 #include <algorithm>
0005 #include <cmath>
0006
0007 #include "Gaudi/Property.h"
0008 #include "GaudiAlg/GaudiTool.h"
0009 #include "GaudiAlg/Transformer.h"
0010 #include "GaudiKernel/PhysicalConstants.h"
0011 #include "GaudiKernel/RndmGenerators.h"
0012
0013 #include "JugBase/DataHandle.h"
0014
0015
0016
0017 #include "edm4hep/MCParticle.h"
0018 #include "edm4hep/SimTrackerHitCollection.h"
0019
0020 #include "edm4eic/RawTrackerHitCollection.h"
0021
0022 namespace Jug::Digi {
0023
0024
0025
0026
0027
0028 class SiliconTrackerDigi : public GaudiAlgorithm {
0029 private:
0030 Gaudi::Property<double> m_timeResolution{this, "timeResolution", 10};
0031 Gaudi::Property<double> m_threshold{this, "threshold", 0. * Gaudi::Units::keV};
0032 Rndm::Numbers m_gaussDist;
0033 DataHandle<edm4hep::SimTrackerHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader,
0034 this};
0035 DataHandle<edm4eic::RawTrackerHitCollection> m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer,
0036 this};
0037
0038 public:
0039
0040 SiliconTrackerDigi(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) {
0041 declareProperty("inputHitCollection", m_inputHitCollection, "");
0042 declareProperty("outputHitCollection", m_outputHitCollection, "");
0043 }
0044 StatusCode initialize() override {
0045 if (GaudiAlgorithm::initialize().isFailure()) {
0046 return StatusCode::FAILURE;
0047 }
0048 IRndmGenSvc* randSvc = svc<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() override {
0056
0057 const auto* const simhits = m_inputHitCollection.get();
0058
0059 auto* rawhits = m_outputHitCollection.createAndPut();
0060
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 edm4eic::RawTrackerHit rawhit(ahit.getCellID(),
0086 ahit.getMCParticle().getTime() * 1e6 + m_gaussDist() * 1e3,
0087 std::llround(ahit.getEDep() * 1e6));
0088 rawhits->push_back(rawhit);
0089 } else {
0090 auto hit = (*rawhits)[cell_hit_map[ahit.getCellID()]];
0091 hit.setTimeStamp(ahit.getMCParticle().getTime() * 1e6 + m_gaussDist() * 1e3);
0092 auto ch = hit.getCharge();
0093 hit.setCharge(ch + std::llround(ahit.getEDep() * 1e6));
0094 }
0095 }
0096 return StatusCode::SUCCESS;
0097 }
0098 };
0099
0100 DECLARE_COMPONENT(SiliconTrackerDigi)
0101
0102 }