From ac36796db4c7fff33ac548c1fbe72711253ea63d Mon Sep 17 00:00:00 2001
From: Tobias Boeckh <tobias.boeckh@cern.ch>
Date: Fri, 18 Nov 2022 23:23:43 +0100
Subject: [PATCH] write out for all truth particles with a momentum larger 50
GeV if they are fiducial in the detector and count number of tracks they are
matched to
---
.../NtupleDumper/src/NtupleDumperAlg.cxx | 84 ++++++++++++++-----
.../NtupleDumper/src/NtupleDumperAlg.h | 9 +-
2 files changed, 69 insertions(+), 24 deletions(-)
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
index dba1c0e68..971118b48 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
@@ -85,9 +85,21 @@ NtupleDumperAlg::Vector3 NtupleDumperAlg::getGlobalPosition(const FaserSiHit &hi
return globalPosition;
}
-StatusCode NtupleDumperAlg::getTruthPositions(int barcode) const {
+std::array<std::array<double, 3>, 4>
+NtupleDumperAlg::getTruthPositions(int barcode) const {
+ // initialize positions as NaN
+ std::array<std::array<double, 3>, 4> positions {};
+ for (auto &row : positions) {
+ for (double &val : row)
+ val = NaN;
+ }
+
+ // get simulated hits
SG::ReadHandle<FaserSiHitCollection> siHitCollection(m_siHitCollectionKey);
- ATH_CHECK(siHitCollection.isValid());
+ if (!siHitCollection.isValid()) {
+ ATH_MSG_WARNING("FaserSiHitCollection not valid.");
+ return positions;
+ }
// create map with truth positions in each station
std::array<std::vector<Vector3>, 4> hitMap {};
@@ -106,35 +118,32 @@ StatusCode NtupleDumperAlg::getTruthPositions(int barcode) const {
m_t_st_y[station].push_back(NaN);
m_t_st_z[station].push_back(NaN);
} else {
+ // count number of FaserSiHits in each station, this should be six
auto const count = static_cast<double>(hits.size());
+ // calculate average position of all FaserSiHits in a station
std::array<double, 3> sums {};
for (const Vector3 &hit : hits) {
sums[0] += hit.x();
sums[1] += hit.y();
sums[2] += hit.z();
}
- m_t_st_x[station].push_back(sums[0] / count);
- m_t_st_y[station].push_back(sums[1] / count);
- m_t_st_z[station].push_back(sums[2] / count);
+ positions[station][0] = sums[0] / count;
+ positions[station][1] = sums[1] / count;
+ positions[station][2] = sums[2] / count;
}
}
- return StatusCode::SUCCESS;
+ return positions;
}
-bool NtupleDumperAlg::isFiducial() const {
+bool NtupleDumperAlg::isFiducial(int barcode) const {
+ auto positions = getTruthPositions(barcode);
bool isFiducial {true};
- for (int station = 0; station < 4; ++station) {
- double st_x {m_t_st_x[station].back()};
- double st_y {m_t_st_y[station].back()};
- if (!std::isnan(st_x) && !std::isnan(st_y)) {
- // distance from center < 100 mm
- if (st_x * st_x + st_y * st_y > 100 * 100)
+ // this does not check if the track is fiducial in the IFT!
+ for (int station = 1; station < 4; ++station) {
+ double x {positions[station][0]};
+ double y {positions[station][1]};
+ if (std::isnan(x) || std::isnan(y) || (x * x + y * y > 100 * 100))
isFiducial = false;
- } else {
- // there have to be simulated hits in stations 1 - 3
- if (station > 0)
- isFiducial = false;
- }
}
return isFiducial;
}
@@ -322,6 +331,10 @@ StatusCode NtupleDumperAlg::initialize()
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");
@@ -619,6 +632,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};
@@ -685,13 +708,19 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
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
+ // 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));
- ATH_CHECK(getTruthPositions(truthParticle->barcode()));
- m_isFiducial.push_back(isFiducial());
+ m_isFiducial.push_back(isFiducial(truthParticle->barcode()));
+ auto positions = getTruthPositions(truthParticle->barcode());
+ for (int station = 0; station < 4; ++station) {
+ m_t_st_x[station].push_back(positions[station][0]);
+ m_t_st_y[station].push_back(positions[station][1]);
+ m_t_st_z[station].push_back(positions[station][2]);
+ }
if (truthParticle->hasProdVtx()) {
m_t_prodVtx_x.push_back(truthParticle->prodVtx()->x());
m_t_prodVtx_y.push_back(truthParticle->prodVtx()->y());
@@ -881,6 +910,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(isFiducial(tp.first));
+ }
+ }
+
/*
// Here we apply the signal selection
// Very simple/unrealistic to start
@@ -1053,6 +1090,9 @@ NtupleDumperAlg::clearTree() const
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;
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
index e35d10f2b..0229bc64a 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
@@ -54,8 +54,8 @@ private:
using Vector3 = HepGeom::Point3D<double>;
Vector3 getGlobalPosition(const FaserSiHit &hit) const;
- StatusCode getTruthPositions(int barcode) const;
- bool isFiducial() const;
+ std::array<std::array<double, 3>, 4> getTruthPositions(int barcode) const;
+ bool isFiducial(int barcode) const;
ServiceHandle <ITHistSvc> m_histSvc;
SG::ReadHandleKey<xAOD::TruthEventContainer> m_truthEventContainer { this, "EventContainer", "TruthEvents", "Truth event container name." };
@@ -95,6 +95,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};
@@ -234,6 +235,10 @@ private:
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;
--
GitLab