Skip to content

Commit a9fd2a2

Browse files
bloodyyugobloodyyugo
committed
rotate momentum and positions from tracking to global frame when making Track objects
1 parent 25e45c0 commit a9fd2a2

File tree

5 files changed

+119
-82
lines changed

5 files changed

+119
-82
lines changed

‎Tracking/include/Tracking/Sim/TrackingUtils.h‎

Lines changed: 22 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -164,7 +164,7 @@ inline Acts::BoundSquareMatrix unpackCov(const std::vector<double>& v_cov) {
164164
return cov;
165165
}
166166

167-
// Rotate to ACTS frame
167+
// Rotate from LDMX --> ACTS frame
168168
// z_->x_, x_->y_, y_->z_
169169

170170
//(0 0 1) x_ = z_
@@ -178,7 +178,28 @@ inline Acts::Vector3 ldmx2Acts(Acts::Vector3 ldmx_v) {
178178

179179
return acts_rot * ldmx_v;
180180
}
181+
// Rotate from ACTs --> LDMX frame
182+
// y_->x_, z_->y_, x_->z_
183+
184+
//(0 1 0) x_ = y_
185+
//(0 0 1) y_ = z_
186+
//(1 0 0) z_ = x_
187+
188+
inline Acts::Vector3 acts2Ldmx(Acts::Vector3 acts_v) {
189+
// TODO::Move it to a static member
190+
Acts::SquareMatrix3 ldmx_rot;
191+
ldmx_rot << 0., 1., 0., 0., 0., 1., 1., 0., 0.;
192+
193+
return ldmx_rot * acts_v;
194+
}
181195

196+
inline std::vector<double> acts2LdmxMomentum(Acts::Vector3 acts_mom) {
197+
Acts::Vector3 ldmx_mom=acts2Ldmx(acts_mom);
198+
return std::vector<double>{ldmx_mom[0]/Acts::UnitConstants::MeV,
199+
ldmx_mom[1]/Acts::UnitConstants::MeV,
200+
ldmx_mom[2]/Acts::UnitConstants::MeV};
201+
}
202+
182203
// Transform position, momentum and charge to free parameters
183204

184205
inline Acts::FreeVector toFreeParameters(Acts::Vector3 pos_, Acts::Vector3 mom,

‎Tracking/python/dqm.py‎

Lines changed: 62 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -33,7 +33,7 @@ def __init__(self, name="TrackingRecoDQM"):
3333
super().__init__(name, 'tracking::dqm::TrackingRecoDQM',
3434
'Tracking')
3535

36-
self.nbins = 1000
36+
self.nbins = 200
3737
self.d0min = -50#-15.
3838
self.d0max = 50# 15.
3939
self.z0min = -200#-45.
@@ -44,8 +44,8 @@ def __init__(self, name="TrackingRecoDQM"):
4444
self.thetamax = 1.7#1.7
4545
self.qopmin = -10#-10
4646
self.qopmax = 10#10
47-
self.pmax = 10.
48-
47+
self.pmax = 10000.#10GeV
48+
self.respmax = 1500.
4949
self.trackStates = ["ecal","target"]
5050
self.doTruth= True
5151

@@ -76,8 +76,8 @@ def buildHistograms(self) :
7676
thetamax=self.thetamax
7777
qopmin =self.qopmin
7878
qopmax =self.qopmax
79-
pmax =self.pmax
80-
79+
pmax =self.pmax
80+
respmax =self.respmax
8181

8282
self.build1DHistogram("N_tracks",
8383
"N tracks",10,0,10)
@@ -96,15 +96,15 @@ def buildHistograms(self) :
9696
"#theta [rad]",nbins,thetamin,thetamax)
9797

9898
self.build1DHistogram("p",
99-
"P [GeV]",nbins,0,pmax)
99+
"P [MeV]",nbins,0,pmax)
100100

101101
self.build1DHistogram("qop",
102102
"qOverP [GeV^{-1}]",nbins,qopmin,qopmax)
103103

104104
self.build1DHistogram("pt_bending",
105-
"pT bending plane [GeV]",nbins,0.,pmax)
105+
"pT bending plane [MeV]",nbins,0.,pmax)
106106
self.build1DHistogram("pt_beam",
107-
"pT beam axis [GeV]",nbins,0,0.5)
107+
"pT beam axis [MeV]",nbins,0,500.)
108108

