From 6cfd0745b7d5a2e003e2b529ab4655fe4b60f7da Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Mon, 29 Sep 2025 16:46:55 -0700 Subject: [PATCH 1/5] rotate momentum and positions from tracking to global frame when making Track objects --- Tracking/include/Tracking/Sim/TrackingUtils.h | 23 +++- Tracking/python/dqm.py | 116 ++++++++++-------- Tracking/src/Tracking/Reco/CKFProcessor.cxx | 15 ++- .../src/Tracking/Reco/TruthSeedProcessor.cxx | 11 +- Tracking/src/Tracking/dqm/TrackingRecoDQM.cxx | 36 +++--- 5 files changed, 119 insertions(+), 82 deletions(-) diff --git a/Tracking/include/Tracking/Sim/TrackingUtils.h b/Tracking/include/Tracking/Sim/TrackingUtils.h index 6acf5b9f99..bdc4339d23 100644 --- a/Tracking/include/Tracking/Sim/TrackingUtils.h +++ b/Tracking/include/Tracking/Sim/TrackingUtils.h @@ -164,7 +164,7 @@ inline Acts::BoundSquareMatrix unpackCov(const std::vector& v_cov) { return cov; } -// Rotate to ACTS frame +// Rotate from LDMX --> ACTS frame // z_->x_, x_->y_, y_->z_ //(0 0 1) x_ = z_ @@ -178,7 +178,28 @@ inline Acts::Vector3 ldmx2Acts(Acts::Vector3 ldmx_v) { return acts_rot * ldmx_v; } +// Rotate from ACTs --> LDMX frame +// y_->x_, z_->y_, x_->z_ + +//(0 1 0) x_ = y_ +//(0 0 1) y_ = z_ +//(1 0 0) z_ = x_ + +inline Acts::Vector3 acts2Ldmx(Acts::Vector3 acts_v) { + // TODO::Move it to a static member + Acts::SquareMatrix3 ldmx_rot; + ldmx_rot << 0., 1., 0., 0., 0., 1., 1., 0., 0.; + + return ldmx_rot * acts_v; +} +inline std::vector acts2LdmxMomentum(Acts::Vector3 acts_mom) { + Acts::Vector3 ldmx_mom=acts2Ldmx(acts_mom); + return std::vector{ldmx_mom[0]/Acts::UnitConstants::MeV, + ldmx_mom[1]/Acts::UnitConstants::MeV, + ldmx_mom[2]/Acts::UnitConstants::MeV}; +} + // Transform position, momentum and charge to free parameters inline Acts::FreeVector toFreeParameters(Acts::Vector3 pos_, Acts::Vector3 mom, diff --git a/Tracking/python/dqm.py b/Tracking/python/dqm.py index 598fd49eaf..37bedab7af 100644 --- a/Tracking/python/dqm.py +++ b/Tracking/python/dqm.py @@ -33,7 +33,7 @@ def __init__(self, name="TrackingRecoDQM"): super().__init__(name, 'tracking::dqm::TrackingRecoDQM', 'Tracking') - self.nbins = 1000 + self.nbins = 200 self.d0min = -50#-15. self.d0max = 50# 15. self.z0min = -200#-45. @@ -44,8 +44,8 @@ def __init__(self, name="TrackingRecoDQM"): self.thetamax = 1.7#1.7 self.qopmin = -10#-10 self.qopmax = 10#10 - self.pmax = 10. - + self.pmax = 10000.#10GeV + self.respmax = 1500. self.trackStates = ["ecal","target"] self.doTruth= True @@ -76,8 +76,8 @@ def buildHistograms(self) : thetamax=self.thetamax qopmin =self.qopmin qopmax =self.qopmax - pmax =self.pmax - + pmax =self.pmax + respmax =self.respmax self.build1DHistogram("N_tracks", "N tracks",10,0,10) @@ -96,15 +96,15 @@ def buildHistograms(self) : "#theta [rad]",nbins,thetamin,thetamax) self.build1DHistogram("p", - "P [GeV]",nbins,0,pmax) + "P [MeV]",nbins,0,pmax) self.build1DHistogram("qop", "qOverP [GeV^{-1}]",nbins,qopmin,qopmax) self.build1DHistogram("pt_bending", - "pT bending plane [GeV]",nbins,0.,pmax) + "pT bending plane [MeV]",nbins,0.,pmax) self.build1DHistogram("pt_beam", - "pT beam axis [GeV]",nbins,0,0.5) + "pT beam axis [MeV]",nbins,0,500.) self.build1DHistogram("nHits", "nHits",15,0,15) @@ -120,12 +120,20 @@ def buildHistograms(self) : "nShared",5,0,5) self.build1DHistogram("nHoles", "nHoles",5,0,5) + +# self.build1DHistogram("px", +# "pX [GeV]",nbins,0.0,pmax) +# self.build1DHistogram("py", +# "pY [GeV]",nbins,-pmax/20.0,pmax/20.0) +# self.build1DHistogram("pz", +# "pZ [GeV]",nbins,-pmax/20.0,pmax/20.0) + self.build1DHistogram("px", - "pX [GeV]",nbins,0.0,pmax) + "pX [MeV]",nbins,-pmax/20.0,pmax/20.0) self.build1DHistogram("py", - "pY [GeV]",nbins,-pmax/20.0,pmax/20.0) + "pY [MeV]",nbins,-pmax/20.0,pmax/20.0) self.build1DHistogram("pz", - "pZ [GeV]",nbins,-pmax/20.0,pmax/20.0) + "pZ [MeV]",nbins,0,pmax) self.build1DHistogram("d0_err", "#sigma_{d0} [mm]",nbins,0,0.2) self.build1DHistogram("z0_err", @@ -137,94 +145,94 @@ def buildHistograms(self) : self.build1DHistogram("qop_err", "#sigma_{qOp} [GeV^{-1}]",nbins,0,1) self.build1DHistogram("p_err", - "#sigma_{p} [GeV]", nbins, 0, 1) + "#sigma_{p} [MeV]", nbins, 0, 1500) self.build2DHistogram("d0_err_vs_p", - "p [GeV]", nbins, 0, pmax, + "p [MeV]", nbins, 0, pmax, "#sigma_{d0} [mm]", nbins, 0,0.2) self.build2DHistogram("z0_err_vs_p", - "p [GeV]", nbins, 0, pmax, + "p [MeV]", nbins, 0, pmax, "#sigma_{z0} [mm]", nbins, 0,0.8) self.build2DHistogram("p_err_vs_p", - "p [GeV]", nbins, 0, pmax, - "#sigma_{p} [GeV]", nbins, 0,1) + "p [MeV]", nbins, 0, pmax, + "#sigma_{p} [MeV]", nbins, 0,1500) self.build2DHistogram("p_err_vs_p_8hits", - "p [GeV]", nbins, 0, pmax, - "#sigma_{p} [GeV]", nbins, 0,1) + "p [MeV]", nbins, 0, pmax, + "#sigma_{p} [MeV]", nbins, 0,1500) self.build2DHistogram("p_err_vs_p_9hits", - "p [GeV]", nbins, 0, pmax, - "#sigma_{p} [GeV]", nbins, 0,1) + "p [MeV]", nbins, 0, pmax, + "#sigma_{p} [MeV]", nbins, 0,1500) self.build2DHistogram("p_err_vs_p_10hits", - "p [GeV]", nbins, 0, pmax, - "#sigma_{p} [GeV]", nbins, 0,1) + "p [MeV]", nbins, 0, pmax, + "#sigma_{p} [MeV]", nbins, 0,1500) self.build2DHistogram("res_p_vs_p", - "p [GeV]", nbins, 0, pmax, - "res_{p} [GeV]", nbins, -3,3) + "p [MeV]", nbins, 0, pmax, + "res_{p} [MeV]", nbins, -respmax,respmax) self.build2DHistogram("res_qop_vs_p", - "p [GeV]", nbins, 0, pmax, + "p [MeV]", nbins, 0, pmax, "res_{qop} [GeV^{-1}]", nbins, -0.5,0.5) self.build2DHistogram("res_d0_vs_p", - "p [GeV]" , nbins, 0, pmax, + "p [MeV]" , nbins, 0, pmax, "res_{d0} [mm]", nbins, -0.05,0.05) self.build2DHistogram("res_z0_vs_p", - "p [GeV]" , nbins, 0, pmax, + "p [MeV]" , nbins, 0, pmax, "res_{z0} [mm]", nbins,-0.5,0.5) self.build2DHistogram("res_phi_vs_p", - "p [GeV]" , nbins, 0, pmax, + "p [MeV]" , nbins, 0, pmax, "res_{#phi}", nbins, -0.005,0.005) self.build2DHistogram("res_theta_vs_p", - "p [GeV]" , nbins, 0, pmax, + "p [MeV]" , nbins, 0, pmax, "res_{#theta}", nbins, -0.01,0.01) self.build2DHistogram("res_p_vs_p_8hits", - "p [GeV]", nbins, 0, pmax, - "res_{p} [GeV]", nbins, -3,3) + "p [MeV]", nbins, 0, pmax, + "res_{p} [MeV]", nbins, -respmax,respmax) self.build2DHistogram("res_p_vs_p_9hits", - "p [GeV]", nbins, 0, pmax, - "res_{p} [GeV]", nbins, -3,3) + "p [MeV]", nbins, 0, pmax, + "res_{p} [MeV]", nbins, -respmax,respmax) self.build2DHistogram("res_p_vs_p_10hits", - "p [GeV]", nbins, 0, pmax, - "res_{p} [GeV]", nbins, -3,3) + "p [MeV]", nbins, 0, pmax, + "res_{p} [MeV]", nbins, -respmax,respmax) self.build2DHistogram("pull_qop_vs_p", - "p [GeV]" , nbins, 0, pmax, + "p [MeV]" , nbins, 0, pmax, "pull_{qop}", nbins, -5,5) self.build2DHistogram("pull_d0_vs_p", - "p [GeV]" , nbins, 0, pmax, + "p [MeV]" , nbins, 0, pmax, "pull_{d0}" , nbins, -5,5) self.build2DHistogram("pull_z0_vs_p", - "p [GeV]" , nbins, 0, pmax, + "p [MeV]" , nbins, 0, pmax, "pull_{z0}" , nbins,-5,5) self.build2DHistogram("pull_phi_vs_p", - "p [GeV]" , nbins, 0, pmax, + "p [MeV]" , nbins, 0, pmax, "pull_{#phi}", nbins, -5,5) self.build2DHistogram("pull_theta_vs_p", - "p [GeV]" , nbins, 0, pmax, + "p [MeV]" , nbins, 0, pmax, "pull_{#theta}", nbins, -5,5) self.build2DHistogram("res_pt_beam_vs_p", - "truth p [GeV]", nbins, 0, pmax, - "res_{pt} beam", nbins, -0.5,0.5) + "truth p [MeV]", nbins, 0, pmax, + "res_{pt} beam", nbins, -respmax,respmax) if self.doTruth: @@ -245,7 +253,7 @@ def buildHistograms(self) : self.build1DHistogram("truth_qop", "truth q/p [GeV^{-1}]",nbins,qopmin,qopmax) self.build1DHistogram("truth_p", - "truth p [GeV]",nbins,0,pmax) + "truth p [MeV]",nbins,0,pmax) self.build1DHistogram("truth_beam_angle", "angle wrt beam axis",20,0,2) self.build1DHistogram("truth_PID", @@ -268,9 +276,9 @@ def buildHistograms(self) : #self.build1DHistogram("truth_pt_bending", - # "pT bending plane [GeV]",100,-pmax,pmax) + # "pT bending plane [MeV]",100,-pmax,pmax) #self.build1DHistogram("truth_pt_beam", - # "pT beam axis [GeV]",100,-pmax,pmax) + # "pT beam axis [MeV]",100,-pmax,pmax) #res self.build1DHistogram("res_d0", @@ -282,10 +290,10 @@ def buildHistograms(self) : self.build1DHistogram("res_theta", "res #theta", nbins, -0.04,0.04) self.build1DHistogram("res_p", - "res p [GeV]",nbins,-1,1) + "res p [MeV]",nbins,-500,500) self.build1DHistogram("res_pt_beam", - "res pt beam [GeV]", nbins, -1,1) + "res pt beam [MeV]", nbins, -500,500) self.build1DHistogram("res_qop", "res q/p [GeV^{-1}]",nbins,-1,1) @@ -319,7 +327,7 @@ def buildHistograms(self) : self.build1DHistogram("match_qop", "truth q/p [GeV^{-1}]",nbins,qopmin,qopmax) self.build1DHistogram("match_p", - "truth p [GeV]", nbins,0,pmax) + "truth p [MeV]", nbins,0,pmax) self.build1DHistogram("match_beam_angle", "angle wrt beam axis", 20, 0, 2) self.build1DHistogram("match_PID", @@ -362,7 +370,7 @@ def buildHistograms(self) : self.build1DHistogram("fake_theta", "fake #theta", nbins, scaling*thetamin,scaling*thetamax) self.build1DHistogram("fake_p", - "fake p [GeV]",nbins,0,pmax) + "fake p [MeV]",nbins,0,pmax) self.build1DHistogram("fake_qop", "fake qOverP [GeV^{-1}]",nbins,-40,40) self.build1DHistogram("fake_nHits", @@ -390,7 +398,7 @@ def buildHistograms(self) : self.build1DHistogram("dup_theta", "dup #theta", 100, thetamin,thetamax) self.build1DHistogram("dup_p", - "dup p [GeV]",100,0,pmax) + "dup p [MeV]",100,0,pmax) self.build1DHistogram("dup_qop", "dup qOverP [GeV^{-1}]",nbins,qopmin,qopmax) self.build1DHistogram("dup_nHits", @@ -424,12 +432,12 @@ def buildHistograms(self) : self.build2DHistogram(trackState+"_res_loc0-vs-N_hits","N_hits", 5,6.5,11.5,trackState+"_res_loc0 [mm]",100,-0.2,0.2) self.build2DHistogram(trackState+"_res_loc1-vs-N_hits","N_hits", 5,6.5,11.5,trackState+"_res_loc1 [mm]",100,-5,5) - self.build2DHistogram(trackState+"_res_loc0-vs-trk_p", "trk_p",200,0,5, trackState+"_res_loc0 [mm]",100,-0.2,0.2) - self.build2DHistogram(trackState+"_res_loc1-vs-trk_p", "trk_p",200,0,5, trackState+"_res_loc1 [mm]",100,-5,5) + self.build2DHistogram(trackState+"_res_loc0-vs-trk_p", "trk_p",200,0,pmax, trackState+"_res_loc0 [mm]",100,-0.2,0.2) + self.build2DHistogram(trackState+"_res_loc1-vs-trk_p", "trk_p",200,0,pmax, trackState+"_res_loc1 [mm]",100,-5,5) self.build2DHistogram(trackState+"_pulls_loc0-vs-N_hits","N_hits",5,6.5,11.5,trackState+"_pulls_loc0 [mm]",100,-3,3) self.build2DHistogram(trackState+"_pulls_loc1-vs-N_hits","N_hits",5,6.5,11.5,trackState+"_pulls_loc1 [mm]",100,-3,3) - self.build2DHistogram(trackState+"_pulls_loc0-vs-trk_p","trk_p",200,0,5, trackState+"_pulls_loc0 [mm]",100,-3,3) - self.build2DHistogram(trackState+"_pulls_loc1-vs-trk_p","trk_p",200,0,5, trackState+"_pulls_loc1 [mm]",100,-3,3) + self.build2DHistogram(trackState+"_pulls_loc0-vs-trk_p","trk_p",200,0,pmax, trackState+"_pulls_loc0 [mm]",100,-3,3) + self.build2DHistogram(trackState+"_pulls_loc1-vs-trk_p","trk_p",200,0,pmax, trackState+"_pulls_loc1 [mm]",100,-3,3) class StraightTracksDQM(ldmxcfg.Analyzer): def __init__(self, name="StraightTracksDQM"): diff --git a/Tracking/src/Tracking/Reco/CKFProcessor.cxx b/Tracking/src/Tracking/Reco/CKFProcessor.cxx index 2c81814b40..ee913a3185 100644 --- a/Tracking/src/Tracking/Reco/CKFProcessor.cxx +++ b/Tracking/src/Tracking/Reco/CKFProcessor.cxx @@ -502,12 +502,17 @@ void CKFProcessor::produce(framework::Event& event) { trk.setNdf(track.nMeasurements() - 5); trk.setNsharedHits(track.nSharedHits()); - ldmx_log(debug) << " Track momentum: px = " << track.momentum()[0] - << " py = " << track.momentum()[1] - << " pz = " << track.momentum()[2]; - trk.setMomentum(track.momentum()[0], track.momentum()[1], - track.momentum()[2]); + Acts::Vector3 trkMomTrking(track.momentum()[0],track.momentum()[1],track.momentum()[2]); + //save the momentum in the global coords and and GeV-->MeV + auto trkMomGbl=tracking::sim::utils::acts2LdmxMomentum(trkMomTrking); + ldmx_log(debug) << " Global Track momentum: px = " << trkMomGbl[0] + << " py = " << trkMomGbl[1] + << " pz = " << trkMomGbl[2]; + + trk.setMomentum(trkMomGbl[0], trkMomGbl[1], trkMomGbl[2]); + + // At least min_hits hits and p > 50 MeV if ((trk.getNhits() <= min_hits_) || (abs(1. / trk.getQoP()) <= 0.05)) { ldmx_log(debug) diff --git a/Tracking/src/Tracking/Reco/TruthSeedProcessor.cxx b/Tracking/src/Tracking/Reco/TruthSeedProcessor.cxx index 0558701122..c7ec04e313 100644 --- a/Tracking/src/Tracking/Reco/TruthSeedProcessor.cxx +++ b/Tracking/src/Tracking/Reco/TruthSeedProcessor.cxx @@ -112,8 +112,8 @@ void TruthSeedProcessor::createTruthTrack( Acts::Vector3 mom{p_vec[0], p_vec[1], p_vec[2]}; // Rotate the position and momentum into the ACTS frame. - pos = tracking::sim::utils::ldmx2Acts(pos); - mom = tracking::sim::utils::ldmx2Acts(mom); + auto posActs = tracking::sim::utils::ldmx2Acts(pos); + auto momActs = tracking::sim::utils::ldmx2Acts(mom); // Get the charge of the particle. // TODO: Add function that uses the PDG ID to calculate this. @@ -125,7 +125,8 @@ void TruthSeedProcessor::createTruthTrack( // BoundTrackState there. // Transform the position, momentum and charge to free parameters. - auto free_params{tracking::sim::utils::toFreeParameters(pos, mom, q)}; + // use position & momentum in ACTS frame + auto free_params{tracking::sim::utils::toFreeParameters(posActs, momActs, q)}; // Create a line surface at the perigee. The perigee position is extracted // from a particle's vertex or the particle's position at a specific @@ -178,8 +179,10 @@ void TruthSeedProcessor::createTruthTrack( trk.setPerigeeParameters( tracking::sim::utils::convertActsToLdmxPars(prop_bound_vec)); + //save the global position & momentum in the ldmx::track object + //save momentum in MeV trk.setPosition(pos(0), pos(1), pos(2)); - trk.setMomentum(mom(0), mom(1), mom(2)); + trk.setMomentum(mom(0)/Acts::UnitConstants::MeV, mom(1)/Acts::UnitConstants::MeV, mom(2)/Acts::UnitConstants::MeV); } // origin_surface is the perigee diff --git a/Tracking/src/Tracking/dqm/TrackingRecoDQM.cxx b/Tracking/src/Tracking/dqm/TrackingRecoDQM.cxx index 6e3bf32229..2e52fe28d7 100644 --- a/Tracking/src/Tracking/dqm/TrackingRecoDQM.cxx +++ b/Tracking/src/Tracking/dqm/TrackingRecoDQM.cxx @@ -154,7 +154,7 @@ void TrackingRecoDQM::efficiencyPlots( auto truth_z0 = truth_trk.getZ0(); auto truth_theta = truth_trk.getTheta(); auto truth_qop = truth_trk.getQoP(); - auto truth_p = 1. / abs(truth_trk.getQoP()); + auto truth_p = 1. / abs(truth_trk.getQoP())*1000.0; auto truth_n_hits = truth_trk.getNhits(); std::vector truth_mom = truth_trk.getMomentum(); @@ -162,9 +162,9 @@ void TrackingRecoDQM::efficiencyPlots( // Polar angle // The momentum in the plane transverse wrt the beam axis auto truth_pt_beam = - std::sqrt(truth_mom[1] * truth_mom[1] + truth_mom[2] * truth_mom[2]); + std::sqrt(truth_mom[1] * truth_mom[1] + truth_mom[0] * truth_mom[0]); - auto truth_beam_angle = std::atan2(truth_pt_beam, truth_mom[0]); + auto truth_beam_angle = std::atan2(truth_pt_beam, truth_mom[2]); histograms_.fill(title + "truth_nHits", truth_n_hits); histograms_.fill(title + "truth_d0", truth_d0); @@ -235,15 +235,15 @@ void TrackingRecoDQM::efficiencyPlots( auto truth_z0 = truth_trk->getZ0(); auto truth_theta = truth_trk->getTheta(); auto truth_qop = truth_trk->getQoP(); - auto truth_p = 1. / abs(truth_trk->getQoP()); + auto truth_p = 1. / abs(truth_trk->getQoP())*1000.0; std::vector truth_mom = truth_trk->getMomentum(); // Polar angle // The momentum in the plane transverse wrt the beam axis auto truth_pt_beam = - std::sqrt(truth_mom[1] * truth_mom[1] + truth_mom[2] * truth_mom[2]); + std::sqrt(truth_mom[1] * truth_mom[1] + truth_mom[0] * truth_mom[0]); - auto truth_beam_angle = std::atan2(truth_pt_beam, truth_mom[0]); + auto truth_beam_angle = std::atan2(truth_pt_beam, truth_mom[2]); // Fill reco plots for efficiencies - numerator. The quantities are truth histograms_.fill(title + "match_prob", track_truth_prob); @@ -310,7 +310,7 @@ void TrackingRecoDQM::trackMonitoring( auto trk_qop = track.getQoP(); auto trk_theta = track.getTheta(); auto trk_phi = track.getPhi(); - auto trk_p = 1. / abs(trk_qop); + auto trk_p = 1. / abs(trk_qop)*1000.0; for (auto track_hit : track.getMeasurementsIdxs()) { histograms_.fill(title + "LayersHit", measurements.at(track_hit).getLayer()); @@ -320,11 +320,11 @@ void TrackingRecoDQM::trackMonitoring( // The transverse momentum in the bending plane double pt_bending = - std::sqrt(trk_mom[0] * trk_mom[0] + trk_mom[1] * trk_mom[1]); + std::sqrt(trk_mom[1] * trk_mom[1] + trk_mom[2] * trk_mom[2]); // The momentum in the plane transverse wrt the beam axis double pt_beam = - std::sqrt(trk_mom[1] * trk_mom[1] + trk_mom[2] * trk_mom[2]); + std::sqrt(trk_mom[0] * trk_mom[0] + trk_mom[1] * trk_mom[1]); // Covariance matrix Acts::BoundSquareMatrix cov = @@ -340,14 +340,14 @@ void TrackingRecoDQM::trackMonitoring( cov(Acts::BoundIndices::eBoundTheta, Acts::BoundIndices::eBoundTheta)); double sigmaqop = sqrt(cov(Acts::BoundIndices::eBoundQOverP, Acts::BoundIndices::eBoundQOverP)); - double sigmap = (1. / trk_qop) * (1. / trk_qop) * sigmaqop; + double sigmap = (1. / trk_qop) * (1. / trk_qop) * sigmaqop*1000;//convert to MeV histograms_.fill(title + "d0", trk_d0); histograms_.fill(title + "z0", trk_z0); histograms_.fill(title + "qop", trk_qop); histograms_.fill(title + "phi", trk_phi); histograms_.fill(title + "theta", trk_theta); - histograms_.fill(title + "p", std::abs(1. / trk_qop)); + histograms_.fill(title + "p", std::abs(1. / trk_qop)*1000); if (doDetail) { histograms_.fill(title + "px", trk_mom[0]); @@ -372,10 +372,10 @@ void TrackingRecoDQM::trackMonitoring( histograms_.fill(title + "p_err", sigmap); // 2D Error plots - double p = std::abs(1. / trk_qop); + double p = std::abs(1. / trk_qop)*1000.0; histograms_.fill(title + "d0_err_vs_p", p, sigmad0); - histograms_.fill(title + "z0_err_vs_p", std::abs(1. / trk_qop), sigmaz0); - histograms_.fill(title + "p_err_vs_p", std::abs(1. / trk_qop), sigmap); + histograms_.fill(title + "z0_err_vs_p", p, sigmaz0); + histograms_.fill(title + "p_err_vs_p", p, sigmap); if (track.getNhits() == 8) histograms_.fill(title + "p_err_vs_p_8hits", p, sigmap); @@ -408,12 +408,12 @@ void TrackingRecoDQM::trackMonitoring( auto truth_phi = truth_trk->getPhi(); auto truth_theta = truth_trk->getTheta(); auto truth_qop = truth_trk->getQoP(); - auto truth_p = 1. / abs(truth_trk->getQoP()); + auto truth_p = 1. / abs(truth_trk->getQoP())*1000.0; std::vector truth_mom = truth_trk->getMomentum(); // Polar angle // The momentum in the plane transverse wrt the beam axis double truth_pt_beam = std::sqrt(truth_mom[1] * truth_mom[1] + - truth_mom[2] * truth_mom[2]); + truth_mom[0] * truth_mom[0]); // histograms_.fill(title+"truth_d0", truth_d0); // histograms_.fill(title+"truth_z0", truth_z0); @@ -475,7 +475,7 @@ void TrackingRecoDQM::trackMonitoring( else if (track.getNhits() == 10) histograms_.fill(title + "res_p_vs_p_10hits", truth_p, res_p); - histograms_.fill(title + "res_pt_beam_vs_p", truth_pt_beam, + histograms_.fill(title + "res_pt_beam_vs_p", truth_p, res_pt_beam); } // found matched track @@ -534,7 +534,7 @@ void TrackingRecoDQM::trackStateMonitoring(const ldmx::Tracks& tracks, Acts::BoundIndices::eBoundQOverP, Acts::BoundIndices::eBoundQOverP)); double trk_qop = track.getQoP(); - double trk_p = 1. / abs(trk_qop); + double trk_p = 1. / abs(trk_qop)*1000.0; double track_state_loc0 = target_state.params_[0]; double track_state_loc1 = target_state.params_[1]; From 8b6337085b9d78858c8059bf2ddddefb7d5ece5a Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Mon, 29 Sep 2025 16:56:20 -0700 Subject: [PATCH 2/5] tidy up --- Tracking/src/Tracking/Reco/CKFProcessor.cxx | 12 ++++++------ Tracking/src/Tracking/Reco/TruthSeedProcessor.cxx | 6 +++--- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/Tracking/src/Tracking/Reco/CKFProcessor.cxx b/Tracking/src/Tracking/Reco/CKFProcessor.cxx index ee913a3185..846c3176aa 100644 --- a/Tracking/src/Tracking/Reco/CKFProcessor.cxx +++ b/Tracking/src/Tracking/Reco/CKFProcessor.cxx @@ -502,15 +502,15 @@ void CKFProcessor::produce(framework::Event& event) { trk.setNdf(track.nMeasurements() - 5); trk.setNsharedHits(track.nSharedHits()); - Acts::Vector3 trkMomTrking(track.momentum()[0],track.momentum()[1],track.momentum()[2]); + Acts::Vector3 trk_mom_trking(track.momentum()[0],track.momentum()[1],track.momentum()[2]); //save the momentum in the global coords and and GeV-->MeV - auto trkMomGbl=tracking::sim::utils::acts2LdmxMomentum(trkMomTrking); + auto trk_mom_gbl=tracking::sim::utils::acts2LdmxMomentum(trk_mom_trking); - ldmx_log(debug) << " Global Track momentum: px = " << trkMomGbl[0] - << " py = " << trkMomGbl[1] - << " pz = " << trkMomGbl[2]; + ldmx_log(debug) << " Global Track momentum: px = " << trk_mom_gbl[0] + << " py = " << trk_mom_gbl[1] + << " pz = " << trk_mom_gbl[2]; - trk.setMomentum(trkMomGbl[0], trkMomGbl[1], trkMomGbl[2]); + trk.setMomentum(trk_mom_gbl[0], trk_mom_gbl[1], trk_mom_gbl[2]); // At least min_hits hits and p > 50 MeV diff --git a/Tracking/src/Tracking/Reco/TruthSeedProcessor.cxx b/Tracking/src/Tracking/Reco/TruthSeedProcessor.cxx index c7ec04e313..b2368c3a99 100644 --- a/Tracking/src/Tracking/Reco/TruthSeedProcessor.cxx +++ b/Tracking/src/Tracking/Reco/TruthSeedProcessor.cxx @@ -112,8 +112,8 @@ void TruthSeedProcessor::createTruthTrack( Acts::Vector3 mom{p_vec[0], p_vec[1], p_vec[2]}; // Rotate the position and momentum into the ACTS frame. - auto posActs = tracking::sim::utils::ldmx2Acts(pos); - auto momActs = tracking::sim::utils::ldmx2Acts(mom); + auto pos_acts = tracking::sim::utils::ldmx2Acts(pos); + auto mom_acts = tracking::sim::utils::ldmx2Acts(mom); // Get the charge of the particle. // TODO: Add function that uses the PDG ID to calculate this. @@ -126,7 +126,7 @@ void TruthSeedProcessor::createTruthTrack( // Transform the position, momentum and charge to free parameters. // use position & momentum in ACTS frame - auto free_params{tracking::sim::utils::toFreeParameters(posActs, momActs, q)}; + auto free_params{tracking::sim::utils::toFreeParameters(pos_acts, mom_acts, q)}; // Create a line surface at the perigee. The perigee position is extracted // from a particle's vertex or the particle's position at a specific From 1fb1fe8b2951837083ae6e04ba93066304d8ae50 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Mon, 29 Sep 2025 17:08:31 -0700 Subject: [PATCH 3/5] rotate momentum from GSF --- Tracking/src/Tracking/Reco/GSFProcessor.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Tracking/src/Tracking/Reco/GSFProcessor.cxx b/Tracking/src/Tracking/Reco/GSFProcessor.cxx index 4e9a86ee52..8547a61993 100644 --- a/Tracking/src/Tracking/Reco/GSFProcessor.cxx +++ b/Tracking/src/Tracking/Reco/GSFProcessor.cxx @@ -439,7 +439,8 @@ void GSFProcessor::produce(framework::Event& event) { tracking::sim::utils::flatCov(trk_cov, v_trk_cov); trk.setPerigeeCov(v_trk_cov); Acts::Vector3 trk_momentum = gsftrk.momentum(); - trk.setMomentum(trk_momentum(0), trk_momentum(1), trk_momentum(2)); + auto trk_mom_gbl=tracking::sim::utils::acts2LdmxMomentum(trk_momentum); + trk.setMomentum(trk_mom_gbl[0], trk_mom_gbl[1], trk_mom_gbl[2]); // truth information trk.setTrackID(track.getTrackID()); From 3bced5c0581dd240fa49524e6bbd84a74c5f08e3 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" Date: Tue, 30 Sep 2025 00:12:02 +0000 Subject: [PATCH 4/5] Apply clang-format --- Tracking/include/Tracking/Sim/TrackingUtils.h | 12 +++++------ Tracking/src/Tracking/Reco/CKFProcessor.cxx | 11 +++++----- Tracking/src/Tracking/Reco/GSFProcessor.cxx | 4 ++-- .../src/Tracking/Reco/TruthSeedProcessor.cxx | 11 ++++++---- Tracking/src/Tracking/dqm/TrackingRecoDQM.cxx | 20 +++++++++---------- 5 files changed, 31 insertions(+), 27 deletions(-) diff --git a/Tracking/include/Tracking/Sim/TrackingUtils.h b/Tracking/include/Tracking/Sim/TrackingUtils.h index bdc4339d23..ca3190e6a4 100644 --- a/Tracking/include/Tracking/Sim/TrackingUtils.h +++ b/Tracking/include/Tracking/Sim/TrackingUtils.h @@ -184,7 +184,7 @@ inline Acts::Vector3 ldmx2Acts(Acts::Vector3 ldmx_v) { //(0 1 0) x_ = y_ //(0 0 1) y_ = z_ //(1 0 0) z_ = x_ - + inline Acts::Vector3 acts2Ldmx(Acts::Vector3 acts_v) { // TODO::Move it to a static member Acts::SquareMatrix3 ldmx_rot; @@ -194,12 +194,12 @@ inline Acts::Vector3 acts2Ldmx(Acts::Vector3 acts_v) { } inline std::vector acts2LdmxMomentum(Acts::Vector3 acts_mom) { - Acts::Vector3 ldmx_mom=acts2Ldmx(acts_mom); - return std::vector{ldmx_mom[0]/Acts::UnitConstants::MeV, - ldmx_mom[1]/Acts::UnitConstants::MeV, - ldmx_mom[2]/Acts::UnitConstants::MeV}; + Acts::Vector3 ldmx_mom = acts2Ldmx(acts_mom); + return std::vector{ldmx_mom[0] / Acts::UnitConstants::MeV, + ldmx_mom[1] / Acts::UnitConstants::MeV, + ldmx_mom[2] / Acts::UnitConstants::MeV}; } - + // Transform position, momentum and charge to free parameters inline Acts::FreeVector toFreeParameters(Acts::Vector3 pos_, Acts::Vector3 mom, diff --git a/Tracking/src/Tracking/Reco/CKFProcessor.cxx b/Tracking/src/Tracking/Reco/CKFProcessor.cxx index 846c3176aa..d56e7690d2 100644 --- a/Tracking/src/Tracking/Reco/CKFProcessor.cxx +++ b/Tracking/src/Tracking/Reco/CKFProcessor.cxx @@ -502,17 +502,18 @@ void CKFProcessor::produce(framework::Event& event) { trk.setNdf(track.nMeasurements() - 5); trk.setNsharedHits(track.nSharedHits()); - Acts::Vector3 trk_mom_trking(track.momentum()[0],track.momentum()[1],track.momentum()[2]); - //save the momentum in the global coords and and GeV-->MeV - auto trk_mom_gbl=tracking::sim::utils::acts2LdmxMomentum(trk_mom_trking); + Acts::Vector3 trk_mom_trking(track.momentum()[0], track.momentum()[1], + track.momentum()[2]); + // save the momentum in the global coords and and GeV-->MeV + auto trk_mom_gbl = + tracking::sim::utils::acts2LdmxMomentum(trk_mom_trking); ldmx_log(debug) << " Global Track momentum: px = " << trk_mom_gbl[0] << " py = " << trk_mom_gbl[1] << " pz = " << trk_mom_gbl[2]; - trk.setMomentum(trk_mom_gbl[0], trk_mom_gbl[1], trk_mom_gbl[2]); + trk.setMomentum(trk_mom_gbl[0], trk_mom_gbl[1], trk_mom_gbl[2]); - // At least min_hits hits and p > 50 MeV if ((trk.getNhits() <= min_hits_) || (abs(1. / trk.getQoP()) <= 0.05)) { ldmx_log(debug) diff --git a/Tracking/src/Tracking/Reco/GSFProcessor.cxx b/Tracking/src/Tracking/Reco/GSFProcessor.cxx index 8547a61993..02c7795b18 100644 --- a/Tracking/src/Tracking/Reco/GSFProcessor.cxx +++ b/Tracking/src/Tracking/Reco/GSFProcessor.cxx @@ -439,8 +439,8 @@ void GSFProcessor::produce(framework::Event& event) { tracking::sim::utils::flatCov(trk_cov, v_trk_cov); trk.setPerigeeCov(v_trk_cov); Acts::Vector3 trk_momentum = gsftrk.momentum(); - auto trk_mom_gbl=tracking::sim::utils::acts2LdmxMomentum(trk_momentum); - trk.setMomentum(trk_mom_gbl[0], trk_mom_gbl[1], trk_mom_gbl[2]); + auto trk_mom_gbl = tracking::sim::utils::acts2LdmxMomentum(trk_momentum); + trk.setMomentum(trk_mom_gbl[0], trk_mom_gbl[1], trk_mom_gbl[2]); // truth information trk.setTrackID(track.getTrackID()); diff --git a/Tracking/src/Tracking/Reco/TruthSeedProcessor.cxx b/Tracking/src/Tracking/Reco/TruthSeedProcessor.cxx index b2368c3a99..74030b6dd3 100644 --- a/Tracking/src/Tracking/Reco/TruthSeedProcessor.cxx +++ b/Tracking/src/Tracking/Reco/TruthSeedProcessor.cxx @@ -126,7 +126,8 @@ void TruthSeedProcessor::createTruthTrack( // Transform the position, momentum and charge to free parameters. // use position & momentum in ACTS frame - auto free_params{tracking::sim::utils::toFreeParameters(pos_acts, mom_acts, q)}; + auto free_params{ + tracking::sim::utils::toFreeParameters(pos_acts, mom_acts, q)}; // Create a line surface at the perigee. The perigee position is extracted // from a particle's vertex or the particle's position at a specific @@ -179,10 +180,12 @@ void TruthSeedProcessor::createTruthTrack( trk.setPerigeeParameters( tracking::sim::utils::convertActsToLdmxPars(prop_bound_vec)); - //save the global position & momentum in the ldmx::track object - //save momentum in MeV + // save the global position & momentum in the ldmx::track object + // save momentum in MeV trk.setPosition(pos(0), pos(1), pos(2)); - trk.setMomentum(mom(0)/Acts::UnitConstants::MeV, mom(1)/Acts::UnitConstants::MeV, mom(2)/Acts::UnitConstants::MeV); + trk.setMomentum(mom(0) / Acts::UnitConstants::MeV, + mom(1) / Acts::UnitConstants::MeV, + mom(2) / Acts::UnitConstants::MeV); } // origin_surface is the perigee diff --git a/Tracking/src/Tracking/dqm/TrackingRecoDQM.cxx b/Tracking/src/Tracking/dqm/TrackingRecoDQM.cxx index 2e52fe28d7..52d4ab4561 100644 --- a/Tracking/src/Tracking/dqm/TrackingRecoDQM.cxx +++ b/Tracking/src/Tracking/dqm/TrackingRecoDQM.cxx @@ -154,7 +154,7 @@ void TrackingRecoDQM::efficiencyPlots( auto truth_z0 = truth_trk.getZ0(); auto truth_theta = truth_trk.getTheta(); auto truth_qop = truth_trk.getQoP(); - auto truth_p = 1. / abs(truth_trk.getQoP())*1000.0; + auto truth_p = 1. / abs(truth_trk.getQoP()) * 1000.0; auto truth_n_hits = truth_trk.getNhits(); std::vector truth_mom = truth_trk.getMomentum(); @@ -235,7 +235,7 @@ void TrackingRecoDQM::efficiencyPlots( auto truth_z0 = truth_trk->getZ0(); auto truth_theta = truth_trk->getTheta(); auto truth_qop = truth_trk->getQoP(); - auto truth_p = 1. / abs(truth_trk->getQoP())*1000.0; + auto truth_p = 1. / abs(truth_trk->getQoP()) * 1000.0; std::vector truth_mom = truth_trk->getMomentum(); // Polar angle @@ -310,7 +310,7 @@ void TrackingRecoDQM::trackMonitoring( auto trk_qop = track.getQoP(); auto trk_theta = track.getTheta(); auto trk_phi = track.getPhi(); - auto trk_p = 1. / abs(trk_qop)*1000.0; + auto trk_p = 1. / abs(trk_qop) * 1000.0; for (auto track_hit : track.getMeasurementsIdxs()) { histograms_.fill(title + "LayersHit", measurements.at(track_hit).getLayer()); @@ -340,14 +340,15 @@ void TrackingRecoDQM::trackMonitoring( cov(Acts::BoundIndices::eBoundTheta, Acts::BoundIndices::eBoundTheta)); double sigmaqop = sqrt(cov(Acts::BoundIndices::eBoundQOverP, Acts::BoundIndices::eBoundQOverP)); - double sigmap = (1. / trk_qop) * (1. / trk_qop) * sigmaqop*1000;//convert to MeV + double sigmap = + (1. / trk_qop) * (1. / trk_qop) * sigmaqop * 1000; // convert to MeV histograms_.fill(title + "d0", trk_d0); histograms_.fill(title + "z0", trk_z0); histograms_.fill(title + "qop", trk_qop); histograms_.fill(title + "phi", trk_phi); histograms_.fill(title + "theta", trk_theta); - histograms_.fill(title + "p", std::abs(1. / trk_qop)*1000); + histograms_.fill(title + "p", std::abs(1. / trk_qop) * 1000); if (doDetail) { histograms_.fill(title + "px", trk_mom[0]); @@ -372,7 +373,7 @@ void TrackingRecoDQM::trackMonitoring( histograms_.fill(title + "p_err", sigmap); // 2D Error plots - double p = std::abs(1. / trk_qop)*1000.0; + double p = std::abs(1. / trk_qop) * 1000.0; histograms_.fill(title + "d0_err_vs_p", p, sigmad0); histograms_.fill(title + "z0_err_vs_p", p, sigmaz0); histograms_.fill(title + "p_err_vs_p", p, sigmap); @@ -408,7 +409,7 @@ void TrackingRecoDQM::trackMonitoring( auto truth_phi = truth_trk->getPhi(); auto truth_theta = truth_trk->getTheta(); auto truth_qop = truth_trk->getQoP(); - auto truth_p = 1. / abs(truth_trk->getQoP())*1000.0; + auto truth_p = 1. / abs(truth_trk->getQoP()) * 1000.0; std::vector truth_mom = truth_trk->getMomentum(); // Polar angle // The momentum in the plane transverse wrt the beam axis @@ -475,8 +476,7 @@ void TrackingRecoDQM::trackMonitoring( else if (track.getNhits() == 10) histograms_.fill(title + "res_p_vs_p_10hits", truth_p, res_p); - histograms_.fill(title + "res_pt_beam_vs_p", truth_p, - res_pt_beam); + histograms_.fill(title + "res_pt_beam_vs_p", truth_p, res_pt_beam); } // found matched track } // do TruthComparison @@ -534,7 +534,7 @@ void TrackingRecoDQM::trackStateMonitoring(const ldmx::Tracks& tracks, Acts::BoundIndices::eBoundQOverP, Acts::BoundIndices::eBoundQOverP)); double trk_qop = track.getQoP(); - double trk_p = 1. / abs(trk_qop)*1000.0; + double trk_p = 1. / abs(trk_qop) * 1000.0; double track_state_loc0 = target_state.params_[0]; double track_state_loc1 = target_state.params_[1]; From 56ad1f9b2ee3152f0258f14940c337a1b3a2b50b Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Tue, 30 Sep 2025 07:37:17 -0700 Subject: [PATCH 5/5] increment class def --- Tracking/include/Tracking/Event/Track.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Tracking/include/Tracking/Event/Track.h b/Tracking/include/Tracking/Event/Track.h index 845673ed31..3ffeb48196 100644 --- a/Tracking/include/Tracking/Event/Track.h +++ b/Tracking/include/Tracking/Event/Track.h @@ -270,7 +270,7 @@ class Track { std::vector track_states_; /// Class declaration needed by the ROOT dictionary. - ClassDef(Track, 4); + ClassDef(Track, 5); }; // Track