diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx index 8a7fabc48085845f921966229c85504e50e70ffa..c06051d2a5d9ee23b527fff81ec78e1ae47a3f20 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 178fd07faa93162f977bbc4565466004730f14f4..9c4a90240c8cbc3fc21681e6b254aecc2064b35d 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 {