109109
self.build1DHistogram("nHits",
110110
"nHits",15,0,15)
@@ -120,12 +120,20 @@ def buildHistograms(self) :
120120
"nShared",5,0,5)
121121
self.build1DHistogram("nHoles",
122122
"nHoles",5,0,5)
123+
124+
# self.build1DHistogram("px",
125+
# "pX [GeV]",nbins,0.0,pmax)
126+
# self.build1DHistogram("py",
127+
# "pY [GeV]",nbins,-pmax/20.0,pmax/20.0)
128+
# self.build1DHistogram("pz",
129+
# "pZ [GeV]",nbins,-pmax/20.0,pmax/20.0)
130+
123131
self.build1DHistogram("px",
124-
"pX [GeV]",nbins,0.0,pmax)
132+
"pX [MeV]",nbins,-pmax/20.0,pmax/20.0)
125133
self.build1DHistogram("py",
126-
"pY [GeV]",nbins,-pmax/20.0,pmax/20.0)
134+
"pY [MeV]",nbins,-pmax/20.0,pmax/20.0)
127135
self.build1DHistogram("pz",
128-
"pZ [GeV]",nbins,-pmax/20.0,pmax/20.0)
136+
"pZ [MeV]",nbins,0,pmax)
129137
self.build1DHistogram("d0_err",
130138
"#sigma_{d0} [mm]",nbins,0,0.2)
131139
self.build1DHistogram("z0_err",
@@ -137,94 +145,94 @@ def buildHistograms(self) :
137145
self.build1DHistogram("qop_err",
138146
"#sigma_{qOp} [GeV^{-1}]",nbins,0,1)
139147
self.build1DHistogram("p_err",
140-
"#sigma_{p} [GeV]", nbins, 0, 1)
148+
"#sigma_{p} [MeV]", nbins, 0, 1500)
141149

142150

143151
self.build2DHistogram("d0_err_vs_p",
144-
"p [GeV]", nbins, 0, pmax,
152+
"p [MeV]", nbins, 0, pmax,
145153
"#sigma_{d0} [mm]", nbins, 0,0.2)
146154

147155
self.build2DHistogram("z0_err_vs_p",
148-
"p [GeV]", nbins, 0, pmax,
156+
"p [MeV]", nbins, 0, pmax,
149157
"#sigma_{z0} [mm]", nbins, 0,0.8)
150158

151159
self.build2DHistogram("p_err_vs_p",
152-
"p [GeV]", nbins, 0, pmax,
153-
"#sigma_{p} [GeV]", nbins, 0,1)
160+
"p [MeV]", nbins, 0, pmax,
161+
"#sigma_{p} [MeV]", nbins, 0,1500)
154162

155163
self.build2DHistogram("p_err_vs_p_8hits",
156-
"p [GeV]", nbins, 0, pmax,
157-
"#sigma_{p} [GeV]", nbins, 0,1)
164+
"p [MeV]", nbins, 0, pmax,
165+
"#sigma_{p} [MeV]", nbins, 0,1500)
158166

159167
self.build2DHistogram("p_err_vs_p_9hits",
160-
"p [GeV]", nbins, 0, pmax,
161-
"#sigma_{p} [GeV]", nbins, 0,1)
168+
"p [MeV]", nbins, 0, pmax,
169+
"#sigma_{p} [MeV]", nbins, 0,1500)
162170

163171
self.build2DHistogram("p_err_vs_p_10hits",
164-
"p [GeV]", nbins, 0, pmax,
165-
"#sigma_{p} [GeV]", nbins, 0,1)
172+
"p [MeV]", nbins, 0, pmax,
173+
"#sigma_{p} [MeV]", nbins, 0,1500)
166174

167175
self.build2DHistogram("res_p_vs_p",
168-
"p [GeV]", nbins, 0, pmax,
169-
"res_{p} [GeV]", nbins, -3,3)
176+
"p [MeV]", nbins, 0, pmax,
177+
"res_{p} [MeV]", nbins, -respmax,respmax)
170178

171179
self.build2DHistogram("res_qop_vs_p",
172-
"p [GeV]", nbins, 0, pmax,
180+
"p [MeV]", nbins, 0, pmax,
173181
"res_{qop} [GeV^{-1}]", nbins, -0.5,0.5)
174182

