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