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
0012 #include <Evaluator/DD4hepUnits.h>
0013 #include <edm4hep/Vector3d.h>
0014 #include <edm4hep/Vector3f.h>
0015 #include <edm4eic/Cov6f.h>
0016 #include <fmt/core.h>
0017 #include <spdlog/common.h>
0018 #include <Eigen/Core>
0019 #include <cmath>
0020 #include <limits>
0021 #include <memory>
0022 #include <utility>
0023
0024 #include "extensions/spdlog/SpdlogFormatters.h" // IWYU pragma: keep
0025
0026
0027 void eicrecon::TrackParamTruthInit::init(std::shared_ptr<const ActsGeometryProvider> geo_svc, const std::shared_ptr<spdlog::logger> logger) {
0028     m_log = logger;
0029     m_geoSvc = geo_svc;
0030 }
0031
0032 std::unique_ptr<edm4eic::TrackParametersCollection>
0033 eicrecon::TrackParamTruthInit::produce(const edm4hep::MCParticleCollection* mcparticles) {
0034     // MCParticles uses numerical values in its specified units,
0035     // while m_cfg is in the DD4hep unit system
0036
0037     // Create output collection
0038     auto track_parameters = std::make_unique<edm4eic::TrackParametersCollection>();
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             m_log->trace("ignoring particle with generatorStatus = {}", mcparticle.getGeneratorStatus());
0046             continue;
0047         }
0048
0049         // require close to interaction vertex
0050         auto v = mcparticle.getVertex();
0051         if (abs(v.x) * dd4hep::mm > m_cfg.maxVertexX ||
0052             abs(v.y) * dd4hep::mm > m_cfg.maxVertexY ||
0053             abs(v.z) * dd4hep::mm > m_cfg.maxVertexZ) {
0054             m_log->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             m_log->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             m_log->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 (abs(charge) < std::numeric_limits<double>::epsilon()) {
0082             m_log->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 * m_normDist(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         {
0108             m_log->error("skipping the track because globaltoLocal function failed");
0109             continue;
0110         }
0111
0112         Acts::Vector2 localpos = local.value();
0113
0114         // Insert into edm4eic::TrackParameters, which uses numerical values in its specified units
0115         auto track_parameter = track_parameters->create();
0116         track_parameter.setType(-1); // type --> seed(-1)
0117         track_parameter.setLoc({static_cast<float>(localpos(0)), 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 (m_log->level() <= spdlog::level::debug) {
0133             m_log->debug("Invoke track finding seeded by truth particle with:");
0134             m_log->debug("   p     = {} GeV (smeared to {} GeV)", pmag / dd4hep::GeV, pinit / dd4hep::GeV);
0135             m_log->debug("   q     = {}", charge);
0136             m_log->debug("   q/p   = {} e/GeV (smeared to {} e/GeV)", charge / (pmag / dd4hep::GeV), charge / (pinit / dd4hep::GeV));
0137             m_log->debug("   theta = {}", theta);
0138             m_log->debug("   phi   = {}", phi);
0139         }
0140     }
0141
0142     return std::move(track_parameters);
0143
0144 }