175183
self.build2DHistogram("res_d0_vs_p",
176-
"p [GeV]" , nbins, 0, pmax,
184+
"p [MeV]" , nbins, 0, pmax,
177185
"res_{d0} [mm]", nbins, -0.05,0.05)
178186

179187
self.build2DHistogram("res_z0_vs_p",
180-
"p [GeV]" , nbins, 0, pmax,
188+
"p [MeV]" , nbins, 0, pmax,
181189
"res_{z0} [mm]", nbins,-0.5,0.5)
182190

183191
self.build2DHistogram("res_phi_vs_p",
184-
"p [GeV]" , nbins, 0, pmax,
192+
"p [MeV]" , nbins, 0, pmax,
185193
"res_{#phi}", nbins, -0.005,0.005)
186194

187195
self.build2DHistogram("res_theta_vs_p",
188-
"p [GeV]" , nbins, 0, pmax,
196+
"p [MeV]" , nbins, 0, pmax,
189197
"res_{#theta}", nbins, -0.01,0.01)
190198

191199
self.build2DHistogram("res_p_vs_p_8hits",
192-
"p [GeV]", nbins, 0, pmax,
193-
"res_{p} [GeV]", nbins, -3,3)
200+
"p [MeV]", nbins, 0, pmax,
201+
"res_{p} [MeV]", nbins, -respmax,respmax)
194202

195203
self.build2DHistogram("res_p_vs_p_9hits",
196-
"p [GeV]", nbins, 0, pmax,
197-
"res_{p} [GeV]", nbins, -3,3)
204+
"p [MeV]", nbins, 0, pmax,
205+
"res_{p} [MeV]", nbins, -respmax,respmax)
198206

199207
self.build2DHistogram("res_p_vs_p_10hits",
200-
"p [GeV]", nbins, 0, pmax,
201-
"res_{p} [GeV]", nbins, -3,3)
208+
"p [MeV]", nbins, 0, pmax,
209+
"res_{p} [MeV]", nbins, -respmax,respmax)
202210

203211

204212
self.build2DHistogram("pull_qop_vs_p",
205-
"p [GeV]" , nbins, 0, pmax,
213+
"p [MeV]" , nbins, 0, pmax,
206214
"pull_{qop}", nbins, -5,5)
207215

208216
self.build2DHistogram("pull_d0_vs_p",
209-
"p [GeV]" , nbins, 0, pmax,
217+
"p [MeV]" , nbins, 0, pmax,
210218
"pull_{d0}" , nbins, -5,5)
211219

212220
self.build2DHistogram("pull_z0_vs_p",
213-
"p [GeV]" , nbins, 0, pmax,
221+
"p [MeV]" , nbins, 0, pmax,
214222
"pull_{z0}" , nbins,-5,5)
215223

216224
self.build2DHistogram("pull_phi_vs_p",
217-
"p [GeV]" , nbins, 0, pmax,
225+
"p [MeV]" , nbins, 0, pmax,
218226
"pull_{#phi}", nbins, -5,5)
219227

220228
self.build2DHistogram("pull_theta_vs_p",
221-
"p [GeV]" , nbins, 0, pmax,
229+
"p [MeV]" , nbins, 0, pmax,
222230
"pull_{#theta}", nbins, -5,5)
223231

224232

225233
self.build2DHistogram("res_pt_beam_vs_p",
226-
"truth p [GeV]", nbins, 0, pmax,
227-
"res_{pt} beam", nbins, -0.5,0.5)
234+
"truth p [MeV]", nbins, 0, pmax,
235+
"res_{pt} beam", nbins, -respmax,respmax)
228236

229237

