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