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