230238
if self.doTruth:
@@ -245,7 +253,7 @@ def buildHistograms(self) :
245253
self.build1DHistogram("truth_qop",
246254
"truth q/p [GeV^{-1}]",nbins,qopmin,qopmax)
247255
self.build1DHistogram("truth_p",
248-
"truth p [GeV]",nbins,0,pmax)
256+
"truth p [MeV]",nbins,0,pmax)
249257
self.build1DHistogram("truth_beam_angle",
250258
"angle wrt beam axis",20,0,2)
251259
self.build1DHistogram("truth_PID",
@@ -268,9 +276,9 @@ def buildHistograms(self) :
268276

269277

270278
#self.build1DHistogram("truth_pt_bending",
271-
# "pT bending plane [GeV]",100,-pmax,pmax)
279+
# "pT bending plane [MeV]",100,-pmax,pmax)
272280
#self.build1DHistogram("truth_pt_beam",
273-
# "pT beam axis [GeV]",100,-pmax,pmax)
281+
# "pT beam axis [MeV]",100,-pmax,pmax)
274282

275283
#res
276284
self.build1DHistogram("res_d0",
@@ -282,10 +290,10 @@ def buildHistograms(self) :
282290
self.build1DHistogram("res_theta",
283291
"res #theta", nbins, -0.04,0.04)
284292
self.build1DHistogram("res_p",
285-
"res p [GeV]",nbins,-1,1)
293+
"res p [MeV]",nbins,-500,500)
286294

287295
self.build1DHistogram("res_pt_beam",
288-
"res pt beam [GeV]", nbins, -1,1)
296+
"res pt beam [MeV]", nbins, -500,500)
289297

