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