Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
79 changes: 54 additions & 25 deletions src/algorithms/fardetectors/FarDetectorLinearTracking.cc
Original file line number Diff line number Diff line change
Expand Up @@ -74,20 +74,22 @@ void FarDetectorLinearTracking::process(const FarDetectorLinearTracking::Input&
return;
}

const bool has_assoc_output = trackLinks != nullptr || assocTracks != nullptr;
// Check if truth associations are possible
const bool do_assoc = hitLinks != nullptr && !hitLinks->empty();
const bool do_assoc = has_assoc_output && ((hitLinks != nullptr && !hitLinks->empty()) ||
(assocHits != nullptr && !assocHits->empty()));
if (!do_assoc) {
debug("Provided MCRecoTrackerHitLink collection is empty. No truth associations "
debug("Truth association inputs or outputs are missing. No truth associations "
"will be performed.");
Comment on lines +82 to 83
}
// Build fast lookup once per event using podio::LinkNavigator
std::optional<podio::LinkNavigator<edm4eic::MCRecoTrackerHitLinkCollection>> link_nav;
if (do_assoc) {
if (hitLinks != nullptr && !hitLinks->empty()) {
link_nav.emplace(*hitLinks);
}
Comment on lines 85 to 89

std::vector<std::vector<Eigen::Vector3d>> convertedHits;
std::vector<std::vector<edm4hep::MCParticle>> assocParts;
std::vector<std::vector<std::optional<edm4hep::MCParticle>>> assocParts;
convertedHits.reserve(m_cfg.n_layer);
assocParts.reserve(m_cfg.n_layer);

Expand All @@ -103,7 +105,8 @@ void FarDetectorLinearTracking::process(const FarDetectorLinearTracking::Input&
trace("No hits in layer");
return;
}
ConvertClusters(*layerHits, *link_nav, *assocHits, convertedHits, assocParts);
ConvertClusters(*layerHits, link_nav ? &(*link_nav) : nullptr, assocHits, convertedHits,
assocParts);
}

// Create a matrix to store the hit positions
Expand All @@ -129,7 +132,7 @@ void FarDetectorLinearTracking::process(const FarDetectorLinearTracking::Input&
if (layer == static_cast<long>(m_cfg.n_layer) - 1) {
// Check the combination, if chi2 limit is passed, add the track to the output
checkHitCombination(&hitMatrix, outputTracks, trackLinks, assocTracks, inputhits,
assocParts, layerHitIndex);
assocParts, layerHitIndex, do_assoc);
} else {
layer++;
continue;
Expand Down Expand Up @@ -162,8 +165,8 @@ void FarDetectorLinearTracking::checkHitCombination(
edm4eic::MCRecoTrackParticleLinkCollection* trackLinks,
edm4eic::MCRecoTrackParticleAssociationCollection* assocTracks,
const std::vector<gsl::not_null<const edm4eic::Measurement2DCollection*>>& inputHits,
const std::vector<std::vector<edm4hep::MCParticle>>& assocParts,
const std::vector<std::size_t>& layerHitIndex) const {
const std::vector<std::vector<std::optional<edm4hep::MCParticle>>>& assocParts,
const std::vector<std::size_t>& layerHitIndex, bool do_assoc) const {

Eigen::Vector3d weightedAnchor = (*hitMatrix) * m_layerWeights / (m_layerWeights.sum());

Expand Down Expand Up @@ -210,20 +213,35 @@ void FarDetectorLinearTracking::checkHitCombination(
std::unordered_map<edm4hep::MCParticle, int> particleCount;
for (std::size_t layer = 0; layer < layerHitIndex.size(); layer++) {
track.addToMeasurements((*inputHits[layer])[layerHitIndex[layer]]);
if (!do_assoc) {
continue;
}
if (layer >= assocParts.size() || layerHitIndex[layer] >= assocParts[layer].size()) {
continue;
}
const auto& assocParticle = assocParts[layer][layerHitIndex[layer]];
particleCount[assocParticle]++;
if (assocParticle.has_value()) {
particleCount[*assocParticle]++;
}
}

// Create track associations for each particle
if (!do_assoc) {
return;
}
for (const auto& [particle, count] : particleCount) {
auto trackLink = trackLinks->create();
trackLink.setFrom(track);
trackLink.setTo(particle);
trackLink.setWeight(count / static_cast<double>(m_cfg.n_layer));
auto trackAssoc = assocTracks->create();
trackAssoc.setRec(track);
trackAssoc.setSim(particle);
trackAssoc.setWeight(count / static_cast<double>(m_cfg.n_layer));
if (trackLinks != nullptr) {
auto trackLink = trackLinks->create();
trackLink.setFrom(track);
trackLink.setTo(particle);
trackLink.setWeight(count / static_cast<double>(m_cfg.n_layer));
}
if (assocTracks != nullptr) {
auto trackAssoc = assocTracks->create();
trackAssoc.setRec(track);
trackAssoc.setSim(particle);
trackAssoc.setWeight(count / static_cast<double>(m_cfg.n_layer));
}
}
}

Expand All @@ -247,17 +265,17 @@ bool FarDetectorLinearTracking::checkHitPair(const Eigen::Vector3d& hit1,
// Convert measurements into global coordinates
void FarDetectorLinearTracking::ConvertClusters(
const edm4eic::Measurement2DCollection& clusters,
const podio::LinkNavigator<edm4eic::MCRecoTrackerHitLinkCollection>& link_nav,
[[maybe_unused]] const edm4eic::MCRecoTrackerHitAssociationCollection& assoc_hits,
const podio::LinkNavigator<edm4eic::MCRecoTrackerHitLinkCollection>* link_nav,
const edm4eic::MCRecoTrackerHitAssociationCollection* assoc_hits,
std::vector<std::vector<Eigen::Vector3d>>& pointPositions,
std::vector<std::vector<edm4hep::MCParticle>>& assoc_parts) const {
std::vector<std::vector<std::optional<edm4hep::MCParticle>>>& assoc_parts) const {

// Get context of first hit
const dd4hep::VolumeManagerContext* context =
m_cellid_converter->findContext(clusters[0].getSurface());

std::vector<Eigen::Vector3d> layerPositions;
std::vector<edm4hep::MCParticle> assocParticles;
std::vector<std::optional<edm4hep::MCParticle>> assocParticles;

for (auto cluster : clusters) {

Expand All @@ -277,17 +295,28 @@ void FarDetectorLinearTracking::ConvertClusters(
}
if (maxIndex == cluster.getWeights().size()) {
// no maximum found (e.g. all weights zero, cluster size zero)
assocParticles.emplace_back(std::nullopt);
continue;
}
auto maxHit = cluster.getHits()[maxIndex];
// Get associated raw hit
auto rawHit = maxHit.getRawHit();

const auto sim_hits = link_nav.getLinked(rawHit);
if (!sim_hits.empty()) {
auto particle = sim_hits[0].o.getParticle();
assocParticles.push_back(particle);
std::optional<edm4hep::MCParticle> assocParticle = std::nullopt;
if (link_nav != nullptr) {
const auto sim_hits = link_nav->getLinked(rawHit);
if (!sim_hits.empty()) {
assocParticle = sim_hits[0].o.getParticle();
}
} else if (assoc_hits != nullptr) {
for (const auto& assoc : *assoc_hits) {
if (assoc.getRawHit() == rawHit) {
assocParticle = assoc.getSimHit().getParticle();
break;
}
}
}
assocParticles.push_back(assocParticle);
}

pointPositions.push_back(layerPositions);
Expand Down
10 changes: 5 additions & 5 deletions src/algorithms/fardetectors/FarDetectorLinearTracking.h
Original file line number Diff line number Diff line change
Expand Up @@ -63,19 +63,19 @@ class FarDetectorLinearTracking : public FarDetectorLinearTrackingAlgorithm,
edm4eic::MCRecoTrackParticleLinkCollection* trackLinks,
edm4eic::MCRecoTrackParticleAssociationCollection* assocTracks,
const std::vector<gsl::not_null<const edm4eic::Measurement2DCollection*>>& inputHits,
const std::vector<std::vector<edm4hep::MCParticle>>& assocParts,
const std::vector<std::size_t>& layerHitIndex) const;
const std::vector<std::vector<std::optional<edm4hep::MCParticle>>>& assocParts,
const std::vector<std::size_t>& layerHitIndex, bool do_assoc) const;

/** Check if the last two hits are within a certain angle of the optimum direction **/
bool checkHitPair(const Eigen::Vector3d& hit1, const Eigen::Vector3d& hit2) const;

/** Convert 2D clusters to 3D coordinates and match associated particle **/
void
ConvertClusters(const edm4eic::Measurement2DCollection& clusters,
const podio::LinkNavigator<edm4eic::MCRecoTrackerHitLinkCollection>& link_nav,
const edm4eic::MCRecoTrackerHitAssociationCollection& assoc_hits,
const podio::LinkNavigator<edm4eic::MCRecoTrackerHitLinkCollection>* link_nav,
const edm4eic::MCRecoTrackerHitAssociationCollection* assoc_hits,
std::vector<std::vector<Eigen::Vector3d>>& pointPositions,
std::vector<std::vector<edm4hep::MCParticle>>& assoc_parts) const;
std::vector<std::vector<std::optional<edm4hep::MCParticle>>>& assoc_parts) const;
};

} // namespace eicrecon
55 changes: 33 additions & 22 deletions src/algorithms/tracking/ActsToTracks.cc
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,9 @@ void ActsToTracks::init() {}
void ActsToTracks::process(const Input& input, const Output& output) const {
const auto [meas2Ds, track_seeds, acts_track_states, acts_tracks, raw_hit_assocs] = input;
auto [trajectories, track_parameters, tracks, tracks_links, tracks_assoc] = output;
const bool do_track_links = tracks_links != nullptr;
const bool do_track_assoc = tracks_assoc != nullptr;
const bool do_raw_assoc = raw_hit_assocs != nullptr && !raw_hit_assocs->empty();
Comment on lines +67 to +69

// Create accessor for seed number dynamic column
Acts::ConstProxyAccessor<unsigned int> seedNumber("seed");
Expand Down Expand Up @@ -233,13 +236,15 @@ void ActsToTracks::process(const Input& input, const Output& output) const {
// Determine track associations if hit associations provided
// FIXME: not able to check whether optional inputs were provided
//if (raw_hit_assocs->has_value()) {
for (const auto& hit : meas2D.getHits()) {
auto raw_hit = hit.getRawHit();
for (const auto raw_hit_assoc : *raw_hit_assocs) {
if (raw_hit_assoc.getRawHit() == raw_hit) {
auto sim_hit = raw_hit_assoc.getSimHit();
auto mc_particle = sim_hit.getParticle();
mcparticle_weight_by_hit_count[mc_particle]++;
if (do_raw_assoc) {
Comment on lines 236 to +239
for (const auto& hit : meas2D.getHits()) {
auto raw_hit = hit.getRawHit();
for (const auto raw_hit_assoc : *raw_hit_assocs) {
if (raw_hit_assoc.getRawHit() == raw_hit) {
auto sim_hit = raw_hit_assoc.getSimHit();
auto mc_particle = sim_hit.getParticle();
mcparticle_weight_by_hit_count[mc_particle]++;
}
}
}
}
Expand All @@ -252,21 +257,27 @@ void ActsToTracks::process(const Input& input, const Output& output) const {
// Store track associations if hit associations provided
// FIXME: not able to check whether optional inputs were provided
//if (raw_hit_assocs->has_value()) {
double total_weight = std::accumulate(
mcparticle_weight_by_hit_count.begin(), mcparticle_weight_by_hit_count.end(), 0,
[](const double sum, const auto& i) { return sum + i.second; });
for (const auto& [mcparticle, weight] : mcparticle_weight_by_hit_count) {
double normalized_weight = weight / total_weight;
auto track_link = tracks_links->create();
track_link.setFrom(track_out);
track_link.setTo(mcparticle);
track_link.setWeight(normalized_weight);
auto track_assoc = tracks_assoc->create();
track_assoc.setRec(track_out);
track_assoc.setSim(mcparticle);
track_assoc.setWeight(normalized_weight);
debug("track {}: mcparticle {} weight {}", track_out.id().index, mcparticle.id().index,
normalized_weight);
if ((do_track_assoc || do_track_links) && !mcparticle_weight_by_hit_count.empty()) {
double total_weight = std::accumulate(
mcparticle_weight_by_hit_count.begin(), mcparticle_weight_by_hit_count.end(), 0.0,
[](const double sum, const auto& i) { return sum + i.second; });
for (const auto& [mcparticle, weight] : mcparticle_weight_by_hit_count) {
double normalized_weight = weight / total_weight;
if (do_track_links) {
auto track_link = tracks_links->create();
track_link.setFrom(track_out);
track_link.setTo(mcparticle);
track_link.setWeight(normalized_weight);
}
if (do_track_assoc) {
auto track_assoc = tracks_assoc->create();
track_assoc.setRec(track_out);
track_assoc.setSim(mcparticle);
track_assoc.setWeight(normalized_weight);
}
debug("track {}: mcparticle {} weight {}", track_out.id().index, mcparticle.id().index,
normalized_weight);
}
}
//}
}
Expand Down
Loading