290298
self.build1DHistogram("res_qop",
291299
"res q/p [GeV^{-1}]",nbins,-1,1)
@@ -319,7 +327,7 @@ def buildHistograms(self) :
319327
self.build1DHistogram("match_qop",
320328
"truth q/p [GeV^{-1}]",nbins,qopmin,qopmax)
321329
self.build1DHistogram("match_p",
322-
"truth p [GeV]", nbins,0,pmax)
330+
"truth p [MeV]", nbins,0,pmax)
323331
self.build1DHistogram("match_beam_angle",
324332
"angle wrt beam axis", 20, 0, 2)
325333
self.build1DHistogram("match_PID",
@@ -362,7 +370,7 @@ def buildHistograms(self) :
362370
self.build1DHistogram("fake_theta",
363371
"fake #theta", nbins, scaling*thetamin,scaling*thetamax)
364372
self.build1DHistogram("fake_p",
365-
"fake p [GeV]",nbins,0,pmax)
373+
"fake p [MeV]",nbins,0,pmax)
366374
self.build1DHistogram("fake_qop",
367375
"fake qOverP [GeV^{-1}]",nbins,-40,40)
368376
self.build1DHistogram("fake_nHits",
@@ -390,7 +398,7 @@ def buildHistograms(self) :
390398
self.build1DHistogram("dup_theta",
391399
"dup #theta", 100, thetamin,thetamax)
392400
self.build1DHistogram("dup_p",
393-
"dup p [GeV]",100,0,pmax)
401+
"dup p [MeV]",100,0,pmax)
394402
self.build1DHistogram("dup_qop",
395403
"dup qOverP [GeV^{-1}]",nbins,qopmin,qopmax)
396404
self.build1DHistogram("dup_nHits",
@@ -424,12 +432,12 @@ def buildHistograms(self) :
424432

425433
self.build2DHistogram(trackState+"_res_loc0-vs-N_hits","N_hits", 5,6.5,11.5,trackState+"_res_loc0 [mm]",100,-0.2,0.2)
426434
self.build2DHistogram(trackState+"_res_loc1-vs-N_hits","N_hits", 5,6.5,11.5,trackState+"_res_loc1 [mm]",100,-5,5)
427-
self.build2DHistogram(trackState+"_res_loc0-vs-trk_p", "trk_p",200,0,5, trackState+"_res_loc0 [mm]",100,-0.2,0.2)
428-
self.build2DHistogram(trackState+"_res_loc1-vs-trk_p", "trk_p",200,0,5, trackState+"_res_loc1 [mm]",100,-5,5)
435+
self.build2DHistogram(trackState+"_res_loc0-vs-trk_p", "trk_p",200,0,pmax, trackState+"_res_loc0 [mm]",100,-0.2,0.2)
436+
self.build2DHistogram(trackState+"_res_loc1-vs-trk_p", "trk_p",200,0,pmax, trackState+"_res_loc1 [mm]",100,-5,5)
429437
self.build2DHistogram(trackState+"_pulls_loc0-vs-N_hits","N_hits",5,6.5,11.5,trackState+"_pulls_loc0 [mm]",100,-3,3)
430438
self.build2DHistogram(trackState+"_pulls_loc1-vs-N_hits","N_hits",5,6.5,11.5,trackState+"_pulls_loc1 [mm]",100,-3,3)
431-
self.build2DHistogram(trackState+"_pulls_loc0-vs-trk_p","trk_p",200,0,5, trackState+"_pulls_loc0 [mm]",100,-3,3)
432-
self.build2DHistogram(trackState+"_pulls_loc1-vs-trk_p","trk_p",200,0,5, trackState+"_pulls_loc1 [mm]",100,-3,3)
439+
self.build2DHistogram(trackState+"_pulls_loc0-vs-trk_p","trk_p",200,0,pmax, trackState+"_pulls_loc0 [mm]",100,-3,3)
440+
self.build2DHistogram(trackState+"_pulls_loc1-vs-trk_p","trk_p",200,0,pmax, trackState+"_pulls_loc1 [mm]",100,-3,3)
433441

434442
class StraightTracksDQM(ldmxcfg.Analyzer):
435443
def __init__(self, name="StraightTracksDQM"):

‎Tracking/src/Tracking/Reco/CKFProcessor.cxx‎

Lines changed: 10 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -502,12 +502,17 @@ void CKFProcessor::produce(framework::Event& event) {
502502
trk.setNdf(track.nMeasurements() - 5);
503503
trk.setNsharedHits(track.nSharedHits());
504504

505-
ldmx_log(debug) << " Track momentum: px = " << track.momentum()[0]
506-
<< " py = " << track.momentum()[1]
507-
<< " pz = " << track.momentum()[2];
508-
trk.setMomentum(track.momentum()[0], track.momentum()[1],
509-
track.momentum()[2]);
505+
Acts::Vector3 trkMomTrking(track.momentum()[0],track.momentum()[1],track.momentum()[2]);
506+
//save the momentum in the global coords and and GeV-->MeV
507+
auto trkMomGbl=tracking::sim::utils::acts2LdmxMomentum(trkMomTrking);
510508

509+
ldmx_log(debug) << " Global Track momentum: px = " << trkMomGbl[0]
510+
<< " py = " << trkMomGbl[1]
511+
<< " pz = " << trkMomGbl[2];
512+
513+
trk.setMomentum(trkMomGbl[0], trkMomGbl[1], trkMomGbl[2]);
514+
515+
511516
// At least min_hits hits and p > 50 MeV
512517
if ((trk.getNhits() <= min_hits_) || (abs(1. / trk.getQoP()) <= 0.05)) {
513518
ldmx_log(debug)

‎Tracking/src/Tracking/Reco/TruthSeedProcessor.cxx‎

Lines changed: 7 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -112,8 +112,8 @@ void TruthSeedProcessor::createTruthTrack(
112112
Acts::Vector3 mom{p_vec[0], p_vec[1], p_vec[2]};
113113

114114
// Rotate the position and momentum into the ACTS frame.
115-
pos = tracking::sim::utils::ldmx2Acts(pos);
116-
mom = tracking::sim::utils::ldmx2Acts(mom);
115+
auto posActs = tracking::sim::utils::ldmx2Acts(pos);
116+
auto momActs = tracking::sim::utils::ldmx2Acts(mom);
117117

118118
// Get the charge of the particle.
119119
// TODO: Add function that uses the PDG ID to calculate this.
@@ -125,7 +125,8 @@ void TruthSeedProcessor::createTruthTrack(
125125
// BoundTrackState there.
126126

127127
// Transform the position, momentum and charge to free parameters.
128-
auto free_params{tracking::sim::utils::toFreeParameters(pos, mom, q)};
128+
// use position & momentum in ACTS frame
129+
auto free_params{tracking::sim::utils::toFreeParameters(posActs, momActs, q)};
129130

130131
// Create a line surface at the perigee. The perigee position is extracted
131132
// from a particle's vertex or the particle's position at a specific
@@ -178,8 +179,10 @@ void TruthSeedProcessor::createTruthTrack(
178179
trk.setPerigeeParameters(
179180
tracking::sim::utils::convertActsToLdmxPars(prop_bound_vec));
180181

182+
//save the global position & momentum in the ldmx::track object
183+
//save momentum in MeV
181184
trk.setPosition(pos(0), pos(1), pos(2));
182-
trk.setMomentum(mom(0), mom(1), mom(2));
185+
trk.setMomentum(mom(0)/Acts::UnitConstants::MeV, mom(1)/Acts::UnitConstants::MeV, mom(2)/Acts::UnitConstants::MeV);
183186
}
184187

185188
// origin_surface is the perigee

0 commit comments

Comments
 (0)