Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
//*************************************************************************
//* =====================
//* TKalDetCradle Class
//* =====================
//*
//* (Description)
//* A singleton class to hold information of detector system
//* used in Kalman filter classes.
//* (Requires)
//* TObjArray
//* TVKalDetector
//* (Provides)
//* class TKalDetCradle
//* (Update Recored)
//* 2005/02/23 A.Yamaguchi Original version.
//* 2005/08/14 K.Fujii Removed CalcTable(), GetMeasLayerTable(),
//* GetPhiTable(), and GetDir() and added
//* Transport() to do their functions.
//* 2010/04/06 K.Fujii Modified Transport() to allow a 1-dim hit,
//* for which pivot is at the expected hit.
//* 2012/11/29 K.Fujii Moved GetEnergyLoss and CalcQms from
//* TKalDetCradle.
//*
//*************************************************************************
#include "TKalDetCradle.h" // from KalTrackLib
#include "TVMeasLayer.h" // from KalTrackLib
#include "TVKalDetector.h" // from KalTrackLib
#include "TKalTrackSite.h" // from KalTrackLib
#include "TKalTrackState.h" // from KalTrackLib
#include "TKalTrack.h" // from KalTrackLib
#include "TVSurface.h" // from GeomLib
#include "TVTrack.h" // from GeomLib
#include "TBField.h" // from Bfield
#include "TRungeKuttaTrack.h"
#include <iostream> // from STL
ClassImp(TKalDetCradle)
Bool_t TKalDetCradle::fUseRKTrack= kFALSE;
//_________________________________________________________________________
// ----------------------------------
// Ctors and Dtor
// ----------------------------------
TKalDetCradle::TKalDetCradle(Int_t n)
: TObjArray(n), fIsMSON(kTRUE), fIsDEDXON(kTRUE),
fDone(kFALSE), fIsClosed(kFALSE)
{
}
TKalDetCradle::~TKalDetCradle()
{
//std::cout << "TKalDetCradle::~TKalDetCradle() " << this << " " << GetEntries() << std::endl;
std::map<TAttElement*, int> det_nelement;
TIter next(this);
TObject *mlp = 0;
while ((mlp = next())) {
TAttElement* det = const_cast<TAttElement*>(&(dynamic_cast<TAttElement *>(mlp)->GetParent(kFALSE)));
if(det_nelement.find(det)!=det_nelement.end()) det_nelement[det]++;
else det_nelement[det] = 1;
}
for (auto it : det_nelement) {
delete it.first;
}
}
//_________________________________________________________________________
// ----------------------------------
// Utility Methods
// ----------------------------------
//_________________________________________________________________________
// -----------------
// Install
// -----------------
// installs a sub-detector into this cradle.
//
void TKalDetCradle::Install(TVKalDetector &det)
{
if (IsClosed()) {
std::cerr << ">>>> Error!! >>>> TKalDetCradle::Install" << std::endl
<< " Cradle already closed. Abort!!" << std::endl;
abort();
}
TIter next(&det);
TObject *mlp = 0; // measment layer pointer
while ((mlp = next())) {
Add(mlp);
dynamic_cast<TAttElement *>(mlp)->SetParentPtr(&det);
det.SetParentPtr(this);
}
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
fDone = kFALSE;
}
void TKalDetCradle::Transport(const TKalTrackSite &from, // site from
TKalTrackSite &to, // site to
TKalMatrix &sv, // state vector
TKalMatrix &F, // propagator matrix
TKalMatrix &Q) // process noise matrix
{
// ---------------------------------------------------------------------
// Sort measurement layers in this cradle if not
// ---------------------------------------------------------------------
if (!fDone) Update();
// ---------------------------------------------------------------------
// Move to site "to"
// ---------------------------------------------------------------------
std::unique_ptr<TVTrack> help(&static_cast<TKalTrackState &>
(from.GetCurState()).CreateTrack()); // tmp track
const TVMeasLayer& ml_to = to.GetHit().GetMeasLayer() ;
TVector3 x0; // local pivot at the "to" site
if(fUseRKTrack) {
this->Transport2(from, ml_to, x0, sv, F, Q, help);
}
else {
this->Transport(from, ml_to, x0, sv, F, Q, help);
}
TVTrack &hel = *help;
// ---------------------------------------------------------------------
// Move pivot from last expected hit to actural hit at site "to"
// ---------------------------------------------------------------------
if (to.GetDimension() > 1) {
double fid = 0.;
Int_t sdim = sv.GetNrows(); // number of track parameters
if(to.GetHit().GetBfield()!=0.) {
TKalMatrix DF(sdim, sdim); // propagator matrix segment
hel.MoveTo(to.GetGlobalPivot(), fid, &DF, 0, kFALSE); // move pivot to actual hit (to)
F = DF * F; // update F accordingly
hel.PutInto(sv); // save updated hel to sv
}
else {
TKalMatrix DF(sdim, sdim);
TStraightTrack strtrk(sv, x0);
strtrk.MoveTo(to.GetPivot(), fid, &DF); // move pivot to actual hit (to)
F = DF * F; // update F accordingly
strtrk.PutInto(sv); // save updated hel to sv
}
} else {
TVector3 x0g = x0;
if(!TBField::IsUsingUniformBfield()) {
x0g = (hel.GetFrame()).Transform(x0, TTrackFrame::kLocalToGlobal);
}
to.SetPivot(x0g); // if it is a 1-dim hit
}
if(!TBField::IsUsingUniformBfield()) {
to.SetFrame(hel.GetFrame()); // set frame
to.SetBfield(hel.GetMagField()); // set B field
}
}
int TKalDetCradle::Transport(const TKalTrackSite &from, // site from
const TVMeasLayer &ml_to, // layer to reach
TVector3 &x0, // pivot for sv
TKalMatrix &sv, // state vector
TKalMatrix &F, // propagator matrix
TKalMatrix &Q) // process noise matrix
{
if (!fDone) Update();
std::unique_ptr<TVTrack> help(&static_cast<TKalTrackState &>
(from.GetCurState()).CreateTrack()); // tmp track
return Transport(from, ml_to, x0, sv, F, Q, help);
}
//
//
//
//_________________________________________________________________________
// -----------------
// Transport
// -----------------
// transports state (sv) from site (from) to layer (ml_to), taking into
// account multiple scattering and energy loss and updates state (sv),
// fills pivot in x0, propagator matrix (F), and process noise matrix (Q).
//
int TKalDetCradle::Transport(const TKalTrackSite &from, // site from
const TVMeasLayer &ml_to, // layer to reach
TVector3 &x0, // pivot for sv
TKalMatrix &sv, // state vector
TKalMatrix &F, // propagator matrix
TKalMatrix &Q, // process noise matrix
std::unique_ptr<TVTrack> &help) // pointer to update track object
{
// ---------------------------------------------------------------------
// Sort measurement layers in this cradle if not
// ---------------------------------------------------------------------
if (!fDone) Update();
// ---------------------------------------------------------------------
// Locate sites from and to in this cradle
// ---------------------------------------------------------------------
Int_t fridx = from.GetHit().GetMeasLayer().GetIndex(); // index of site from
Int_t toidx = ml_to.GetIndex(); // index of layer to
Int_t di = fridx > toidx ? -1 : 1; // layer increment
TVTrack &hel = *help;
//=====================
// FIXME
//=====================
TVector3 xfrom = from.GetGlobalPivot(); // get the referenece point
TVector3 xto; // reference point at destination to be returned by CalcXingPointWith
Double_t fito = 0; // deflection angle to destination to be returned by CalcXingPointWith
const TVSurface *sfp = dynamic_cast<const TVSurface *>(&ml_to);// surface at destination
double eps = 1.e-8;
if(!TBField::IsUsingUniformBfield()) {
eps = 1.e-5;
}
sfp->CalcXingPointWith(hel, xto, fito, 0, eps);
// as mode is 0 here the closest point crossing point is taken
// this means that if we are at the top of a looping track
// and the point to which we want to move is on the other side of
// the loop but has a lower radius the transport will move down
// through all layers and segfault on reaching index -1
// if( does_cross < 1 ) return does_cross ;
TMatrixD dxdphi = hel.CalcDxDphi(fito); // tangent vector at destination surface
TVector3 dxdphiv(dxdphi(0,0),dxdphi(1,0),dxdphi(2,0)); // convert matirix diagonal to vector
// Double_t cpa = hel.GetKappa(); // get pt
Bool_t isout = -fito*dxdphiv.Dot(sfp->GetOutwardNormal(xto)) < 0 ? kTRUE : kFALSE; // out-going or in-coming at the destination surface
//=====================
// ENDFIXME
//=====================
TVector3 xx; // expected hit position vector
Double_t fid = 0.; // deflection angle from the last hit
Int_t sdim = sv.GetNrows(); // number of track parameters
F.UnitMatrix(); // set the propagator matrix to the unit matrix
Q.Zero(); // zero the noise matrix
TKalMatrix DF(sdim, sdim); // propagator matrix segment
// ---------------------------------------------------------------------
// Loop over layers and transport sv, F, and Q step by step
// ---------------------------------------------------------------------
Int_t ifr = fridx; // set index to the index of the intitial starting layer
// here we make first make sure that the helix is at the crossing point of the current surface.
// this is necessary to ensure that the material is only accounted for between fridx and toidx
// otherwise it is possible to have inconsistencies with material treatment.
// loop until we reach the index toidx, which is the surface we need to reach
for (Int_t ito=fridx; (di>0 && ito<=toidx)||(di<0 && ito>=toidx); ito += di) {
Double_t fid_temp = fid; // deflection angle from the last layer crossing
int mode = ito!=fridx ? di : 0; // need to move to the from site as the helix may not be on the crossing point yet, meaning that the eloss and ms will be incorrectely attributed ...
if (static_cast<TVSurface *>(At(ito))->CalcXingPointWith(hel, xx, fid, mode, eps)) { // if we have a crossing point at this surface, note di specifies if we are moving forwards or backwards
//=====================
// FIXME
//=====================
static const Double_t kMergin = 1.0;
// if the distance from the current crossing point to the starting point - kMergin(1mm) is greater than the distance from the destination to the starting point
// this is needed to skip crossing points which come from the far side of the IP, for a cylinder this would not be a problem
// but for the bounded planes it is perfectly posible due to the sorting in R
// reset the deflection angle and skip this layer
// this would at stop layers being added which are too far away but I am not sure how this will work with the problem described above.
if( (xx-xfrom).Mag() - kMergin > (xto-xfrom).Mag() ){
fid = fid_temp;
continue ;
}
//=====================
// ENDFIXME
//=====================
const TVMeasLayer &ml = *dynamic_cast<TVMeasLayer *>(At(ifr)); // get the last layer
TKalMatrix Qms(sdim, sdim);
if (IsMSOn()&& ito!=fridx ){
ml.CalcQms(isout, hel, fid, Qms); // Qms for this step, using the fact that the material was found to be outgoing or incomming above, and the distance from the last layer
}
hel.MoveTo(xx, fid, &DF); // move the helix to the present crossing point, DF will simply have its values overwritten so it could be explicitly set to unity here
if (sdim == 6) DF(5, 5) = 1.; // t0 stays the same
F = DF * F; // update F
TKalMatrix DFt = TKalMatrix(TMatrixD::kTransposed, DF);
Q = DF * (Q + Qms) * DFt; // transport Q to the present crossing point
if (IsDEDXOn() && ito!=fridx) {
hel.PutInto(sv); // copy hel to sv
// whether the helix is moving forwards or backwards is calculated using the sign of the charge and the sign of the deflection angle
// Bool_t isfwd = ((cpa > 0 && df < 0) || (cpa <= 0 && df > 0)) ? kForward : kBackward; // taken from TVMeasurmentLayer::GetEnergyLoss not df = fid
sv(2,0) += ml.GetEnergyLoss(isout, hel, fid); // correct for dE/dx, returns delta kappa i.e. the change in pt
hel.SetTo(sv, hel.GetPivot()); // save sv back to hel
}
ifr = ito; // for the next iteration set the "previous" layer to the current layer moved to
//fg: need to set the deflection angle to 0. as we moved the helix to a new position
// as this fid is used in the next call to TVSurface::CalcXingPointWith(), i.e. to
// compute the start position of the newtonian solver (which should be the reference point)
fid = 0. ;
} else {
fid = fid_temp;
}
} // end of loop over surfaces
// // ---------------------------------------------------------------------
// // Move pivot to crossing point with layer to move to
// // ---------------------------------------------------------------------
// dynamic_cast<const TVSurface *>(&ml_to)->CalcXingPointWith(hel, xx, fid);
// hel.MoveTo(xx, fid, &DF); // move pivot to expected hit, DF will simply have its values overwritten so it could be explicitly set to unity here
// F = DF * F; // update F accordingly
x0 = hel.GetPivot(); // local pivot corresponding to sv
hel.PutInto(sv); // save updated hel to sv
return 0;
}
int TKalDetCradle::Transport2(const TKalTrackSite &from, // site from
const TVMeasLayer &ml_to, // layer to reach
TVector3 &x0, // pivot for sv
TKalMatrix &sv, // state vector
TKalMatrix &F, // propagator matrix
TKalMatrix &Q, // process noise matrix
std::unique_ptr<TVTrack> &help) // pointer to update track object
{
// ---------------------------------------------------------------------
// Sort measurement layers in this cradle if not
// ---------------------------------------------------------------------
if (!fDone) Update();
// ---------------------------------------------------------------------
// Locate sites from and to in this cradle
// ---------------------------------------------------------------------
Int_t fridx = from.GetHit().GetMeasLayer().GetIndex(); // index of site from
Int_t toidx = ml_to.GetIndex(); // index of layer to
Int_t di = fridx > toidx ? -1 : 1; // layer increment
// FIXME
THelicalTrack& hel = dynamic_cast<THelicalTrack&>(*help);
TVector3 xx; // expected hit position vector
Double_t fid = 0.; // deflection angle from the last hit
Int_t sdim = sv.GetNrows(); // number of track parameters
F.UnitMatrix(); // set the propagator matrix to the unit matrix
Q.Zero(); // zero the noise matrix
TKalMatrix DF(sdim, sdim); // propagator matrix segment
// ---------------------------------------------------------------------
// Loop over layers and transport sv, F, and Q step by step
// ---------------------------------------------------------------------
Int_t ifr = fridx; // set index to the index of the intitial starting layer
// here we make first make sure that the helix is at the crossing point of the current surface.
// this is necessary to ensure that the material is only accounted for between fridx and toidx
// otherwise it is possible to have inconsistencies with material treatment.
// loop until we reach the index toidx, which is the surface we need to reach
TRungeKuttaTrack rk;
TVTrack* trackPtr = 0;
for (Int_t ito=fridx; (di>0 && ito<=toidx)||(di<0 && ito>=toidx); ito += di) {
Bool_t isout = di>0 ? kTRUE: kFALSE;
// need to move to the from site as the helix may not be on the crossing point yet,
// meaning that the eloss and ms will be incorrectely attributed
// the current TCylinder uses the Newtonian method to calculate crossing point,
// in which the mode is always 0.
int mode = ito!=fridx ? di : 0;
const TVMeasLayer &ml = *dynamic_cast<TVMeasLayer *>(At(ifr)); // get the last layer
//The crossing point between a TRungeKuttaTrack and a TVSurface
TVector3 rkxx;
//The initial step of Runge-Kutta algorithm
Double_t step = 0.01;
if(ito==fridx) {
// If ito==fridx, it means the track will move from the current pivot (hit point)
// to the crossing point at the SAME layer.
// Because the distance between the two points is very small, so the difference of b field
// of two points are also small, and we use the helical track model.
// helix is defined by local track parameters.
// xx is a global coordinate.
// angle fid is difference of direction for pivot and crossing point in a helix.
dynamic_cast<TVSurface *>(At(ito))->CalcXingPointWith(hel, xx, fid, mode);
}
else {
// If ito!=fridx, it means the track will move from the currnt pivot (crossing point)
// to the next crossing point of track and next layer.
// Considering the non-uniformity of the magnetic field, we use the Runge-Kutta track model here,
// and we create a runge-kutta track to calculate the crossing point.
rk.SetFromTrack(hel);
dynamic_cast<TVSurface *>(At(ito))->CalcXingPointWith(rk, rkxx, step, mode);
}
TKalMatrix Qms(sdim, sdim);
if (IsMSOn()&& ito!=fridx ){
ml.CalcQms(isout, hel, fid, Qms);
// Qms for this step, using the fact that the material was found to be outgoing
// or incomming above, and the distance from the last layer
}
if(ito==fridx) {
// Move the helix to the present crossing point.
// Frame rotation is not needed.
hel.MoveTo(xx, fid, &DF, 0, kFALSE);
trackPtr = &hel;
}
else {
// Move the track to the new layer.
// A rotation is done in TRungeKuttaTrack::MoveTo.
rk.MoveTo(rkxx, step, DF);
trackPtr = &rk;
rk.SetToTrack(hel);
}
if (sdim == 6) DF(5, 5) = 1.; // t0 stays the same
F = DF * F; // update F
TKalMatrix DFt = TKalMatrix(TMatrixD::kTransposed, DF);
Q = DF * (Q + Qms) * DFt; // transport Q to the present crossing point
if (IsDEDXOn() && ito!=fridx) {
hel.PutInto(sv); // copy hel to sv
// whether the helix is moving forwards or backwards is calculated using
// the sign of the charge and the sign of the deflection angle
// Bool_t isfwd = ((cpa > 0 && df < 0) || (cpa <= 0 && df > 0)) ? kForward : kBackward;
// taken from TVMeasurmentLayer::GetEnergyLoss not df = fid
sv(2,0) += ml.GetEnergyLoss(isout, hel, fid);
// correct for dE/dx, returns delta kappa i.e. the change in pt
hel.SetTo(sv, hel.GetPivot()); // save sv back to hel
}
ifr = ito; // for the next iteration set the "previous" layer to the current layer moved to
} // end of loop over surfaces
//hel or rk
x0 = trackPtr->GetPivot();
trackPtr->PutInto(sv); // save updated hel to sv
return 0;
}
//_________________________________________________________________________
// -----------------
// Update
// -----------------
// sorts meaurement layers according to layer's sorting policy
// and puts index to layers from inside to outside.
//
void TKalDetCradle::Update()
{
fDone = kTRUE;
UnSort(); // unsort
Sort(); // sort layers according to sorting policy
TIter next(this);
TVMeasLayer *mlp = 0;
Int_t i = 0;
while ((mlp = dynamic_cast<TVMeasLayer *>(next()))) {
mlp->SetIndex(i++);
}
}