From a9081b5dfb55b19f34a92eb48d6a76b34cfa1bb0 Mon Sep 17 00:00:00 2001
From: Deion Elgin Fellers <dfellers@lxplus7106.cern.ch>
Date: Tue, 21 Feb 2023 10:19:54 +0100
Subject: [PATCH] update Tobias' code that stores track position at maximum
track radius
---
.../NtupleDumper/src/NtupleDumperAlg.cxx | 74 ++++---------------
.../NtupleDumper/src/NtupleDumperAlg.h | 6 --
2 files changed, 14 insertions(+), 66 deletions(-)
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
index 8a7fabc4..c06051d2 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
@@ -390,16 +390,6 @@ StatusCode NtupleDumperAlg::initialize()
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge13", m_HistRandomCharge[13]));
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge14", m_HistRandomCharge[14]));
- // get average position of each tracking station
- for (const TrackerDD::SiDetectorElement *elem : *m_detMgr->getDetectorElementCollection()) {
- int station = m_sctHelper->station(elem->identify());
- m_tracking_station_positions[station] += elem->center().z();
- }
- for (size_t i = 0; i < m_tracking_station_positions.size(); ++i) {
- m_tracking_station_positions[i] /= 48; // each station has 8 * 3 * 2 = 48 wafers
- ATH_MSG_VERBOSE("Station " << i << " : " << m_tracking_station_positions[i]);
- }
-
if (m_onlyBlinded){
ATH_MSG_INFO("Only events that would be blinded are saved in ntuple");
m_doBlinding = false;
@@ -861,6 +851,20 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_pzdown.push_back(downstreamParameters->momentum().z());
m_pdown.push_back(sqrt( pow(downstreamParameters->momentum().x(),2) + pow(downstreamParameters->momentum().y(),2) + pow(downstreamParameters->momentum().z(),2) ));
+ // find and store the position of the measurement on track with the largest radius
+ double maxRadius = 0;
+ Amg::Vector3D positionAtMaxRadius {};
+ for (const Trk::TrackParameters* trackParameters : *track->trackParameters()) {
+ if (radius(trackParameters->position()) > maxRadius) {
+ maxRadius = radius(trackParameters->position());
+ positionAtMaxRadius = trackParameters->position();
+ }
+ }
+ m_r_atMaxRadius.push_back(maxRadius);
+ m_x_atMaxRadius.push_back(positionAtMaxRadius.x());
+ m_y_atMaxRadius.push_back(positionAtMaxRadius.y());
+ m_z_atMaxRadius.push_back(positionAtMaxRadius.z());
+
if (isMC) { // if simulation, store track truth info as well
auto [truthParticle, hitCount] = m_trackTruthMatchingTool->getTruthParticle(track);
if (truthParticle != nullptr) {
@@ -940,12 +944,6 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_yCalo.push_back(999999);
m_thetaxCalo.push_back(999999);
m_thetayCalo.push_back(999999);
- m_x_atMaxRadius.push_back(999999);
- m_y_atMaxRadius.push_back(999999);
- m_z_atMaxRadius.push_back(999999);
- m_r_atMaxRadius.push_back(999999);
-
- std::optional<Acts::Vector3> propagationResult0, propagationResult1;
// Define paramters for track extrapolation from most upstream measurement
Amg::Vector3D position_up = upstreamParameters->position();
@@ -1065,20 +1063,6 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
ATH_MSG_INFO("Calo null targetParameters");
}
- // get the track position with the maximum radius
- propagationResult0 = positionAtMaxRadius(ctx, startParameters_up, m_tracking_station_positions[2], Acts::forward); // get track radius at maximum position between station 1 and 2
- propagationResult1 = positionAtMaxRadius(ctx, startParameters_down, m_tracking_station_positions[2], Acts::backward); // get track radius at maximum position between station 3 and 2
- if (propagationResult0.has_value() && propagationResult1.has_value()) {
- auto position0 = propagationResult0.value();
- auto position1 = propagationResult1.value();
- auto position = radius(position0) > radius(position1) ? position0 : position1;
-
- m_x_atMaxRadius[m_longTracks] = position.x();
- m_y_atMaxRadius[m_longTracks] = position.y();
- m_z_atMaxRadius[m_longTracks] = position.z();
- m_r_atMaxRadius[m_longTracks] = radius(position);
- }
-
m_longTracks++;
}
@@ -1421,33 +1405,3 @@ void NtupleDumperAlg::clearTrackTruth() const {
double NtupleDumperAlg::radius(const Acts::Vector3 &position) const {
return std::sqrt(position.x() * position.x() + position.y() * position.y());
}
-
-std::optional<Acts::Vector3> NtupleDumperAlg::positionAtMaxRadius(
- const EventContext &ctx, const Acts::BoundTrackParameters &startParameters,
- double targetPosition, Acts::NavigationDirection navDir) const {
- // create target surface from translation and rotation
- Acts::AngleAxis3 rotZ(-M_PI / 2, Acts::Vector3(0., 0., 1.));
- Acts::Translation3 targetTranslation{0., 0., targetPosition};
- auto targetTransform = Acts::Transform3(targetTranslation * rotZ);
- auto targetSurface = Acts::Surface::makeShared<Acts::PlaneSurface>(targetTransform);
- // propagate start parameters to target surface
- ActsPropagationOutput propagationOutput = m_extrapolationTool->propagationSteps(
- ctx, startParameters, *targetSurface, navDir);
- std::vector<Acts::detail::Step> propagationSteps = propagationOutput.first;
- if (propagationSteps.size() > 0) {
- for (const auto &step : propagationSteps) {
- ATH_MSG_VERBOSE("step: " << step.position.x() << ", " << step.position.y() << ", " << step.position.z());
- }
- // find step with maximum radius
- auto maxStep = std::max_element(propagationSteps.begin(), propagationSteps.end(),
- [&](const auto &lhs, const auto &rhs) {
- return radius(lhs.position) < radius(rhs.position);
- });
- return maxStep->position;
- } else {
- ATH_MSG_WARNING("Extrapolation failed");
- return {};
- }
-}
-
-
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
index 178fd07f..9c4a9024 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
@@ -59,10 +59,6 @@ private:
void FillWaveBranches(const xAOD::WaveformHitContainer &wave) const;
void addCalibratedBranches(const std::string &name, int nchannels, int first);
double radius(const Acts::Vector3 &position) const;
- std::optional<Acts::Vector3> positionAtMaxRadius(
- const EventContext &ctx, const Acts::BoundTrackParameters &startParameters,
- double targetPosition, Acts::NavigationDirection navDir = Acts::forward) const;
-
ServiceHandle <ITHistSvc> m_histSvc;
@@ -347,8 +343,6 @@ private:
mutable int m_eventsPassed = 0;
- mutable std::array<double, 4> m_tracking_station_positions {}; // get average position of each tracking station for track extrapolation
-
};
inline const ServiceHandle <ITHistSvc> &NtupleDumperAlg::histSvc() const {
--
GitLab