Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-02-24 09:23:26

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/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 // Event Model related classes
0016 // edm4hep's tracker hit is the input collectiopn
0017 #include "edm4hep/MCParticle.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 GaudiAlgorithm {
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   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   //  ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
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     // 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         edm4eic::RawTrackerHit rawhit(ahit.getCellID(),
0086                                    ahit.getMCParticle().getTime() * 1e6 + m_gaussDist() * 1e3, // ns->fs
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 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0100 DECLARE_COMPONENT(SiliconTrackerDigi)
0101 
0102 } // namespace Jug::Digi