acts/src/TruthTrackSeeding.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Whitney Armstrong, Wouter Deconinck, Sylvester Joosten
0003
0004 #include <cmath>
0005 // Gaudi
0006 #include "GaudiAlg/GaudiAlgorithm.h"
0007 #include "GaudiKernel/ToolHandle.h"
0008 #include "GaudiAlg/Transformer.h"
0009 #include "GaudiAlg/GaudiTool.h"
0010 #include "GaudiKernel/RndmGenerators.h"
0011 #include "Gaudi/Property.h"
0012
0013 #include "JugBase/DataHandle.h"
0014 #include "JugBase/IGeoSvc.h"
0015 #include "JugBase/IParticleSvc.h"
0016 #include "JugTrack/Track.hpp"
0017
0018 //#include "Acts/Definitions/Units.hpp"
0019 //#include "Acts/Definitions/Common.hpp"
0020 //#include "Acts/EventData/Charge.hpp"
0021
0022 #include "edm4hep/MCParticleCollection.h"
0023 #include "edm4eic/TrackParametersCollection.h"
0024 #include "Math/Vector3D.h"
0025
0026
0027 namespace Jug::Reco {
0028
0029   /** Track seeding using MC truth.
0030    *
0031    *  \note "Seeding" algorithms are required to output a edm4eic::TrackParametersCollection, as opposed to the legacy "init"
0032    *  algorithms, such as  TrackParamTruthInit.
0033    *
0034    *  \ingroup tracking
0035    */
0036   class TruthTrackSeeding : public GaudiAlgorithm {
0037   private:
0038     DataHandle<edm4hep::MCParticleCollection> m_inputMCParticles{"inputMCParticles", Gaudi::DataHandle::Reader,
0039                                                                     this};
0040     DataHandle<edm4eic::TrackParametersCollection> m_outputTrackParameters{"outputTrackParameters",
0041                                                                        Gaudi::DataHandle::Writer, this};
0042     SmartIF<IParticleSvc> m_pidSvc;
0043
0044   public:
0045     TruthTrackSeeding(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) {
0046       declareProperty("inputMCParticles", m_inputMCParticles, "mcparticle truth data from npsim");
0047       declareProperty("outputTrackParameters", m_outputTrackParameters, "Output initial track parameters");
0048     }
0049
0050     StatusCode initialize() override {
0051       if (GaudiAlgorithm::initialize().isFailure()) {
0052         return StatusCode::FAILURE;
0053       }
0054       IRndmGenSvc* randSvc = svc<IRndmGenSvc>("RndmGenSvc", true);
0055       if (randSvc == nullptr) {
0056         return StatusCode::FAILURE;
0057       }
0058       m_pidSvc = service("ParticleSvc");
0059       if (!m_pidSvc) {
0060         error() << "Unable to locate Particle Service. "
0061                 << "Make sure you have ParticleSvc in the configuration."
0062                 << endmsg;
0063         return StatusCode::FAILURE;
0064       }
0065       return StatusCode::SUCCESS;
0066     }
0067
0068     StatusCode execute() override {
0069       // input collection
0070       const auto* const mcparts = m_inputMCParticles.get();
0071       // Create output collections
0072       auto* init_trk_params = m_outputTrackParameters.createAndPut();
0073
0074       for(const auto& part : *mcparts) {
0075
0076         // getGeneratorStatus = 1 means thrown G4Primary
0077         if(part.getGeneratorStatus() != 1 ) {
0078           continue;
0079         }
0080
0081         const auto& pvec = part.getMomentum();
0082         const auto p = std::hypot(pvec.x, pvec.y, pvec.z);
0083         const auto phi = std::atan2(pvec.x, pvec.y);
0084         const auto theta = std::atan2(std::hypot(pvec.x, pvec.y), pvec.z);
0085
0086         // get the particle charge
0087         // note that we cannot trust the mcparticles charge, as DD4hep
0088         // sets this value to zero! let's lookup by PDGID instead
0089         const auto charge = static_cast<float>(m_pidSvc->particle(part.getPDG()).charge);
0090         if (abs(charge) < std::numeric_limits<double>::epsilon()) {
0091           continue;
0092         }
0093
0094         const auto q_over_p = charge / p;
0095
0096         edm4eic::TrackParameters params{-1,                // type --> seed (-1)
0097                                      {0.0F, 0.0F},      // location on surface
0098                                      {0.1, 0.1, 0.1},   // Covariance on location
0099                                      theta,             // theta (rad)
0100                                      phi,               // phi  (rad)
0101                                      q_over_p * .05F,   // Q/P (e/GeV)
0102                                      {0.1, 0.1, 0.1},   // Covariance on theta/phi/Q/P
0103                                      part.getTime(),    // Time (ns)
0104                                      0.1,               // Error on time
0105                                      charge};           // Charge
0106
0107         ////// Construct a perigee surface as the target surface
0108         //auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(
0109         //    Acts::Vector3{part.getVertex().x * mm, part.getVertex().y * mm, part.getVertex().z * mm});
0110
0111         init_trk_params->push_back(params);
0112
0113         if (msgLevel(MSG::DEBUG)) {
0114           debug() << "Invoke track finding seeded by truth particle with p = " << p  << " GeV" << endmsg;
0115           debug() << "                                              charge = " << charge << endmsg;
0116           debug() << "                                                 q/p = " << charge / p << endmsg;
0117         }
0118       }
0119       return StatusCode::SUCCESS;
0120     }
0121   };
0122   // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0123   DECLARE_COMPONENT(TruthTrackSeeding)
0124
0125 } // namespace Jug::reco