algorithms/tracking/TrackParamTruthInit.cc

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