Skip to content
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Tracking/include/Tracking/Event/Track.h
Original file line number Diff line number Diff line change
Expand Up @@ -270,7 +270,7 @@ class Track {
std::vector<TrackState> track_states_;

/// Class declaration needed by the ROOT dictionary.
ClassDef(Track, 4);
ClassDef(Track, 5);

}; // Track

Expand Down
23 changes: 22 additions & 1 deletion Tracking/include/Tracking/Sim/TrackingUtils.h
Original file line number Diff line number Diff line change
Expand Up @@ -164,7 +164,7 @@ inline Acts::BoundSquareMatrix unpackCov(const std::vector<double>& v_cov) {
return cov;
}

// Rotate to ACTS frame
// Rotate from LDMX --> ACTS frame
// z_->x_, x_->y_, y_->z_

//(0 0 1) x_ = z_
Expand All @@ -178,6 +178,27 @@ 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<double> acts2LdmxMomentum(Acts::Vector3 acts_mom) {
Acts::Vector3 ldmx_mom = acts2Ldmx(acts_mom);
return std::vector<double>{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

Expand Down
116 changes: 62 additions & 54 deletions Tracking/python/dqm.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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

Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand All @@ -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",
Expand All @@ -137,94 +145,94 @@ def buildHistograms(self) :
self.build1DHistogram("qop_err",
"#sigma_{qOp} [GeV^{-1}]",nbins,0,1)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
"#sigma_{qOp} [GeV^{-1}]",nbins,0,1)
"#sigma_{qOp} [MeV^{-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:
Expand All @@ -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",
Expand All @@ -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",
Expand All @@ -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)
Expand Down Expand Up @@ -319,7 +327,7 @@ def buildHistograms(self) :
self.build1DHistogram("match_qop",
"truth q/p [GeV^{-1}]",nbins,qopmin,qopmax)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
"truth q/p [GeV^{-1}]",nbins,qopmin,qopmax)
"truth q/p [MeV^{-1}]",nbins,qopmin,qopmax)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm keeping q/p (which is a track fit parameter) in GeV so that it aligns with ACTs.

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",
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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"):
Expand Down
16 changes: 11 additions & 5 deletions Tracking/src/Tracking/Reco/CKFProcessor.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -502,11 +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 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]);

// At least min_hits hits and p > 50 MeV
if ((trk.getNhits() <= min_hits_) || (abs(1. / trk.getQoP()) <= 0.05)) {
Expand Down
3 changes: 2 additions & 1 deletion Tracking/src/Tracking/Reco/GSFProcessor.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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());
Expand Down
14 changes: 10 additions & 4 deletions Tracking/src/Tracking/Reco/TruthSeedProcessor.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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 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.
Expand All @@ -125,7 +125,9 @@ 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(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
Expand Down Expand Up @@ -178,8 +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
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
Expand Down
Loading