diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
index 7b3a631a763f4767601482dd8fb1ab19892314b2..468a6739bdd9606cc8f34c4ff0996cbda393273d 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
@@ -335,6 +335,18 @@ StatusCode NtupleDumperAlg::initialize()
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("t_st0_px", &m_t_st_px[0]);
+ m_tree->Branch("t_st0_py", &m_t_st_py[0]);
+ m_tree->Branch("t_st0_pz", &m_t_st_pz[0]);
+ m_tree->Branch("t_st1_px", &m_t_st_px[1]);
+ m_tree->Branch("t_st1_py", &m_t_st_py[1]);
+ m_tree->Branch("t_st1_pz", &m_t_st_pz[1]);
+ m_tree->Branch("t_st2_px", &m_t_st_px[2]);
+ m_tree->Branch("t_st2_py", &m_t_st_py[2]);
+ m_tree->Branch("t_st2_pz", &m_t_st_pz[2]);
+ m_tree->Branch("t_st3_px", &m_t_st_px[3]);
+ m_tree->Branch("t_st3_py", &m_t_st_py[3]);
+ m_tree->Branch("t_st3_pz", &m_t_st_pz[3]);
m_tree->Branch("isFiducial", &m_isFiducial);
m_tree->Branch("truthParticleBarcode", &m_truthParticleBarcode);
@@ -1054,6 +1066,14 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_t_st_y[station].push_back(positions[station].y());
m_t_st_z[station].push_back(positions[station].z());
}
+
+ auto momenta = m_fiducialParticleTool->getTruthMomenta(truthParticle->barcode());
+ for (int station = 0; station < 4; ++station) {
+ m_t_st_px[station].push_back(momenta[station].x());
+ m_t_st_py[station].push_back(momenta[station].y());
+ m_t_st_pz[station].push_back(momenta[station].z());
+ }
+
if (truthParticle->hasProdVtx()) {
m_t_prodVtx_x.push_back(truthParticle->prodVtx()->x());
m_t_prodVtx_y.push_back(truthParticle->prodVtx()->y());
@@ -1550,6 +1570,11 @@ NtupleDumperAlg::clearTree() const
m_t_st_y[station].clear();
m_t_st_z[station].clear();
}
+ for (int station = 0; station < 4; ++station) {
+ m_t_st_px[station].clear();
+ m_t_st_py[station].clear();
+ m_t_st_pz[station].clear();
+ }
m_truthParticleBarcode.clear();
m_truthParticleMatchedTracks.clear();
m_truthParticleIsFiducial.clear();
@@ -1634,6 +1659,11 @@ void NtupleDumperAlg::clearTrackTruth() const {
m_t_st_y[station].push_back(999999);
m_t_st_z[station].push_back(999999);
}
+ for (int station = 0; station < 4; ++station) {
+ m_t_st_px[station].push_back(999999);
+ m_t_st_py[station].push_back(999999);
+ m_t_st_pz[station].push_back(999999);
+ }
m_isFiducial.push_back(false);
}
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
index d636d8ef70dcdd7fbcb90fd44eddccd8c2ac36ac..d4fc22b9754617ea8c68aba3c365ef334abd0d5d 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
@@ -317,6 +317,9 @@ private:
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::array<std::vector<double>, 4> m_t_st_px; // x components of the true momentum at each station
+ mutable std::array<std::vector<double>, 4> m_t_st_py; // y components of the true momentum at each station
+ mutable std::array<std::vector<double>, 4> m_t_st_pz; // z components of the true momentum at 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
diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/IFiducialParticleTool.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/IFiducialParticleTool.h
index 1df9aeab2ce991ee5d683d92fa5efb72eca7a5dd..081e6dc27913107fbd0cf4d4cb7a1646316eff49 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/IFiducialParticleTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/IFiducialParticleTool.h
@@ -23,6 +23,12 @@ public:
*/
virtual std::array<HepGeom::Point3D<double>, 4>
getTruthPositions(int barcode) const = 0;
+
+ /** Return average truth momentum in each station.
+ * @param barcode of a xAOD::TruthParticle
+ */
+ virtual std::array<HepGeom::Point3D<double>, 4>
+ getTruthMomenta(int barcode) const = 0;
};
#endif // FASERACTSKALMANFILTER_IFIDUCIALEVENTSELECTIONTOOL_H
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/FiducialParticleTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/FiducialParticleTool.cxx
index 101ddfdf99a645aaef69f6d824e782ff7525928f..71a469ba4cb1b20e682f610d0be13238394e9278 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/FiducialParticleTool.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/FiducialParticleTool.cxx
@@ -6,8 +6,7 @@
constexpr float NaN = std::numeric_limits<double>::quiet_NaN();
-FiducialParticleTool::FiducialParticleTool(const std::string &type,
- const std::string &name,
+FiducialParticleTool::FiducialParticleTool(const std::string &type, const std::string &name,
const IInterface *parent)
: base_class(type, name, parent) {}
@@ -32,23 +31,32 @@ bool FiducialParticleTool::isFiducial(int barcode) const {
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());
+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;
+ 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 {
+HepGeom::Point3D<double> FiducialParticleTool::getMomentum(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();
+ const TrackerDD::SiDetectorElement *element = m_detMgr->getDetectorElement(waferId);
+ auto globalStartPosition = Amg::EigenTransformToCLHEP(element->transformHit()) * localStartPos;
+ auto globalEndPosition = Amg::EigenTransformToCLHEP(element->transformHit()) * localEndPos;
+ double p = hit.particleLink()->momentum().rho();
+ HepGeom::Point3D<double> globalDirection = globalEndPosition - globalStartPosition;
+ HepGeom::Point3D<double> momentum = globalDirection / globalDirection.mag() * p;
+ return momentum;
+}
+
+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) {
@@ -89,3 +97,46 @@ FiducialParticleTool::getTruthPositions(int barcode) const {
}
return positions;
}
+
+std::array<HepGeom::Point3D<double>, 4> FiducialParticleTool::getTruthMomenta(int barcode) const {
+ // initialize positions as NaN
+ std::array<HepGeom::Point3D<double>, 4> momenta{};
+ for (auto &station : momenta) {
+ 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 momenta;
+ }
+
+ // create map with truth momenta in each station
+ std::array<std::vector<HepGeom::Point3D<double>>, 4> hitMap{};
+ for (const FaserSiHit &hit : *siHitCollection) {
+ if (hit.trackNumber() == barcode) {
+ HepGeom::Point3D<double> momentum = getMomentum(hit);
+ hitMap[hit.getStation()].push_back(momentum);
+ }
+ }
+
+ // calculate average momentum 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 momentum 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;
+ }
+ momenta[station] = sums / count;
+ }
+ }
+ return momenta;
+}
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/FiducialParticleTool.h b/Tracking/Acts/FaserActsKalmanFilter/src/FiducialParticleTool.h
index 8d1435734c6982e9dcbd0abd52acb07003914366..ecba3a8318f3730127ddf8786722a3509b96b4ac 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/FiducialParticleTool.h
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/FiducialParticleTool.h
@@ -23,6 +23,10 @@ public:
std::array<HepGeom::Point3D<double>, 4>
getTruthPositions(int barcode) const override;
+ std::array<HepGeom::Point3D<double>, 4>
+ getTruthMomenta(int barcode) const override;
+ HepGeom::Point3D<double>
+ getMomentum(const FaserSiHit &hit) const;
private:
HepGeom::Point3D<double> getGlobalPosition(const FaserSiHit &hit) const;