File indexing completed on 2025-10-23 08:23:30
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "ActsExamples/Io/Root/VertexNTupleWriter.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/EventData/GenericBoundTrackParameters.hpp"
0014 #include "Acts/EventData/TrackParameters.hpp"
0015 #include "Acts/Propagator/Propagator.hpp"
0016 #include "Acts/Propagator/SympyStepper.hpp"
0017 #include "Acts/Surfaces/PerigeeSurface.hpp"
0018 #include "Acts/Surfaces/Surface.hpp"
0019 #include "Acts/Utilities/Intersection.hpp"
0020 #include "Acts/Utilities/Logger.hpp"
0021 #include "Acts/Utilities/UnitVectors.hpp"
0022 #include "Acts/Vertexing/TrackAtVertex.hpp"
0023 #include "ActsExamples/EventData/SimParticle.hpp"
0024 #include "ActsExamples/EventData/SimVertex.hpp"
0025 #include "ActsExamples/EventData/Track.hpp"
0026 #include "ActsExamples/EventData/TruthMatching.hpp"
0027 #include "ActsFatras/EventData/Barcode.hpp"
0028
0029 #include <algorithm>
0030 #include <cmath>
0031 #include <cstddef>
0032 #include <cstdint>
0033 #include <ios>
0034 #include <limits>
0035 #include <map>
0036 #include <memory>
0037 #include <numbers>
0038 #include <optional>
0039 #include <ostream>
0040 #include <set>
0041 #include <stdexcept>
0042 #include <utility>
0043
0044 #include <TFile.h>
0045 #include <TTree.h>
0046
0047 using Acts::VectorHelpers::eta;
0048 using Acts::VectorHelpers::perp;
0049 using Acts::VectorHelpers::phi;
0050 using Acts::VectorHelpers::theta;
0051
0052 namespace ActsExamples {
0053
0054 namespace {
0055
0056 constexpr double nan = std::numeric_limits<double>::quiet_NaN();
0057
0058 std::uint32_t getNumberOfReconstructableVertices(
0059 const SimParticleContainer& collection) {
0060
0061 std::map<std::uint32_t, std::uint32_t> fmap;
0062
0063 std::vector<std::uint32_t> reconstructableTruthVertices;
0064
0065
0066 for (const SimParticle& p : collection) {
0067 std::uint32_t generation = p.particleId().generation();
0068 if (generation > 0) {
0069
0070 continue;
0071 }
0072 std::uint32_t priVtxId = p.particleId().vertexPrimary();
0073 fmap[priVtxId]++;
0074 }
0075
0076
0077 for (const auto& [priVtxId, occurrence] : fmap) {
0078
0079 if (occurrence > 1) {
0080 reconstructableTruthVertices.push_back(priVtxId);
0081 }
0082 }
0083
0084 return reconstructableTruthVertices.size();
0085 }
0086
0087 std::uint32_t getNumberOfTruePriVertices(
0088 const SimParticleContainer& collection) {
0089
0090 std::set<std::uint32_t> allPriVtxIds;
0091 for (const SimParticle& p : collection) {
0092 std::uint32_t priVtxId = p.particleId().vertexPrimary();
0093 std::uint32_t generation = p.particleId().generation();
0094 if (generation > 0) {
0095
0096 continue;
0097 }
0098
0099 allPriVtxIds.insert(priVtxId);
0100 }
0101
0102 return allPriVtxIds.size();
0103 }
0104
0105 double calcSumPt2(const VertexNTupleWriter::Config& config,
0106 const Acts::Vertex& vtx) {
0107 double sumPt2 = 0;
0108 for (const auto& trk : vtx.tracks()) {
0109 if (trk.trackWeight > config.minTrkWeight) {
0110 double pt = trk.originalParams.as<Acts::BoundTrackParameters>()
0111 ->transverseMomentum();
0112 sumPt2 += pt * pt;
0113 }
0114 }
0115 return sumPt2;
0116 }
0117
0118
0119 double pull(double diff, double variance, const std::string& variableStr,
0120 bool afterFit, const Acts::Logger& logger) {
0121 if (variance <= 0) {
0122 std::string tempStr;
0123 if (afterFit) {
0124 tempStr = "after";
0125 } else {
0126 tempStr = "before";
0127 }
0128 ACTS_WARNING("Nonpositive variance " << tempStr << " vertex fit: Var("
0129 << variableStr << ") = " << variance
0130 << " <= 0.");
0131 return nan;
0132 }
0133 double std = std::sqrt(variance);
0134 return diff / std;
0135 }
0136
0137 double calculateTruthPrimaryVertexDensity(
0138 const VertexNTupleWriter::Config& config,
0139 const SimVertexContainer& truthVertices, const Acts::Vertex& vtx) {
0140 double z = vtx.fullPosition()[Acts::CoordinateIndices::eZ];
0141 int count = 0;
0142 for (const SimVertex& truthVertex : truthVertices) {
0143 if (truthVertex.vertexId().vertexSecondary() != 0) {
0144 continue;
0145 }
0146 double zTruth = truthVertex.position4[Acts::CoordinateIndices::eZ];
0147 if (std::abs(z - zTruth) <= config.vertexDensityWindow) {
0148 ++count;
0149 }
0150 }
0151 return count / (2 * config.vertexDensityWindow);
0152 }
0153
0154 const SimParticle* findParticle(
0155 const SimParticleContainer& particles,
0156 const TrackParticleMatching& trackParticleMatching, ConstTrackProxy track,
0157 const Acts::Logger& logger) {
0158
0159 auto imatched = trackParticleMatching.find(track.index());
0160 if (imatched == trackParticleMatching.end() ||
0161 !imatched->second.particle.has_value()) {
0162 ACTS_DEBUG("No truth particle associated with this track, index = "
0163 << track.index() << " tip index = " << track.tipIndex());
0164 return nullptr;
0165 }
0166
0167 const TrackMatchEntry& particleMatch = imatched->second;
0168
0169 auto iparticle = particles.find(SimBarcode{particleMatch.particle.value()});
0170 if (iparticle == particles.end()) {
0171 ACTS_DEBUG(
0172 "Truth particle found but not monitored with this track, index = "
0173 << track.index() << " tip index = " << track.tipIndex()
0174 << " and this barcode = " << particleMatch.particle.value());
0175 return {};
0176 }
0177
0178 return &(*iparticle);
0179 }
0180
0181 std::optional<ConstTrackProxy> findTrack(const ConstTrackContainer& tracks,
0182 const Acts::TrackAtVertex& trkAtVtx) {
0183
0184 const Acts::BoundTrackParameters& origTrack =
0185 *trkAtVtx.originalParams.as<Acts::BoundTrackParameters>();
0186
0187
0188
0189
0190
0191 for (ConstTrackProxy inputTrk : tracks) {
0192 const auto& params = inputTrk.parameters();
0193
0194 if (origTrack.parameters() == params) {
0195 return inputTrk;
0196 }
0197 }
0198
0199 return std::nullopt;
0200 }
0201
0202 }
0203
0204 VertexNTupleWriter::VertexNTupleWriter(const VertexNTupleWriter::Config& config,
0205 Acts::Logging::Level level)
0206 : WriterT(config.inputVertices, "VertexNTupleWriter", level),
0207 m_cfg(config) {
0208 if (m_cfg.filePath.empty()) {
0209 throw std::invalid_argument("Missing output filename");
0210 }
0211 if (m_cfg.treeName.empty()) {
0212 throw std::invalid_argument("Missing tree name");
0213 }
0214 if (m_cfg.inputTruthVertices.empty()) {
0215 throw std::invalid_argument("Collection with truth vertices missing");
0216 }
0217 if (m_cfg.inputParticles.empty()) {
0218 throw std::invalid_argument("Collection with particles missing");
0219 }
0220 if (m_cfg.inputSelectedParticles.empty()) {
0221 throw std::invalid_argument("Collection with selected particles missing");
0222 }
0223
0224 m_inputTracks.maybeInitialize(m_cfg.inputTracks);
0225 m_inputTruthVertices.initialize(m_cfg.inputTruthVertices);
0226 m_inputParticles.initialize(m_cfg.inputParticles);
0227 m_inputSelectedParticles.initialize(m_cfg.inputSelectedParticles);
0228 m_inputTrackParticleMatching.maybeInitialize(
0229 m_cfg.inputTrackParticleMatching);
0230
0231 if (m_cfg.writeTrackInfo && !m_inputTracks.isInitialized()) {
0232 throw std::invalid_argument("Missing input tracks collection");
0233 }
0234 if (m_cfg.writeTrackInfo && !m_inputTrackParticleMatching.isInitialized()) {
0235 throw std::invalid_argument("Missing input track particles matching");
0236 }
0237
0238
0239 m_outputFile = TFile::Open(m_cfg.filePath.c_str(), m_cfg.fileMode.c_str());
0240 if (m_outputFile == nullptr) {
0241 throw std::ios_base::failure("Could not open '" + m_cfg.filePath + "'");
0242 }
0243 m_outputFile->cd();
0244 m_outputTree = new TTree(m_cfg.treeName.c_str(), m_cfg.treeName.c_str());
0245 if (m_outputTree == nullptr) {
0246 throw std::bad_alloc();
0247 }
0248
0249 m_outputTree->Branch("event_nr", &m_eventNr);
0250
0251 m_outputTree->Branch("nRecoVtx", &m_nRecoVtx);
0252 m_outputTree->Branch("nTrueVtx", &m_nTrueVtx);
0253 m_outputTree->Branch("nCleanVtx", &m_nCleanVtx);
0254 m_outputTree->Branch("nMergedVtx", &m_nMergedVtx);
0255 m_outputTree->Branch("nSplitVtx", &m_nSplitVtx);
0256 m_outputTree->Branch("nVtxDetectorAcceptance", &m_nVtxDetAcceptance);
0257 m_outputTree->Branch("nVtxReconstructable", &m_nVtxReconstructable);
0258
0259 m_outputTree->Branch("nTracksRecoVtx", &m_nTracksOnRecoVertex);
0260
0261 m_outputTree->Branch("recoVertexTrackWeights", &m_recoVertexTrackWeights);
0262
0263 m_outputTree->Branch("sumPt2", &m_sumPt2);
0264
0265 m_outputTree->Branch("recoX", &m_recoX);
0266 m_outputTree->Branch("recoY", &m_recoY);
0267 m_outputTree->Branch("recoZ", &m_recoZ);
0268 m_outputTree->Branch("recoT", &m_recoT);
0269
0270 m_outputTree->Branch("covXX", &m_covXX);
0271 m_outputTree->Branch("covYY", &m_covYY);
0272 m_outputTree->Branch("covZZ", &m_covZZ);
0273 m_outputTree->Branch("covTT", &m_covTT);
0274 m_outputTree->Branch("covXY", &m_covXY);
0275 m_outputTree->Branch("covXZ", &m_covXZ);
0276 m_outputTree->Branch("covXT", &m_covXT);
0277 m_outputTree->Branch("covYZ", &m_covYZ);
0278 m_outputTree->Branch("covYT", &m_covYT);
0279 m_outputTree->Branch("covZT", &m_covZT);
0280
0281 m_outputTree->Branch("seedX", &m_seedX);
0282 m_outputTree->Branch("seedY", &m_seedY);
0283 m_outputTree->Branch("seedZ", &m_seedZ);
0284 m_outputTree->Branch("seedT", &m_seedT);
0285
0286 m_outputTree->Branch("vertex_primary", &m_vertexPrimary);
0287 m_outputTree->Branch("vertex_secondary", &m_vertexSecondary);
0288
0289 m_outputTree->Branch("nTracksTruthVtx", &m_nTracksOnTruthVertex);
0290
0291 m_outputTree->Branch("truthPrimaryVertexDensity",
0292 &m_truthPrimaryVertexDensity);
0293
0294 m_outputTree->Branch("truthVertexTrackWeights", &m_truthVertexTrackWeights);
0295 m_outputTree->Branch("truthVertexMatchRatio", &m_truthVertexMatchRatio);
0296 m_outputTree->Branch("recoVertexContamination", &m_recoVertexContamination);
0297
0298 m_outputTree->Branch("recoVertexClassification", &m_recoVertexClassification);
0299
0300 m_outputTree->Branch("truthX", &m_truthX);
0301 m_outputTree->Branch("truthY", &m_truthY);
0302 m_outputTree->Branch("truthZ", &m_truthZ);
0303 m_outputTree->Branch("truthT", &m_truthT);
0304
0305 m_outputTree->Branch("resX", &m_resX);
0306 m_outputTree->Branch("resY", &m_resY);
0307 m_outputTree->Branch("resZ", &m_resZ);
0308 m_outputTree->Branch("resT", &m_resT);
0309
0310 m_outputTree->Branch("resSeedZ", &m_resSeedZ);
0311 m_outputTree->Branch("resSeedT", &m_resSeedT);
0312
0313 m_outputTree->Branch("pullX", &m_pullX);
0314 m_outputTree->Branch("pullY", &m_pullY);
0315 m_outputTree->Branch("pullZ", &m_pullZ);
0316 m_outputTree->Branch("pullT", &m_pullT);
0317
0318 if (m_cfg.writeTrackInfo) {
0319 m_outputTree->Branch("trk_weight", &m_trkWeight);
0320
0321 m_outputTree->Branch("trk_recoPhi", &m_recoPhi);
0322 m_outputTree->Branch("trk_recoTheta", &m_recoTheta);
0323 m_outputTree->Branch("trk_recoQOverP", &m_recoQOverP);
0324
0325 m_outputTree->Branch("trk_recoPhiFitted", &m_recoPhiFitted);
0326 m_outputTree->Branch("trk_recoThetaFitted", &m_recoThetaFitted);
0327 m_outputTree->Branch("trk_recoQOverPFitted", &m_recoQOverPFitted);
0328
0329 m_outputTree->Branch("trk_particleId", &m_trkParticleId);
0330
0331 m_outputTree->Branch("trk_truthPhi", &m_truthPhi);
0332 m_outputTree->Branch("trk_truthTheta", &m_truthTheta);
0333 m_outputTree->Branch("trk_truthQOverP", &m_truthQOverP);
0334
0335 m_outputTree->Branch("trk_resPhi", &m_resPhi);
0336 m_outputTree->Branch("trk_resTheta", &m_resTheta);
0337 m_outputTree->Branch("trk_resQOverP", &m_resQOverP);
0338 m_outputTree->Branch("trk_momOverlap", &m_momOverlap);
0339
0340 m_outputTree->Branch("trk_resPhiFitted", &m_resPhiFitted);
0341 m_outputTree->Branch("trk_resThetaFitted", &m_resThetaFitted);
0342 m_outputTree->Branch("trk_resQOverPFitted", &m_resQOverPFitted);
0343 m_outputTree->Branch("trk_momOverlapFitted", &m_momOverlapFitted);
0344
0345 m_outputTree->Branch("trk_pullPhi", &m_pullPhi);
0346 m_outputTree->Branch("trk_pullTheta", &m_pullTheta);
0347 m_outputTree->Branch("trk_pullQOverP", &m_pullQOverP);
0348
0349 m_outputTree->Branch("trk_pullPhiFitted", &m_pullPhiFitted);
0350 m_outputTree->Branch("trk_pullThetaFitted", &m_pullThetaFitted);
0351 m_outputTree->Branch("trk_pullQOverPFitted", &m_pullQOverPFitted);
0352 }
0353 }
0354
0355 VertexNTupleWriter::~VertexNTupleWriter() {
0356 if (m_outputFile != nullptr) {
0357 m_outputFile->Close();
0358 }
0359 }
0360
0361 ProcessCode VertexNTupleWriter::finalize() {
0362 m_outputFile->cd();
0363 m_outputTree->Write();
0364 m_outputFile->Close();
0365
0366 return ProcessCode::SUCCESS;
0367 }
0368
0369 ProcessCode VertexNTupleWriter::writeT(
0370 const AlgorithmContext& ctx, const std::vector<Acts::Vertex>& vertices) {
0371
0372
0373 const static TrackParticleMatching emptyTrackParticleMatching;
0374 const static Acts::ConstVectorTrackContainer emptyConstTrackContainer;
0375 const static Acts::ConstVectorMultiTrajectory emptyConstTrackStateContainer;
0376 const static ConstTrackContainer emptyTracks(
0377 std::make_shared<Acts::ConstVectorTrackContainer>(
0378 emptyConstTrackContainer),
0379 std::make_shared<Acts::ConstVectorMultiTrajectory>(
0380 emptyConstTrackStateContainer));
0381
0382
0383 const SimVertexContainer& truthVertices = m_inputTruthVertices(ctx);
0384
0385 const SimParticleContainer& particles = m_inputParticles(ctx);
0386 const SimParticleContainer& selectedParticles = m_inputSelectedParticles(ctx);
0387 const TrackParticleMatching& trackParticleMatching =
0388 (m_inputTrackParticleMatching.isInitialized()
0389 ? m_inputTrackParticleMatching(ctx)
0390 : emptyTrackParticleMatching);
0391 const ConstTrackContainer& tracks =
0392 (m_inputTracks.isInitialized() ? m_inputTracks(ctx) : emptyTracks);
0393 SimParticleContainer recoParticles;
0394
0395 for (ConstTrackProxy track : tracks) {
0396 if (!track.hasReferenceSurface()) {
0397 ACTS_DEBUG("No reference surface on this track, index = "
0398 << track.index() << " tip index = " << track.tipIndex());
0399 continue;
0400 }
0401
0402 if (const SimParticle* particle =
0403 findParticle(particles, trackParticleMatching, track, logger());
0404 particle != nullptr) {
0405 recoParticles.insert(*particle);
0406 }
0407 }
0408 if (tracks.size() == 0) {
0409
0410
0411 recoParticles = particles;
0412 }
0413
0414
0415 std::lock_guard<std::mutex> lock(m_writeMutex);
0416
0417 m_nRecoVtx = vertices.size();
0418 m_nCleanVtx = 0;
0419 m_nMergedVtx = 0;
0420 m_nSplitVtx = 0;
0421
0422 ACTS_DEBUG("Number of reco vertices in event: " << m_nRecoVtx);
0423
0424
0425 m_nTrueVtx = getNumberOfTruePriVertices(particles);
0426
0427 m_nVtxDetAcceptance = getNumberOfTruePriVertices(selectedParticles);
0428
0429 ACTS_DEBUG("Number of truth particles in event : " << particles.size());
0430 ACTS_DEBUG("Number of truth primary vertices : " << m_nTrueVtx);
0431 ACTS_DEBUG("Number of detector-accepted truth primary vertices : "
0432 << m_nVtxDetAcceptance);
0433
0434
0435 m_eventNr = ctx.eventNumber;
0436
0437
0438 m_nVtxReconstructable = getNumberOfReconstructableVertices(recoParticles);
0439
0440 ACTS_DEBUG("Number of reconstructed tracks : " << tracks.size());
0441 ACTS_DEBUG("Number of reco track-associated truth particles in event : "
0442 << recoParticles.size());
0443 ACTS_DEBUG("Maximum number of reconstructible primary vertices : "
0444 << m_nVtxReconstructable);
0445
0446 struct ToTruthMatching {
0447 std::optional<SimVertexBarcode> vertexId;
0448 double totalTrackWeight{};
0449 double truthMajorityTrackWeights{};
0450 double matchFraction{};
0451
0452 RecoVertexClassification classification{RecoVertexClassification::Unknown};
0453 };
0454 struct ToRecoMatching {
0455 std::size_t recoIndex{};
0456
0457 double recoSumPt2{};
0458 };
0459
0460 std::vector<ToTruthMatching> recoToTruthMatching;
0461 std::map<SimVertexBarcode, ToRecoMatching> truthToRecoMatching;
0462
0463
0464 for (const auto& [vtxIndex, vtx] : Acts::enumerate(vertices)) {
0465
0466 const std::vector<Acts::TrackAtVertex>& tracksAtVtx = vtx.tracks();
0467
0468
0469
0470 std::vector<std::pair<SimVertexBarcode, double>> contributingTruthVertices;
0471
0472 double totalTrackWeight = 0;
0473 if (!tracksAtVtx.empty()) {
0474 for (const Acts::TrackAtVertex& trk : tracksAtVtx) {
0475 if (trk.trackWeight < m_cfg.minTrkWeight) {
0476 continue;
0477 }
0478
0479 totalTrackWeight += trk.trackWeight;
0480
0481 std::optional<ConstTrackProxy> trackOpt = findTrack(tracks, trk);
0482 if (!trackOpt.has_value()) {
0483 ACTS_DEBUG("Track has no matching input track.");
0484 continue;
0485 }
0486 const ConstTrackProxy& inputTrk = *trackOpt;
0487 const SimParticle* particle =
0488 findParticle(particles, trackParticleMatching, inputTrk, logger());
0489 if (particle == nullptr) {
0490 ACTS_VERBOSE("Track has no matching truth particle.");
0491 } else {
0492 contributingTruthVertices.emplace_back(
0493 SimBarcode{particle->particleId()}.vertexId(), trk.trackWeight);
0494 }
0495 }
0496 } else {
0497
0498 for (auto& particle : recoParticles) {
0499 contributingTruthVertices.emplace_back(
0500 SimBarcode{particle.particleId()}.vertexId(), 1.);
0501 }
0502 }
0503
0504
0505 std::map<SimVertexBarcode, std::pair<int, double>> fmap;
0506 for (const auto& [vtxId, weight] : contributingTruthVertices) {
0507 ++fmap[vtxId].first;
0508 fmap[vtxId].second += weight;
0509 }
0510 double truthMajorityVertexTrackWeights = 0;
0511 std::optional<SimVertexBarcode> truthMajorityVertexId = std::nullopt;
0512 for (const auto& [vtxId, counter] : fmap) {
0513 if (counter.second > truthMajorityVertexTrackWeights) {
0514 truthMajorityVertexId = vtxId;
0515 truthMajorityVertexTrackWeights = counter.second;
0516 }
0517 }
0518
0519 double sumPt2 = calcSumPt2(m_cfg, vtx);
0520
0521 double vertexMatchFraction =
0522 (totalTrackWeight > 0
0523 ? truthMajorityVertexTrackWeights / totalTrackWeight
0524 : 0.);
0525 RecoVertexClassification recoVertexClassification =
0526 RecoVertexClassification::Unknown;
0527
0528 if (vertexMatchFraction >= m_cfg.vertexMatchThreshold) {
0529 recoVertexClassification = RecoVertexClassification::Clean;
0530 } else {
0531 recoVertexClassification = RecoVertexClassification::Merged;
0532 }
0533
0534 recoToTruthMatching.push_back({truthMajorityVertexId, totalTrackWeight,
0535 truthMajorityVertexTrackWeights,
0536 vertexMatchFraction,
0537 recoVertexClassification});
0538
0539 auto& recoToTruth = recoToTruthMatching.back();
0540
0541
0542 if (!truthMajorityVertexId.has_value()) {
0543
0544 ACTS_DEBUG("No truth vertex matched to this reconstructed vertex.");
0545 continue;
0546 }
0547
0548 if (auto it = truthToRecoMatching.find(truthMajorityVertexId.value());
0549 it != truthToRecoMatching.end()) {
0550
0551
0552
0553
0554
0555 if (sumPt2 <= it->second.recoSumPt2) {
0556
0557
0558 recoToTruth.classification = RecoVertexClassification::Split;
0559
0560
0561 } else {
0562
0563
0564 auto& otherRecoToTruth = recoToTruthMatching.at(it->second.recoIndex);
0565
0566 recoToTruth.classification = otherRecoToTruth.classification;
0567 otherRecoToTruth.classification = RecoVertexClassification::Split;
0568
0569
0570 it->second = {vtxIndex, sumPt2};
0571 }
0572 } else {
0573 truthToRecoMatching[truthMajorityVertexId.value()] = {vtxIndex, sumPt2};
0574 }
0575 }
0576
0577
0578
0579 for (const auto& [vtxIndex, vtx] : Acts::enumerate(vertices)) {
0580
0581 const auto& tracksAtVtx = vtx.tracks();
0582
0583 std::vector<std::uint32_t> trackIndices;
0584
0585 const auto& toTruthMatching = recoToTruthMatching[vtxIndex];
0586
0587 m_recoX.push_back(vtx.fullPosition()[Acts::CoordinateIndices::eX]);
0588 m_recoY.push_back(vtx.fullPosition()[Acts::CoordinateIndices::eY]);
0589 m_recoZ.push_back(vtx.fullPosition()[Acts::CoordinateIndices::eZ]);
0590 m_recoT.push_back(vtx.fullPosition()[Acts::CoordinateIndices::eTime]);
0591
0592 double varX = vtx.fullCovariance()(Acts::CoordinateIndices::eX,
0593 Acts::CoordinateIndices::eX);
0594 double varY = vtx.fullCovariance()(Acts::CoordinateIndices::eY,
0595 Acts::CoordinateIndices::eY);
0596 double varZ = vtx.fullCovariance()(Acts::CoordinateIndices::eZ,
0597 Acts::CoordinateIndices::eZ);
0598 double varTime = vtx.fullCovariance()(Acts::CoordinateIndices::eTime,
0599 Acts::CoordinateIndices::eTime);
0600
0601 m_covXX.push_back(varX);
0602 m_covYY.push_back(varY);
0603 m_covZZ.push_back(varZ);
0604 m_covTT.push_back(varTime);
0605 m_covXY.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eX,
0606 Acts::CoordinateIndices::eY));
0607 m_covXZ.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eX,
0608 Acts::CoordinateIndices::eZ));
0609 m_covXT.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eX,
0610 Acts::CoordinateIndices::eTime));
0611 m_covYZ.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eY,
0612 Acts::CoordinateIndices::eZ));
0613 m_covYT.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eY,
0614 Acts::CoordinateIndices::eTime));
0615 m_covZT.push_back(vtx.fullCovariance()(Acts::CoordinateIndices::eZ,
0616 Acts::CoordinateIndices::eTime));
0617
0618 double sumPt2 = calcSumPt2(m_cfg, vtx);
0619 m_sumPt2.push_back(sumPt2);
0620
0621 double recoVertexTrackWeights = 0;
0622 for (const Acts::TrackAtVertex& trk : tracksAtVtx) {
0623 if (trk.trackWeight < m_cfg.minTrkWeight) {
0624 continue;
0625 }
0626 recoVertexTrackWeights += trk.trackWeight;
0627 }
0628 m_recoVertexTrackWeights.push_back(recoVertexTrackWeights);
0629
0630 unsigned int nTracksOnRecoVertex = std::count_if(
0631 tracksAtVtx.begin(), tracksAtVtx.end(), [this](const auto& trkAtVtx) {
0632 return trkAtVtx.trackWeight > m_cfg.minTrkWeight;
0633 });
0634 m_nTracksOnRecoVertex.push_back(nTracksOnRecoVertex);
0635
0636
0637 bool truthInfoWritten = false;
0638 std::optional<Acts::Vector4> truthPos;
0639 if (toTruthMatching.vertexId.has_value()) {
0640 auto iTruthVertex = truthVertices.find(toTruthMatching.vertexId.value());
0641 if (iTruthVertex == truthVertices.end()) {
0642 ACTS_ERROR("Truth vertex not found for id: "
0643 << toTruthMatching.vertexId.value());
0644 continue;
0645 }
0646 const SimVertex& truthVertex = *iTruthVertex;
0647
0648 m_vertexPrimary.push_back(truthVertex.vertexId().vertexPrimary());
0649 m_vertexSecondary.push_back(truthVertex.vertexId().vertexSecondary());
0650
0651
0652 int nTracksOnTruthVertex = 0;
0653 for (const auto& particle : selectedParticles) {
0654 if (static_cast<SimVertexBarcode>(particle.particleId().vertexId()) ==
0655 truthVertex.vertexId()) {
0656 ++nTracksOnTruthVertex;
0657 }
0658 }
0659 m_nTracksOnTruthVertex.push_back(nTracksOnTruthVertex);
0660
0661 double truthPrimaryVertexDensity =
0662 calculateTruthPrimaryVertexDensity(m_cfg, truthVertices, vtx);
0663 m_truthPrimaryVertexDensity.push_back(truthPrimaryVertexDensity);
0664
0665 double truthVertexTrackWeights =
0666 toTruthMatching.truthMajorityTrackWeights;
0667 m_truthVertexTrackWeights.push_back(truthVertexTrackWeights);
0668
0669 double truthVertexMatchRatio = toTruthMatching.matchFraction;
0670 m_truthVertexMatchRatio.push_back(truthVertexMatchRatio);
0671
0672 double recoVertexContamination = 1 - truthVertexMatchRatio;
0673 m_recoVertexContamination.push_back(recoVertexContamination);
0674
0675 RecoVertexClassification recoVertexClassification =
0676 toTruthMatching.classification;
0677 m_recoVertexClassification.push_back(
0678 static_cast<int>(recoVertexClassification));
0679
0680 if (recoVertexClassification == RecoVertexClassification::Clean) {
0681 ++m_nCleanVtx;
0682 } else if (recoVertexClassification == RecoVertexClassification::Merged) {
0683 ++m_nMergedVtx;
0684 } else if (recoVertexClassification == RecoVertexClassification::Split) {
0685 ++m_nSplitVtx;
0686 }
0687
0688 const Acts::Vector4& truePos = truthVertex.position4;
0689 truthPos = truePos;
0690 m_truthX.push_back(truePos[Acts::CoordinateIndices::eX]);
0691 m_truthY.push_back(truePos[Acts::CoordinateIndices::eY]);
0692 m_truthZ.push_back(truePos[Acts::CoordinateIndices::eZ]);
0693 m_truthT.push_back(truePos[Acts::CoordinateIndices::eTime]);
0694
0695 const Acts::Vector4 diffPos = vtx.fullPosition() - truePos;
0696 m_resX.push_back(diffPos[Acts::CoordinateIndices::eX]);
0697 m_resY.push_back(diffPos[Acts::CoordinateIndices::eY]);
0698 m_resZ.push_back(diffPos[Acts::CoordinateIndices::eZ]);
0699 m_resT.push_back(diffPos[Acts::CoordinateIndices::eTime]);
0700
0701 m_pullX.push_back(pull(diffPos[Acts::CoordinateIndices::eX], varX, "X",
0702 true, logger()));
0703 m_pullY.push_back(pull(diffPos[Acts::CoordinateIndices::eY], varY, "Y",
0704 true, logger()));
0705 m_pullZ.push_back(pull(diffPos[Acts::CoordinateIndices::eZ], varZ, "Z",
0706 true, logger()));
0707 m_pullT.push_back(pull(diffPos[Acts::CoordinateIndices::eTime], varTime,
0708 "T", true, logger()));
0709
0710 truthInfoWritten = true;
0711 }
0712 if (!truthInfoWritten) {
0713 m_vertexPrimary.push_back(-1);
0714 m_vertexSecondary.push_back(-1);
0715
0716 m_nTracksOnTruthVertex.push_back(-1);
0717
0718 m_truthPrimaryVertexDensity.push_back(nan);
0719
0720 m_truthVertexTrackWeights.push_back(nan);
0721 m_truthVertexMatchRatio.push_back(nan);
0722 m_recoVertexContamination.push_back(nan);
0723
0724 m_recoVertexClassification.push_back(
0725 static_cast<int>(RecoVertexClassification::Unknown));
0726
0727 m_truthX.push_back(nan);
0728 m_truthY.push_back(nan);
0729 m_truthZ.push_back(nan);
0730 m_truthT.push_back(nan);
0731
0732 m_resX.push_back(nan);
0733 m_resY.push_back(nan);
0734 m_resZ.push_back(nan);
0735 m_resT.push_back(nan);
0736
0737 m_pullX.push_back(nan);
0738 m_pullY.push_back(nan);
0739 m_pullZ.push_back(nan);
0740 m_pullT.push_back(nan);
0741 }
0742
0743 if (m_cfg.writeTrackInfo) {
0744 writeTrackInfo(ctx, particles, tracks, trackParticleMatching, truthPos,
0745 tracksAtVtx);
0746 }
0747 }
0748
0749
0750 m_outputTree->Fill();
0751
0752 m_nTracksOnRecoVertex.clear();
0753 m_recoVertexTrackWeights.clear();
0754 m_recoX.clear();
0755 m_recoY.clear();
0756 m_recoZ.clear();
0757 m_recoT.clear();
0758 m_covXX.clear();
0759 m_covYY.clear();
0760 m_covZZ.clear();
0761 m_covTT.clear();
0762 m_covXY.clear();
0763 m_covXZ.clear();
0764 m_covXT.clear();
0765 m_covYZ.clear();
0766 m_covYT.clear();
0767 m_covZT.clear();
0768 m_seedX.clear();
0769 m_seedY.clear();
0770 m_seedZ.clear();
0771 m_seedT.clear();
0772 m_vertexPrimary.clear();
0773 m_vertexSecondary.clear();
0774 m_nTracksOnTruthVertex.clear();
0775 m_truthPrimaryVertexDensity.clear();
0776 m_truthVertexTrackWeights.clear();
0777 m_truthVertexMatchRatio.clear();
0778 m_recoVertexContamination.clear();
0779 m_recoVertexClassification.clear();
0780 m_truthX.clear();
0781 m_truthY.clear();
0782 m_truthZ.clear();
0783 m_truthT.clear();
0784 m_resX.clear();
0785 m_resY.clear();
0786 m_resZ.clear();
0787 m_resT.clear();
0788 m_resSeedZ.clear();
0789 m_resSeedT.clear();
0790 m_pullX.clear();
0791 m_pullY.clear();
0792 m_pullZ.clear();
0793 m_pullT.clear();
0794 m_sumPt2.clear();
0795 m_trkWeight.clear();
0796 m_recoPhi.clear();
0797 m_recoTheta.clear();
0798 m_recoQOverP.clear();
0799 m_recoPhiFitted.clear();
0800 m_recoThetaFitted.clear();
0801 m_recoQOverPFitted.clear();
0802 m_trkParticleId.clear();
0803 m_truthPhi.clear();
0804 m_truthTheta.clear();
0805 m_truthQOverP.clear();
0806 m_resPhi.clear();
0807 m_resTheta.clear();
0808 m_resQOverP.clear();
0809 m_momOverlap.clear();
0810 m_resPhiFitted.clear();
0811 m_resThetaFitted.clear();
0812 m_resQOverPFitted.clear();
0813 m_momOverlapFitted.clear();
0814 m_pullPhi.clear();
0815 m_pullTheta.clear();
0816 m_pullQOverP.clear();
0817 m_pullPhiFitted.clear();
0818 m_pullThetaFitted.clear();
0819 m_pullQOverPFitted.clear();
0820
0821 return ProcessCode::SUCCESS;
0822 }
0823
0824 void VertexNTupleWriter::writeTrackInfo(
0825 const AlgorithmContext& ctx, const SimParticleContainer& particles,
0826 const ConstTrackContainer& tracks,
0827 const TrackParticleMatching& trackParticleMatching,
0828 const std::optional<Acts::Vector4>& truthPos,
0829 const std::vector<Acts::TrackAtVertex>& tracksAtVtx) {
0830
0831
0832
0833 Acts::SympyStepper stepper(m_cfg.bField);
0834 using Propagator = Acts::Propagator<Acts::SympyStepper>;
0835 auto propagator = std::make_shared<Propagator>(stepper);
0836
0837
0838
0839
0840
0841
0842
0843 auto& innerTrkWeight = m_trkWeight.emplace_back();
0844
0845 auto& innerRecoPhi = m_recoPhi.emplace_back();
0846 auto& innerRecoTheta = m_recoTheta.emplace_back();
0847 auto& innerRecoQOverP = m_recoQOverP.emplace_back();
0848
0849 auto& innerRecoPhiFitted = m_recoPhiFitted.emplace_back();
0850 auto& innerRecoThetaFitted = m_recoThetaFitted.emplace_back();
0851 auto& innerRecoQOverPFitted = m_recoQOverPFitted.emplace_back();
0852
0853 auto& innerTrkParticleId = m_trkParticleId.emplace_back();
0854
0855 auto& innerTruthPhi = m_truthPhi.emplace_back();
0856 auto& innerTruthTheta = m_truthTheta.emplace_back();
0857 auto& innerTruthQOverP = m_truthQOverP.emplace_back();
0858
0859 auto& innerResPhi = m_resPhi.emplace_back();
0860 auto& innerResTheta = m_resTheta.emplace_back();
0861 auto& innerResQOverP = m_resQOverP.emplace_back();
0862
0863 auto& innerResPhiFitted = m_resPhiFitted.emplace_back();
0864 auto& innerResThetaFitted = m_resThetaFitted.emplace_back();
0865 auto& innerResQOverPFitted = m_resQOverPFitted.emplace_back();
0866
0867 auto& innerMomOverlap = m_momOverlap.emplace_back();
0868 auto& innerMomOverlapFitted = m_momOverlapFitted.emplace_back();
0869
0870 auto& innerPullPhi = m_pullPhi.emplace_back();
0871 auto& innerPullTheta = m_pullTheta.emplace_back();
0872 auto& innerPullQOverP = m_pullQOverP.emplace_back();
0873
0874 auto& innerPullPhiFitted = m_pullPhiFitted.emplace_back();
0875 auto& innerPullThetaFitted = m_pullThetaFitted.emplace_back();
0876 auto& innerPullQOverPFitted = m_pullQOverPFitted.emplace_back();
0877
0878
0879 std::shared_ptr<Acts::PerigeeSurface> perigeeSurface;
0880 if (truthPos.has_value()) {
0881 perigeeSurface =
0882 Acts::Surface::makeShared<Acts::PerigeeSurface>(truthPos->head<3>());
0883 }
0884
0885 auto propagateToVtx =
0886 [&](const auto& params) -> std::optional<Acts::BoundTrackParameters> {
0887 if (!perigeeSurface) {
0888 return std::nullopt;
0889 }
0890
0891 Acts::Intersection3D intersection =
0892 perigeeSurface
0893 ->intersect(ctx.geoContext, params.position(ctx.geoContext),
0894 params.direction(), Acts::BoundaryTolerance::Infinite())
0895 .closest();
0896
0897
0898 using PropagatorOptions = Propagator::Options<>;
0899 PropagatorOptions pOptions(ctx.geoContext, ctx.magFieldContext);
0900 pOptions.direction =
0901 Acts::Direction::fromScalarZeroAsPositive(intersection.pathLength());
0902
0903 auto result = propagator->propagate(params, *perigeeSurface, pOptions);
0904 if (!result.ok()) {
0905 ACTS_ERROR("Propagation to true vertex position failed.");
0906 return std::nullopt;
0907 }
0908 auto& paramsAtVtx = *result->endParameters;
0909 return std::make_optional(paramsAtVtx);
0910 };
0911
0912 for (const Acts::TrackAtVertex& trk : tracksAtVtx) {
0913 if (trk.trackWeight < m_cfg.minTrkWeight) {
0914 continue;
0915 }
0916
0917 innerTrkWeight.push_back(trk.trackWeight);
0918
0919 Acts::Vector3 trueUnitDir = Acts::Vector3::Zero();
0920 Acts::Vector3 trueMom = Acts::Vector3::Zero();
0921
0922 const SimParticle* particle = nullptr;
0923 std::optional<ConstTrackProxy> trackOpt = findTrack(tracks, trk);
0924 if (trackOpt.has_value()) {
0925 particle =
0926 findParticle(particles, trackParticleMatching, *trackOpt, logger());
0927 }
0928
0929 if (particle != nullptr) {
0930 innerTrkParticleId.push_back(particle->particleId().asVector());
0931
0932 trueUnitDir = particle->direction();
0933 trueMom.head<2>() = Acts::makePhiThetaFromDirection(trueUnitDir);
0934 trueMom[2] = particle->qOverP();
0935
0936 innerTruthPhi.push_back(trueMom[0]);
0937 innerTruthTheta.push_back(trueMom[1]);
0938 innerTruthQOverP.push_back(trueMom[2]);
0939 } else {
0940 ACTS_VERBOSE("Track has no matching truth particle.");
0941
0942 innerTrkParticleId.push_back(ActsFatras::Barcode::Invalid().asVector());
0943
0944 innerTruthPhi.push_back(nan);
0945 innerTruthTheta.push_back(nan);
0946 innerTruthQOverP.push_back(nan);
0947 }
0948
0949
0950 const auto paramsAtVtx =
0951 propagateToVtx(*(trk.originalParams.as<Acts::BoundTrackParameters>()));
0952 if (paramsAtVtx.has_value()) {
0953 Acts::Vector3 recoMom =
0954 paramsAtVtx->parameters().segment(Acts::eBoundPhi, 3);
0955 const Acts::SquareMatrix3& momCov =
0956 paramsAtVtx->covariance()->template block<3, 3>(Acts::eBoundPhi,
0957 Acts::eBoundPhi);
0958 innerRecoPhi.push_back(recoMom[0]);
0959 innerRecoTheta.push_back(recoMom[1]);
0960 innerRecoQOverP.push_back(recoMom[2]);
0961
0962 if (particle != nullptr) {
0963 Acts::Vector3 diffMom = recoMom - trueMom;
0964
0965
0966 diffMom[0] = Acts::detail::difference_periodic(recoMom(0), trueMom(0),
0967 2 * std::numbers::pi);
0968 innerResPhi.push_back(diffMom[0]);
0969 innerResTheta.push_back(diffMom[1]);
0970 innerResQOverP.push_back(diffMom[2]);
0971
0972 innerPullPhi.push_back(
0973 pull(diffMom[0], momCov(0, 0), "phi", false, logger()));
0974 innerPullTheta.push_back(
0975 pull(diffMom[1], momCov(1, 1), "theta", false, logger()));
0976 innerPullQOverP.push_back(
0977 pull(diffMom[2], momCov(2, 2), "q/p", false, logger()));
0978
0979 const auto& recoUnitDir = paramsAtVtx->direction();
0980 double overlap = trueUnitDir.dot(recoUnitDir);
0981 innerMomOverlap.push_back(overlap);
0982 } else {
0983 innerResPhi.push_back(nan);
0984 innerResTheta.push_back(nan);
0985 innerResQOverP.push_back(nan);
0986 innerPullPhi.push_back(nan);
0987 innerPullTheta.push_back(nan);
0988 innerPullQOverP.push_back(nan);
0989 innerMomOverlap.push_back(nan);
0990 }
0991 } else {
0992 innerRecoPhi.push_back(nan);
0993 innerRecoTheta.push_back(nan);
0994 innerRecoQOverP.push_back(nan);
0995 innerResPhi.push_back(nan);
0996 innerResTheta.push_back(nan);
0997 innerResQOverP.push_back(nan);
0998 innerPullPhi.push_back(nan);
0999 innerPullTheta.push_back(nan);
1000 innerPullQOverP.push_back(nan);
1001 innerMomOverlap.push_back(nan);
1002 }
1003
1004
1005 const auto paramsAtVtxFitted = propagateToVtx(trk.fittedParams);
1006 if (paramsAtVtxFitted.has_value()) {
1007 Acts::Vector3 recoMomFitted =
1008 paramsAtVtxFitted->parameters().segment(Acts::eBoundPhi, 3);
1009 const Acts::SquareMatrix3& momCovFitted =
1010 paramsAtVtxFitted->covariance()->block<3, 3>(Acts::eBoundPhi,
1011 Acts::eBoundPhi);
1012 innerRecoPhiFitted.push_back(recoMomFitted[0]);
1013 innerRecoThetaFitted.push_back(recoMomFitted[1]);
1014 innerRecoQOverPFitted.push_back(recoMomFitted[2]);
1015
1016 if (particle != nullptr) {
1017 Acts::Vector3 diffMomFitted = recoMomFitted - trueMom;
1018
1019
1020 diffMomFitted[0] = Acts::detail::difference_periodic(
1021 recoMomFitted(0), trueMom(0), 2 * std::numbers::pi);
1022 innerResPhiFitted.push_back(diffMomFitted[0]);
1023 innerResThetaFitted.push_back(diffMomFitted[1]);
1024 innerResQOverPFitted.push_back(diffMomFitted[2]);
1025
1026 innerPullPhiFitted.push_back(
1027 pull(diffMomFitted[0], momCovFitted(0, 0), "phi", true, logger()));
1028 innerPullThetaFitted.push_back(pull(
1029 diffMomFitted[1], momCovFitted(1, 1), "theta", true, logger()));
1030 innerPullQOverPFitted.push_back(
1031 pull(diffMomFitted[2], momCovFitted(2, 2), "q/p", true, logger()));
1032
1033 const auto& recoUnitDirFitted = paramsAtVtxFitted->direction();
1034 double overlapFitted = trueUnitDir.dot(recoUnitDirFitted);
1035 innerMomOverlapFitted.push_back(overlapFitted);
1036 } else {
1037 innerResPhiFitted.push_back(nan);
1038 innerResThetaFitted.push_back(nan);
1039 innerResQOverPFitted.push_back(nan);
1040 innerPullPhiFitted.push_back(nan);
1041 innerPullThetaFitted.push_back(nan);
1042 innerPullQOverPFitted.push_back(nan);
1043 innerMomOverlapFitted.push_back(nan);
1044 }
1045 } else {
1046 innerRecoPhiFitted.push_back(nan);
1047 innerRecoThetaFitted.push_back(nan);
1048 innerRecoQOverPFitted.push_back(nan);
1049 innerResPhiFitted.push_back(nan);
1050 innerResThetaFitted.push_back(nan);
1051 innerResQOverPFitted.push_back(nan);
1052 innerPullPhiFitted.push_back(nan);
1053 innerPullThetaFitted.push_back(nan);
1054 innerPullQOverPFitted.push_back(nan);
1055 innerMomOverlapFitted.push_back(nan);
1056 }
1057 }
1058 }
1059
1060 }