From e36b3989457d65296e28b50e59c752b13d7361fc Mon Sep 17 00:00:00 2001
From: Tobias Bockh <tobias.boeckh@cern.ch>
Date: Wed, 23 Nov 2022 19:27:59 +0100
Subject: [PATCH] track truth matching
---
PhysicsAnalysis/NtupleDumper/CMakeLists.txt | 2 +-
.../NtupleDumper/src/NtupleDumperAlg.cxx | 180 +++++++++++++++++-
.../NtupleDumper/src/NtupleDumperAlg.h | 34 ++++
.../Acts/FaserActsKalmanFilter/CMakeLists.txt | 6 +
.../IFiducialParticleTool.h | 28 +++
.../ITrackTruthMatchingTool.h | 21 ++
.../src/FiducialParticleTool.cxx | 91 +++++++++
.../src/FiducialParticleTool.h | 37 ++++
.../src/TrackTruthMatchingTool.cxx | 103 ++++++++++
.../src/TrackTruthMatchingTool.h | 40 ++++
.../FaserActsKalmanFilter_entries.cxx | 4 +
11 files changed, 536 insertions(+), 10 deletions(-)
create mode 100644 Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/IFiducialParticleTool.h
create mode 100644 Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/ITrackTruthMatchingTool.h
create mode 100644 Tracking/Acts/FaserActsKalmanFilter/src/FiducialParticleTool.cxx
create mode 100644 Tracking/Acts/FaserActsKalmanFilter/src/FiducialParticleTool.h
create mode 100644 Tracking/Acts/FaserActsKalmanFilter/src/TrackTruthMatchingTool.cxx
create mode 100644 Tracking/Acts/FaserActsKalmanFilter/src/TrackTruthMatchingTool.h
diff --git a/PhysicsAnalysis/NtupleDumper/CMakeLists.txt b/PhysicsAnalysis/NtupleDumper/CMakeLists.txt
index 944ed968..b81b9818 100644
--- a/PhysicsAnalysis/NtupleDumper/CMakeLists.txt
+++ b/PhysicsAnalysis/NtupleDumper/CMakeLists.txt
@@ -5,7 +5,7 @@ atlas_add_component(
src/NtupleDumperAlg.h
src/NtupleDumperAlg.cxx
src/component/NtupleDumper_entries.cxx
- LINK_LIBRARIES AthenaBaseComps StoreGateLib xAODFaserWaveform xAODFaserTrigger ScintIdentifier FaserCaloIdentifier GeneratorObjects FaserActsGeometryLib TrackerSimEvent TrackerSimData TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives TrackerRIO_OnTrack TrackerSpacePoint
+ LINK_LIBRARIES AthenaBaseComps StoreGateLib xAODFaserWaveform xAODFaserTrigger ScintIdentifier FaserCaloIdentifier GeneratorObjects FaserActsGeometryLib TrackerSimEvent TrackerSimData TrackerIdentifier TrackerReadoutGeometry TrkTrack GeoPrimitives TrackerRIO_OnTrack TrackerSpacePoint FaserActsKalmanFilterLib
)
atlas_install_python_modules(python/*.py)
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
index d1a07f42..587a894b 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
@@ -16,6 +16,9 @@
#include "xAODTruth/TruthParticle.h"
#include <cmath>
#include <TH1F.h>
+#include <numeric>
+
+constexpr float NaN = std::numeric_limits<double>::quiet_NaN();
NtupleDumperAlg::NtupleDumperAlg(const std::string &name,
@@ -70,7 +73,7 @@ void NtupleDumperAlg::FillWaveBranches(const xAOD::WaveformHitContainer &wave) c
}
}
-StatusCode NtupleDumperAlg::initialize()
+StatusCode NtupleDumperAlg::initialize()
{
ATH_CHECK(m_truthEventContainer.initialize());
ATH_CHECK(m_truthParticleContainer.initialize());
@@ -85,6 +88,7 @@ StatusCode NtupleDumperAlg::initialize()
ATH_CHECK(m_simDataCollection.initialize());
ATH_CHECK(m_FaserTriggerData.initialize());
ATH_CHECK(m_ClockWaveformContainer.initialize());
+ ATH_CHECK(m_siHitCollectionKey.initialize());
ATH_CHECK(detStore()->retrieve(m_sctHelper, "FaserSCT_ID"));
ATH_CHECK(detStore()->retrieve(m_vetoNuHelper, "VetoNuID"));
@@ -96,6 +100,8 @@ StatusCode NtupleDumperAlg::initialize()
ATH_CHECK(detStore()->retrieve(m_detMgr, "SCT"));
ATH_CHECK(m_extrapolationTool.retrieve());
ATH_CHECK(m_trackingGeometryTool.retrieve());
+ ATH_CHECK(m_trackTruthMatchingTool.retrieve());
+ ATH_CHECK(m_fiducialParticleTool.retrieve());
ATH_CHECK(m_spacePointContainerKey.initialize());
@@ -220,6 +226,41 @@ StatusCode NtupleDumperAlg::initialize()
m_tree->Branch("Track_ThetaX_atCalo", &m_thetaxCalo);
m_tree->Branch("Track_ThetaY_atCalo", &m_thetayCalo);
+ m_tree->Branch("t_pdg", &m_t_pdg);
+ m_tree->Branch("t_barcode", &m_t_barcode);
+ m_tree->Branch("t_truthHitRatio", &m_t_truthHitRatio);
+ m_tree->Branch("t_prodVtx_x", &m_t_prodVtx_x);
+ m_tree->Branch("t_prodVtx_y", &m_t_prodVtx_y);
+ m_tree->Branch("t_prodVtx_z", &m_t_prodVtx_z);
+ m_tree->Branch("t_decayVtx_x", &m_t_decayVtx_x);
+ m_tree->Branch("t_decayVtx_y", &m_t_decayVtx_y);
+ m_tree->Branch("t_decayVtx_z", &m_t_decayVtx_z);
+ m_tree->Branch("t_px", &m_t_px);
+ m_tree->Branch("t_py", &m_t_py);
+ m_tree->Branch("t_pz", &m_t_pz);
+ m_tree->Branch("t_theta", &m_t_theta);
+ m_tree->Branch("t_phi", &m_t_phi);
+ m_tree->Branch("t_p", &m_t_p);
+ m_tree->Branch("t_pT", &m_t_pT);
+ m_tree->Branch("t_eta", &m_t_eta);
+ m_tree->Branch("t_st0_x", &m_t_st_x[0]);
+ m_tree->Branch("t_st0_y", &m_t_st_y[0]);
+ m_tree->Branch("t_st0_z", &m_t_st_z[0]);
+ m_tree->Branch("t_st1_x", &m_t_st_x[1]);
+ m_tree->Branch("t_st1_y", &m_t_st_y[1]);
+ m_tree->Branch("t_st1_z", &m_t_st_z[1]);
+ m_tree->Branch("t_st2_x", &m_t_st_x[2]);
+ m_tree->Branch("t_st2_y", &m_t_st_y[2]);
+ m_tree->Branch("t_st2_z", &m_t_st_z[2]);
+ m_tree->Branch("t_st3_x", &m_t_st_x[3]);
+ m_tree->Branch("t_st3_y", &m_t_st_y[3]);
+ m_tree->Branch("t_st3_z", &m_t_st_z[3]);
+ m_tree->Branch("isFiducial", &m_isFiducial);
+
+ m_tree->Branch("truthParticleBarcode", &m_truthParticleBarcode);
+ m_tree->Branch("truthParticleMatchedTracks", &m_truthParticleMatchedTracks);
+ m_tree->Branch("truthParticleIsFiducial", &m_truthParticleIsFiducial);
+
m_tree->Branch("pTruthLepton", &m_truthLeptonMomentum, "pTruthLepton/D");
m_tree->Branch("truthBarcode", &m_truthBarcode, "truthBarcode/I");
m_tree->Branch("truthPdg", &m_truthPdg, "truthPdg/I");
@@ -273,7 +314,7 @@ StatusCode NtupleDumperAlg::initialize()
}
-StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
+StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
{
clearTree();
@@ -339,7 +380,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
return StatusCode::SUCCESS; // finished with this event
- } else if ( ((m_tap&8)==0) && (((m_tap&4)==0)||((m_tap&2)==0)) && (((m_tap&4)==0)||((m_tap&1)==0)) && (((m_tap&2)==0)||((m_tap&1)==0)) ) { // don't process events that don't trigger coincidence triggers: 1=calo, 2=veotnu|neto1|preshower, 4=TimingLayer, 8=(VetoNu|Veto2)&Preshower
+ } else if ( ((m_tap&8)==0) && (((m_tap&4)==0)||((m_tap&2)==0)) && (((m_tap&4)==0)||((m_tap&1)==0)) && (((m_tap&2)==0)||((m_tap&1)==0)) ) { // don't process events that don't trigger coincidence triggers: 1=calo, 2=veotnu|neto1|preshower, 4=TimingLayer, 8=(VetoNu|Veto2)&Preshower
return StatusCode::SUCCESS;
}
m_tbp=triggerData->tbp();
@@ -420,7 +461,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
FillWaveBranches(*triggerContainer);
FillWaveBranches(*preshowerContainer);
FillWaveBranches(*ecalContainer);
-
+
m_calo_total=m_wave_charge[0]+m_wave_charge[1]+m_wave_charge[2]+m_wave_charge[3];
m_calo_rawtotal=m_wave_raw_charge[0]+m_wave_raw_charge[1]+m_wave_raw_charge[2]+m_wave_raw_charge[3];
@@ -430,13 +471,13 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_Calo1_Edep = (m_wave_charge[1] / 24.333) * m_MIP_sim_Edep_calo;
m_Calo2_Edep = (m_wave_charge[2] / 24.409) * m_MIP_sim_Edep_calo;
m_Calo3_Edep = (m_wave_charge[3] / 25.555) * m_MIP_sim_Edep_calo;
- } else if (m_CaloConfig == "Low_gain") { // assume low gain calo
+ } else if (m_CaloConfig == "Low_gain") { // assume low gain calo
m_Calo0_Edep = (m_wave_charge[0] / 0.7909) * m_MIP_sim_Edep_calo;
m_Calo1_Edep = (m_wave_charge[1] / 0.8197) * m_MIP_sim_Edep_calo;
m_Calo2_Edep = (m_wave_charge[2] / 0.8256) * m_MIP_sim_Edep_calo;
m_Calo3_Edep = (m_wave_charge[3] / 0.8821) * m_MIP_sim_Edep_calo;
} else {
- ATH_MSG_WARNING("Run config is neither High_gain nor Low_gain, it is " << m_CaloConfig << ", calo calibration will be zero");
+ ATH_MSG_WARNING("Run config is neither High_gain nor Low_gain, it is " << m_CaloConfig << ", calo calibration will be zero");
}
m_Calo_Total_Edep = m_Calo0_Edep + m_Calo1_Edep + m_Calo2_Edep + m_Calo3_Edep;
@@ -517,6 +558,16 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_trackseg_pz.push_back(SegMomentum.z());
}
+ std::map<int, size_t> truthParticleCount {};
+ if (!realData) {
+ SG::ReadHandle<xAOD::TruthParticleContainer> truthParticleContainer { m_truthParticleContainer, ctx };
+ ATH_CHECK(truthParticleContainer.isValid() && truthParticleContainer->size() > 0);
+ for (const xAOD::TruthParticle *tp : *truthParticleContainer) {
+ if (tp->p4().P() > m_minMomentum)
+ truthParticleCount[tp->barcode()] = 0;
+ }
+ }
+
SG::ReadHandle<TrackCollection> trackCollection {m_trackCollection, ctx};
ATH_CHECK(trackCollection.isValid());
const Trk::TrackParameters* candidateParameters {nullptr};
@@ -580,6 +631,56 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_pzdown.push_back(candidateDownParameters->momentum().z());
m_pdown.push_back(sqrt( pow(candidateDownParameters->momentum().x(),2) + pow(candidateDownParameters->momentum().y(),2) + pow(candidateDownParameters->momentum().z(),2) ));
+ if (!realData) {
+ auto [truthParticle, hitCount] = m_trackTruthMatchingTool->getTruthParticle(track);
+ if (truthParticle != nullptr) {
+ if (truthParticleCount.count(truthParticle->barcode()) > 0)
+ truthParticleCount[truthParticle->barcode()] += 1;
+ m_t_pdg.push_back(truthParticle->pdgId());
+ m_t_barcode.push_back(truthParticle->barcode());
+ // the track fit eats up 5 degrees of freedom, thus the number of hits on track is m_DoF + 5
+ m_t_truthHitRatio.push_back(hitCount / (m_DoF.back() + 5));
+ m_isFiducial.push_back(m_fiducialParticleTool->isFiducial(truthParticle->barcode()));
+ auto positions = m_fiducialParticleTool->getTruthPositions(truthParticle->barcode());
+ for (int station = 0; station < 4; ++station) {
+ m_t_st_x[station].push_back(positions[station].x());
+ m_t_st_y[station].push_back(positions[station].y());
+ m_t_st_z[station].push_back(positions[station].z());
+ }
+ if (truthParticle->hasProdVtx()) {
+ m_t_prodVtx_x.push_back(truthParticle->prodVtx()->x());
+ m_t_prodVtx_y.push_back(truthParticle->prodVtx()->y());
+ m_t_prodVtx_z.push_back(truthParticle->prodVtx()->z());
+ } else {
+ m_t_prodVtx_x.push_back(NaN);
+ m_t_prodVtx_y.push_back(NaN);
+ m_t_prodVtx_z.push_back(NaN);
+ }
+ if (truthParticle->hasDecayVtx()) {
+ m_t_decayVtx_x.push_back(truthParticle->decayVtx()->x());
+ m_t_decayVtx_y.push_back(truthParticle->decayVtx()->y());
+ m_t_decayVtx_z.push_back(truthParticle->decayVtx()->z());
+ } else {
+ m_t_decayVtx_x.push_back(NaN);
+ m_t_decayVtx_y.push_back(NaN);
+ m_t_decayVtx_z.push_back(NaN);
+ }
+ m_t_px.push_back(truthParticle->px());
+ m_t_py.push_back(truthParticle->py());
+ m_t_pz.push_back(truthParticle->pz());
+ m_t_theta.push_back(truthParticle->p4().Theta());
+ m_t_phi.push_back(truthParticle->p4().Phi());
+ m_t_p.push_back(truthParticle->p4().P());
+ m_t_pT.push_back(truthParticle->p4().Pt());
+ m_t_eta.push_back(truthParticle->p4().Eta());
+ } else {
+ setNaN();
+ }
+ } else {
+ ATH_MSG_WARNING("Can not find truthParticle.");
+ setNaN();
+ }
+
// fill extrapolation vectors with filler values that get changed iif the track extrapolation succeeds
m_xVetoNu.push_back(-10000);
m_yVetoNu.push_back(-10000);
@@ -676,7 +777,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
ATH_MSG_INFO("Trig null targetParameters");
}
- }
+ }
// extrapolate track from tracking station 3
if (stationMap.count(3) > 0) { // extrapolation crashes if the track does not end in the Station 3, as it is too far away to extrapolate
@@ -735,6 +836,14 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_longTracks++;
}
+ if (!realData) {
+ for (auto &tp : truthParticleCount) {
+ m_truthParticleBarcode.push_back(tp.first);
+ m_truthParticleMatchedTracks.push_back(tp.second);
+ m_truthParticleIsFiducial.push_back(m_fiducialParticleTool->isFiducial(tp.first));
+ }
+ }
+
/*
// Here we apply the signal selection
// Very simple/unrealistic to start
@@ -751,7 +860,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
}
-StatusCode NtupleDumperAlg::finalize()
+StatusCode NtupleDumperAlg::finalize()
{
return StatusCode::SUCCESS;
}
@@ -843,7 +952,7 @@ NtupleDumperAlg::clearTree() const
m_charge.clear();
m_nLayers.clear();
m_longTracks = 0;
-
+
m_nHit0.clear();
m_nHit1.clear();
m_nHit2.clear();
@@ -884,7 +993,60 @@ NtupleDumperAlg::clearTree() const
m_thetaxCalo.clear();
m_thetayCalo.clear();
+ m_t_pdg.clear();
+ m_t_barcode.clear();
+ m_t_truthHitRatio.clear();
+ m_t_prodVtx_x.clear();
+ m_t_prodVtx_y.clear();
+ m_t_prodVtx_z.clear();
+ m_t_decayVtx_x.clear();
+ m_t_decayVtx_y.clear();
+ m_t_decayVtx_z.clear();
+ m_t_px.clear();
+ m_t_py.clear();
+ m_t_pz.clear();
+ m_t_theta.clear();
+ m_t_phi.clear();
+ m_t_p.clear();
+ m_t_pT.clear();
+ m_t_eta.clear();
+ m_isFiducial.clear();
+ for (int station = 0; station < 4; ++station) {
+ m_t_st_x[station].clear();
+ m_t_st_y[station].clear();
+ m_t_st_z[station].clear();
+ }
+ m_truthParticleBarcode.clear();
+ m_truthParticleMatchedTracks.clear();
+ m_truthParticleIsFiducial.clear();
+
m_truthLeptonMomentum = 0;
m_truthBarcode = 0;
m_truthPdg = 0;
}
+
+void NtupleDumperAlg::setNaN() const {
+ m_t_pdg.push_back(0);
+ m_t_barcode.push_back(-1);
+ m_t_truthHitRatio.push_back(NaN);
+ m_t_prodVtx_x.push_back(NaN);
+ m_t_prodVtx_y.push_back(NaN);
+ m_t_prodVtx_z.push_back(NaN);
+ m_t_decayVtx_x.push_back(NaN);
+ m_t_decayVtx_y.push_back(NaN);
+ m_t_decayVtx_z.push_back(NaN);
+ m_t_px.push_back(NaN);
+ m_t_py.push_back(NaN);
+ m_t_pz.push_back(NaN);
+ m_t_theta.push_back(NaN);
+ m_t_phi.push_back(NaN);
+ m_t_p.push_back(NaN);
+ m_t_pT.push_back(NaN);
+ m_t_eta.push_back(NaN);
+ for (int station = 0; station < 4; ++station) {
+ m_t_st_x[station].push_back(NaN);
+ m_t_st_y[station].push_back(NaN);
+ m_t_st_z[station].push_back(NaN);
+ }
+ m_isFiducial.push_back(false);
+}
\ No newline at end of file
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
index 1b26caad..fb805829 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
@@ -15,6 +15,9 @@
#include "TrackerSimData/TrackerSimDataCollection.h"
#include "FaserActsGeometryInterfaces/IFaserActsExtrapolationTool.h"
#include "FaserActsGeometryInterfaces/IFaserActsTrackingGeometryTool.h"
+#include "FaserActsKalmanFilter/IFiducialParticleTool.h"
+#include "FaserActsKalmanFilter/ITrackTruthMatchingTool.h"
+#include "TrackerSimEvent/FaserSiHitCollection.h"
#include <vector>
@@ -44,6 +47,7 @@ private:
bool waveformHitOK(const xAOD::WaveformHit* hit) const;
void clearTree() const;
+ void setNaN() const;
void addBranch(const std::string &name,float* var);
void addBranch(const std::string &name,unsigned int* var);
void addWaveBranches(const std::string &name, int nchannels, int first);
@@ -64,11 +68,14 @@ private:
SG::ReadHandleKey<xAOD::WaveformHitContainer> m_ecalContainer { this, "EcalContainer", "CaloWaveformHits", "Ecal hit container name" };
SG::ReadHandleKey<Tracker::FaserSCT_ClusterContainer> m_clusterContainer { this, "ClusterContainer", "SCT_ClusterContainer", "Tracker cluster container name" };
SG::ReadHandleKey<FaserSCT_SpacePointContainer> m_spacePointContainerKey { this, "SpacePoints", "SCT_SpacePointContainer", "space point container"};
+ SG::ReadHandleKey<FaserSiHitCollection> m_siHitCollectionKey{this, "FaserSiHitCollection", "SCT_Hits"};
SG::ReadHandleKey<xAOD::FaserTriggerData> m_FaserTriggerData { this, "FaserTriggerDataKey", "FaserTriggerData", "ReadHandleKey for xAOD::FaserTriggerData"};
SG::ReadHandleKey<xAOD::WaveformClock> m_ClockWaveformContainer { this, "WaveformClockKey", "WaveformClock", "ReadHandleKey for ClockWaveforms Container"};
ToolHandle<IFaserActsExtrapolationTool> m_extrapolationTool { this, "ExtrapolationTool", "FaserActsExtrapolationTool" };
ToolHandle<IFaserActsTrackingGeometryTool> m_trackingGeometryTool {this, "TrackingGeometryTool", "FaserActsTrackingGeometryTool"};
+ ToolHandle<ITrackTruthMatchingTool> m_trackTruthMatchingTool {this, "TrackTruthMatchingTool", "TrackTruthMatchingTool"};
+ ToolHandle<IFiducialParticleTool> m_fiducialParticleTool {this, "FiducialParticleTool", "FiducialParticleTool"};
const TrackerDD::SCT_DetectorManager* m_detMgr {nullptr};
@@ -86,6 +93,7 @@ private:
IntegerProperty m_flukaCollisions { this, "FlukaCollisions", 137130000, "Number of proton-proton collisions in FLUKA sample." };
DoubleProperty m_flukaCrossSection { this, "FlukaCrossSection", 80.0, "Fluka p-p inelastic cross-section in millibarns." };
DoubleProperty m_genieLuminosity { this, "GenieLuminosity", 150.0, "Genie luminosity in inverse fb." };
+ DoubleProperty m_minMomentum { this, "MinMomentum", 50000.0, "Write out all truth particles with a momentum larger MinMomentum"};
double m_baseEventCrossSection {1.0};
const double kfemtoBarnsPerMilliBarn {1.0e12};
@@ -203,6 +211,32 @@ private:
mutable std::vector<double> m_thetaxCalo;
mutable std::vector<double> m_thetayCalo;
+ mutable std::vector<int> m_t_pdg; // pdg code of the truth matched particle
+ mutable std::vector<int> m_t_barcode; // barcode of the truth matched particle
+ mutable std::vector<double> m_t_truthHitRatio; // ratio of hits on track matched to the truth particle over all hits on track
+ mutable std::vector<double> m_t_prodVtx_x; // x component of the production vertex in mm
+ mutable std::vector<double> m_t_prodVtx_y; // y component of the production vertex in mm
+ mutable std::vector<double> m_t_prodVtx_z; // z component of the production vertex in mm
+ mutable std::vector<double> m_t_decayVtx_x; // x component of the decay vertex in mm
+ mutable std::vector<double> m_t_decayVtx_y; // y component of the decay vertex in mm
+ mutable std::vector<double> m_t_decayVtx_z; // z component of the decay vertex in mm
+ mutable std::vector<double> m_t_px; // truth momentum px in MeV
+ mutable std::vector<double> m_t_py; // truth momentum py in MeV
+ mutable std::vector<double> m_t_pz; // truth momentum pz in MeV
+ mutable std::vector<double> m_t_theta; // angle of truth particle with respsect to the beam axis in rad, theta = arctan(sqrt(px * px + py * py) / pz)
+ mutable std::vector<double> m_t_phi; // polar angle of truth particle in rad, phi = arctan(py / px)
+ mutable std::vector<double> m_t_p; // truth momentum p in MeV
+ mutable std::vector<double> m_t_pT; // transverse truth momentum pT in MeV
+ mutable std::vector<double> m_t_eta; // eta of truth particle
+ mutable std::array<std::vector<double>, 4> m_t_st_x; // vector of the x components of the simulated hits of the truth particle for each station
+ mutable std::array<std::vector<double>, 4> m_t_st_y; // vector of the y components of the simulated hits of the truth particle for each station
+ mutable std::array<std::vector<double>, 4> m_t_st_z; // vector of the z components of the simulated hits of the truth particle for each station
+ mutable std::vector<bool> m_isFiducial; // track is fiducial if there are simulated hits for stations 1 - 3 and the distance from the center is smaller than 100 mm
+
+ mutable std::vector<int> m_truthParticleBarcode; // vector of barcodes of all truth particles with a momentum larger 50 GeV
+ mutable std::vector<int> m_truthParticleMatchedTracks; // vector of number of tracks to which a truth particle is matched to
+ mutable std::vector<bool> m_truthParticleIsFiducial; // vector of boolean showing whether a truth particle is fiducial
+
mutable double m_truthLeptonMomentum;
mutable int m_truthBarcode;
mutable int m_truthPdg;
diff --git a/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt b/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt
index 9377c471..f93c0942 100755
--- a/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt
+++ b/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt
@@ -39,8 +39,10 @@ atlas_add_component(FaserActsKalmanFilter
MyTrackSeedTool.h
FaserActsKalmanFilter/IdentifierLink.h
FaserActsKalmanFilter/IndexSourceLink.h
+ FaserActsKalmanFilter/IFiducialParticleTool.h
FaserActsKalmanFilter/ITrackFinderTool.h
FaserActsKalmanFilter/ITrackSeedTool.h
+ FaserActsKalmanFilter/ITrackTruthMatchingTool.h
KalmanFitterTool.h
LinearFit.h
# ClusterTrackSeedTool.h
@@ -77,6 +79,8 @@ atlas_add_component(FaserActsKalmanFilter
src/CombinatorialKalmanFilterAlg.cxx
src/EffPlotTool.cxx
src/FaserActsKalmanFilterAlg.cxx
+ src/FiducialParticleTool.h
+ src/FiducialParticleTool.cxx
src/GhostBusters.cxx
src/MyTrackSeedTool.cxx
src/KalmanFitterTool.cxx
@@ -101,6 +105,8 @@ atlas_add_component(FaserActsKalmanFilter
src/TrackClassification.cxx
src/TrackSeedWriterTool.cxx
src/TrackSelection.cxx
+ src/TrackTruthMatchingTool.h
+ src/TrackTruthMatchingTool.cxx
# src/TruthTrackFinderTool.cxx
# src/TruthSeededTrackFinderTool.cxx
src/ThreeStationTrackSeedTool.cxx
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/IFiducialParticleTool.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/IFiducialParticleTool.h
new file mode 100644
index 00000000..1df9aeab
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/IFiducialParticleTool.h
@@ -0,0 +1,28 @@
+/*
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS and FASER
+ collaborations
+*/
+
+#ifndef FASERACTSKALMANFILTER_IFIDUCIALEVENTSELECTIONTOOL_H
+#define FASERACTSKALMANFILTER_IFIDUCIALEVENTSELECTIONTOOL_H
+
+#include "CLHEP/Geometry/Point3D.h"
+
+class IFiducialParticleTool : virtual public IAlgTool {
+public:
+ DeclareInterfaceID(IFiducialParticleTool, 1, 0);
+
+ /** Check if truth hits for the given barcode exist and are in the
+ * fiducial volume (r < 100mm) in stations 1 to 3.
+ * @param barcode of a xAOD::TruthParticle
+ */
+ virtual bool isFiducial(int barcode) const = 0;
+
+ /** Return average truth position in each station.
+ * @param barcode of a xAOD::TruthParticle
+ */
+ virtual std::array<HepGeom::Point3D<double>, 4>
+ getTruthPositions(int barcode) const = 0;
+};
+
+#endif // FASERACTSKALMANFILTER_IFIDUCIALEVENTSELECTIONTOOL_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/ITrackTruthMatchingTool.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/ITrackTruthMatchingTool.h
new file mode 100644
index 00000000..045e597a
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/ITrackTruthMatchingTool.h
@@ -0,0 +1,21 @@
+/*
+ Copyright (C) 2002-2022 CERN for the benefit of the ATLAS and FASER
+ collaborations
+*/
+
+#ifndef FASERACTSKALMANFILTER_ITRACKTRUTHMATCHINGTOOL_H
+#define FASERACTSKALMANFILTER_ITRACKTRUTHMATCHINGTOOL_H
+
+#include "GaudiKernel/IAlgTool.h"
+#include "TrkTrack/Track.h"
+#include "xAODTruth/TruthParticle.h"
+
+class ITrackTruthMatchingTool : virtual public IAlgTool {
+public:
+ DeclareInterfaceID(ITrackTruthMatchingTool, 1, 0);
+
+ virtual std::pair<const xAOD::TruthParticle*, int>
+ getTruthParticle(const Trk::Track *track) const = 0;
+};
+
+#endif /* FASERACTSKALMANFILTER_ITRACKTRUTHMATCHINGTOOL_H */
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/FiducialParticleTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/FiducialParticleTool.cxx
new file mode 100644
index 00000000..101ddfdf
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/FiducialParticleTool.cxx
@@ -0,0 +1,91 @@
+#include "FiducialParticleTool.h"
+#include "GeoPrimitives/CLHEPtoEigenConverter.h"
+#include "TrackerIdentifier/FaserSCT_ID.h"
+#include "TrackerReadoutGeometry/SCT_DetectorManager.h"
+#include "TrackerReadoutGeometry/SiDetectorElement.h"
+
+constexpr float NaN = std::numeric_limits<double>::quiet_NaN();
+
+FiducialParticleTool::FiducialParticleTool(const std::string &type,
+ const std::string &name,
+ const IInterface *parent)
+ : base_class(type, name, parent) {}
+
+StatusCode FiducialParticleTool::initialize() {
+ ATH_CHECK(m_siHitCollectionKey.initialize());
+ ATH_CHECK(detStore()->retrieve(m_sctHelper, "FaserSCT_ID"));
+ ATH_CHECK(detStore()->retrieve(m_detMgr, "SCT"));
+ return StatusCode::SUCCESS;
+}
+
+StatusCode FiducialParticleTool::finalize() { return StatusCode::SUCCESS; }
+
+bool FiducialParticleTool::isFiducial(int barcode) const {
+ auto positions = getTruthPositions(barcode);
+ bool isFiducial{true};
+ for (int station = 1; station < 4; ++station) {
+ double x{positions[station].x()};
+ double y{positions[station].y()};
+ if (std::isnan(x) || std::isnan(y) || (x * x + y * y > 100 * 100))
+ isFiducial = false;
+ }
+ return isFiducial;
+}
+
+HepGeom::Point3D<double>
+FiducialParticleTool::getGlobalPosition(const FaserSiHit &hit) const {
+ Identifier waferId =
+ m_sctHelper->wafer_id(hit.getStation(), hit.getPlane(), hit.getRow(),
+ hit.getModule(), hit.getSensor());
+ auto localStartPos = hit.localStartPosition();
+ auto localEndPos = hit.localEndPosition();
+ HepGeom::Point3D<double> localPos = 0.5 * (localEndPos + localStartPos);
+ const TrackerDD::SiDetectorElement *element =
+ m_detMgr->getDetectorElement(waferId);
+ auto globalPosition =
+ Amg::EigenTransformToCLHEP(element->transformHit()) * localPos;
+ return globalPosition;
+}
+
+std::array<HepGeom::Point3D<double>, 4>
+FiducialParticleTool::getTruthPositions(int barcode) const {
+ // initialize positions as NaN
+ std::array<HepGeom::Point3D<double>, 4> positions{};
+ for (auto &station : positions) {
+ station.setX(NaN);
+ station.setY(NaN);
+ station.setZ(NaN);
+ }
+
+ // get simulated hits
+ SG::ReadHandle<FaserSiHitCollection> siHitCollection(m_siHitCollectionKey);
+ if (!siHitCollection.isValid()) {
+ ATH_MSG_WARNING("FaserSiHitCollection not valid.");
+ return positions;
+ }
+
+ // create map with truth positions in each station
+ std::array<std::vector<HepGeom::Point3D<double>>, 4> hitMap{};
+ for (const FaserSiHit &hit : *siHitCollection) {
+ if (hit.trackNumber() == barcode) {
+ auto position = getGlobalPosition(hit);
+ hitMap[hit.getStation()].push_back(position);
+ }
+ }
+
+ // calculate average position in each station
+ for (int station = 0; station < 4; ++station) {
+ std::vector<HepGeom::Point3D<double>> &hits{hitMap[station]};
+ if (hits.empty()) {
+ continue;
+ } else {
+ // calculate average position of all FaserSiHits in a station
+ auto const count = static_cast<double>(hits.size());
+ HepGeom::Point3D<double> sums{};
+ for (const HepGeom::Point3D<double> &hit : hits)
+ sums += hit;
+ positions[station] = sums / count;
+ }
+ }
+ return positions;
+}
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/FiducialParticleTool.h b/Tracking/Acts/FaserActsKalmanFilter/src/FiducialParticleTool.h
new file mode 100644
index 00000000..8d143573
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/FiducialParticleTool.h
@@ -0,0 +1,37 @@
+#ifndef FASERACTSKALMANFILTER_FIDUCIALEVENTSELECTIONTOOL_H
+#define FASERACTSKALMANFILTER_FIDUCIALEVENTSELECTIONTOOL_H
+
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "FaserActsKalmanFilter/IFiducialParticleTool.h"
+#include "GeoPrimitives/GeoPrimitives.h"
+#include "TrackerSimEvent/FaserSiHitCollection.h"
+
+class FaserSCT_ID;
+namespace TrackerDD {
+class SCT_DetectorManager;
+}
+
+class FiducialParticleTool : public extends<AthAlgTool, IFiducialParticleTool> {
+public:
+ FiducialParticleTool(const std::string &type, const std::string &name,
+ const IInterface *parent);
+ virtual ~FiducialParticleTool() = default;
+ virtual StatusCode initialize() override;
+ virtual StatusCode finalize() override;
+
+ bool isFiducial(int barcode) const override;
+
+ std::array<HepGeom::Point3D<double>, 4>
+ getTruthPositions(int barcode) const override;
+
+private:
+ HepGeom::Point3D<double> getGlobalPosition(const FaserSiHit &hit) const;
+
+ SG::ReadHandleKey<FaserSiHitCollection> m_siHitCollectionKey{
+ this, "FaserSiHitCollection", "SCT_Hits"};
+
+ const FaserSCT_ID *m_sctHelper{nullptr};
+ const TrackerDD::SCT_DetectorManager *m_detMgr{nullptr};
+};
+
+#endif // FASERACTSKALMANFILTER_FIDUCIALEVENTSELECTIONTOOL_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/TrackTruthMatchingTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/TrackTruthMatchingTool.cxx
new file mode 100644
index 00000000..502faa3d
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/TrackTruthMatchingTool.cxx
@@ -0,0 +1,103 @@
+#include "TrackTruthMatchingTool.h"
+#include "TrackerPrepRawData/FaserSCT_Cluster.h"
+#include "TrackerRIO_OnTrack/FaserSCT_ClusterOnTrack.h"
+
+TrackTruthMatchingTool::TrackTruthMatchingTool(const std::string &type,
+ const std::string &name,
+ const IInterface *parent)
+ : base_class(type, name, parent) {}
+
+StatusCode TrackTruthMatchingTool::initialize() {
+ ATH_CHECK(m_simDataCollectionKey.initialize());
+ ATH_CHECK(m_truthParticleContainerKey.initialize());
+ return StatusCode::SUCCESS;
+}
+
+std::pair<const xAOD::TruthParticle *, int>
+TrackTruthMatchingTool::getTruthParticle(const Trk::Track *track) const {
+ const xAOD::TruthParticle *truthParticle = nullptr;
+ const EventContext &ctx = Gaudi::Hive::currentContext();
+ SG::ReadHandle<xAOD::TruthParticleContainer> truthParticleContainer{
+ m_truthParticleContainerKey, ctx};
+ if (!truthParticleContainer.isValid()) {
+ ATH_MSG_WARNING("xAOD::TruthParticleContainer is not valid.");
+ return {truthParticle, -1};
+ }
+ SG::ReadHandle<TrackerSimDataCollection> simDataCollection{
+ m_simDataCollectionKey, ctx};
+ if (!simDataCollection.isValid()) {
+ ATH_MSG_WARNING("TrackerSimDataCollection is not valid.");
+ return {truthParticle, -1};
+ }
+ std::vector<ParticleHitCount> particleHitCounts{};
+ identifyContributingParticles(*track, *simDataCollection, particleHitCounts);
+ if (particleHitCounts.empty()) {
+ ATH_MSG_WARNING("Cannot find any truth particle matched to the track.");
+ return {truthParticle, -1};
+ }
+ int barcode = particleHitCounts.front().barcode;
+ int hitCount = particleHitCounts.front().hitCount;
+ auto it = std::find_if(truthParticleContainer->begin(),
+ truthParticleContainer->end(),
+ [barcode](const xAOD::TruthParticle_v1 *particle) {
+ return particle->barcode() == barcode;
+ });
+ if (it == truthParticleContainer->end()) {
+ ATH_MSG_WARNING("Cannot find particle with barcode "
+ << barcode << " in truth particle container.");
+ return {truthParticle, -1};
+ }
+ truthParticle = *it;
+ return {truthParticle, hitCount};
+}
+
+StatusCode TrackTruthMatchingTool::finalize() { return StatusCode::SUCCESS; }
+
+void TrackTruthMatchingTool::increaseHitCount(
+ std::vector<ParticleHitCount> &particleHitCounts, int barcode) {
+ auto it = std::find_if(
+ particleHitCounts.begin(), particleHitCounts.end(),
+ [=](const ParticleHitCount &phc) { return (phc.barcode == barcode); });
+ // either increase count if we saw the particle before or add it
+ if (it != particleHitCounts.end()) {
+ it->hitCount += 1u;
+ } else {
+ particleHitCounts.push_back({barcode, 1u});
+ }
+}
+
+void TrackTruthMatchingTool::sortHitCount(
+ std::vector<ParticleHitCount> &particleHitCounts) {
+ std::sort(particleHitCounts.begin(), particleHitCounts.end(),
+ [](const ParticleHitCount &lhs, const ParticleHitCount &rhs) {
+ return (lhs.hitCount > rhs.hitCount);
+ });
+}
+
+void TrackTruthMatchingTool::identifyContributingParticles(
+ const Trk::Track &track, const TrackerSimDataCollection &simDataCollection,
+ std::vector<ParticleHitCount> &particleHitCounts) {
+ for (const Trk::MeasurementBase *meas : *track.measurementsOnTrack()) {
+ const auto *clusterOnTrack =
+ dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack *>(meas);
+ if (!clusterOnTrack)
+ continue;
+ std::vector<int> barcodes{};
+ const Tracker::FaserSCT_Cluster *cluster = clusterOnTrack->prepRawData();
+ for (Identifier id : cluster->rdoList()) {
+ if (simDataCollection.count(id) == 0)
+ continue;
+ const auto &deposits = simDataCollection.at(id).getdeposits();
+ for (const TrackerSimData::Deposit &deposit : deposits) {
+ int barcode = deposit.first->barcode();
+ // count each barcode only once for a wafer
+ if (std::find(barcodes.begin(), barcodes.end(), barcode) ==
+ barcodes.end()) {
+ barcodes.push_back(barcode);
+ increaseHitCount(particleHitCounts, barcode);
+ }
+ }
+ }
+ }
+ sortHitCount(particleHitCounts);
+}
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/TrackTruthMatchingTool.h b/Tracking/Acts/FaserActsKalmanFilter/src/TrackTruthMatchingTool.h
new file mode 100644
index 00000000..94ab4405
--- /dev/null
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/TrackTruthMatchingTool.h
@@ -0,0 +1,40 @@
+#ifndef FASERACTSKALMANFILTER_TRACKTRUTHMATCHINGTOOL_H
+#define FASERACTSKALMANFILTER_TRACKTRUTHMATCHINGTOOL_H
+
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "FaserActsKalmanFilter/ITrackTruthMatchingTool.h"
+#include "TrackerSimData/TrackerSimDataCollection.h"
+#include "TrkTrack/Track.h"
+#include "xAODTruth/TruthParticleContainer.h"
+
+class TrackTruthMatchingTool
+ : public extends<AthAlgTool, ITrackTruthMatchingTool> {
+public:
+ TrackTruthMatchingTool(const std::string &type, const std::string &name,
+ const IInterface *parent);
+ virtual ~TrackTruthMatchingTool() = default;
+ virtual StatusCode initialize() override;
+ virtual StatusCode finalize() override;
+
+ std::pair<const xAOD::TruthParticle*, int>
+ getTruthParticle(const Trk::Track *track) const;
+
+private:
+ struct ParticleHitCount {
+ int barcode;
+ size_t hitCount;
+ };
+ static void increaseHitCount(std::vector<ParticleHitCount> &particleHitCounts,
+ int particleId);
+ static void sortHitCount(std::vector<ParticleHitCount> &particleHitCounts);
+ static void identifyContributingParticles(
+ const Trk::Track &track, const TrackerSimDataCollection &simDataCollection,
+ std::vector<ParticleHitCount> &particleHitCounts);
+
+ SG::ReadHandleKey<TrackerSimDataCollection> m_simDataCollectionKey{
+ this, "TrackerSimDataCollection", "SCT_SDO_Map"};
+ SG::ReadHandleKey<xAOD::TruthParticleContainer> m_truthParticleContainerKey{
+ this, "ParticleContainer", "TruthParticles"};
+};
+
+#endif /* FASERACTSKALMANFILTER_TRACKTRUTHMATCHINGTOOL_H */
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/components/FaserActsKalmanFilter_entries.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/components/FaserActsKalmanFilter_entries.cxx
index 967f8064..d9fefbdd 100755
--- a/Tracking/Acts/FaserActsKalmanFilter/src/components/FaserActsKalmanFilter_entries.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/components/FaserActsKalmanFilter_entries.cxx
@@ -29,6 +29,8 @@
#include "../CircleFitTrackSeedTool.h"
#include "../GhostBusters.h"
#include "../CreateTrkTrackTool.h"
+#include "../TrackTruthMatchingTool.h"
+#include "../FiducialParticleTool.h"
DECLARE_COMPONENT(FaserActsKalmanFilterAlg)
DECLARE_COMPONENT(CombinatorialKalmanFilterAlg)
@@ -57,3 +59,5 @@ DECLARE_COMPONENT(SeedingAlg)
DECLARE_COMPONENT(CircleFitTrackSeedTool)
DECLARE_COMPONENT(GhostBusters)
DECLARE_COMPONENT(CreateTrkTrackTool)
+DECLARE_COMPONENT(TrackTruthMatchingTool)
+DECLARE_COMPONENT(FiducialParticleTool)
--
GitLab