algorithms/tracking/TrackParamTruthInit.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 - 2025 Whitney Armstrong, Wouter Deconinck, Sylvester Joosten, Dmitry Romanov
0003
0004 #include "TrackParamTruthInit.h"
0005
0006 #include <Acts/Definitions/Algebra.hpp>
0007 #include <Acts/Surfaces/PerigeeSurface.hpp>
0008 #include <Acts/Surfaces/Surface.hpp>
0009 #include <Acts/Utilities/Result.hpp>
0010 #include <Evaluator/DD4hepUnits.h>
0011 #include <algorithms/logger.h>
0012 #include <edm4eic/Cov6f.h>
0013 #include <edm4hep/Vector2f.h>
0014 #include <edm4hep/Vector3d.h>
0015 #include <fmt/core.h>
0016 #include <Eigen/Core>
0017 #include <cmath>
0018 #include <cstdlib>
0019 #include <gsl/pointers>
0020 #include <limits>
0021 #include <memory>
0022 #include <random>
0023
0024 #include "extensions/spdlog/SpdlogFormatters.h" // IWYU pragma: keep
0025
0026 namespace eicrecon {
0027
0028 void TrackParamTruthInit::process(const Input& input, const Output& output) const {
0029   // MCParticles uses numerical values in its specified units,
0030   // while m_cfg is in the DD4hep unit system
0031
0032   const auto [headers, mcparticles] = input;
0033   auto [track_parameters]           = output;
0034
0035   // local random generator
0036   auto seed = m_uid.getUniqueID(*headers, "TrackParamTruthInit");
0037   std::default_random_engine generator(seed);
0038   std::normal_distribution<double> gaussian;
0039
0040   // Loop over input particles
0041   for (const auto& mcparticle : *mcparticles) {
0042
0043     // require generatorStatus == 1 for stable generated particles in HepMC3 and DDSim gun
0044     if (mcparticle.getGeneratorStatus() != 1) {
0045       trace("ignoring particle with generatorStatus = {}", mcparticle.getGeneratorStatus());
0046       continue;
0047     }
0048
0049     // require close to interaction vertex
0050     auto v = mcparticle.getVertex();
0051     if (std::abs(v.x) * dd4hep::mm > m_cfg.maxVertexX ||
0052         std::abs(v.y) * dd4hep::mm > m_cfg.maxVertexY ||
0053         std::abs(v.z) * dd4hep::mm > m_cfg.maxVertexZ) {
0054       trace("ignoring particle with vs = {} [mm]", v);
0055       continue;
0056     }
0057
0058     // require minimum momentum
0059     const auto& p   = mcparticle.getMomentum();
0060     const auto pmag = std::hypot(p.x, p.y, p.z);
0061     if (pmag * dd4hep::GeV < m_cfg.minMomentum) {
0062       trace("ignoring particle with p = {} GeV ", pmag);
0063       continue;
0064     }
0065
0066     // require minimum pseudorapidity
0067     const auto phi   = std::atan2(p.y, p.x);
0068     const auto theta = std::atan2(std::hypot(p.x, p.y), p.z);
0069     const auto eta   = -std::log(std::tan(theta / 2));
0070     if (eta > m_cfg.maxEtaForward || eta < -std::abs(m_cfg.maxEtaBackward)) {
0071       trace("ignoring particle with Eta = {}", eta);
0072       continue;
0073     }
0074
0075     // get the particle charge
0076     // note that we cannot trust the mcparticles charge, as DD4hep
0077     // sets this value to zero! let's lookup by PDGID instead
0078     const auto pdg      = mcparticle.getPDG();
0079     const auto particle = m_particleSvc.particle(pdg);
0080     double charge       = std::copysign(1.0, particle.charge);
0081     if (std::abs(particle.charge) < std::numeric_limits<double>::epsilon()) {
0082       trace("ignoring neutral particle");
0083       continue;
0084     }
0085
0086     // modify initial momentum to avoid bleeding truth to results when fit fails
0087     const auto pinit = pmag * (1.0 + m_cfg.momentumSmear * gaussian(generator));
0088
0089     // define line surface for local position values
0090     auto perigee = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3(0, 0, 0));
0091
0092     // track particle back to transverse point-of-closest approach
0093     // with respect to the defined line surface
0094     auto linesurface_parameter = -(v.x * p.x + v.y * p.y) / (p.x * p.x + p.y * p.y);
0095
0096     auto xpca = v.x + linesurface_parameter * p.x;
0097     auto ypca = v.y + linesurface_parameter * p.y;
0098     auto zpca = v.z + linesurface_parameter * p.z;
0099
0100     Acts::Vector3 global(xpca, ypca, zpca);
0101     Acts::Vector3 direction(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
0102
0103     // convert from global to local coordinates using the defined line surface
0104     auto local = perigee->globalToLocal(m_geoSvc->getActsGeometryContext(), global, direction);
0105
0106     if (!local.ok()) {
0107       error("skipping the track because globaltoLocal function failed");
0108       continue;
0109     }
0110
0111     Acts::Vector2 localpos = local.value();
0112
0113     // Insert into edm4eic::TrackParameters, which uses numerical values in its specified units
0114     auto track_parameter = track_parameters->create();
0115     track_parameter.setType(-1); // type --> seed(-1)
0116     track_parameter.setLoc({static_cast<float>(localpos(0)),
0117                             static_cast<float>(localpos(1))}); // 2d location on surface [mm]
0118     track_parameter.setPhi(phi);                               // phi [rad]
0119     track_parameter.setTheta(theta);                           // theta [rad]
0120     track_parameter.setQOverP(charge / (pinit / dd4hep::GeV)); // Q/p [e/GeV]
0121     track_parameter.setTime(mcparticle.getTime());             // time [ns]
0122     edm4eic::Cov6f cov;
0123     cov(0, 0) = 1.0;  // loc0
0124     cov(1, 1) = 1.0;  // loc1
0125     cov(2, 2) = 0.05; // phi
0126     cov(3, 3) = 0.01; // theta
0127     cov(4, 4) = 0.1;  // qOverP
0128     cov(5, 5) = 10e9; // time
0129     track_parameter.setCovariance(cov);
0130
0131     // Debug output
0132     if (level() <= algorithms::LogLevel::kDebug) {
0133       debug("Invoke track finding seeded by truth particle with:");
0134       debug("   p     = {} GeV (smeared to {} GeV)", pmag / dd4hep::GeV, pinit / dd4hep::GeV);
0135       debug("   q     = {}", charge);
0136       debug("   q/p   = {} e/GeV (smeared to {} e/GeV)", charge / (pmag / dd4hep::GeV),
0137             charge / (pinit / dd4hep::GeV));
0138       debug("   theta = {}", theta);
0139       debug("   phi   = {}", phi);
0140     }
0141   }
0142 }
0143
0144 } // namespace eicrecon