Skip to content

Commit

Permalink
Merge pull request #220 from eic/217-save-associations
Browse files Browse the repository at this point in the history 
217 save associations
  • Loading branch information
@DraTeots
DraTeots authored Oct 15, 2022
2 parents d4483c2 + c59bd1f commit 8930fe9
Show file tree
Hide file tree
Showing 21 changed files with 497 additions and 87 deletions.
147 changes: 103 additions & 44 deletions src/algorithms/tracking/ParticlesWithTruthPID.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,31 +1,33 @@
// Original licence: SPDX-License-Identifier: LGPL-3.0-or-later
// Copyright (C) 2022 Sylvester Joosten, Wouter Deconinck, Dmitry Romanov

#include <algorithm>
#include <cmath>
#include "ParticlesWithTruthPID.h"

#include <fmt/format.h>
// Event Model related classes
#include <edm4hep/MutableMCParticle.h>
#include <edm4eic/ReconstructedParticle.h>
#include <edm4eic/MutableReconstructedParticle.h>
#include <edm4eic/MutableMCRecoParticleAssociation.h>
#include <edm4eic/vector_utils.h>

#include <spdlog/spdlog.h>
// Getting access to PODIO indexes
#include <extensions/podio_access/accessor.h>
#include <edm4eic/ReconstructedParticleObj.h>

// Event Model related classes
#include "edm4hep/MCParticleCollection.h"
#include "edm4eic/MCRecoParticleAssociationCollection.h"
#include "edm4eic/ReconstructedParticleCollection.h"
#include "edm4eic/TrackParametersCollection.h"
#include "edm4eic/vector_utils.h"
using ReconstructedParticleObjPtr = edm4eic::ReconstructedParticleObj*;
ALLOW_ACCESS(edm4eic::MutableReconstructedParticle, m_obj, ReconstructedParticleObjPtr);

#include "ParticlesWithTruthPIDConfig.h"
#include "ParticlesWithTruthPID.h"
#include "algorithms/reco/ParticlesWithAssociation.h"
#include <algorithms/interfaces/WithPodConfig.h>

