src/Utilities/VertexTruthUtility.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/Utilities/VertexTruthUtility.hpp"
0010
0011 #include "ActsExamples/EventData/SimParticle.hpp"
0012 #include "ActsExamples/EventData/Track.hpp"
0013
0014 #include <map>
0015 #include <set>
0016 #include <vector>
0017
0018 std::uint32_t ActsExamples::getNumberOfReconstructableVertices(
0019     const SimParticleContainer& collection) {
0020   // map for finding frequency
0021   std::map<std::uint32_t, std::uint32_t> fmap;
0022
0023   std::vector<std::uint32_t> reconstructableTruthVertices;
0024
0025   // traverse the array for frequency
0026   for (const SimParticle& p : collection) {
0027     std::uint32_t generation = p.particleId().generation();
0028     if (generation > 0) {
0029       // truthparticle from secondary vtx
0030       continue;
0031     }
0032     std::uint32_t priVtxId = p.particleId().vertexPrimary();
0033     fmap[priVtxId]++;
0034   }
0035
0036   // iterate over the map
0037   for (const auto& [priVtxId, occurrence] : fmap) {
0038     // Require at least 2 tracks
0039     if (occurrence > 1) {
0040       reconstructableTruthVertices.push_back(priVtxId);
0041     }
0042   }
0043
0044   return reconstructableTruthVertices.size();
0045 }
0046
0047 std::uint32_t ActsExamples::getNumberOfTruePriVertices(
0048     const SimParticleContainer& collection) {
0049   // Vector to store indices of all primary vertices
0050   std::set<std::uint32_t> allPriVtxIds;
0051   for (const SimParticle& p : collection) {
0052     std::uint32_t priVtxId = p.particleId().vertexPrimary();
0053     std::uint32_t generation = p.particleId().generation();
0054     if (generation > 0) {
0055       // truthparticle from secondary vtx
0056       continue;
0057     }
0058     // Insert to set, removing duplicates
0059     allPriVtxIds.insert(priVtxId);
0060   }
0061   // Size of set corresponds to total number of primary vertices
0062   return allPriVtxIds.size();
0063 }
0064
0065 double ActsExamples::calcSumPt2(const Acts::Vertex& vtx, double minTrkWeight) {
0066   double sumPt2 = 0;
0067   for (const auto& trk : vtx.tracks()) {
0068     if (trk.trackWeight > minTrkWeight) {
0069       double pt = trk.originalParams.as<Acts::BoundTrackParameters>()
0070                       ->transverseMomentum();
0071       sumPt2 += pt * pt;
0072     }
0073   }
0074   return sumPt2;
0075 }
0076
0077 double ActsExamples::calculateTruthPrimaryVertexDensity(
0078     const SimVertexContainer& truthVertices, const Acts::Vertex& vtx,
0079     double vertexDensityWindow) {
0080   double z = vtx.fullPosition()[Acts::CoordinateIndices::eZ];
0081   int count = 0;
0082   for (const SimVertex& truthVertex : truthVertices) {
0083     if (truthVertex.vertexId().vertexSecondary() != 0) {
0084       continue;
0085     }
0086     double zTruth = truthVertex.position4[Acts::CoordinateIndices::eZ];
0087     if (std::abs(z - zTruth) <= vertexDensityWindow) {
0088       ++count;
0089     }
0090   }
0091   return count / (2 * vertexDensityWindow);
0092 }
0093
0094 const ActsExamples::SimParticle* ActsExamples::findParticle(
0095     const SimParticleContainer& particles,
0096     const TrackParticleMatching& trackParticleMatching, ConstTrackProxy track,
0097     const Acts::Logger& logger) {
0098   // Get the truth-matched particle
0099   auto imatched = trackParticleMatching.find(track.index());
0100   if (imatched == trackParticleMatching.end() ||
0101       !imatched->second.particle.has_value()) {
0102     ACTS_DEBUG("No truth particle associated with this track, index = "
0103                << track.index() << " tip index = " << track.tipIndex());
0104     return nullptr;
0105   }
0106
0107   const TrackMatchEntry& particleMatch = imatched->second;
0108
0109   auto iparticle = particles.find(SimBarcode{particleMatch.particle.value()});
0110   if (iparticle == particles.end()) {
0111     ACTS_DEBUG(
0112         "Truth particle found but not monitored with this track, index = "
0113         << track.index() << " tip index = " << track.tipIndex()
0114         << " and this barcode = " << particleMatch.particle.value());
0115     return {};
0116   }
0117
0118   return &(*iparticle);
0119 }
0120
0121 std::optional<ActsExamples::ConstTrackProxy> ActsExamples::findTrack(
0122     const ConstTrackContainer& tracks, const Acts::TrackAtVertex& trkAtVtx) {
0123   // Track parameters before the vertex fit
0124   const Acts::BoundTrackParameters& origTrack =
0125       *trkAtVtx.originalParams.as<Acts::BoundTrackParameters>();
0126
0127   // Finding the matching parameters in the container of all track
0128   // parameters. This allows us to identify the corresponding particle.
0129   // TODO this should not be necessary if the tracks at vertex would keep
0130   // this information
0131   for (ConstTrackProxy inputTrk : tracks) {
0132     const auto& params = inputTrk.parameters();
0133
0134     if (origTrack.parameters() == params) {
0135       return inputTrk;
0136     }
0137   }
0138
0139   return std::nullopt;
0140 }