From de635472dda42f7558e6455adbf0a802a85b5889 Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Sat, 22 Jun 2024 15:14:10 -0700
Subject: [PATCH] Working alma9 NtupleDumper
---
.../NtupleDumper/src/NtupleDumperAlg.cxx | 579 +++++++++---------
1 file changed, 291 insertions(+), 288 deletions(-)
diff --git a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
index 77043fa6..edcfbb03 100644
--- a/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
+++ b/PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
@@ -1110,309 +1110,312 @@ 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 {};
- for (const Trk::Track* track : *trackCollection)
+ if (trackCollection.isValid())
{
- 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));
+ 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();
}
- }
- 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());
- // }
-
- // 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();
- // }
- } else {
- clearTrackTruth();
- }
+ 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());
+ // }
- // 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 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 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");
- }
+ // 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();
+ }
- 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");
- }
+ // 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_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_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");
+ }
- 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_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_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_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_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_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");
+ }
- m_longTracks++;
+ 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_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");
+ }
+
+ m_longTracks++;
+ }
}
if (m_runVertexing && goodTracks.size() >= 2) {
--
GitLab