namespace eicrecon {

void ParticlesWithTruthPID::init(std::shared_ptr<spdlog::logger> logger) {
m_log = logger;

ParticlesWithAssociation * ParticlesWithTruthPID::execute(
}

ParticlesWithAssociation *ParticlesWithTruthPID::process(
const std::vector<const edm4hep::MCParticle *> &mc_particles,
const std::vector<const edm4eic::TrackParameters*> &track_params) {
const std::vector<const edm4eic::TrackParameters *> &track_params) {
// input collection

/// Resulting reconstructed particles
Expand All @@ -34,23 +36,55 @@ namespace eicrecon {
/// Resulting associations
std::vector<edm4eic::MCRecoParticleAssociation *> associations;


const double sinPhiOver2Tolerance = sin(0.5 * m_cfg.phiTolerance);
std::vector<bool> consumed(mc_particles.size(), false);
tracePhiToleranceOnce(sinPhiOver2Tolerance, m_cfg.phiTolerance);

std::vector<bool> mc_prt_is_consumed(mc_particles.size(), false); // MCParticle is already consumed flag

for (const auto &trk: track_params) {
const auto mom = edm4eic::sphericalToVector(1.0 / std::abs(trk->getQOverP()), trk->getTheta(), trk->getPhi());
const auto mom = edm4eic::sphericalToVector(1.0 / std::abs(trk->getQOverP()), trk->getTheta(),
trk->getPhi());
const auto charge_rec = trk->getCharge();


m_log->debug("Match: [id] [mom] [theta] [phi] [charge] [PID]");
m_log->debug(" Track : {:<4} {:<8.3f} {:<8.3f} {:<8.2f} {:<4}",
trk->getObjectID().index, edm4eic::magnitude(mom), edm4eic::anglePolar(mom), edm4eic::angleAzimuthal(mom), charge_rec);

// utility variables for matching
int best_match = -1;
double best_delta = std::numeric_limits<double>::max();
for (size_t ip = 0; ip < mc_particles.size(); ++ip) {
const auto &mc_part = mc_particles[ip];
if (consumed[ip] || mc_part->getGeneratorStatus() > 1 || mc_part->getCharge() == 0 ||
mc_part->getCharge() * charge_rec < 0) {
m_log->debug("ignoring non-primary/neutral/opposite charge particle");
const auto &p = mc_part->getMomentum();

m_log->trace(" MCParticle with id={:<4} mom={:<8.3f} charge={}", mc_part->getObjectID().index,
edm4eic::magnitude(p), mc_part->getCharge());

// Check if used
if (mc_prt_is_consumed[ip]) {
m_log->trace(" Ignoring. Particle is already used");
continue;
}
const auto &p = mc_part->getMomentum();

// Check if non-primary
if (mc_part->getGeneratorStatus() > 1) {
m_log->trace(" Ignoring. GeneratorStatus > 1 => Non-primary particle");
continue;
}

// Check if neutral
if (mc_part->getCharge() == 0) {
m_log->trace(" Ignoring. Neutral particle");
continue;
}

// Check opposite charge
if (mc_part->getCharge() * charge_rec < 0) {
m_log->trace(" Ignoring. Opposite charge particle");
continue;
}

const auto p_mag = std::hypot(p.x, p.y, p.z);
const auto p_phi = std::atan2(p.y, p.x);
const auto p_eta = std::atanh(p.z / p_mag);
Expand All @@ -60,13 +94,24 @@ namespace eicrecon {
const double dsphi = std::abs(sin(0.5 * (edm4eic::angleAzimuthal(mom) - p_phi)));
const double deta = std::abs((edm4eic::eta(mom) - p_eta));

if (dp_rel < m_cfg.pRelativeTolerance && deta < m_cfg.etaTolerance && dsphi < sinPhiOver2Tolerance) {
bool is_matching = dp_rel < m_cfg.momentumRelativeTolerance &&
deta < m_cfg.etaTolerance &&
dsphi < sinPhiOver2Tolerance;

m_log->trace(" Decision: {} dp: {:.4f} < {} && d_eta: {:.6f} < {} && d_sin_phi: {:.4e} < {:.4e} ",
is_matching? "Matching":"Ignoring",
dp_rel, m_cfg.momentumRelativeTolerance,
deta, m_cfg.etaTolerance,
dsphi, sinPhiOver2Tolerance);

if (is_matching) {
const double delta =
std::hypot(dp_rel / m_cfg.pRelativeTolerance, deta / m_cfg.etaTolerance,
std::hypot(dp_rel / m_cfg.momentumRelativeTolerance, deta / m_cfg.etaTolerance,
dsphi / sinPhiOver2Tolerance);
if (delta < best_delta) {
best_match = ip;
best_delta = delta;
m_log->trace(" Is the best match now");
}
}
}
Expand All @@ -76,17 +121,26 @@ namespace eicrecon {
// float time = 0;
float mass = 0;
if (best_match >= 0) {
consumed[best_match] = true;
m_log->trace("Best match is found and is: {}", best_match);
mc_prt_is_consumed[best_match] = true;
const auto &best_mc_part = mc_particles[best_match];
best_pid = best_mc_part->getPDG();
referencePoint = {
static_cast<float>(best_mc_part->getVertex().x), static_cast<float>(best_mc_part->getVertex().y),
static_cast<float>(best_mc_part->getVertex().x),
static_cast<float>(best_mc_part->getVertex().y),
static_cast<float>(best_mc_part->getVertex().z)}; // @TODO: not sure if vertex/reference point makes sense here
// time = mcpart.getTime();
mass = best_mc_part->getMass();
}

// Getting access to private PODIO member to set the index in the collection
// Evil, Evil, Evil, Evil is here TODO change asap
// Nothing is private in this world anymore!
auto obj = ACCESS(rec_part, m_obj);
obj->id.index = (int) reco_particles.size(); // Set the index in the collection

rec_part.setType(static_cast<int16_t>(best_match >= 0 ? 0 : -1)); // @TODO: determine type codes
rec_part.setEnergy(std::hypot(edm4eic::magnitude(mom), mass));
rec_part.setEnergy((float) std::hypot(edm4eic::magnitude(mom), mass));
rec_part.setMomentum(mom);
rec_part.setReferencePoint(referencePoint);
rec_part.setCharge(charge_rec);
Expand All @@ -101,30 +155,28 @@ namespace eicrecon {
rec_assoc.setSimID(mc_particles[best_match]->getObjectID().index);
rec_assoc.setWeight(1);
rec_assoc.setRec(rec_part);
//rec_assoc.setSim(mc[best_match]);
auto sim = mc_particles[best_match]->clone();
rec_assoc.setSim(sim);
associations.emplace_back(new edm4eic::MCRecoParticleAssociation(rec_assoc));
}
if (m_log->level() <= spdlog::level::debug) {
if (best_match > 0) {

if (m_log->level() <= spdlog::level::debug) {

const auto &mcpart = mc_particles[best_match];
m_log->debug("Matched track with MC particle {}\n", best_match);
m_log->debug(" - Track: (mom: {}, theta: {}, phi: {}, charge: {})",
edm4eic::magnitude(mom),
edm4eic::anglePolar(mom), edm4eic::angleAzimuthal(mom), charge_rec);
const auto &p = mcpart->getMomentum();
const auto p_mag = edm4eic::magnitude(p);
const auto p_phi = edm4eic::angleAzimuthal(p);
const auto p_theta = edm4eic::anglePolar(p);
m_log->debug(" - MC particle: (mom: {}, theta: {}, phi: {}, charge: {}, type: {}",
p_mag, p_theta,
p_phi, mcpart->getCharge(), mcpart->getPDG());
} else {
m_log->debug("Did not find a good match for track \n");
m_log->debug(" - Track: (mom: {}, theta: {}, phi: {}, charge: {})",
edm4eic::magnitude(mom),
edm4eic::anglePolar(mom), edm4eic::angleAzimuthal(mom), charge_rec);
m_log->debug(" MCPart: {:<4} {:<8.3f} {:<8.3f} {:<8.2f} {:<6}",
mcpart->getObjectID().index, p_mag, p_theta, p_phi, mcpart->getCharge(),
mcpart->getPDG());

m_log->debug(" Assoc: id={} SimId={} RecId={}",
rec_assoc.getObjectID().index, rec_assoc.getSimID(), rec_assoc.getSimID());
}
}
else {
m_log->debug(" MCPart: Did not find a good match");
}

// Add particle to the output vector
reco_particles.push_back(new edm4eic::ReconstructedParticle(rec_part));
Expand All @@ -134,4 +186,11 @@ namespace eicrecon {
return new ParticlesWithAssociation(std::move(reco_particles), std::move(associations));
}

void ParticlesWithTruthPID::tracePhiToleranceOnce(const double sinPhiOver2Tolerance, double phiTolerance) {
// This function is called once to print tolerances useful for tracing
static std::once_flag do_it_once;
std::call_once(do_it_once, [this, sinPhiOver2Tolerance, phiTolerance]() {
m_log->trace("m_cfg.phiTolerance: {:<8.4f} => sinPhiOver2Tolerance: {:<8.4f}", sinPhiOver2Tolerance, phiTolerance);
});
}
}
23 changes: 9 additions & 14 deletions src/algorithms/tracking/ParticlesWithTruthPID.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,38 +8,33 @@
#include <algorithm>
#include <cmath>

#include <fmt/format.h>

#include <spdlog/spdlog.h>

// Event Model related classes
#include "edm4hep/MCParticleCollection.h"
#include "edm4eic/MCRecoParticleAssociationCollection.h"
#include "edm4eic/ReconstructedParticleCollection.h"
#include "edm4eic/TrackParametersCollection.h"
#include "edm4eic/vector_utils.h"
#include <edm4eic/TrackParameters.h>
#include <edm4hep/MCParticle.h>

#include "ParticlesWithTruthPIDConfig.h"
#include "algorithms/reco/ParticlesWithAssociation.h"
#include <algorithms/reco/ParticlesWithAssociation.h>
#include <algorithms/interfaces/WithPodConfig.h>
#include "ParticlesWithTruthPIDConfig.h"



namespace eicrecon {
class ParticlesWithTruthPID : public WithPodConfig<ParticlesWithTruthPIDConfig> {

public:

void init(std::shared_ptr<spdlog::logger> logger) {
m_log = logger;
}
void init(std::shared_ptr<spdlog::logger> logger);

ParticlesWithAssociation *execute(
ParticlesWithAssociation *process(
const std::vector<const edm4hep::MCParticle *> &mc_particles,
const std::vector<const edm4eic::TrackParameters*> &track_params);

private:

std::shared_ptr<spdlog::logger> m_log;

void tracePhiToleranceOnce(const double sinPhiOver2Tolerance, double phiTolerance);
};
}

Expand Down
6 changes: 3 additions & 3 deletions src/algorithms/tracking/ParticlesWithTruthPIDConfig.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,11 +8,11 @@
namespace eicrecon {

struct ParticlesWithTruthPIDConfig {
double pRelativeTolerance = 0.1; /// Matching momentum tolerance requires 10% by default;
double momentumRelativeTolerance = 0.1; /// Matching momentum tolerance requires 10% by default;

double phiTolerance = 0.030; /// Matching phi tolerance of 10 mrad
double phiTolerance = 0.030; /// Matching phi tolerance [mrad]

double etaTolerance = 0.2; /// Matching eta tolerance of 0.1
double etaTolerance = 0.2; /// Matching eta tolerance of 0.2
};

} // eicrecon
Expand Down
1 change: 0 additions & 1 deletion src/benchmarks/reconstruction/tracking_efficiency/TrackingEfficiency_processor.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,6 @@
#include <services/rootfile/RootFile_service.h>
#include <extensions/spdlog/SpdlogExtensions.h>


//--------------------------------
// OccupancyAnalysis (Constructor)
//--------------------------------
Expand Down
Loading

0 comments on commit 8930fe9

Please sign in to comment.