src/components/TrackParamTruthInit.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/Property.h"
0007 #include "Gaudi/Algorithm.h"
0008 #include "GaudiKernel/PhysicalConstants.h"
0009 #include "GaudiKernel/RndmGenerators.h"
0010 #include "GaudiKernel/ToolHandle.h"
0011
0012 #include <k4FWCore/DataHandle.h>
0013 #include <k4Interface/IGeoSvc.h>
0014 #include "JugBase/IParticleSvc.h"
0015 #include "ActsExamples/EventData/Track.hpp"
0016 #include "Acts/Definitions/Units.hpp"
0017 #include "Acts/Definitions/Common.hpp"
0018
0019 #include "edm4eic/TrackerHitCollection.h"
0020 #include "edm4hep/MCParticleCollection.h"
0021 #include "Math/Vector3D.h"
0022 #include "Acts/Surfaces/PerigeeSurface.hpp"
0023
0024
0025   ///// (Reconstructed) track parameters e.g. close to the vertex.
0026   //using TrackParameters = Acts::CurvilinearTrackParameters;
0027
0028   ///// Container of reconstructed track states for multiple tracks.
0029   //using TrackParametersContainer = std::vector<TrackParameters>;
0030
0031   ///// MultiTrajectory definition
0032   //using Trajectory = Acts::MultiTrajectory<SourceLink>;
0033
0034   ///// Container for the truth fitting/finding track(s)
0035   //using TrajectoryContainer = std::vector<SimMultiTrajectory>;
0036
0037 namespace Jug::Reco {
0038
0039   /** Initial Track parameters from MC truth.
0040    *
0041    *  TrackParmetersContainer
0042    *  \ingroup tracking
0043    */
0044   class TrackParamTruthInit : public Gaudi::Algorithm {
0045   private:
0046     mutable DataHandle<edm4hep::MCParticleCollection> m_inputMCParticles{"inputMCParticles", Gaudi::DataHandle::Reader,
0047                                                                     this};
0048     mutable DataHandle<ActsExamples::TrackParametersContainer>         m_outputInitialTrackParameters{"outputInitialTrackParameters",
0049                                                                         Gaudi::DataHandle::Writer, this};
0050
0051     // selection settings
0052     Gaudi::Property<double> m_maxVertexX{this, "maxVertexX", 80. * Gaudi::Units::mm};
0053     Gaudi::Property<double> m_maxVertexY{this, "maxVertexY", 80. * Gaudi::Units::mm};
0054     Gaudi::Property<double> m_maxVertexZ{this, "maxVertexZ", 200. * Gaudi::Units::mm};
0055     Gaudi::Property<double> m_minMomentum{this, "minMomentum", 100. * Gaudi::Units::MeV};
0056     Gaudi::Property<double> m_maxEtaForward{this, "maxEtaForward", 4.0};
0057     Gaudi::Property<double> m_maxEtaBackward{this, "maxEtaBackward", 4.1};
0058
0059     SmartIF<IParticleSvc> m_pidSvc;
0060
0061   public:
0062     TrackParamTruthInit(const std::string& name, ISvcLocator* svcLoc)
0063         : Gaudi::Algorithm(name, svcLoc) {
0064       declareProperty("inputMCParticles", m_inputMCParticles, "");
0065       declareProperty("outputInitialTrackParameters", m_outputInitialTrackParameters, "");
0066     }
0067
0068     StatusCode initialize() override {
0069       if (Gaudi::Algorithm::initialize().isFailure()) {
0070         return StatusCode::FAILURE;
0071       }
0072       IRndmGenSvc* randSvc = Gaudi::svcLocator()->service<IRndmGenSvc>("RndmGenSvc", true);
0073       if (randSvc == nullptr) {
0074         return StatusCode::FAILURE;
0075       }
0076       m_pidSvc = service("ParticleSvc");
0077       if (!m_pidSvc) {
0078         error() << "Unable to locate Particle Service. "
0079                 << "Make sure you have ParticleSvc in the configuration."
0080                 << endmsg;
0081         return StatusCode::FAILURE;
0082       }
0083       return StatusCode::SUCCESS;
0084     }
0085
0086     StatusCode execute(const EventContext&) const override {
0087       // input collection
0088       const auto* const mcparts = m_inputMCParticles.get();
0089       // Create output collections
0090       auto* init_trk_params = m_outputInitialTrackParameters.createAndPut();
0091
0092       for(const auto& part : *mcparts) {
0093
0094         // getGeneratorStatus = 1 means thrown G4Primary, but dd4gun uses getGeneratorStatus == 0
0095         if (part.getGeneratorStatus() > 1 ) {
0096           if (msgLevel(MSG::DEBUG)) {
0097             debug() << "ignoring particle with generatorStatus = " << part.getGeneratorStatus() << endmsg;
0098           }
0099           continue;
0100         }
0101
0102         // require close to interaction vertex
0103         if (abs(part.getVertex().x) * Gaudi::Units::mm > m_maxVertexX
0104          || abs(part.getVertex().y) * Gaudi::Units::mm > m_maxVertexY
0105          || abs(part.getVertex().z) * Gaudi::Units::mm > m_maxVertexZ) {
0106           if (msgLevel(MSG::DEBUG)) {
0107             debug() << "ignoring particle with vs = " << part.getVertex() << " mm" << endmsg;
0108           }
0109           continue;
0110         }
0111
0112         // require minimum momentum
0113         const auto& p = part.getMomentum();
0114         const auto pmag = std::hypot(p.x, p.y, p.z);
0115         if (pmag * Gaudi::Units::GeV < m_minMomentum) {
0116           if (msgLevel(MSG::DEBUG)) {
0117             debug() << "ignoring particle with p = " << pmag << " GeV" << endmsg;
0118           }
0119           continue;
0120         }
0121
0122         // require minimum pseudorapidity
0123         const auto phi   = std::atan2(p.y, p.x);
0124         const auto theta = std::atan2(std::hypot(p.x, p.y), p.z);
0125         const auto eta   = -std::log(std::tan(theta/2));
0126         if (eta > m_maxEtaForward || eta < -std::abs(m_maxEtaBackward)) {
0127           if (msgLevel(MSG::DEBUG)) {
0128             debug() << "ignoring particle with Eta = " << eta << endmsg;
0129           }
0130           continue;
0131         }
0132
0133         // get the particle charge
0134         // note that we cannot trust the mcparticles charge, as DD4hep
0135         // sets this value to zero! let's lookup by PDGID instead
0136         const double charge = m_pidSvc->particle(part.getPDG()).charge;
0137         if (abs(charge) < std::numeric_limits<double>::epsilon()) {
0138           if (msgLevel(MSG::DEBUG)) {
0139             debug() << "ignoring neutral particle" << endmsg;
0140           }
0141           continue;
0142         }
0143
0144         using Acts::UnitConstants::GeV;
0145         using Acts::UnitConstants::MeV;
0146         using Acts::UnitConstants::mm;
0147         using Acts::UnitConstants::um;
0148         using Acts::UnitConstants::ns;
0149
0150         // build some track cov matrix
0151         Acts::BoundSquareMatrix cov                    = Acts::BoundSquareMatrix::Zero();
0152         cov(Acts::eBoundLoc0, Acts::eBoundLoc0)     = 1000*um*1000*um;
0153         cov(Acts::eBoundLoc1, Acts::eBoundLoc1)     = 1000*um*1000*um;
0154         cov(Acts::eBoundPhi, Acts::eBoundPhi)       = 0.05*0.05;
0155         cov(Acts::eBoundTheta, Acts::eBoundTheta)   = 0.01*0.01;
0156         cov(Acts::eBoundQOverP, Acts::eBoundQOverP) = (0.1*0.1) / (GeV*GeV);
0157         cov(Acts::eBoundTime, Acts::eBoundTime)     = 10.0e9*ns*10.0e9*ns;
0158
0159         Acts::BoundVector  params;
0160         params(Acts::eBoundLoc0)   = 0.0 * mm ;  // cylinder radius
0161         params(Acts::eBoundLoc1)   = 0.0 * mm ; // cylinder length
0162         params(Acts::eBoundPhi)    = phi;
0163         params(Acts::eBoundTheta)  = theta;
0164         params(Acts::eBoundQOverP) = charge / (pmag * GeV);
0165         params(Acts::eBoundTime)   = part.getTime() * ns;
0166
0167         //// Construct a perigee surface as the target surface
0168         auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(
0169             Acts::Vector3{part.getVertex().x * mm, part.getVertex().y * mm, part.getVertex().z * mm});
0170
0171         //params(Acts::eBoundQOverP) = charge/p;
0172         init_trk_params->push_back({pSurface, params, cov, Acts::ParticleHypothesis::pion()});
0173         // std::make_optional(std::move(cov))
0174
0175         if (msgLevel(MSG::DEBUG)) {
0176           debug() << "Invoke track finding seeded by truth particle with p = " << pmag << " GeV" << endmsg;
0177           debug() << "                                              charge = " << charge << endmsg;
0178           debug() << "                                                 q/p = " << charge / pmag << " e/GeV" << endmsg;
0179         }
0180         //Acts::BoundMatrix cov           = Acts::BoundMatrix::Zero();
0181         //cov(Acts::eLOC_0, Acts::eLOC_0) = ahit.covMatrix(0)*ahit.covMatrix(0);
0182         //cov(Acts::eLOC_1, Acts::eLOC_1) = ahit.covMatrix(1)*ahit.covMatrix(1);
0183       }
0184       return StatusCode::SUCCESS;
0185     }
0186   };
0187   // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0188   DECLARE_COMPONENT(TrackParamTruthInit)
0189
0190 } // namespace Jug::reco