Project
Loading...
Searching...
No Matches
GPUTPCGMPropagator.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
14
15#include "GPUTPCGMPropagator.h"
16#include "GPUTPCGMPhysicalTrackModel.h"
17#include "GPUParam.h"
18#include "GPUTPCGMMergedTrackHit.h"
19#include "GPUO2DataTypes.h"
20#include "GPUParam.inc"
21#include "GPUTPCGMMergerTypes.h"
22#include "GPUDebugStreamer.h"
23#include "GPUTPCGMMerger.h"
24
25using namespace o2::gpu;
26
27GPUd() void GPUTPCGMPropagator::GetBxByBzBase(float cosAlpha, float sinAlpha, float X, float Y, float Z, float B[3]) const
28{
29 // get global coordinates
30
31 float gx = getGlobalX(cosAlpha, sinAlpha, X, Y);
32 float gy = getGlobalY(cosAlpha, sinAlpha, X, Y);
33
34 float bb[3];
35 switch (mFieldRegion) {
36 case ITS:
37 mField->GetFieldIts(gx, gy, Z, bb);
38 break;
39 case TRD:
40 mField->GetFieldTrd(gx, gy, Z, bb);
41 break;
42 case TPC:
43 default:
44 mField->GetField(gx, gy, Z, bb);
45 }
46
47 // rotate field to local coordinates
48
49 B[0] = bb[0] * cosAlpha + bb[1] * sinAlpha;
50 B[1] = -bb[0] * sinAlpha + bb[1] * cosAlpha;
51 B[2] = bb[2];
52}
53
54GPUd() float GPUTPCGMPropagator::GetBzBase(float cosAlpha, float sinAlpha, float X, float Y, float Z) const
55{
56 float gx = getGlobalX(cosAlpha, sinAlpha, X, Y);
57 float gy = getGlobalY(cosAlpha, sinAlpha, X, Y);
58
59 switch (mFieldRegion) {
60 case ITS:
61 return mField->GetFieldItsBz(gx, gy, Z);
62 case TRD:
63 return mField->GetFieldTrdBz(gx, gy, Z);
64 case TPC:
65 default:
66 return mField->GetFieldBz(gx, gy, Z);
67 }
68}
69
70GPUd() int32_t GPUTPCGMPropagator::RotateToAlpha(float newAlpha)
71{
72 //
73 // Rotate the track coordinate system in XY to the angle newAlpha
74 // return value is error code (0==no error)
75 //
76
77 float newCosAlpha, newSinAlpha;
78 CAMath::SinCos(newAlpha, newSinAlpha, newCosAlpha);
79
80 float cc = newCosAlpha * mCosAlpha + newSinAlpha * mSinAlpha; // CAMath::Cos(newAlpha - mAlpha);
81 float ss = newSinAlpha * mCosAlpha - newCosAlpha * mSinAlpha; // CAMath::Sin(newAlpha - mAlpha);
82
83 GPUTPCGMPhysicalTrackModel t0 = mT0;
84
85 float x0 = mT0.X();
86 float y0 = mT0.Y();
87 float px0 = mT0.Px();
88 float py0 = mT0.Py();
89 // float pt0 = mT0.GetPt();
90
91 if (CAMath::Abs(mT->GetSinPhi()) >= mMaxSinPhi || CAMath::Abs(px0) < (1 - mMaxSinPhi)) {
92 return -1;
93 }
94
95 // rotate t0 track
96 float px1 = px0 * cc + py0 * ss;
97 float py1 = -px0 * ss + py0 * cc;
98
99 {
100 t0.X() = x0 * cc + y0 * ss;
101 t0.Y() = -x0 * ss + y0 * cc;
102 t0.Px() = px1;
103 t0.Py() = py1;
104 t0.UpdateValues();
105 }
106
107 if (CAMath::Abs(py1) > mMaxSinPhi * mT0.GetPt() || CAMath::Abs(px1) < (1 - mMaxSinPhi)) {
108 return -1;
109 }
110
111 // calculate X of rotated track:
112 float trackX = x0 * cc + ss * mT->Y();
113
114 // transport t0 to trackX
115 float dLp = 0;
116 float Bz;
117 if (mPropagateBzOnly) {
118 Bz = GetBzBase(newCosAlpha, newSinAlpha, t0.X(), t0.Y(), t0.Z());
119 if (t0.PropagateToXBzLight(trackX, Bz, dLp)) {
120 return -1;
121 }
122 } else {
123 float B[3];
124 GetBxByBz(newAlpha, t0.X(), t0.Y(), t0.Z(), B);
125 Bz = B[2];
126 if (t0.PropagateToXBxByBz(trackX, B[0], B[1], B[2], dLp)) {
127 return -1;
128 }
129 }
130
131 if (CAMath::Abs(t0.SinPhi()) >= mMaxSinPhi) {
132 return -1;
133 }
134
135 // now t0 is rotated and propagated, all checks are passed
136
137 // Rotate track using mT0 for linearisation. After rotation X is not fixed, but has a covariance
138
139 // Y Z Sin DzDs q/p
140 // Jacobian J0 = { { j0, 0, 0, 0, 0 }, // Y
141 // { 0, 1, 0, 0, 0 }, // Z
142 // { 0, 0, j1, 0, 0 }, // SinPhi
143 // { 0, 0, 0, 1, 0 }, // DzDs
144 // { 0, 0, 0, 0, 1 }, // q/p
145 // { j2, 0, 0, 0, 0 } }// X (rotated )
146
147 float j0 = cc;
148 float j1 = px1 / px0;
149 float j2 = ss;
150 // float dy = mT->Y() - y0;
151 // float ds = mT->SinPhi() - mT0.SinPhi();
152
153 mT->X() = trackX; // == x0*cc + ss*mT->Y() == t0.X() + j0*dy;
154 mT->Y() = -x0 * ss + cc * mT->Y(); //== t0.Y() + j0*dy;
155 // mT->SinPhi() = py1/pt0 + j1*ds; // == t0.SinPhi() + j1*ds; // use py1, since t0.SinPhi can have different sign
156 mT->SinPhi() = -CAMath::Sqrt(1.f - mT->SinPhi() * mT->SinPhi()) * ss + mT->SinPhi() * cc;
157
158 // Rotate cov. matrix Cr = J0 x C x J0T. Cr has one more row+column for X:
159 float* c = mT->Cov();
160
161 float c15 = c[0] * j0 * j2;
162 float c16 = c[1] * j2;
163 float c17 = c[3] * j1 * j2;
164 float c18 = c[6] * j2;
165 float c19 = c[10] * j2;
166 float c20 = c[0] * j2 * j2;
167
168 c[0] *= j0 * j0;
169 c[3] *= j0;
170 c[10] *= j0;
171
172 c[3] *= j1;
173 c[5] *= j1 * j1;
174 c[12] *= j1;
175
176 if (!mFitInProjections && mT->NDF() > 0) {
177 c[1] *= j0;
178 c[6] *= j0;
179 c[4] *= j1;
180 c[8] *= j1;
181 }
182
183 if (t0.SetDirectionAlongX()) { // change direction if Px < 0
184 mT->SinPhi() = -mT->SinPhi();
185 mT->DzDs() = -mT->DzDs();
186 mT->QPt() = -mT->QPt();
187 c[3] = -c[3]; // covariances with SinPhi
188 c[4] = -c[4];
189 c17 = -c17;
190 c[6] = -c[6]; // covariances with DzDs
191 c[7] = -c[7];
192 c18 = -c18;
193 c[10] = -c[10]; // covariances with QPt
194 c[11] = -c[11];
195 c19 = -c19;
196 }
197
198 // Now fix the X coordinate: so to say, transport track T to fixed X = mT->X().
199 // only covariance changes. Use rotated and transported t0 for linearisation
200 float j3 = -t0.Py() / t0.Px();
201 float j4 = -t0.Pz() / t0.Px();
202 float j5 = t0.QPt() * Bz;
203
204 // Y Z Sin DzDs q/p X
205 // Jacobian J1 = { { 1, 0, 0, 0, 0, j3 }, // Y
206 // { 0, 1, 0, 0, 0, j4 }, // Z
207 // { 0, 0, 1, 0, 0, j5 }, // SinPhi
208 // { 0, 0, 0, 1, 0, 0 }, // DzDs
209 // { 0, 0, 0, 0, 1, 0 } }; // q/p
210
211 float h15 = c15 + c20 * j3;
212 float h16 = c16 + c20 * j4;
213 float h17 = c17 + c20 * j5;
214
215 c[0] += j3 * (c15 + h15);
216
217 c[2] += j4 * (c16 + h16);
218
219 c[3] += c17 * j3 + h15 * j5;
220 c[5] += j5 * (c17 + h17);
221
222 c[7] += c18 * j4;
223 // c[ 9] = c[ 9];
224
225 c[10] += c19 * j3;
226 c[12] += c19 * j5;
227 // c[14] = c[14];
228
229 if (!mFitInProjections && mT->NDF() > 0) {
230 c[1] += c16 * j3 + h15 * j4;
231 c[4] += c17 * j4 + h16 * j5;
232 c[6] += c18 * j3;
233 c[8] += c18 * j5;
234 c[11] += c19 * j4;
235 // c[13] = c[13];
236 }
237
238 mAlpha = newAlpha;
239 mCosAlpha = newCosAlpha;
240 mSinAlpha = newSinAlpha;
241 mT0 = t0;
242
243 return 0;
244}
245
246GPUd() int32_t GPUTPCGMPropagator::PropagateToXAlpha(float posX, float posAlpha, bool inFlyDirection)
247{
248 if (mPropagateBzOnly) {
249 return PropagateToXAlphaBz(posX, posAlpha, inFlyDirection);
250 }
251 if (CAMath::Abs(posAlpha - mAlpha) > 1.e-4f) {
252 if (RotateToAlpha(posAlpha) != 0) {
253 return -2;
254 }
255 }
256 if (CAMath::Abs(posX - mT->X()) < 1.e-7f) {
257 mT->SetX(posX);
258 return 0;
259 }
260
261 float B[3];
262 GetBxByBz(mT0.X(), mT0.Y(), mT0.Z(), B);
263
264 // propagate mT0 to t0e
265
266 GPUTPCGMPhysicalTrackModel t0e(mT0);
267 float dLp = 0;
268 if (t0e.PropagateToXBxByBz(posX, B[0], B[1], B[2], dLp) && t0e.PropagateToXBzLight(posX, B[2], dLp)) {
269 return 1;
270 }
271
272 if (CAMath::Abs(t0e.SinPhi()) >= mMaxSinPhi) {
273 return -3;
274 }
275
276 return FollowLinearization(t0e, B[2], dLp, inFlyDirection);
277}
278
279GPUd() int32_t GPUTPCGMPropagator::PropagateToXAlphaBz(float posX, float posAlpha, bool inFlyDirection)
280{
281 if (CAMath::Abs(posAlpha - mAlpha) > 1.e-4f) {
282 if (RotateToAlpha(posAlpha) != 0) {
283 return -2;
284 }
285 }
286 if (CAMath::Abs(posX - mT->X()) < 1.e-7f) {
287 mT->SetX(posX);
288 return 0;
289 }
290
291 float Bz = GetBz(mT0.X(), mT0.Y(), mT0.Z());
292
293 // propagate mT0 to t0e
294
295 GPUTPCGMPhysicalTrackModel t0e(mT0);
296 float dLp = 0;
297 if (t0e.PropagateToXBzLight(posX, Bz, dLp)) {
298 return 1;
299 }
300
301 if (CAMath::Abs(t0e.SinPhi()) >= mMaxSinPhi) {
302 return -3;
303 }
304
305 return FollowLinearization(t0e, Bz, dLp, inFlyDirection);
306}
307
308GPUd() int32_t GPUTPCGMPropagator::FollowLinearization(const GPUTPCGMPhysicalTrackModel& GPUrestrict() t0e, float Bz, float dLp, bool inFlyDirection)
309{
310 // t0e is alrerady extrapolated t0
311
312 // propagate track and cov matrix with derivatives for (0,0,Bz) field
313
314 float dS = dLp * t0e.GetPt();
315 float dL = CAMath::Abs(dLp * t0e.GetP());
316
317 if (inFlyDirection) {
318 dL = -dL;
319 }
320
321 float ey = mT0.SinPhi();
322 float ex = mT0.CosPhi();
323 float exi = mT0.SecPhi();
324 float ey1 = t0e.GetSinPhi();
325 float ex1 = t0e.GetCosPhi();
326 float ex1i = t0e.GetSecPhi();
327
328 float k = -mT0.QPt() * Bz;
329 float dx = t0e.GetX() - mT0.X();
330 float kdx = k * dx;
331 float cc = ex + ex1;
332 float cci = 1.f / cc;
333
334 float dxcci = dx * cci;
335 float hh = dxcci * ex1i * (1.f + ex * ex1 + ey * ey1);
336
337 float j02 = exi * hh;
338 float j04 = -Bz * dxcci * hh;
339 float j13 = dS;
340 float j24 = -dx * Bz;
341
342 float* p = mT->Par();
343
344 float d0 = p[0] - mT0.Y();
345 float d1 = p[1] - mT0.Z();
346 float d2 = p[2] - mT0.SinPhi();
347 float d3 = p[3] - mT0.DzDs();
348 float d4 = p[4] - mT0.QPt();
349
350 float newSinPhi = ey1 + d2 + j24 * d4;
351 if (mT->NDF() >= 15 && CAMath::Abs(newSinPhi) > GPUCA_MAX_SIN_PHI) {
352 return -4;
353 }
354
355 if (mMatLUT) {
356 UpdateMaterial(t0e);
357 }
358
359 mT0 = t0e;
360 mT->X() = t0e.GetX();
361 p[0] = t0e.GetY() + d0 + j02 * d2 + j04 * d4;
362 p[1] = t0e.GetZ() + d1 + j13 * d3;
363 p[2] = newSinPhi;
364 p[3] = t0e.GetDzDs() + d3;
365 p[4] = t0e.GetQPt() + d4;
366
367 float* c = mT->Cov();
368
369 float c20 = c[3];
370 float c21 = c[4];
371 float c22 = c[5];
372
373 float c30 = c[6];
374 float c31 = c[7];
375 float c32 = c[8];
376 float c33 = c[9];
377
378 float c40 = c[10];
379 float c41 = c[11];
380 float c42 = c[12];
381 float c43 = c[13];
382 float c44 = c[14];
383
384 if (mFitInProjections || mT->NDF() <= 0) {
385 float c20ph04c42 = c20 + j04 * c42;
386 float j02c22 = j02 * c22;
387 float j04c44 = j04 * c44;
388
389 float n6 = c30 + j02 * c32 + j04 * c43;
390 float n7 = c31 + j13 * c33;
391 float n10 = c40 + j02 * c42 + j04c44;
392 float n11 = c41 + j13 * c43;
393 float n12 = c42 + j24 * c44;
394
395 c[0] += j02 * j02c22 + j04 * j04c44 + 2.f * (j02 * c20ph04c42 + j04 * c40);
396 c[1] += j02 * c21 + j04 * c41 + j13 * n6;
397 c[2] += j13 * (c31 + n7);
398 c[3] = c20ph04c42 + j02c22 + j24 * n10;
399 c[4] = c21 + j13 * c32 + j24 * n11;
400 c[5] = c22 + j24 * (c42 + n12);
401 c[6] = n6;
402 c[7] = n7;
403 c[8] = c32 + c43 * j24;
404 c[10] = n10;
405 c[11] = n11;
406 c[12] = n12;
407 } else {
408 float c00 = c[0];
409 float c10 = c[1];
410 float c11 = c[2];
411
412 float ss = ey + ey1;
413 float tg = ss * cci;
414 float xx = 1.f - 0.25f * kdx * kdx * (1.f + tg * tg);
415 if (xx < 1.e-8f) {
416 return -1;
417 }
418 xx = CAMath::Sqrt(xx);
419 float yy = CAMath::Sqrt(ss * ss + cc * cc);
420
421 float j12 = dx * mT0.DzDs() * tg * (2.f + tg * (ey * exi + ey1 * ex1i)) / (xx * yy);
422 float j14 = 0;
423 if (CAMath::Abs(mT0.QPt()) > 1.e-6f) {
424 j14 = (2.f * xx * ex1i * dx / yy - dS) * mT0.DzDs() / mT0.QPt();
425 } else {
426 j14 = -mT0.DzDs() * Bz * dx * dx * exi * exi * exi * (0.5f * ey + (1.f / 3.f) * kdx * (1 + 2.f * ey * ey) * exi * exi);
427 }
428
429 p[1] += j12 * d2 + j14 * d4;
430
431 float h00 = c00 + c20 * j02 + c40 * j04;
432 // float h01 = c10 + c21*j02 + c41*j04;
433 float h02 = c20 + c22 * j02 + c42 * j04;
434 // float h03 = c30 + c32*j02 + c43*j04;
435 float h04 = c40 + c42 * j02 + c44 * j04;
436
437 float h10 = c10 + c20 * j12 + c30 * j13 + c40 * j14;
438 float h11 = c11 + c21 * j12 + c31 * j13 + c41 * j14;
439 float h12 = c21 + c22 * j12 + c32 * j13 + c42 * j14;
440 float h13 = c31 + c32 * j12 + c33 * j13 + c43 * j14;
441 float h14 = c41 + c42 * j12 + c43 * j13 + c44 * j14;
442
443 float h20 = c20 + c40 * j24;
444 float h21 = c21 + c41 * j24;
445 float h22 = c22 + c42 * j24;
446 float h23 = c32 + c43 * j24;
447 float h24 = c42 + c44 * j24;
448
449 c[0] = h00 + h02 * j02 + h04 * j04;
450
451 c[1] = h10 + h12 * j02 + h14 * j04;
452 c[2] = h11 + h12 * j12 + h13 * j13 + h14 * j14;
453
454 c[3] = h20 + h22 * j02 + h24 * j04;
455 c[4] = h21 + h22 * j12 + h23 * j13 + h24 * j14;
456 c[5] = h22 + h24 * j24;
457
458 c[6] = c30 + c32 * j02 + c43 * j04;
459 c[7] = c31 + c32 * j12 + c33 * j13 + c43 * j14;
460 c[8] = c32 + c43 * j24;
461 // c[ 9] = c33;
462
463 c[10] = c40 + c42 * j02 + c44 * j04;
464 c[11] = c41 + c42 * j12 + c43 * j13 + c44 * j14;
465 c[12] = c42 + c44 * j24;
466 // c[13] = c43;
467 // c[14] = c44;
468 }
469
470 float& mC22 = c[5];
471 float& mC33 = c[9];
472 float& mC40 = c[10];
473 float& mC41 = c[11];
474 float& mC42 = c[12];
475 float& mC43 = c[13];
476 float& mC44 = c[14];
477
478 float dLmask = 0.f;
479 bool maskMS = (CAMath::Abs(dL) < mMaterial.DLMax);
480 if (maskMS) {
481 dLmask = dL;
482 }
483 float dLabs = CAMath::Abs(dLmask);
484
485 // Energy Loss
486 if (true) {
487 // std::cout<<"APPLY ENERGY LOSS!!!"<<std::endl;
488 float corr = 1.f - mMaterial.EP2 * dLmask;
489 float corrInv = 1.f / corr;
490 mT0.Px() *= corrInv;
491 mT0.Py() *= corrInv;
492 mT0.Pz() *= corrInv;
493 mT0.Pt() *= corrInv;
494 mT0.P() *= corrInv;
495 mT0.QPt() *= corr;
496
497 p[4] *= corr;
498
499 mC40 *= corr;
500 mC41 *= corr;
501 mC42 *= corr;
502 mC43 *= corr;
503 mC44 = mC44 * corr * corr + dLabs * mMaterial.sigmadE2;
504 } else {
505 // std::cout<<"DONT APPLY ENERGY LOSS!!!"<<std::endl;
506 }
507
508 // Multiple Scattering
509 if (true) {
510 mC22 += dLabs * mMaterial.k22 * mT0.CosPhi() * mT0.CosPhi();
511 mC33 += dLabs * mMaterial.k33;
512 mC43 += dLabs * mMaterial.k43;
513 mC44 += dLabs * mMaterial.k44;
514 }
515 return 0;
516}
517
518GPUd() int32_t GPUTPCGMPropagator::GetPropagatedYZ(float x, float& GPUrestrict() projY, float& GPUrestrict() projZ)
519{
520 float bz = GetBz(mT->X(), mT->Y(), mT->Z());
521 float k = mT0.QPt() * bz;
522 float dx = x - mT->X();
523 float ex = mT0.CosPhi();
524 float ey = mT0.SinPhi();
525 float ey1 = ey - k * dx;
526 if (CAMath::Abs(ey1) > GPUCA_MAX_SIN_PHI) {
527 return 1;
528 }
529 float ss = ey + ey1;
530 float ex1 = CAMath::Sqrt(1.f - ey1 * ey1);
531 if (ex < 0) {
532 ex1 = -ex1;
533 }
534 float cc = ex + ex1;
535 float dxcci = dx / cc;
536 float dy = dxcci * ss;
537 float norm2 = 1.f + ey * ey1 + ex * ex1;
538 float dl = dxcci * CAMath::Sqrt(norm2 + norm2);
539 float dS;
540 {
541 float dSin = 0.5f * k * dl;
542 float a = dSin * dSin;
543 const float k2 = 1.f / 6.f;
544 const float k4 = 3.f / 40.f;
545 dS = dl + dl * a * (k2 + a * (k4));
546 }
547 float dz = dS * mT0.DzDs();
548 projY = mT->Y() + dy;
549 projZ = mT->Z() + dz;
550 return 0;
551}
552
553GPUd() void GPUTPCGMPropagator::GetErr2(float& GPUrestrict() err2Y, float& GPUrestrict() err2Z, const GPUParam& GPUrestrict() param, float posZ, int32_t iRow, int16_t clusterState, int8_t sector, float time, float avgCharge, float charge) const
554{
555 GetErr2(err2Y, err2Z, param, mT0.GetSinPhi(), mT0.DzDs(), posZ, mT->GetX(), mT->GetY(), iRow, clusterState, sector, time, avgCharge, charge, mSeedingErrors);
556}
557
558GPUd() void GPUTPCGMPropagator::GetErr2(float& GPUrestrict() err2Y, float& GPUrestrict() err2Z, const GPUParam& GPUrestrict() param, float snp, float tgl, float posZ, float trackX, float trackY, int32_t iRow, int16_t clusterState, int8_t sector, float time, float avgCharge, float charge, bool seedingErrors)
559{
560#ifndef GPUCA_TPC_GEOMETRY_O2
561 if (seedingErrors) {
562 param.GetClusterErrorsSeeding2(sector, iRow, posZ, snp, tgl, time, err2Y, err2Z);
563 } else
564#endif
565 {
566 param.GetClusterErrors2(sector, iRow, posZ, snp, tgl, time, avgCharge, charge, err2Y, err2Z);
567 }
568 param.UpdateClusterError2ByState(clusterState, err2Y, err2Z);
569 float statErr2 = param.GetSystematicClusterErrorIFC2(trackX, trackY, posZ, sector >= (GPUCA_NSECTORS / 2));
570 if (sector >= GPUCA_NSECTORS / 2 + 1 && sector <= GPUCA_NSECTORS / 2 + 2) {
571 statErr2 += param.GetSystematicClusterErrorC122(trackX, trackY, sector);
572 }
573 err2Y += statErr2;
574 err2Z += statErr2;
575}
576
577GPUd() float GPUTPCGMPropagator::PredictChi2(float posY, float posZ, int32_t iRow, const GPUParam& GPUrestrict() param, int16_t clusterState, int8_t sector, float time, float avgCharge, float charge) const
578{
579 float err2Y, err2Z;
580 GetErr2(err2Y, err2Z, param, posZ, iRow, clusterState, sector, time, avgCharge, charge);
581 return PredictChi2(posY, posZ, err2Y, err2Z);
582}
583
584GPUd() float GPUTPCGMPropagator::PredictChi2(float posY, float posZ, float err2Y, float err2Z) const
585{
586 const float* mC = mT->Cov();
587 const float* mP = mT->Par();
588 const float z0 = posY - mP[0];
589 const float z1 = posZ - mP[1];
590
591 if (!mFitInProjections || mT->NDF() <= 0) {
592 const float w0 = 1.f / (err2Y + mC[0]);
593 const float w2 = 1.f / (err2Z + mC[2]);
594 return w0 * z0 * z0 + w2 * z1 * z1;
595 } else {
596 float w0 = mC[2] + err2Z, w1 = mC[1], w2 = mC[0] + err2Y;
597 { // Invert symmetric matrix
598 float det = w0 * w2 - w1 * w1;
599 if (CAMath::Abs(det) < 1.e-10f) {
600 det = 1.e-10f;
601 }
602 det = 1.f / det;
603 w0 = w0 * det;
604 w1 = -w1 * det;
605 w2 = w2 * det;
606 }
607 return CAMath::Abs((w0 * z0 + w1 * z1) * z0) + CAMath::Abs((w1 * z0 + w2 * z1) * z1);
608 }
609}
610
611GPUd() int32_t GPUTPCGMPropagator::Update(float posY, float posZ, int32_t iRow, const GPUParam& GPUrestrict() param, int16_t clusterState, int8_t rejectChi2, gputpcgmmergertypes::InterpolationErrorHit* inter, bool refit, int8_t sector, float time, float avgInvCharge, float invCharge GPUCA_DEBUG_STREAMER_CHECK(, DebugStreamerVals* debugVals))
612{
613 float err2Y, err2Z;
614 GetErr2(err2Y, err2Z, param, posZ, iRow, clusterState, sector, time, avgInvCharge, invCharge);
615 GPUCA_DEBUG_STREAMER_CHECK(if (debugVals) { debugVals->err2Y = err2Y; debugVals->err2Z = err2Z; });
616
617 if (rejectChi2 >= rejectInterFill) {
618 if (rejectChi2 == rejectInterReject && inter->errorY < (GPUCA_PAR_MERGER_INTERPOLATION_ERROR_TYPE_A)0) {
619 rejectChi2 = rejectDirect;
620 } else {
621 int32_t retVal = InterpolateReject(param, posY, posZ, clusterState, rejectChi2, inter, err2Y, err2Z);
622 GPUCA_DEBUG_STREAMER_CHECK(if (debugVals) { debugVals->retVal = retVal; });
623 if (retVal) {
624 return retVal;
625 }
626 }
627 }
628
629 if (mT->NDF() == -5) { // first measurement: no need to filter, as the result is known in advance. just set it.
630 mT->ResetCovariance();
631 float* mC = mT->Cov();
632 float* mP = mT->Par();
633 if (refit) {
634 mC[14] = CAMath::Max(0.5f, CAMath::Abs(mP[4]));
635 mC[5] = CAMath::Max(0.2f, CAMath::Abs(mP[2]) / 2);
636 mC[9] = CAMath::Max(0.5f, CAMath::Abs(mP[3]) / 2);
637 }
638 mP[0] = posY;
639 mP[1] = posZ;
640 mC[0] = err2Y;
641 mC[2] = err2Z;
642 mT->NDF() = -3;
643 return 0;
644 }
645
646 return Update(posY, posZ, clusterState, rejectChi2 == rejectDirect || (param.rec.tpc.mergerInterpolateRejectAlsoOnCurrentPosition && rejectChi2 == rejectInterReject), err2Y, err2Z, &param);
647}
648
649GPUd() int32_t GPUTPCGMPropagator::InterpolateReject(const GPUParam& GPUrestrict() param, float posY, float posZ, int16_t clusterState, int8_t rejectChi2, gputpcgmmergertypes::InterpolationErrorHit* inter, float err2Y, float err2Z)
650{
651 float* GPUrestrict() mC = mT->Cov();
652 float* GPUrestrict() mP = mT->Par();
653 if (rejectChi2 == rejectInterFill) {
654 inter->posY = mP[0];
655 inter->posZ = mP[1];
656 inter->errorY = mC[0];
657 inter->errorZ = mC[2];
658 } else if (rejectChi2 == rejectInterReject) {
659 float chi2Y, chi2Z;
660 if (mFitInProjections || mT->NDF() <= 0) {
661 const float Iz0 = inter->posY - mP[0];
662 const float Iz1 = inter->posZ - mP[1];
663 const float Iw0 = 1.f / (mC[0] + (float)inter->errorY);
664 const float Iw2 = 1.f / (mC[2] + (float)inter->errorZ);
665 const float Ik00 = mC[0] * Iw0;
666 const float Ik11 = mC[2] * Iw2;
667 const float ImP0 = mP[0] + Ik00 * Iz0;
668 const float ImP1 = mP[1] + Ik11 * Iz1;
669 const float ImC0 = mC[0] - Ik00 * mC[0];
670 const float ImC2 = mC[2] - Ik11 * mC[2];
671 // printf("\t%21sInterpo ----- abde artaf%16s Y %8.3f, Z %8.3f (Errors %f <-- (%f, %f) %f <-- (%f, %f))\n", "", "", ImP0, ImP1, sqrtf(ImC0), sqrtf(mC[0]), sqrtf(inter->errorY), sqrtf(ImC2), sqrtf(mC[2]), sqrtf(inter->errorZ));
672 const float Jz0 = posY - ImP0;
673 const float Jz1 = posZ - ImP1;
674 const float Jw0 = 1.f / (ImC0 + err2Y);
675 const float Jw2 = 1.f / (ImC2 + err2Z);
676 chi2Y = Jw0 * Jz0 * Jz0;
677 chi2Z = Jw2 * Jz1 * Jz1;
678 } else {
679 const float Iz0 = inter->posY - mP[0];
680 const float Iz1 = inter->posZ - mP[1];
681 float Iw0 = mC[2] + (float)inter->errorZ;
682 float Iw2 = mC[0] + (float)inter->errorY;
683 float Idet = CAMath::Max(1e-10f, Iw0 * Iw2 - mC[1] * mC[1]);
684 Idet = 1.f / Idet;
685 Iw0 *= Idet;
686 const float Iw1 = mC[1] * Idet;
687 Iw2 *= Idet;
688 const float Ik00 = mC[0] * Iw0 + mC[1] * Iw1;
689 const float Ik01 = mC[0] * Iw1 + mC[1] * Iw2;
690 const float Ik10 = mC[1] * Iw0 + mC[2] * Iw1;
691 const float Ik11 = mC[1] * Iw1 + mC[2] * Iw2;
692 const float ImP0 = mP[0] + Ik00 * Iz0 + Ik01 * Iz1;
693 const float ImP1 = mP[1] + Ik10 * Iz0 + Ik11 * Iz1;
694 const float ImC0 = mC[0] - Ik00 * mC[0] + Ik01 * mC[1];
695 const float ImC1 = mC[1] - Ik10 * mC[0] + Ik11 * mC[1];
696 const float ImC2 = mC[2] - Ik10 * mC[1] + Ik11 * mC[2];
697 const float Jz0 = posY - ImP0;
698 const float Jz1 = posZ - ImP1;
699 float Jw0 = ImC2 + err2Z;
700 float Jw2 = ImC0 + err2Y;
701 float Jdet = CAMath::Max(1e-10f, Jw0 * Jw2 - ImC1 * ImC1);
702 Jdet = 1.f / Jdet;
703 Jw0 *= Jdet;
704 const float Jw1 = ImC1 * Jdet;
705 Jw2 *= Jdet;
706 chi2Y = CAMath::Abs((Jw0 * Jz0 + Jw1 * Jz1) * Jz0);
707 chi2Z = CAMath::Abs((Jw1 * Jz0 + Jw2 * Jz1) * Jz1);
708 }
709 if (RejectCluster(chi2Y * param.rec.tpc.clusterRejectChi2TolleranceY, chi2Z * param.rec.tpc.clusterRejectChi2TolleranceZ, clusterState)) { // TODO: Relative Pt resolution decreases slightly, why?
710 return updateErrorClusterRejectedInInterpolation;
711 }
712 }
713 return 0;
714}
715
716GPUd() int32_t GPUTPCGMPropagator::Update(float posY, float posZ, int16_t clusterState, bool rejectChi2, float err2Y, float err2Z, const GPUParam* GPUrestrict() param)
717{
718 float* GPUrestrict() mC = mT->Cov();
719 float* GPUrestrict() mP = mT->Par();
720
721 const float d00 = mC[0], d01 = mC[1], d02 = mC[3], d03 = mC[6], d04 = mC[10];
722 const float d10 = mC[1], d11 = mC[2], d12 = mC[4], d13 = mC[7], d14 = mC[11];
723
724 const float z0 = posY - mP[0];
725 const float z1 = posZ - mP[1];
726 float w0, w1, w2, chi2Y, chi2Z;
727 if (mFitInProjections || mT->NDF() <= 0) {
728 w0 = 1.f / (err2Y + d00);
729 w1 = 0;
730 w2 = 1.f / (err2Z + d11);
731 chi2Y = w0 * z0 * z0;
732 chi2Z = w2 * z1 * z1;
733 } else {
734 w0 = d11 + err2Z, w1 = d10, w2 = d00 + err2Y;
735 { // Invert symmetric matrix
736 float det = w0 * w2 - w1 * w1;
737 if (CAMath::Abs(det) < 1.e-10f) {
738 return updateErrorFitFailed;
739 }
740 det = 1.f / det;
741 w0 = w0 * det;
742 w1 = -w1 * det;
743 w2 = w2 * det;
744 }
745 chi2Y = CAMath::Abs((w0 * z0 + w1 * z1) * z0);
746 chi2Z = CAMath::Abs((w1 * z0 + w2 * z1) * z1);
747 }
748 float dChi2 = chi2Y + chi2Z;
749 // GPUInfo("hits %d chi2 %f, new %f %f (dy %f dz %f)", N, mChi2, chi2Y, chi2Z, z0, z1);
750 if (rejectChi2 && RejectCluster(chi2Y * param->rec.tpc.clusterRejectChi2TolleranceY, chi2Z * param->rec.tpc.clusterRejectChi2TolleranceZ, clusterState)) {
751 return updateErrorClusterRejectedInUpdate;
752 }
753 mT->Chi2() += dChi2;
754 mT->NDF() += 2;
755
756 if (mFitInProjections || mT->NDF() <= 0) {
757 const float k00 = d00 * w0;
758 const float k20 = d02 * w0;
759 const float k40 = d04 * w0;
760 const float k11 = d11 * w2;
761 const float k31 = d13 * w2;
762 mP[0] += k00 * z0;
763 mP[1] += k11 * z1;
764 mP[2] += k20 * z0;
765 mP[3] += k31 * z1;
766 mP[4] += k40 * z0;
767
768 mC[0] -= k00 * d00;
769 mC[2] -= k11 * d11;
770 mC[3] -= k20 * d00;
771 mC[5] -= k20 * d02;
772 mC[7] -= k31 * d11;
773 mC[9] -= k31 * d13;
774 mC[10] -= k00 * d04;
775 mC[12] -= k40 * d02;
776 mC[14] -= k40 * d04;
777 } else {
778 const float k00 = d00 * w0 + d01 * w1;
779 const float k01 = d00 * w1 + d10 * w2;
780 const float k10 = d01 * w0 + d11 * w1;
781 const float k11 = d01 * w1 + d11 * w2;
782 const float k20 = d02 * w0 + d12 * w1;
783 const float k21 = d02 * w1 + d12 * w2;
784 const float k30 = d03 * w0 + d13 * w1;
785 const float k31 = d03 * w1 + d13 * w2;
786 const float k40 = d04 * w0 + d14 * w1;
787 const float k41 = d04 * w1 + d14 * w2;
788
789 mP[0] += k00 * z0 + k01 * z1;
790 mP[1] += k10 * z0 + k11 * z1;
791 mP[2] += k20 * z0 + k21 * z1;
792 mP[3] += k30 * z0 + k31 * z1;
793 mP[4] += k40 * z0 + k41 * z1;
794
795 mC[0] -= k00 * d00 + k01 * d10;
796
797 mC[2] -= k10 * d01 + k11 * d11;
798
799 mC[3] -= k20 * d00 + k21 * d10;
800 mC[5] -= k20 * d02 + k21 * d12;
801
802 mC[7] -= k30 * d01 + k31 * d11;
803 mC[9] -= k30 * d03 + k31 * d13;
804
805 mC[10] -= k40 * d00 + k41 * d10;
806 mC[12] -= k40 * d02 + k41 * d12;
807 mC[14] -= k40 * d04 + k41 * d14;
808
809 mC[1] -= k10 * d00 + k11 * d10;
810 mC[4] -= k20 * d01 + k21 * d11;
811 mC[6] -= k30 * d00 + k31 * d10;
812 mC[8] -= k30 * d02 + k31 * d12;
813 mC[11] -= k40 * d01 + k41 * d11;
814 mC[13] -= k40 * d03 + k41 * d13;
815 }
816 return 0;
817}
818
819//*
820//* Multiple scattering and energy losses
821//*
822
823GPUd() float GPUTPCGMPropagator::ApproximateBetheBloch(float beta2)
824{
825 //------------------------------------------------------------------
826 // This is an approximation of the Bethe-Bloch formula with
827 // the density effect taken into account at beta*gamma > 3.5
828 // (the approximation is reasonable only for solid materials)
829 //------------------------------------------------------------------
830
831 const float log0 = CAMath::Log(5940.f);
832 const float log1 = CAMath::Log(3.5f * 5940.f);
833
834 bool bad = (beta2 >= .999f) || (beta2 < 1.e-8f);
835
836 if (bad) {
837 beta2 = 0.5f;
838 }
839
840 float a = beta2 / (1.f - beta2);
841 float b = 0.5f * CAMath::Log(a);
842 float d = 0.153e-3f / beta2;
843 float c = b - beta2;
844
845 float ret = d * (log0 + b + c);
846 float case1 = d * (log1 + c);
847
848 if (a > 3.5f * 3.5f) {
849 ret = case1;
850 }
851 if (bad) {
852 ret = 0.f;
853 }
854
855 return ret;
856}
857
858GPUd() void GPUTPCGMPropagator::CalculateMaterialCorrection()
859{
860 //*!
861
862 const float mass = 0.13957f;
863
864 float qpt = mT0.GetQPt();
865 if (CAMath::Abs(qpt) > 20) {
866 qpt = 20;
867 }
868
869 float w2 = (1.f + mT0.GetDzDs() * mT0.GetDzDs()); //==(P/pt)2
870 float pti2 = qpt * qpt;
871 if (pti2 < 1.e-4f) {
872 pti2 = 1.e-4f;
873 }
874
875 float mass2 = mass * mass;
876 float beta2 = w2 / (w2 + mass2 * pti2);
877
878 float p2 = w2 / pti2; // impuls 2
879 float betheRho = ApproximateBetheBloch(p2 / mass2) * mMaterial.rho;
880 float E = CAMath::Sqrt(p2 + mass2);
881 float theta2 = (14.1f * 14.1f / 1.e6f) / (beta2 * p2) * mMaterial.radLenInv;
882
883 mMaterial.EP2 = E / p2;
884
885 // Approximate energy loss fluctuation (M.Ivanov)
886
887 const float knst = 0.0007f; // To be tuned.
888 mMaterial.sigmadE2 = knst * mMaterial.EP2 * qpt;
889 mMaterial.sigmadE2 = mMaterial.sigmadE2 * mMaterial.sigmadE2;
890
891 mMaterial.k22 = theta2 * w2;
892 mMaterial.k33 = mMaterial.k22 * w2;
893 mMaterial.k43 = 0.f;
894 mMaterial.k44 = theta2 * mT0.GetDzDs() * mT0.GetDzDs() * pti2;
895
896 float br = (betheRho > 1.e-8f) ? betheRho : 1.e-8f;
897 mMaterial.DLMax = 0.3f * E / br;
898 mMaterial.EP2 *= betheRho;
899 mMaterial.sigmadE2 = mMaterial.sigmadE2 * betheRho; // + mMaterial.fK44;
900}
901
902GPUd() void GPUTPCGMPropagator::Rotate180()
903{
904 mT0.X() = -mT0.X();
905 mT0.Y() = -mT0.Y();
906 mT0.Q() = -mT0.Q();
907 mT0.Pz() = -mT0.Pz();
908 mT0.UpdateValues();
909
910 mT->X() = -mT->X();
911 mT->Y() = -mT->Y();
912 mT->QPt() = -mT->QPt();
913 mT->DzDs() = -mT->DzDs();
914
915 mAlpha = mAlpha + CAMath::Pi();
916 while (mAlpha >= CAMath::Pi()) {
917 mAlpha -= CAMath::TwoPi();
918 }
919 while (mAlpha < -CAMath::Pi()) {
920 mAlpha += CAMath::TwoPi();
921 }
922 mCosAlpha = -mCosAlpha;
923 mSinAlpha = -mSinAlpha;
924
925 float* c = mT->Cov();
926 c[6] = -c[6];
927 c[7] = -c[7];
928 c[8] = -c[8];
929 c[10] = -c[10];
930 c[11] = -c[11];
931 c[12] = -c[12];
932}
933
934GPUd() void GPUTPCGMPropagator::ChangeDirection()
935{
936 mT0.Py() = -mT0.Py();
937 mT0.Pz() = -mT0.Pz();
938 mT0.Q() = -mT0.Q();
939 mT->SinPhi() = -mT->SinPhi();
940 mT->DzDs() = -mT->DzDs();
941 mT->QPt() = -mT->QPt();
942 mT0.UpdateValues();
943
944 float* c = mT->Cov();
945 c[3] = -c[3];
946 c[4] = -c[4];
947 c[6] = -c[6];
948 c[7] = -c[7];
949 c[10] = -c[10];
950 c[11] = -c[11];
951}
952
953GPUd() void GPUTPCGMPropagator::Mirror(bool inFlyDirection)
954{
955 // mirror the track and the track approximation to the point which has the same X, but located on the other side of trajectory
956 float B[3];
957 GetBxByBz(mT0.X(), mT0.Y(), mT0.Z(), B);
958 float Bz = B[2];
959 if (CAMath::Abs(Bz) < 1.e-8f) {
960 Bz = 1.e-8f;
961 }
962
963 float dy = -2.f * mT0.Q() * mT0.Px() / Bz;
964 float dS; // path in XY
965 {
966 float chord = dy; // chord to the extrapolated point == |dy|*sign(x direction)
967 float sa = -mT0.CosPhi(); // sin( half of the rotation angle ) == (chord/2) / radius
968
969 // dS = (Pt/b)*2*arcsin( sa )
970 // = (Pt/b)*2*sa*(1 + 1/6 sa^2 + 3/40 sa^4 + 5/112 sa^6 +... )
971 // = chord*(1 + 1/6 sa^2 + 3/40 sa^4 + 5/112 sa^6 +... )
972
973 float sa2 = sa * sa;
974 const float k2 = 1.f / 6.f;
975 const float k4 = 3.f / 40.f;
976 // const float k6 = 5.f/112.f;
977 dS = chord + chord * sa2 * (k2 + k4 * sa2);
978 // dS = CAMath::Sqrt(pt2)/b*2.f*CAMath::ASin( sa );
979 }
980
981 if (mT0.SinPhi() < 0.f) {
982 dS = -dS;
983 }
984
985 mT0.Y() = mT0.Y() + dy;
986 mT0.Z() = mT0.Z() + mT0.DzDs() * dS;
987 mT0.Px() = mT0.Px(); // should be positive
988 mT->Y() = mT->Y() + dy;
989 mT->Z() = mT->Z() + mT0.DzDs() * dS;
990 ChangeDirection();
991
992 // Energy Loss
993 if (true) {
994 // std::cout<<"MIRROR: APPLY ENERGY LOSS!!!"<<std::endl;
995
996 float dL = CAMath::Abs(dS * mT0.GetDlDs());
997
998 if (inFlyDirection) {
999 dL = -dL;
1000 }
1001
1002 float* c = mT->Cov();
1003 float& mC40 = c[10];
1004 float& mC41 = c[11];
1005 float& mC42 = c[12];
1006 float& mC43 = c[13];
1007 float& mC44 = c[14];
1008
1009 float dLmask = 0.f;
1010 bool maskMS = (CAMath::Abs(dL) < mMaterial.DLMax);
1011 if (maskMS) {
1012 dLmask = dL;
1013 }
1014 float dLabs = CAMath::Abs(dLmask);
1015 float corr = 1.f - mMaterial.EP2 * dLmask;
1016
1017 float corrInv = 1.f / corr;
1018 mT0.Px() *= corrInv;
1019 mT0.Py() *= corrInv;
1020 mT0.Pz() *= corrInv;
1021 mT0.Pt() *= corrInv;
1022 mT0.P() *= corrInv;
1023 mT0.QPt() *= corr;
1024
1025 mT->QPt() *= corr;
1026
1027 mC40 *= corr;
1028 mC41 *= corr;
1029 mC42 *= corr;
1030 mC43 *= corr;
1031 mC44 = mC44 * corr * corr + dLabs * mMaterial.sigmadE2;
1032 } else {
1033 // std::cout<<"MIRROR: DONT APPLY ENERGY LOSS!!!"<<std::endl;
1034 }
1035}
1036
1037GPUd() o2::base::MatBudget GPUTPCGMPropagator::getMatBudget(const float* p1, const float* p2)
1038{
1039 return mMatLUT->getMatBudget(p1[0], p1[1], p1[2], p2[0], p2[1], p2[2]);
1040}
1041
1042GPUdic(0, 1) void GPUTPCGMPropagator::UpdateMaterial(const GPUTPCGMPhysicalTrackModel& GPUrestrict() t0e)
1043{
1044 float xyz1[3] = {getGlobalX(mT0.GetX(), mT0.GetY()), getGlobalY(mT0.GetX(), mT0.GetY()), mT0.GetZ()};
1045 float xyz2[3] = {getGlobalX(t0e.GetX(), t0e.GetY()), getGlobalY(t0e.GetX(), t0e.GetY()), t0e.GetZ()};
1046 o2::base::MatBudget mat = getMatBudget(xyz1, xyz2);
1047 if (mat.meanX2X0 > 1.e-8) {
1048 SetMaterial(mat.length / mat.meanX2X0, mat.meanRho);
1049 } else {
1050 SetMaterialTPC();
1051 }
1052}
charge
int16_t charge
Definition RawEventData.h:5
time
int16_t time
Definition RawEventData.h:4
GPUdic
#define GPUdic(...)
Definition GPUCommonDefAPI.h:188
GPUrestrict
#define GPUrestrict()
Definition GPUCommonDefAPI.h:212
GPUCA_DEBUG_STREAMER_CHECK
#define GPUCA_DEBUG_STREAMER_CHECK(...)
Definition GPUCommonDef.h:68
GPUCA_MAX_SIN_PHI
#define GPUCA_MAX_SIN_PHI
Definition GPUDefConstantsAndSettings.h:42
retVal
int32_t retVal
Definition GPUReconstruction.cxx:60
bb
const int16_t bb
Definition GPUTPCGMMerger.cxx:1479
GPUd
GPUd() void GPUTPCGMPropagator
Definition GPUTPCGMPropagator.cxx:27
GPUCA_NSECTORS
#define GPUCA_NSECTORS
Definition GPUTPCGeometry.h:22
p2
constexpr int p2()
Definition HmpidCoder2.cxx:90
p1
constexpr int p1()
constexpr to accelerate the coordinates changing
Definition HmpidCoder2.cxx:89
c
uint32_t c
Definition RawData.h:2
B
Definition B.h:16
x
GLint GLenum GLint x
Definition glcorearb.h:403
b
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
void
typedef void(APIENTRYP PFNGLCULLFACEPROC)(GLenum mode)
x0
GLuint GLfloat x0
Definition glcorearb.h:5034
t0
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t0
Definition glcorearb.h:5034
param
GLenum GLfloat param
Definition glcorearb.h:271
a
GLboolean GLboolean GLboolean GLboolean a
Definition glcorearb.h:1233
y0
GLuint GLfloat GLfloat y0
Definition glcorearb.h:5034
o2::event_visualisation::TRD
@ TRD
Definition VisualisationConstants.h:31
o2::event_visualisation::ITS
@ ITS
Definition VisualisationConstants.h:29
o2::its::math_utils::p
float p
Definition MathUtils.h:99
o2::its::math_utils::k2
const float k2
Definition MathUtils.h:72
o2::its::math_utils::z1
float float float float z1
Definition MathUtils.h:79
o2::track::d2
value_T d2
Definition TrackUtils.h:135
o2::track::beta2
kp1 *kp2 *value_T beta2
Definition TrackUtils.h:131
o2::trigger::TPC
constexpr uint32_t TPC
Definition Triggers.h:39
o2
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
Definition BitstreamReader.h:24
o2::base::MatBudget::length
float length
length in material
Definition MatCell.h:55
o2::base::MatCell::meanRho
float meanRho
mean density, g/cm^3
Definition MatCell.h:30
o2::base::MatCell::meanX2X0
float meanX2X0
fraction of radiaton lenght
Definition MatCell.h:31
cc
std::vector< o2::mch::ChannelCode > cc
Definition testStatusMap.cxx:37