Skip to content

Commit

Permalink
refactor: simplify RootVertexPerformanceWriter (#1567)
Browse files Browse the repository at this point in the history 
- deduplicate code
- improve doc
- make track selector optional

cc @timadye @paulgessinger
  • Loading branch information
@andiwand
andiwand committed Oct 4, 2022
1 parent 2c40516 commit 9b793b2
Show file tree
Hide file tree
Showing 3 changed files with 106 additions and 195 deletions.
8 changes: 5 additions & 3 deletions Examples/Io/Root/include/ActsExamples/Io/Root/RootVertexPerformanceWriter.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,11 +39,13 @@ class RootVertexPerformanceWriter final
std::string inputSelectedTruthParticles;
/// Truth particles associated to fitted tracks
std::string inputAssociatedTruthParticles;
/// All event fitted tracks
/// Selected fitted tracks
std::string inputFittedTracks;
/// All event fitted tracks indices
/// Selected fitted tracks indices (points from `inputSelectedFittedTracks`
/// to `inputAllFittedTracks`). If empty all we assume selected == all.
std::string inputFittedTracksIndices;
/// All event fitted tracks tips
/// All event fitted tracks tips (points from `inputAllFittedTracks` to
/// `inputTrajectories`)
std::string inputAllFittedTracksTips;
/// Trajectories object from track finidng
std::string inputTrajectories;
Expand Down
290 changes: 99 additions & 191 deletions Examples/Io/Root/src/RootVertexPerformanceWriter.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,8 +53,7 @@ ActsExamples::RootVertexPerformanceWriter::RootVertexPerformanceWriter(
"Collection with selected truth particles missing");
}
if (m_cfg.inputAssociatedTruthParticles.empty() &&
(m_cfg.inputFittedTracksIndices.empty() ||
m_cfg.inputAllFittedTracksTips.empty() ||
(m_cfg.inputAllFittedTracksTips.empty() ||
m_cfg.inputTrajectories.empty())) {
throw std::invalid_argument(
"You need to either provide collection of truth particles matching 1:1 "
Expand Down Expand Up @@ -209,21 +208,12 @@ ActsExamples::ProcessCode ActsExamples::RootVertexPerformanceWriter::writeT(
ACTS_VERBOSE(
"Total number of reconstructed tracks : " << inputFittedTracks.size());

SimParticleContainer associatedTruthParticles;

if (!m_cfg.inputAssociatedTruthParticles.empty()) {
// Read track-associated truth particle input collection
const auto& associatedTruthParticles =
ctx.eventStore.get<SimParticleContainer>(
m_cfg.inputAssociatedTruthParticles);
// Get number of track-associated true primary vertices
m_nVtxReconstructable =
getNumberOfReconstructableVertices(associatedTruthParticles);

ACTS_VERBOSE(
"Total number of reco track-associated truth particles in event : "
<< associatedTruthParticles.size());
ACTS_VERBOSE(
"Total number of reco track-associated truth primary vertices : "
<< m_nVtxReconstructable);
associatedTruthParticles = ctx.eventStore.get<SimParticleContainer>(
m_cfg.inputAssociatedTruthParticles);

/***************** Start x,y,z resolution plots here *****************/
// Matching tracks at vertex to fitted tracks that are in turn matched
Expand All @@ -236,9 +226,7 @@ ActsExamples::ProcessCode ActsExamples::RootVertexPerformanceWriter::writeT(
<< inputFittedTracks.size()
<< " != " << associatedTruthParticles.size()
<< ") Not able to match fitted tracks at reconstructed "
"vertex to "
"truth "
"vertex."
"vertex to truth vertex."
<< extra);
};

Expand All @@ -247,87 +235,6 @@ ActsExamples::ProcessCode ActsExamples::RootVertexPerformanceWriter::writeT(
} else if (associatedTruthParticles.size() > inputFittedTracks.size()) {
mismatchMsg(Acts::Logging::INFO,
" This is likely due to track efficiency < 1");
} else {
// Loop over all reco vertices and find associated truth particles
std::vector<SimParticleContainer> truthParticlesAtVtxContainer;
for (const auto& vtx : vertices) {
const auto tracks = vtx.tracks();
// Store all associated truth particles to current vtx
SimParticleContainer particleAtVtx;

std::vector<int> contributingTruthVertices;

for (const auto& trk : tracks) {
Acts::BoundTrackParameters origTrack = *(trk.originalParams);

// Find associated truth particle now
int idx = 0;
for (const auto& particle : associatedTruthParticles) {
if (origTrack.parameters() == inputFittedTracks[idx].parameters()) {
particleAtVtx.insert(particleAtVtx.end(), particle);

int priVtxId = particle.particleId().vertexPrimary();
contributingTruthVertices.push_back(priVtxId);
}
idx++;
}
} // end loop tracks

// Now find true vtx with most matching tracks at reco vtx
// and check if it contributes more than 50% of all tracks
std::map<int, int> fmap;
for (int priVtxId : contributingTruthVertices) {
fmap[priVtxId]++;
}
int maxOccurrenceId = -1;
int maxOccurrence = -1;
for (auto it : fmap) {
if (it.second > maxOccurrence) {
maxOccurrence = it.second;
maxOccurrenceId = it.first;
}
}

// Match reco to truth vertex if at least 50% of tracks match
double trackVtxMatchFraction =
(double)fmap[maxOccurrenceId] / tracks.size();
if (trackVtxMatchFraction > m_cfg.minTrackVtxMatchFraction) {
for (const auto& particle : associatedTruthParticles) {
int priVtxId = particle.particleId().vertexPrimary();
int secVtxId = particle.particleId().vertexSecondary();

if (secVtxId != 0) {
// truthparticle from secondary vtx
continue;
}

if (priVtxId == maxOccurrenceId) {
// Vertex found, fill varibles
const auto& truePos = particle.position();

m_diffx.push_back(vtx.position()[0] - truePos[0]);
m_diffy.push_back(vtx.position()[1] - truePos[1]);
m_diffz.push_back(vtx.position()[2] - truePos[2]);

m_truthX.push_back(truePos[0]);
m_truthY.push_back(truePos[1]);
m_truthZ.push_back(truePos[2]);

m_recoX.push_back(vtx.position()[0]);
m_recoY.push_back(vtx.position()[1]);
m_recoZ.push_back(vtx.position()[2]);

m_covXX.push_back(vtx.covariance()(0, 0));
m_covYY.push_back(vtx.covariance()(1, 1));
m_covXY.push_back(vtx.covariance()(0, 1));
m_covYX.push_back(vtx.covariance()(1, 0));
m_trackVtxMatchFraction.push_back(trackVtxMatchFraction);
// Next vertex now
break;
}
}
}
} // end loop vertices
}
} else {
// get active tips
Expand All @@ -336,23 +243,25 @@ ActsExamples::ProcessCode ActsExamples::RootVertexPerformanceWriter::writeT(
const auto& allTracksTips =
ctx.eventStore.get<std::vector<std::pair<size_t, size_t>>>(
m_cfg.inputAllFittedTracksTips);
const auto& trackIndices = ctx.eventStore.get<std::vector<uint32_t>>(
m_cfg.inputFittedTracksIndices);
const std::vector<uint32_t>* trackIndices = nullptr;

if (ctx.eventStore.exists(m_cfg.inputFittedTracksIndices)) {
trackIndices = &ctx.eventStore.get<std::vector<uint32_t>>(
m_cfg.inputFittedTracksIndices);

throw_assert(
trackIndices.size() == inputFittedTracks.size(),
"Selected track indices count does not match fitted tracks count");
throw_assert(
trackIndices->size() == inputFittedTracks.size(),
"Selected track indices count does not match fitted tracks count");
}

std::vector<ParticleHitCount> particleHitCounts;

using HitParticlesMap = IndexMultimap<ActsFatras::Barcode>;
const auto& hitParticlesMap =
ctx.eventStore.get<HitParticlesMap>(m_cfg.inputMeasurementParticlesMap);

SimParticleContainer associatedTruthParticles;

for (size_t i = 0; i < inputFittedTracks.size(); i++) {
auto fittedTrackIndex = trackIndices[i];
auto fittedTrackIndex = trackIndices != nullptr ? (*trackIndices)[i] : i;
auto& [iTraj, tip] = allTracksTips[fittedTrackIndex];
const auto& traj = trajectories[iTraj];
identifyContributingParticles(hitParticlesMap, traj, tip,
Expand Down Expand Up @@ -382,98 +291,97 @@ ActsExamples::ProcessCode ActsExamples::RootVertexPerformanceWriter::writeT(
associatedTruthParticles.emplace_hint(associatedTruthParticles.end(),
majorityParticle);
}
}

// Get number of track-associated true primary vertices
m_nVtxReconstructable =
getNumberOfReconstructableVertices(associatedTruthParticles);

ACTS_INFO("Total number of reco track-associated truth particles in event : "
<< associatedTruthParticles.size());
ACTS_INFO("Total number of reco track-associated truth primary vertices : "
<< m_nVtxReconstructable);

// Loop over all reco vertices and find associated truth particles
std::vector<SimParticleContainer> truthParticlesAtVtxContainer;
for (const auto& vtx : vertices) {
const auto tracks = vtx.tracks();
// Store all associated truth particles to current vtx
SimParticleContainer particleAtVtx;

// Get number of track-associated true primary vertices
m_nVtxReconstructable =
getNumberOfReconstructableVertices(associatedTruthParticles);

ACTS_INFO(
"Total number of reco track-associated truth particles in event : "
<< associatedTruthParticles.size());
ACTS_INFO("Total number of reco track-associated truth primary vertices : "
<< m_nVtxReconstructable);

// Loop over all reco vertices and find associated truth particles
std::vector<SimParticleContainer> truthParticlesAtVtxContainer;
for (const auto& vtx : vertices) {
const auto tracks = vtx.tracks();
// Store all associated truth particles to current vtx
SimParticleContainer particleAtVtx;

std::vector<int> contributingTruthVertices;

for (const auto& trk : tracks) {
Acts::BoundTrackParameters origTrack = *(trk.originalParams);

// Find associated truth particle now
int idx = 0;
for (const auto& particle : associatedTruthParticles) {
if (origTrack.parameters() == inputFittedTracks[idx].parameters()) {
particleAtVtx.insert(particleAtVtx.end(), particle);

int priVtxId = particle.particleId().vertexPrimary();
contributingTruthVertices.push_back(priVtxId);
}
idx++;
std::vector<int> contributingTruthVertices;

for (const auto& trk : tracks) {
Acts::BoundTrackParameters origTrack = *(trk.originalParams);

// Find associated truth particle now
int idx = 0;
for (const auto& particle : associatedTruthParticles) {
if (origTrack.parameters() == inputFittedTracks[idx].parameters()) {
particleAtVtx.insert(particleAtVtx.end(), particle);

int priVtxId = particle.particleId().vertexPrimary();
contributingTruthVertices.push_back(priVtxId);
}
} // end loop tracks
idx++;
}
} // end loop tracks

// Now find true vtx with most matching tracks at reco vtx
// and check if it contributes more than 50 of all tracks
std::map<int, int> fmap;
for (int priVtxId : contributingTruthVertices) {
fmap[priVtxId]++;
// Now find true vtx with most matching tracks at reco vtx
// and check if it contributes more than 50 of all tracks
std::map<int, int> fmap;
for (int priVtxId : contributingTruthVertices) {
fmap[priVtxId]++;
}
int maxOccurrenceId = -1;
int maxOccurence = -1;
for (auto it : fmap) {
if (it.second > maxOccurence) {
maxOccurence = it.second;
maxOccurrenceId = it.first;
}
int maxOccurrenceId = -1;
int maxOccurence = -1;
for (auto it : fmap) {
if (it.second > maxOccurence) {
maxOccurence = it.second;
maxOccurrenceId = it.first;
}

// Match reco to truth vertex if at least 50% of tracks match
double trackVtxMatchFraction =
(double)fmap[maxOccurrenceId] / tracks.size();
if (trackVtxMatchFraction > m_cfg.minTrackVtxMatchFraction) {
for (const auto& particle : associatedTruthParticles) {
int priVtxId = particle.particleId().vertexPrimary();
int secVtxId = particle.particleId().vertexSecondary();

if (secVtxId != 0) {
// truthparticle from secondary vtx
continue;
}
}

// Match reco to truth vertex if at least 50% of tracks match
double trackVtxMatchFraction =
(double)fmap[maxOccurrenceId] / tracks.size();
if (trackVtxMatchFraction > m_cfg.minTrackVtxMatchFraction) {
for (const auto& particle : associatedTruthParticles) {
int priVtxId = particle.particleId().vertexPrimary();
int secVtxId = particle.particleId().vertexSecondary();

if (secVtxId != 0) {
// truthparticle from secondary vtx
continue;
}

if (priVtxId == maxOccurrenceId) {
// Vertex found, fill varibles
const auto& truePos = particle.position();

m_diffx.push_back(vtx.position()[0] - truePos[0]);
m_diffy.push_back(vtx.position()[1] - truePos[1]);
m_diffz.push_back(vtx.position()[2] - truePos[2]);

m_truthX.push_back(truePos[0]);
m_truthY.push_back(truePos[1]);
m_truthZ.push_back(truePos[2]);

m_recoX.push_back(vtx.position()[0]);
m_recoY.push_back(vtx.position()[1]);
m_recoZ.push_back(vtx.position()[2]);

m_covXX.push_back(vtx.covariance()(0, 0));
m_covYY.push_back(vtx.covariance()(1, 1));
m_covXY.push_back(vtx.covariance()(0, 1));
m_covYX.push_back(vtx.covariance()(1, 0));
m_trackVtxMatchFraction.push_back(trackVtxMatchFraction);
// Next vertex now
break;
}
if (priVtxId == maxOccurrenceId) {
// Vertex found, fill varibles
const auto& truePos = particle.position();

m_diffx.push_back(vtx.position()[0] - truePos[0]);
m_diffy.push_back(vtx.position()[1] - truePos[1]);
m_diffz.push_back(vtx.position()[2] - truePos[2]);

m_truthX.push_back(truePos[0]);
m_truthY.push_back(truePos[1]);
m_truthZ.push_back(truePos[2]);

m_recoX.push_back(vtx.position()[0]);
m_recoY.push_back(vtx.position()[1]);
m_recoZ.push_back(vtx.position()[2]);

m_covXX.push_back(vtx.covariance()(0, 0));
m_covYY.push_back(vtx.covariance()(1, 1));
m_covXY.push_back(vtx.covariance()(0, 1));
m_covYX.push_back(vtx.covariance()(1, 0));
m_trackVtxMatchFraction.push_back(trackVtxMatchFraction);
// Next vertex now
break;
}
}
} // end loop vertices
}
}
} // end loop vertices

// Retrieve and set reconstruction time
if (!m_cfg.inputTime.empty()) {
Expand Down
3 changes: 2 additions & 1 deletion Examples/Python/python/acts/examples/reconstruction.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -908,7 +908,6 @@ def addVertexFitting(
associatedParticles: str = "particles_input",
trajectories: Optional[str] = None,
trackParameters: str = "fittedTrackParameters",
trackIndices: str = "outputTrackIndices",
vertexFinder: VertexFinder = VertexFinder.Truth,
trackSelectorRanges: TrackSelectorRanges = TrackSelectorRanges(),
logLevel: Optional[acts.logging.Level] = None,
Expand Down Expand Up @@ -944,6 +943,7 @@ def addVertexFitting(

if trackSelectorRanges is not None:
selectedTrackParameters = "trackparameters"
trackIndices = "outputTrackIndices"

trackSelector = acts.examples.TrackSelector(
level=customLogLevel(),
Expand Down Expand Up @@ -972,6 +972,7 @@ def addVertexFitting(
s.addAlgorithm(trackSelector)
else:
selectedTrackParameters = trackParameters
trackIndices = ""

inputParticles = "particles_input"
outputVertices = "fittedVertices"
Expand Down

0 comments on commit 9b793b2

Please sign in to comment.