diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
index edcfbb032088a07937e17c20a59532b3e58c5701..6880f5dce596a18be3f51aaf014f548ce409bb54 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
@@ -15,6 +15,14 @@
#include "TrackerPrepRawData/FaserSCT_Cluster.h"
#include "xAODTruth/TruthParticle.h"
#include "AtlasHepMC/GenEvent.h"
+
+#include "Acts/Definitions/TrackParametrization.hpp"
+#include "Acts/EventData/Measurement.hpp"
+#include "Acts/Geometry/GeometryContext.hpp"
+#include "Acts/Geometry/GeometryIdentifier.hpp"
+#include "Acts/Geometry/TrackingGeometry.hpp"
+#include "Acts/Surfaces/Surface.hpp"
+
#include <cmath>
#include <TH1F.h>
#include <numeric>
@@ -979,24 +987,24 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_truthd1_px.push_back(particle->p4().X());
m_truthd1_py.push_back(particle->p4().Y());
m_truthd1_pz.push_back(particle->p4().Z());
- m_truthd1_isFiducial.push_back(m_fiducialParticleTool->isFiducial(particle->barcode()));
+ m_truthd1_isFiducial.push_back(m_fiducialParticleTool->isFiducial(particle->barcode()));
- // if ( particle->hasProdVtx()) {
- // m_truthd1_x.push_back(particle->prodVtx()->x());
- // m_truthd1_y.push_back(particle->prodVtx()->y());
- // m_truthd1_z.push_back(particle->prodVtx()->z());
- // } else {
- // m_truthd1_x.push_back(999999);
- // m_truthd1_y.push_back(999999);
- // m_truthd1_z.push_back(999999);
- // }
- // for (int station = 1; station < 4; ++station) {
- // m_truthd1_x.push_back(positions[station].x());
- // m_truthd1_y.push_back(positions[station].y());
- // m_truthd1_z.push_back(positions[station].z());
- // }
- // }
- // }
+ if ( particle->hasProdVtx()) {
+ m_truthd1_x.push_back(particle->prodVtx()->x());
+ m_truthd1_y.push_back(particle->prodVtx()->y());
+ m_truthd1_z.push_back(particle->prodVtx()->z());
+ } else {
+ m_truthd1_x.push_back(999999);
+ m_truthd1_y.push_back(999999);
+ m_truthd1_z.push_back(999999);
+ }
+ for (int station = 1; station < 4; ++station) {
+ m_truthd1_x.push_back(positions[station].x());
+ m_truthd1_y.push_back(positions[station].y());
+ m_truthd1_z.push_back(positions[station].z());
+ }
+ }
+ }
}
}
}
@@ -1110,320 +1118,332 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
SG::ReadHandle<TrackCollection> tc {m_trackCollectionWithoutIFT, ctx}; // use track collection that excludes IFT
trackCollection = tc;
}
- // ATH_CHECK(trackCollection.isValid());
+ ATH_CHECK(trackCollection.isValid());
// create list of good tracks in fiducial region (r < 95mm) and chi2/nDoF < 25, #Layers >= 7, #Hits >= 12
std::vector<const Trk::Track*> goodTracks {};
- if (trackCollection.isValid())
+ for (const Trk::Track* track : *trackCollection)
{
- for (const Trk::Track* track : *trackCollection)
- {
- if (track == nullptr) continue;
-
- std::set<std::pair<int, int>> layerMap;
- std::set<int> stationMap;
- // Check for hit in the three downstream stations
- HitSet hitSet {24};
- for (auto measurement : *(track->measurementsOnTrack())) {
- const Tracker::FaserSCT_ClusterOnTrack* cluster = dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*>(measurement);
- if (cluster != nullptr) {
- Identifier id = cluster->identify();
- int station = m_sctHelper->station(id);
- int layer = m_sctHelper->layer(id);
- int side = m_sctHelper->side(id);
- stationMap.emplace(station);
- layerMap.emplace(station, layer);
- hitSet.set(23 - (station * 6 + layer * 2 + side));
- }
- }
- m_hitSet.push_back(hitSet.to_ulong());
- if (stationMap.count(1) == 0 || stationMap.count(2) == 0 || stationMap.count(3) == 0) continue;
-
- const Trk::TrackParameters* upstreamParameters = track->trackParameters()->front();
- const Trk::TrackParameters* downstreamParameters = track->trackParameters()->back();
-
- if ((upstreamParameters == nullptr) || (downstreamParameters == nullptr)) continue;
-
- m_nLayers.push_back(layerMap.size());
-
- m_Chi2.push_back(track->fitQuality()->chiSquared());
- m_DoF.push_back(track->fitQuality()->numberDoF());
- const Trk::MeasurementBase *measurement = track->measurementsOnTrack()->front();
- const Tracker::FaserSCT_ClusterOnTrack* cluster =
- dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*>(measurement);
- Identifier id = cluster->identify();
- m_module_eta0.push_back(m_sctHelper->eta_module(id));
- m_module_phi0.push_back(m_sctHelper->phi_module(id));
-
- m_nHit0.push_back(stationMap.count(0));
- m_nHit1.push_back(stationMap.count(1));
- m_nHit2.push_back(stationMap.count(2));
- m_nHit3.push_back(stationMap.count(3));
-
- m_charge.push_back( (int) upstreamParameters->charge() );
-
- m_xup.push_back(upstreamParameters->position().x());
- m_yup.push_back(upstreamParameters->position().y());
- m_zup.push_back(upstreamParameters->position().z());
- m_pxup.push_back(upstreamParameters->momentum().x());
- m_pyup.push_back(upstreamParameters->momentum().y());
- m_pzup.push_back(upstreamParameters->momentum().z());
- m_pup.push_back(sqrt( pow(upstreamParameters->momentum().x(),2) + pow(upstreamParameters->momentum().y(),2) + pow(upstreamParameters->momentum().z(),2) ));
-
- m_xdown.push_back(downstreamParameters->position().x());
- m_ydown.push_back(downstreamParameters->position().y());
- m_zdown.push_back(downstreamParameters->position().z());
- m_pxdown.push_back(downstreamParameters->momentum().x());
- m_pydown.push_back(downstreamParameters->momentum().y());
- 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();
+ if (track == nullptr) continue;
+
+ std::set<std::pair<int, int>> layerMap;
+ std::set<int> stationMap;
+ // Check for hit in the three downstream stations
+ HitSet hitSet {24};
+ for (auto measurement : *(track->measurementsOnTrack())) {
+ const Tracker::FaserSCT_ClusterOnTrack* cluster = dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*>(measurement);
+ if (cluster != nullptr) {
+ Identifier id = cluster->identify();
+ int station = m_sctHelper->station(id);
+ int layer = m_sctHelper->layer(id);
+ int side = m_sctHelper->side(id);
+ stationMap.emplace(station);
+ layerMap.emplace(station, layer);
+ hitSet.set(23 - (station * 6 + layer * 2 + side));
}
+ }
+ m_hitSet.push_back(hitSet.to_ulong());
+ if (stationMap.count(1) == 0 || stationMap.count(2) == 0 || stationMap.count(3) == 0) continue;
+
+ const Trk::TrackParameters* upstreamParameters = track->trackParameters()->front();
+ const Trk::TrackParameters* downstreamParameters = track->trackParameters()->back();
+
+ if ((upstreamParameters == nullptr) || (downstreamParameters == nullptr)) continue;
+
+ m_nLayers.push_back(layerMap.size());
+
+ m_Chi2.push_back(track->fitQuality()->chiSquared());
+ m_DoF.push_back(track->fitQuality()->numberDoF());
+ const Trk::MeasurementBase *measurement = track->measurementsOnTrack()->front();
+ const Tracker::FaserSCT_ClusterOnTrack* cluster =
+ dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*>(measurement);
+ Identifier id = cluster->identify();
+ m_module_eta0.push_back(m_sctHelper->eta_module(id));
+ m_module_phi0.push_back(m_sctHelper->phi_module(id));
+
+ m_nHit0.push_back(stationMap.count(0));
+ m_nHit1.push_back(stationMap.count(1));
+ m_nHit2.push_back(stationMap.count(2));
+ m_nHit3.push_back(stationMap.count(3));
+
+ m_charge.push_back( (int) upstreamParameters->charge() );
+
+ m_xup.push_back(upstreamParameters->position().x());
+ m_yup.push_back(upstreamParameters->position().y());
+ m_zup.push_back(upstreamParameters->position().z());
+ m_pxup.push_back(upstreamParameters->momentum().x());
+ m_pyup.push_back(upstreamParameters->momentum().y());
+ m_pzup.push_back(upstreamParameters->momentum().z());
+ m_pup.push_back(sqrt( pow(upstreamParameters->momentum().x(),2) + pow(upstreamParameters->momentum().y(),2) + pow(upstreamParameters->momentum().z(),2) ));
+
+ m_xdown.push_back(downstreamParameters->position().x());
+ m_ydown.push_back(downstreamParameters->position().y());
+ m_zdown.push_back(downstreamParameters->position().z());
+ m_pxdown.push_back(downstreamParameters->momentum().x());
+ m_pydown.push_back(downstreamParameters->momentum().y());
+ 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 ((m_Chi2.back() / m_DoF.back() < m_maxChi2NDOF) &&
- (m_nLayers.back() >= m_minLayers) &&
- (m_DoF.back() + 5 >= m_minHits) &&
- (m_r_atMaxRadius.back() < 95)) {
- goodTracks.push_back(track);
- }
+ }
+ 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 ((m_Chi2.back() / m_DoF.back() < m_maxChi2NDOF) &&
+ (m_nLayers.back() >= m_minLayers) &&
+ (m_DoF.back() + 5 >= m_minHits) &&
+ (m_r_atMaxRadius.back() < 95)) {
+ goodTracks.push_back(track);
+ }
- if (isMC) { // if simulation, store track truth info as well
- // 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());
- // if (truthParticle->nParents() > 0){
- // m_t_pdg_parent.push_back(truthParticle->parent(0)->pdgId());
- // m_t_barcode_parent.push_back(truthParticle->parent(0)->barcode());
- // } else {
- // m_t_pdg_parent.push_back(0);
- // m_t_barcode_parent.push_back(-1);
- // }
- // // 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 (isMC) { // if simulation, store track truth info as well
+ 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());
+ if (truthParticle->nParents() > 0){
+ m_t_pdg_parent.push_back(truthParticle->parent(0)->pdgId());
+ m_t_barcode_parent.push_back(truthParticle->parent(0)->barcode());
+ } else {
+ m_t_pdg_parent.push_back(0);
+ m_t_barcode_parent.push_back(-1);
+ }
+ // 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());
+ }
- // 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());
- // }
+ 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());
- // m_t_prodVtx_z.push_back(truthParticle->prodVtx()->z());
- // } else {
- // m_t_prodVtx_x.push_back(999999);
- // m_t_prodVtx_y.push_back(999999);
- // m_t_prodVtx_z.push_back(999999);
- // }
- // 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(999999);
- // m_t_decayVtx_y.push_back(999999);
- // m_t_decayVtx_z.push_back(999999);
- // }
- // 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 {
- // ATH_MSG_WARNING("Can not find truthParticle.");
- // clearTrackTruth();
- // }
+ 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(999999);
+ m_t_prodVtx_y.push_back(999999);
+ m_t_prodVtx_z.push_back(999999);
+ }
+ 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(999999);
+ m_t_decayVtx_y.push_back(999999);
+ m_t_decayVtx_z.push_back(999999);
+ }
+ 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 {
+ ATH_MSG_WARNING("Can not find truthParticle.");
clearTrackTruth();
}
+ } else {
+ clearTrackTruth();
+ }
- // fill extrapolation vectors with dummy values, will get set to real number if the track extrapolation succeeds
- m_xVetoNu.push_back(999999);
- m_yVetoNu.push_back(999999);
- m_thetaxVetoNu.push_back(999999);
- m_thetayVetoNu.push_back(999999);
- m_xVetoStation1.push_back(999999);
- m_yVetoStation1.push_back(999999);
- m_thetaxVetoStation1.push_back(999999);
- m_thetayVetoStation1.push_back(999999);
- m_xVetoStation2.push_back(999999);
- m_yVetoStation2.push_back(999999);
- m_thetaxVetoStation2.push_back(999999);
- m_thetayVetoStation2.push_back(999999);
- m_xTrig.push_back(999999);
- m_yTrig.push_back(999999);
- m_thetaxTrig.push_back(999999);
- m_thetayTrig.push_back(999999);
- m_xPreshower1.push_back(999999);
- m_yPreshower1.push_back(999999);
- m_thetaxPreshower1.push_back(999999);
- m_thetayPreshower1.push_back(999999);
- m_xPreshower2.push_back(999999);
- m_yPreshower2.push_back(999999);
- m_thetaxPreshower2.push_back(999999);
- m_thetayPreshower2.push_back(999999);
- m_xCalo.push_back(999999);
- m_yCalo.push_back(999999);
- m_thetaxCalo.push_back(999999);
- m_thetayCalo.push_back(999999);
-
- // Define paramters for track extrapolation from most upstream measurement
- Amg::Vector3D position_up = upstreamParameters->position();
- Amg::Vector3D momentum_up = upstreamParameters->momentum();
- Acts::BoundVector params_up = Acts::BoundVector::Zero();
- params_up[Acts::eBoundLoc0] = -position_up.y();
- params_up[Acts::eBoundLoc1] = position_up.x();
- params_up[Acts::eBoundPhi] = momentum_up.phi();
- params_up[Acts::eBoundTheta] = momentum_up.theta();
- params_up[Acts::eBoundQOverP] = upstreamParameters->charge() / (momentum_up.mag() * 1_MeV);
- params_up[Acts::eBoundTime] = 0;
- std::optional<Acts::BoundSquareMatrix> cov_up;
- auto startSurface_up = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, position_up.z()), Acts::Vector3(0, 0, 1));
- Acts::BoundTrackParameters startParameters_up(startSurface_up, params_up, cov_up, Acts::ParticleHypothesis::muon());
-
- // Define paramters for track extrapolation from most downstream measurement
- Amg::Vector3D position_down = downstreamParameters->position();
- Amg::Vector3D momentum_down = downstreamParameters->momentum();
- Acts::BoundVector params_down = Acts::BoundVector::Zero();
- params_down[Acts::eBoundLoc0] = -position_down.y();
- params_down[Acts::eBoundLoc1] = position_down.x();
- params_down[Acts::eBoundPhi] = momentum_down.phi();
- params_down[Acts::eBoundTheta] = momentum_down.theta();
- params_down[Acts::eBoundQOverP] = downstreamParameters->charge() / (momentum_down.mag() * 1_MeV);
- params_down[Acts::eBoundTime] = 0;
- std::optional<Acts::BoundSquareMatrix> cov_down;
- auto startSurface_down = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, position_down.z()), Acts::Vector3(0, 0, 1));
- Acts::BoundTrackParameters startParameters_down(startSurface_down, params_down, cov_down, Acts::ParticleHypothesis::muon());
-
- // Extrapolate track to scintillators
- auto targetSurface_VetoNu = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, -3112.0), Acts::Vector3(0, 0, 1)); // -3112 mm is z position of VetoNu planes touching
- std::optional<const Acts::BoundTrackParameters> targetParameters_VetoNu =m_extrapolationTool->propagate(ctx, startParameters_up, *targetSurface_VetoNu, Acts::Direction::Backward);
- if (targetParameters_VetoNu.has_value()) {
- auto targetPosition_VetoNu = targetParameters_VetoNu->position(gctx);
- auto targetMomentum_VetoNu = targetParameters_VetoNu->momentum();
- m_xVetoNu[m_longTracks] = targetPosition_VetoNu.x();
- m_yVetoNu[m_longTracks] = targetPosition_VetoNu.y();
- m_thetaxVetoNu[m_longTracks] = atan(targetMomentum_VetoNu[0]/targetMomentum_VetoNu[2]);
- m_thetayVetoNu[m_longTracks] = atan(targetMomentum_VetoNu[1]/targetMomentum_VetoNu[2]);
- } else {
- ATH_MSG_INFO("vetoNu null targetParameters");
- }
-
- auto targetSurface_Veto1 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, -1769.65), Acts::Vector3(0, 0, 1)); // -1769.65 mm is z position of center of operational layer in Veto station 1
- std::optional<const Acts::BoundTrackParameters> targetParameters_Veto1 =m_extrapolationTool->propagate(ctx, startParameters_up, *targetSurface_Veto1, Acts::Direction::Backward);
- if (targetParameters_Veto1.has_value()) {
- auto targetPosition_Veto1 = targetParameters_Veto1->position(gctx);
- auto targetMomentum_Veto1 = targetParameters_Veto1->momentum();
- m_xVetoStation1[m_longTracks] = targetPosition_Veto1.x();
- m_yVetoStation1[m_longTracks] = targetPosition_Veto1.y();
- m_thetaxVetoStation1[m_longTracks] = atan(targetMomentum_Veto1[0]/targetMomentum_Veto1[2]);
- m_thetayVetoStation1[m_longTracks] = atan(targetMomentum_Veto1[1]/targetMomentum_Veto1[2]);
- } else {
- ATH_MSG_INFO("veto1 null targetParameters");
- }
+ // fill extrapolation vectors with dummy values, will get set to real number if the track extrapolation succeeds
+ m_xVetoNu.push_back(999999);
+ m_yVetoNu.push_back(999999);
+ m_thetaxVetoNu.push_back(999999);
+ m_thetayVetoNu.push_back(999999);
+ m_xVetoStation1.push_back(999999);
+ m_yVetoStation1.push_back(999999);
+ m_thetaxVetoStation1.push_back(999999);
+ m_thetayVetoStation1.push_back(999999);
+ m_xVetoStation2.push_back(999999);
+ m_yVetoStation2.push_back(999999);
+ m_thetaxVetoStation2.push_back(999999);
+ m_thetayVetoStation2.push_back(999999);
+ m_xTrig.push_back(999999);
+ m_yTrig.push_back(999999);
+ m_thetaxTrig.push_back(999999);
+ m_thetayTrig.push_back(999999);
+ m_xPreshower1.push_back(999999);
+ m_yPreshower1.push_back(999999);
+ m_thetaxPreshower1.push_back(999999);
+ m_thetayPreshower1.push_back(999999);
+ m_xPreshower2.push_back(999999);
+ m_yPreshower2.push_back(999999);
+ m_thetaxPreshower2.push_back(999999);
+ m_thetayPreshower2.push_back(999999);
+ m_xCalo.push_back(999999);
+ m_yCalo.push_back(999999);
+ m_thetaxCalo.push_back(999999);
+ m_thetayCalo.push_back(999999);
+
+ // Define paramters for track extrapolation from most upstream measurement
+ Amg::Vector3D position_up = upstreamParameters->position();
+ Amg::Vector3D momentum_up = upstreamParameters->momentum();
+ Acts::BoundVector params_up = Acts::BoundVector::Zero();
+ params_up[Acts::eBoundLoc0] = -position_up.y();
+ params_up[Acts::eBoundLoc1] = position_up.x();
+ params_up[Acts::eBoundPhi] = momentum_up.phi();
+ params_up[Acts::eBoundTheta] = momentum_up.theta();
+ params_up[Acts::eBoundQOverP] = upstreamParameters->charge() / (momentum_up.mag() * 1_MeV);
+ params_up[Acts::eBoundTime] = 0;
+
+ Acts::BoundSquareMatrix cov = Acts::BoundSquareMatrix::Zero();
+ cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = 1;
+ cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = 1;
+ cov(Acts::eBoundPhi, Acts::eBoundPhi) = 1;
+ cov(Acts::eBoundTheta, Acts::eBoundTheta) = 1;
+ cov(Acts::eBoundQOverP, Acts::eBoundQOverP) = 1;
+ cov(Acts::eBoundTime, Acts::eBoundTime) = 1;
+
+
+ auto startSurface_up = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, position_up.z()), Acts::Vector3(0, 0, 1));
+ Acts::BoundTrackParameters startParameters_up(std::move(startSurface_up), params_up, cov, Acts::ParticleHypothesis::muon());
+
+ // Define paramters for track extrapolation from most downstream measurement
+ Amg::Vector3D position_down = downstreamParameters->position();
+ Amg::Vector3D momentum_down = downstreamParameters->momentum();
+ Acts::BoundVector params_down = Acts::BoundVector::Zero();
+ params_down[Acts::eBoundLoc0] = -position_down.y();
+ params_down[Acts::eBoundLoc1] = position_down.x();
+ params_down[Acts::eBoundPhi] = momentum_down.phi();
+ params_down[Acts::eBoundTheta] = momentum_down.theta();
+ params_down[Acts::eBoundQOverP] = downstreamParameters->charge() / (momentum_down.mag() * 1_MeV);
+ params_down[Acts::eBoundTime] = 0;
+ auto startSurface_down = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, position_down.z()), Acts::Vector3(0, 0, 1));
+ Acts::BoundTrackParameters startParameters_down(std::move(startSurface_down), params_down, cov, Acts::ParticleHypothesis::muon());
+
+ // Extrapolate track to scintillators
+ auto targetSurface_VetoNu = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, -3112.0), Acts::Vector3(0, 0, 1)); // -3112 mm is z position of VetoNu planes touching
+ std::optional<const Acts::BoundTrackParameters> targetParameters_VetoNu =m_extrapolationTool->propagate(ctx, startParameters_up, *targetSurface_VetoNu, Acts::Direction::Backward);
+
+ if (targetParameters_VetoNu) {
+ auto targetPosition_VetoNu = targetParameters_VetoNu->position(gctx);
+ auto targetMomentum_VetoNu = targetParameters_VetoNu->momentum();
+ m_xVetoNu[m_longTracks] = targetPosition_VetoNu.x();
+ m_yVetoNu[m_longTracks] = targetPosition_VetoNu.y();
+ m_thetaxVetoNu[m_longTracks] = atan(targetMomentum_VetoNu[0]/targetMomentum_VetoNu[2]);
+ m_thetayVetoNu[m_longTracks] = atan(targetMomentum_VetoNu[1]/targetMomentum_VetoNu[2]);
+ } else {
+ ATH_MSG_INFO("vetoNu null targetParameters");
+ }
- auto targetSurface_Veto2 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, -1609.65), Acts::Vector3(0, 0, 1)); // -1609.65 mm is z position of where planes touch in Veto station 2
- std::optional<const Acts::BoundTrackParameters> targetParameters_Veto2 =m_extrapolationTool->propagate(ctx, startParameters_up, *targetSurface_Veto2, Acts::Direction::Backward);
- if (targetParameters_Veto2.has_value()) {
- auto targetPosition_Veto2 = targetParameters_Veto2->position(gctx);
- auto targetMomentum_Veto2 = targetParameters_Veto2->momentum();
- m_xVetoStation2[m_longTracks] = targetPosition_Veto2.x();
- m_yVetoStation2[m_longTracks] = targetPosition_Veto2.y();
- m_thetaxVetoStation2[m_longTracks] = atan(targetMomentum_Veto2[0]/targetMomentum_Veto2[2]);
- m_thetayVetoStation2[m_longTracks] = atan(targetMomentum_Veto2[1]/targetMomentum_Veto2[2]);
- } else {
- ATH_MSG_INFO("veto2 null targetParameters");
- }
+ auto targetSurface_Veto1 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, -1769.65), Acts::Vector3(0, 0, 1)); // -1769.65 mm is z position of center of operational layer in Veto station 1
+ std::optional<const Acts::BoundTrackParameters> targetParameters_Veto1 =m_extrapolationTool->propagate(ctx, startParameters_up, *targetSurface_Veto1, Acts::Direction::Backward);
+
+ if (targetParameters_Veto1) {
+ auto targetPosition_Veto1 = targetParameters_Veto1->position(gctx);
+ auto targetMomentum_Veto1 = targetParameters_Veto1->momentum();
+ m_xVetoStation1[m_longTracks] = targetPosition_Veto1.x();
+ m_yVetoStation1[m_longTracks] = targetPosition_Veto1.y();
+ m_thetaxVetoStation1[m_longTracks] = atan(targetMomentum_Veto1[0]/targetMomentum_Veto1[2]);
+ m_thetayVetoStation1[m_longTracks] = atan(targetMomentum_Veto1[1]/targetMomentum_Veto1[2]);
+ } else {
+ ATH_MSG_INFO("veto1 null targetParameters");
+ }
- auto targetSurface_Trig = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 0.0), Acts::Vector3(0, 0, 1)); // 0 mm is z position of Trig planes overlapping
- std::optional<const Acts::BoundTrackParameters> targetParameters_Trig =m_extrapolationTool->propagate(ctx, startParameters_up, *targetSurface_Trig, Acts::Direction::Backward); // must extrapolate backsards to trig plane if track starts in station 1
- if (targetParameters_Trig.has_value()) {
- auto targetPosition_Trig = targetParameters_Trig->position(gctx);
- auto targetMomentum_Trig = targetParameters_Trig->momentum();
- m_xTrig[m_longTracks] = targetPosition_Trig.x();
- m_yTrig[m_longTracks] = targetPosition_Trig.y();
- m_thetaxTrig[m_longTracks] = atan(targetMomentum_Trig[0]/targetMomentum_Trig[2]);
- m_thetayTrig[m_longTracks] = atan(targetMomentum_Trig[1]/targetMomentum_Trig[2]);
- } else {
- ATH_MSG_INFO("Trig null targetParameters");
- }
+ auto targetSurface_Veto2 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, -1609.65), Acts::Vector3(0, 0, 1)); // -1609.65 mm is z position of where planes touch in Veto station 2
+ std::optional<const Acts::BoundTrackParameters> targetParameters_Veto2 =m_extrapolationTool->propagate(ctx, startParameters_up, *targetSurface_Veto2, Acts::Direction::Backward);
+
+ if (targetParameters_Veto2) {
+ auto targetPosition_Veto2 = targetParameters_Veto2->position(gctx);
+ auto targetMomentum_Veto2 = targetParameters_Veto2->momentum();
+ m_xVetoStation2[m_longTracks] = targetPosition_Veto2.x();
+ m_yVetoStation2[m_longTracks] = targetPosition_Veto2.y();
+ m_thetaxVetoStation2[m_longTracks] = atan(targetMomentum_Veto2[0]/targetMomentum_Veto2[2]);
+ m_thetayVetoStation2[m_longTracks] = atan(targetMomentum_Veto2[1]/targetMomentum_Veto2[2]);
+ } else {
+ ATH_MSG_INFO("veto2 null targetParameters");
+ }
- auto targetSurface_Preshower1 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 2582.68), Acts::Vector3(0, 0, 1)); // 2582.68 mm is z position of center of upstream preshower layer
- std::optional<const Acts::BoundTrackParameters> targetParameters_Preshower1 =m_extrapolationTool->propagate(ctx, startParameters_down, *targetSurface_Preshower1, Acts::Direction::Forward);
- if (targetParameters_Preshower1.has_value()) {
- auto targetPosition_Preshower1 = targetParameters_Preshower1->position(gctx);
- auto targetMomentum_Preshower1 = targetParameters_Preshower1->momentum();
- m_xPreshower1[m_longTracks] = targetPosition_Preshower1.x();
- m_yPreshower1[m_longTracks] = targetPosition_Preshower1.y();
- m_thetaxPreshower1[m_longTracks] = atan(targetMomentum_Preshower1[0]/targetMomentum_Preshower1[2]);
- m_thetayPreshower1[m_longTracks] = atan(targetMomentum_Preshower1[1]/targetMomentum_Preshower1[2]);
- } else {
- ATH_MSG_INFO("Preshower1 null targetParameters");
- }
+ auto targetSurface_Trig = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 0.0), Acts::Vector3(0, 0, 1)); // 0 mm is z position of Trig planes overlapping
+ std::optional<const Acts::BoundTrackParameters> targetParameters_Trig =m_extrapolationTool->propagate(ctx, startParameters_up, *targetSurface_Trig, Acts::Direction::Backward); // must extrapolate backsards to trig plane if track starts in station 1
+
+ if (targetParameters_Trig) {
+ auto targetPosition_Trig = targetParameters_Trig->position(gctx);
+ auto targetMomentum_Trig = targetParameters_Trig->momentum();
+ m_xTrig[m_longTracks] = targetPosition_Trig.x();
+ m_yTrig[m_longTracks] = targetPosition_Trig.y();
+ m_thetaxTrig[m_longTracks] = atan(targetMomentum_Trig[0]/targetMomentum_Trig[2]);
+ m_thetayTrig[m_longTracks] = atan(targetMomentum_Trig[1]/targetMomentum_Trig[2]);
+ } else {
+ ATH_MSG_INFO("Trig null targetParameters");
+ }
- auto targetSurface_Preshower2 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 2657.68), Acts::Vector3(0, 0, 1)); // 2657.68 mm is z position of center of downstream preshower layer
- std::optional<const Acts::BoundTrackParameters> targetParameters_Preshower2 =m_extrapolationTool->propagate(ctx, startParameters_down, *targetSurface_Preshower2, Acts::Direction::Forward);
- if (targetParameters_Preshower2.has_value()) {
- auto targetPosition_Preshower2 = targetParameters_Preshower2->position(gctx);
- auto targetMomentum_Preshower2 = targetParameters_Preshower2->momentum();
- m_xPreshower2[m_longTracks] = targetPosition_Preshower2.x();
- m_yPreshower2[m_longTracks] = targetPosition_Preshower2.y();
- m_thetaxPreshower2[m_longTracks] = atan(targetMomentum_Preshower2[0]/targetMomentum_Preshower2[2]);
- m_thetayPreshower2[m_longTracks] = atan(targetMomentum_Preshower2[1]/targetMomentum_Preshower2[2]);
- } else {
- ATH_MSG_INFO("Preshower2 null targetParameters");
- }
+ auto targetSurface_Preshower1 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 2582.68), Acts::Vector3(0, 0, 1)); // 2582.68 mm is z position of center of upstream preshower layer
+ std::optional<const Acts::BoundTrackParameters> targetParameters_Preshower1 =m_extrapolationTool->propagate(ctx, startParameters_down, *targetSurface_Preshower1, Acts::Direction::Forward);
+
+ if (targetParameters_Preshower1) {
+ auto targetPosition_Preshower1 = targetParameters_Preshower1->position(gctx);
+ auto targetMomentum_Preshower1 = targetParameters_Preshower1->momentum();
+ m_xPreshower1[m_longTracks] = targetPosition_Preshower1.x();
+ m_yPreshower1[m_longTracks] = targetPosition_Preshower1.y();
+ m_thetaxPreshower1[m_longTracks] = atan(targetMomentum_Preshower1[0]/targetMomentum_Preshower1[2]);
+ m_thetayPreshower1[m_longTracks] = atan(targetMomentum_Preshower1[1]/targetMomentum_Preshower1[2]);
+ } else {
+ ATH_MSG_INFO("Preshower1 null targetParameters");
+ }
- auto targetSurface_Calo = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 2760.0), Acts::Vector3(0, 0, 1)); // 2760 mm is estimated z position of calorimeter face
- std::optional<const Acts::BoundTrackParameters> targetParameters_Calo =m_extrapolationTool->propagate(ctx, startParameters_down, *targetSurface_Calo, Acts::Direction::Forward);
- if (targetParameters_Calo.has_value()) {
- auto targetPosition_Calo = targetParameters_Calo->position(gctx);
- auto targetMomentum_Calo = targetParameters_Calo->momentum();
- m_xCalo[m_longTracks] = targetPosition_Calo.x();
- m_yCalo[m_longTracks] = targetPosition_Calo.y();
- m_thetaxCalo[m_longTracks] = atan(targetMomentum_Calo[0]/targetMomentum_Calo[2]) ;
- m_thetayCalo[m_longTracks] = atan(targetMomentum_Calo[1]/targetMomentum_Calo[2]) ;
- } else {
- ATH_MSG_INFO("Calo null targetParameters");
- }
+ auto targetSurface_Preshower2 = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 2657.68), Acts::Vector3(0, 0, 1)); // 2657.68 mm is z position of center of downstream preshower layer
+ std::optional<const Acts::BoundTrackParameters> targetParameters_Preshower2 =m_extrapolationTool->propagate(ctx, startParameters_down, *targetSurface_Preshower2, Acts::Direction::Forward);
+
+ if (targetParameters_Preshower2) {
+ auto targetPosition_Preshower2 = targetParameters_Preshower2->position(gctx);
+ auto targetMomentum_Preshower2 = targetParameters_Preshower2->momentum();
+ m_xPreshower2[m_longTracks] = targetPosition_Preshower2.x();
+ m_yPreshower2[m_longTracks] = targetPosition_Preshower2.y();
+ m_thetaxPreshower2[m_longTracks] = atan(targetMomentum_Preshower2[0]/targetMomentum_Preshower2[2]);
+ m_thetayPreshower2[m_longTracks] = atan(targetMomentum_Preshower2[1]/targetMomentum_Preshower2[2]);
+ } else {
+ ATH_MSG_INFO("Preshower2 null targetParameters");
+ }
- m_longTracks++;
+ auto targetSurface_Calo = Acts::Surface::makeShared<Acts::PlaneSurface>(Acts::Vector3(0, 0, 2760.0), Acts::Vector3(0, 0, 1)); // 2760 mm is estimated z position of calorimeter face
+ std::optional<const Acts::BoundTrackParameters> targetParameters_Calo =m_extrapolationTool->propagate(ctx, startParameters_down, *targetSurface_Calo, Acts::Direction::Forward);
+
+ if (targetParameters_Calo) {
+ auto targetPosition_Calo = targetParameters_Calo->position(gctx);
+ auto targetMomentum_Calo = targetParameters_Calo->momentum();
+ m_xCalo[m_longTracks] = targetPosition_Calo.x();
+ m_yCalo[m_longTracks] = targetPosition_Calo.y();
+ m_thetaxCalo[m_longTracks] = atan(targetMomentum_Calo[0]/targetMomentum_Calo[2]) ;
+ m_thetayCalo[m_longTracks] = atan(targetMomentum_Calo[1]/targetMomentum_Calo[2]) ;
+ } else {
+ ATH_MSG_INFO("Calo null targetParameters");
}
+
+ m_longTracks++;
}
- if (m_runVertexing && goodTracks.size() >= 2) {
- // sort tracks my momentum and reconstruct vertex
- // std::sort(goodTracks.begin(), goodTracks.end(), [](const Trk::Track *lhs, const Trk::Track *rhs){
- // return lhs->trackParameters()->front()->momentum().z() < rhs->trackParameters()->front()->momentum().z();
- // });
- // std::vector<const Trk::Track*> lowMomentumTracks {goodTracks[0], goodTracks[1]};
+ // if (m_runVertexing && goodTracks.size() >= 2) {
+ // // sort tracks my momentum and reconstruct vertex
+ // std::sort(goodTracks.begin(), goodTracks.end(), [](const Trk::Track *lhs, const Trk::Track *rhs){
+ // return lhs->trackParameters()->front()->momentum().z() < rhs->trackParameters()->front()->momentum().z();
+ // });
+ // std::vector<const Trk::Track*> lowMomentumTracks {goodTracks[0], goodTracks[1]};
// std::optional<FaserTracking::POCA> vertex = m_vertexingTool->getVertex(lowMomentumTracks);
// if (vertex) {
@@ -1451,7 +1471,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
// m_vertex_p1 = std::sqrt(m_vertex_px1*m_vertex_px1 + m_vertex_py1*m_vertex_py1 + m_vertex_pz1*m_vertex_pz1);
// }
// }
- }
+ //}
if (!isMC) {
if (m_doTrackFilter) { // filter events: colliding bunches have at least one long track, non-colliding bunches have at least one long track or one calo module with raw_peak > 3 mV
@@ -1469,7 +1489,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
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));
+ m_truthParticleIsFiducial.push_back(m_fiducialParticleTool->isFiducial(tp.first));
}
}
@@ -1873,6 +1893,6 @@ void NtupleDumperAlg::clearTrackTruth() const {
m_isFiducial.push_back(false);
}
-double NtupleDumperAlg::radius(const Amg::Vector3D &position) const {
+double NtupleDumperAlg::radius(const Acts::Vector3 &position) const {
return std::sqrt(position.x() * position.x() + position.y() * position.y());
}
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
index a539248e89bdc5374956af0e151c278d94c7d480..c5051ae797e6eee967abf9de69394b88f86d6650 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
@@ -18,8 +18,8 @@
#include "TrackerSimData/TrackerSimDataCollection.h"
#include "FaserActsGeometryInterfaces/IFaserActsExtrapolationTool.h"
#include "FaserActsGeometryInterfaces/IFaserActsTrackingGeometryTool.h"
-// #include "FaserActsKalmanFilter/IFiducialParticleTool.h"
-// #include "FaserActsKalmanFilter/ITrackTruthMatchingTool.h"
+#include "FaserActsKalmanFilter/IFiducialParticleTool.h"
+#include "FaserActsKalmanFilter/ITrackTruthMatchingTool.h"
#include "TrackerSimEvent/FaserSiHitCollection.h"
#include "xAODEventInfo/EventInfo.h"
#include "FaserActsVertexing/IVertexingTool.h"
@@ -69,7 +69,7 @@ private:
void addWaveBranches(const std::string &name, std::list<int> channel_list) const;
void FillWaveBranches(const xAOD::WaveformHitContainer &wave, bool isMC) const;
void addCalibratedBranches(const std::string &name, int nchannels, int first);
- double radius(const Amg::Vector3D &position) const;
+ double radius(const Acts::Vector3 &position) const;
ServiceHandle <ITHistSvc> m_histSvc;
diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
index bf405a6692704773345842131dfb4dec45e902c9..687908a4956dcc545bc92d2806712d080e8f0fac 100644
--- a/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
+++ b/Tracking/Acts/FaserActsKalmanFilter/src/CKF2.cxx
@@ -87,7 +87,7 @@ StatusCode CKF2::execute() {
const EventContext& ctx = Gaudi::Hive::currentContext();
m_numberOfEvents++;
- SG::WriteDecorHandle<xAOD::EventInfo,uint32_t> eventInfo (m_eventInfoKey, ctx);
+ SG::ReadHandle<xAOD::EventInfo> eventInfo { m_eventInfoKey, ctx };
SG::WriteHandle trackContainer{m_trackCollection, ctx};
std::unique_ptr<TrackCollection> outputTracks = std::make_unique<TrackCollection>();
diff --git a/package_filters.txt b/package_filters.txt
index 82549fe7d8a7e3869d9271fb2801cf75104f8e4f..c3167d0ab1df2c713b94f6d221dd8d7f1b027414 100644
--- a/package_filters.txt
+++ b/package_filters.txt
@@ -142,6 +142,7 @@
+Tracking/Acts/FaserActsGeometryInterfaces
+Tracking/Acts/ActsInterop
+Tracking/Acts/FaserActsKalmanFilter
+#+Tracking/Acts/FaserActsKalmanFilter/FaserActsVertexing
+Tracking/TrkEventCnv/TrkEventAthenaPool
+Tracking/TrkEventCnv/TrkEventCnvTools
+Tracking/TrkEventCnv/TrkEventTopLevelCnv