22 #ifndef PNDCATRACKPARAM_H 23 #define PNDCATRACKPARAM_H 46 for (
int j = 0; j < 5; ++j)
47 fP[j] = v.
Par()[j][
i];
48 for (
int j = 0; j < 15; ++j)
49 fC[j] = v.
Cov()[j][
i];
64 float X()
const {
return fX; }
65 float Y()
const {
return fP[0]; }
66 float Z()
const {
return fP[1]; }
67 float SinPhi()
const {
return fP[2]; }
68 float DzDs()
const {
return fP[3]; }
69 float QPt()
const {
return fP[4]; }
77 float Chi2()
const {
return fChi2; }
78 int NDF()
const {
return fNDF; }
80 float Err2Y()
const {
return fC[0]; }
81 float Err2Z()
const {
return fC[2]; }
86 float Angle()
const {
return fAlpha; }
87 int ISec()
const {
return fISec; }
89 float GetX()
const {
return fX; }
90 float GetY()
const {
return fP[0]; }
91 float GetZ()
const {
return fP[1]; }
94 float GetQPt()
const {
return fP[4]; }
99 float GetKappa(
float Bz)
const {
return fP[4] * Bz; }
113 const float *
Par()
const {
return fP; }
114 const float *
Cov()
const {
return fC; }
116 const float *
GetPar()
const {
return fP; }
117 const float *
GetCov()
const {
return fC; }
139 const float r =
sqrt(r0 * r0 + r1 * r1);
148 float GetS(
float x,
float y,
float Bz)
const;
150 void GetDCAPoint(
float x,
float y,
float z,
float &px,
float &py,
float &
pz,
float Bz)
const;
152 bool TransportToX(
float x,
float Bz,
float maxSinPhi = .999);
157 bool TransportToX(
float x,
float sinPhi0,
float cosPhi0,
float Bz,
float maxSinPhi = .999);
164 static float BetheBlochGeant(
float bg,
float kp0 = 2.33,
float kp1 = 0.20,
float kp2 = 3.00,
float kp3 = 173e-9,
float kp4 = 0.49848);
174 void RotateXY(
float alpha,
float &x,
float &y,
float &
sin)
const;
175 bool Filter(
float y,
float z,
float err2Y,
float errYZ,
float err2Z,
float maxSinPhi = .999);
204 for (
int i = 0;
i < 5;
i++)
206 for (
int i = 0;
i < 15;
i++)
238 x = (
X() * cA +
Y() * sA);
239 y = (-
X() * sA +
Y() * cA);
245 assert(maxSinPhi > 0.
f);
248 const float c00 = fC[0];
249 const float c10 = fC[1];
250 const float c11 = fC[2];
251 const float c20 = fC[3];
252 const float c21 = fC[4];
254 const float c30 = fC[6];
255 const float c31 = fC[7];
258 const float c40 = fC[10];
259 const float c41 = fC[11];
264 float d = 1.f / (err2Y * err2Z + err2Y * c11 + err2Z * c00 + c00 * c11 - c10 * c10 - 2 * errYZ * c10 - errYZ * errYZ);
269 const float z0 = y - fP[0], z1 = z - fP[1];
274 const float mS0 = err2Z * d;
275 const float mS1 = -errYZ * d;
276 const float mS2 = err2Y * d;
280 const float k00 = c00 * mS0 + c10 * mS1, k01 = c00 * mS1 + c10 * mS2, k10 = c10 * mS0 + c11 * mS1, k11 = c10 * mS1 + c11 * mS2, k20 = c20 * mS0 + c21 * mS1,
281 k21 = c20 * mS1 + c21 * mS2, k30 = c30 * mS0 + c31 * mS1, k31 = c30 * mS1 + c31 * mS2, k40 = c40 * mS0 + c41 * mS1, k41 = c40 * mS1 + c41 * mS2;
283 const float sinPhi = fP[2] + k20 * z0 + k21 * z1;
289 fChi2 += mS0 * z0 * z0 + mS2 * z1 * z1 + 2 * z0 * z1 * mS1;
291 fP[0] += k00 * z0 + k01 * z1;
292 fP[1] += k10 * z0 + k11 * z1;
294 fP[3] += k30 * z0 + k31 * z1;
295 fP[4] += k40 * z0 + k41 * z1;
297 fC[0] -= (k00 * c00 + k01 * c10);
299 fC[1] -= (k10 * c00 + k11 * c10);
300 fC[2] -= (k10 * c10 + k11 * c11);
302 fC[3] -= (k20 * c00 + k21 * c10);
303 fC[4] -= (k20 * c10 + k21 * c11);
304 fC[5] -= (k20 * c20 + k21 * c21);
306 fC[6] -= (k30 * c00 + k31 * c10);
307 fC[7] -= (k30 * c10 + k31 * c11);
308 fC[8] -= (k30 * c20 + k31 * c21);
309 fC[9] -= (k30 * c30 + k31 * c31);
311 fC[10] -= (k40 * c00 + k41 * c10);
312 fC[11] -= (k40 * c10 + k41 * c11);
313 fC[12] -= (k40 * c20 + k41 * c21);
314 fC[13] -= (k40 * c30 + k41 * c31);
315 fC[14] -= (k40 * c40 + k41 * c41);
330 const float beta2_beta21i = beta2 / (1 - beta2);
331 if (beta2_beta21i > 12.25)
332 return 0.153e-3 / beta2 * (9.94223 + 0.5 *
log(beta2_beta21i) - beta2);
334 return 0.153e-3 / beta2 * (8.6895 +
log(beta2_beta21i) - beta2);
340 const float p2 = (1. + fP[3] * fP[3]);
341 const float k2 = fP[4] * fP[4];
342 const float mass2 = mass * mass;
344 const float beta2 = p2 / (p2 + mass2 * k2);
346 const float pp2 = (k2 > 1.e-8) ? p2 / k2 : 10000;
351 par.
fTheta2 = 198.81e-6 / (beta2 * pp2);
356 const float knst = 0.07;
362 par.
fK43 = fP[3] * fP[4] * par.
fK22;
363 par.
fK44 = (p2 - 1.f) * k2;
377 const float ex = t0.
CosPhi();
378 const float ey = t0.
SinPhi();
379 const float k = t0.
QPt() * Bz;
380 const float dx = x -
X();
382 const float ey1 = k * dx + ey;
393 const float dx2 = dx * dx;
394 const float ss = ey + ey1;
395 const float cc = ex + ex1;
400 const float cci = 1.f / cc;
401 const float exi = 1.f / ex;
402 const float ex1i = 1.f / ex1;
404 const float tg = ss * cci;
406 const float dy = dx * tg;
411 float dSin = dl * k * 0.5;
417 const float dz = dS * t0.
DzDs();
422 const float d[3] = {fP[2] - t0.
SinPhi(), fP[3] - t0.
DzDs(), fP[4] - t0.
QPt()};
430 const float h2 = dx * (1.f + ey * ey1 + ex * ex1) * exi * ex1i * cci;
431 const float h4 = dx2 * (cc + ss * ey1 * ex1i) * cci * cci * Bz;
432 const float dxBz = dx * Bz;
438 fP[0] =
Y() + dy + h2 * d[0] + h4 * d[2];
439 fP[1] =
Z() + dz + dS * d[1];
440 fP[2] = t0.
SinPhi() + d[0] + dxBz * d[2];
445 const float c00 = fC[0];
446 const float c10 = fC[1];
447 const float c11 = fC[2];
448 const float c20 = fC[3];
449 const float c21 = fC[4];
450 const float c22 = fC[5];
451 const float c30 = fC[6];
452 const float c31 = fC[7];
453 const float c32 = fC[8];
454 const float c33 = fC[9];
455 const float c40 = fC[10];
456 const float c41 = fC[11];
457 const float c42 = fC[12];
458 const float c43 = fC[13];
459 const float c44 = fC[14];
461 fC[0] = (c00 + h2 * h2 * c22 + h4 * h4 * c44 + 2.f * (h2 * c20 + h4 * c40 + h2 * h4 * c42));
463 fC[1] = c10 + h2 * c21 + h4 * c41 + dS * (c30 + h2 * c32 + h4 * c43);
464 fC[2] = c11 + 2.f * dS * c31 + dS * dS * c33;
466 fC[3] = c20 + h2 * c22 + h4 * c42 + dxBz * (c40 + h2 * c42 + h4 * c44);
467 fC[4] = c21 + dS * c32 + dxBz * (c41 + dS * c43);
468 fC[5] = c22 + 2.f * dxBz * c42 + dxBz * dxBz * c44;
470 fC[6] = c30 + h2 * c32 + h4 * c43;
471 fC[7] = c31 + dS * c33;
472 fC[8] = c32 + dxBz * c43;
475 fC[10] = c40 + h2 * c42 + h4 * c44;
476 fC[11] = c41 + dS * c43;
477 fC[12] = c42 + dxBz * c44;
481 fC[0] = (fC[0] + h2 * h2 * fC[5] + h4 * h4 * fC[14] + 2 * (h2 * fC[3] + h4 * fC[10] + h2 * h4 * fC[12]));
483 fC[1] = fC[1] + h2 * fC[4] + h4 * fC[11] + dS * (fC[6] + h2 * fC[8] + h4 * fC[13]);
484 fC[2] = fC[2] + 2 * dS * fC[7] + dS * dS * fC[9];
486 fC[3] = fC[3] + h2 * fC[5] + h4 * fC[12] + dxBz * (fC[10] + h2 * fC[12] + h4 * fC[14]);
487 fC[4] = fC[4] + dS * fC[8] + dxBz * (fC[11] + dS * fC[13]);
488 fC[5] = fC[5] + 2 * dxBz * fC[12] + dxBz * dxBz * fC[14];
490 fC[6] = fC[6] + h2 * fC[8] + h4 * fC[13];
491 fC[7] = fC[7] + dS * fC[9];
492 fC[8] = fC[8] + dxBz * fC[13];
495 fC[10] = fC[10] + h2 * fC[12] + h4 * fC[14];
496 fC[11] = fC[11] + dS * fC[13];
497 fC[12] = fC[12] + dxBz * fC[14];
524 const float dE = par.
fBethe * xTimesRho;
527 const float corr = (1. - par.
fEP2 * dE);
528 if (corr < 0.3 || corr > 1.3)
536 fC[14] *= corr * corr;
543 fC[5] += theta2 * par.
fK22 * (1. - fP[2] * fP[2]);
544 fC[9] += theta2 * par.
fK33;
545 fC[13] += theta2 * par.
fK43;
546 fC[14] += theta2 * par.
fK44;
555 const float kRho = 1.54e-3;
558 const float kRhoOverRadLen = 7.68e-5;
564 UNUSED_PARAM3(kRho, kRhoOverRadLen, par);
583 const bool rotated =
Rotate(-fAlpha + hit.
Angle(), tR, .999f);
586 return rotated & transported;
static T ASin(const T &x)
bool CorrectForMeanMaterial(float xOverX0, float xTimesRho, const PndCATrackFitParam &par)
float GetDist2(const PndCATrackParam &t) const
bool Filter(float y, float z, float err2Y, float errYZ, float err2Z, float maxSinPhi=.999)
float GetS(float x, float y, float Bz) const
const float * GetCov() const
const float * Cov() const
void SetCov(int i, float v)
bool Transport(const PndCAHit &hit, float Bz)
bool Rotate(float alpha, float maxSinPhi=.999)
friend F32vec4 sqrt(const F32vec4 &a)
static float BetheBlochGeant(float bg, float kp0=2.33, float kp1=0.20, float kp2=3.00, float kp3=173e-9, float kp4=0.49848)
static T Sqrt(const T &x)
friend F32vec4 sin(const F32vec4 &a)
PndCATrackParam TrackParam
void SetSignCosPhi(float v)
friend F32vec4 log(const F32vec4 &a)
void GetDCAPoint(float x, float y, float z, float &px, float &py, float &pz, float Bz) const
const float * Par() const
bool TransportToXWithMaterial(float x, float Bz, float maxSinPhi=.999)
const float * GetPar() const
void InitDirection(float r0, float r1, float r2)
float GetErr2SinPhi() const
static float BetheBlochSolid(float bg)
void SetPar(int i, float v)
float GetSignCosPhi() const
float GetDistXZ2(const PndCATrackParam &t) const
float Err2QMomentum() const
float GetKappa(float Bz) const
bool TransportToX(float x, float Bz, float maxSinPhi=.999)
static float BetheBlochGas(float bg)
PndCATrackParam GetGlobalParam(float alpha) const
static float ApproximateBetheBloch(float beta2)
float GetErr2DzDs() const
const float_v & Cov(int i) const
void CalculateFitParameters(PndCATrackFitParam &par, float mass=0.13957)
float GetCosPhiPositive() const
PndCATrackParam(const TrackParamVector &v, int i)
const float_v & Par(int i) const
void RotateXY(float alpha, float &x, float &y, float &sin) const