11 #ifndef PNDCATRACKPARAM_H
12 #define PNDCATRACKPARAM_H
43 for (
int j = 0; j < 5; ++j )
fP[j] = v.
Par()[j][
i];
44 for (
int j = 0; j < 15; ++j )
fC[j] = v.
Cov()[j][
i];
59 float X()
const {
return fX; }
60 float Y()
const {
return fP[0]; }
61 float Z()
const {
return fP[1]; }
64 float QPt()
const {
return fP[4]; }
108 const float *
Par()
const {
return fP; }
109 const float *
Cov()
const {
return fC; }
134 const float r =
sqrt( r0*r0+r1*r1 );
144 float GetS(
float x,
float y,
float Bz )
const;
147 float &px,
float &py,
float &
pz,
float Bz )
const;
150 bool TransportToX(
float x,
float Bz,
float maxSinPhi = .999 );
154 float Bz,
float maxSinPhi = .999,
float *DL = 0 );
156 bool TransportToX(
float x,
float sinPhi0,
float cosPhi0,
float Bz,
float maxSinPhi = .999 );
160 PndCATrackFitParam &
par,
float Bz,
float maxSinPhi = .999 );
163 PndCATrackFitParam &
par,
float Bz,
float maxSinPhi = .999 );
182 bool Rotate(
float alpha,
float maxSinPhi = .999 );
186 bool Filter(
float y,
float z,
float err2Y,
float errYZ,
float err2Z,
float maxSinPhi = .999 );
194 fC[1] = 0.f;
fC[2] = 10.f;
195 fC[3] = 0.f;
fC[4] = 0.f;
fC[5] = 1.f;
196 fC[6] = 0.f;
fC[7] = 0.f;
fC[8] = 0.f;
fC[9] = 1.f;
197 fC[10] = 0.f;
fC[11] = 0.f;
fC[12] = 0.f;
fC[13] = 0.f;
fC[14] = 10.f;
231 x = (
X()*cA +
Y()*sA );
232 y = ( -
X()*sA +
Y()*cA );
239 assert( maxSinPhi > 0.
f );
242 const float c00 =
fC[0];
243 const float c10 =
fC[1];
244 const float c11 =
fC[2];
245 const float c20 =
fC[3];
246 const float c21 =
fC[4];
248 const float c30 =
fC[6];
249 const float c31 =
fC[7];
252 const float c40 =
fC[10];
253 const float c41 =
fC[11];
258 float d = 1.f / ( err2Y*err2Z + err2Y*c11 + err2Z*c00 + c00*c11 - c10*c10 - 2*errYZ*c10 - errYZ*errYZ );
269 const float mS0 = err2Z*
d;
270 const float mS1 = -errYZ*
d;
271 const float mS2 = err2Y*
d;
276 k00 = c00 * mS0 + c10*mS1, k01 = c00 * mS1 + c10*mS2,
277 k10 = c10 * mS0 + c11*mS1, k11 = c10 * mS1 + c11*mS2,
278 k20 = c20 * mS0 + c21*mS1, k21 = c20 * mS1 + c21*mS2,
279 k30 = c30 * mS0 + c31*mS1, k31 = c30 * mS1 + c31*mS2,
280 k40 = c40 * mS0 + c41*mS1, k41 = c40 * mS1 + c41*mS2;
282 const float sinPhi = fP[2] + k20 * z0 + k21 * z1;
287 fChi2 += mS0 * z0 * z0 + mS2 * z1 * z1 + 2 * z0 * z1 * mS1;
289 fP[ 0] += k00 * z0 + k01 * z1;
290 fP[ 1] += k10 * z0 + k11 * z1;
292 fP[ 3] += k30 * z0 + k31 * z1;
293 fP[ 4] += k40 * z0 + k41 * z1;
295 fC[ 0] -= (k00 * c00 + k01 *
c10);
297 fC[ 1] -= (k10 * c00 + k11 *
c10);
298 fC[ 2] -= (k10 * c10 + k11 *
c11);
300 fC[ 3] -= (k20 * c00 + k21 *
c10);
301 fC[ 4] -= (k20 * c10 + k21 *
c11);
302 fC[ 5] -= (k20 * c20 + k21 *
c21);
304 fC[ 6] -= (k30 * c00 + k31 *
c10);
305 fC[ 7] -= (k30 * c10 + k31 *
c11);
306 fC[ 8] -= (k30 * c20 + k31 *
c21);
307 fC[ 9] -= (k30 * c30 + k31 * c31);
309 fC[10] -= (k40 * c00 + k41 *
c10);
310 fC[11] -= (k40 * c10 + k41 *
c11);
311 fC[12] -= (k40 * c20 + k41 *
c21);
312 fC[13] -= (k40 * c30 + k41 * c31);
313 fC[14] -= (k40 * c40 + k41 * c41);
328 const float beta2_beta21i = beta2 / ( 1 - beta2 );
329 if ( beta2_beta21i > 12.25 )
330 return 0.153e-3 / beta2 * ( 9.94223 + 0.5 *
log( beta2_beta21i ) - beta2 );
332 return 0.153e-3 / beta2 * ( 8.6895 +
log( beta2_beta21i ) - beta2 );
338 const float p2 = ( 1. +
fP[3] *
fP[3] );
339 const float k2 =
fP[4] *
fP[4];
340 const float mass2 = mass * mass;
342 const float beta2 = p2 / ( p2 + mass2 * k2 );
344 const float pp2 = ( k2 > 1.e-8 ) ? p2 / k2 : 10000;
349 par.
fTheta2 = 198.81e-6 / ( beta2 * pp2 );
354 const float knst = 0.07;
360 par.
fK43 = fP[3] * fP[4] * par.
fK22;
361 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;
389 if ( ex < 0 ) ex1 = -ex1;
391 const float dx2 = dx *
dx;
392 const float ss = ey + ey1;
393 const float cc = ex + ex1;
397 const float cci = 1.f / cc;
398 const float exi = 1.f / ex;
399 const float ex1i = 1.f / ex1;
401 const float tg = ss * cci;
403 const float dy = dx * tg;
406 if ( cc < 0 ) dl = -dl;
407 float dSin = dl * k * 0.5;
408 if ( dSin > 1.
f ) dSin = 1.f;
409 if ( dSin < -1.
f ) dSin = -1.f;
411 const float dz = dS * t0.
DzDs();
424 const float h2 = dx * ( 1.f + ey * ey1 + ex * ex1 ) * exi * ex1i * cci;
425 const float h4 = dx2 * ( cc + ss * ey1 * ex1i ) * cci * cci * Bz;
426 const float dxBz = dx * Bz;
432 fP[0] =
Y() + dy + h2 * d[0] + h4 * d[2];
433 fP[1] =
Z() + dz + dS * d[1];
434 fP[2] = t0.
SinPhi() + d[0] + dxBz * d[2];
438 const float c00 =
fC[0];
439 const float c10 =
fC[1];
440 const float c11 =
fC[2];
441 const float c20 =
fC[3];
442 const float c21 =
fC[4];
443 const float c22 =
fC[5];
444 const float c30 =
fC[6];
445 const float c31 =
fC[7];
446 const float c32 =
fC[8];
447 const float c33 =
fC[9];
448 const float c40 =
fC[10];
449 const float c41 =
fC[11];
450 const float c42 =
fC[12];
451 const float c43 =
fC[13];
452 const float c44 =
fC[14];
454 fC[0] = ( c00 + h2 * h2 * c22 + h4 * h4 * c44
455 + 2.f * ( h2 * c20 + h4 * c40 + h2 * h4 * c42 ) );
457 fC[1] = c10 + h2 * c21 + h4 * c41 + dS * ( c30 + h2 * c32 + h4 * c43 );
458 fC[2] = c11 + 2.f * dS * c31 + dS * dS * c33;
460 fC[3] = c20 + h2 * c22 + h4 * c42 + dxBz * ( c40 + h2 * c42 + h4 * c44 );
461 fC[4] = c21 + dS * c32 + dxBz * ( c41 + dS * c43 );
462 fC[5] = c22 + 2.f * dxBz * c42 + dxBz * dxBz * c44;
464 fC[6] = c30 + h2 * c32 + h4 * c43;
465 fC[7] = c31 + dS * c33;
466 fC[8] = c32 + dxBz * c43;
469 fC[10] = c40 + h2 * c42 + h4 * c44;
470 fC[11] = c41 + dS * c43;
471 fC[12] = c42 + dxBz * c44;
475 fC[0] = ( fC[0] + h2 * h2 * fC[5] + h4 * h4 * fC[14]
476 + 2 * ( h2 * fC[3] + h4 * fC[10] + h2 * h4 * fC[12] ) );
478 fC[1] = fC[1] + h2 * fC[4] + h4 * fC[11] + dS * ( fC[6] + h2 * fC[8] + h4 * fC[13] );
479 fC[2] = fC[2] + 2 * dS * fC[7] + dS * dS * fC[9];
481 fC[3] = fC[3] + h2 * fC[5] + h4 * fC[12] + dxBz * ( fC[10] + h2 * fC[12] + h4 * fC[14] );
482 fC[4] = fC[4] + dS * fC[8] + dxBz * ( fC[11] + dS * fC[13] );
483 fC[5] = fC[5] + 2 * dxBz * fC[12] + dxBz * dxBz * fC[14];
485 fC[6] = fC[6] + h2 * fC[8] + h4 * fC[13];
486 fC[7] = fC[7] + dS * fC[9];
487 fC[8] = fC[8] + dxBz * fC[13];
490 fC[10] = fC[10] + h2 * fC[12] + h4 * fC[14];
491 fC[11] = fC[11] + dS * fC[13];
492 fC[12] = fC[12] + dxBz * fC[14];
518 const float dE = par.
fBethe * xTimesRho;
521 if ( corr < 0.3 || corr > 1.3 )
return 0;
535 fC[5] += theta2 * par.
fK22 * ( 1. -
fP[2] *
fP[2] );
536 fC[9] += theta2 * par.
fK33;
537 fC[13] += theta2 * par.
fK43;
538 fC[14] += theta2 * par.
fK44;
547 const float kRho = 1.54e-3;
550 const float kRhoOverRadLen = 7.68e-5;
553 if ( !
TransportToX( x, t0, Bz, maxSinPhi, &dl ) )
return 0;
555 UNUSED_PARAM3(kRho, kRhoOverRadLen, par);
578 return rotated & transported;
static T ASin(const T &x)
bool CorrectForMeanMaterial(float xOverX0, float xTimesRho, const PndCATrackFitParam &par)
void GetDCAPoint(float x, float y, float z, float &px, float &py, float &pz, float Bz) const
float GetErr2SinPhi() const
static float BetheBlochSolid(float bg)
bool Filter(float y, float z, float err2Y, float errYZ, float err2Z, float maxSinPhi=.999)
void SetCov(int i, float v)
bool Transport(const PndCAHit &hit, float Bz)
const float_v & Cov(int i) const
float GetKappa(float Bz) const
bool Rotate(float alpha, float maxSinPhi=.999)
friend F32vec4 sqrt(const F32vec4 &a)
static T Sqrt(const T &x)
const float * Par() const
PndCATrackParam GetGlobalParam(float alpha) const
friend F32vec4 sin(const F32vec4 &a)
PndCATrackParam TrackParam
float GetSignCosPhi() const
void SetSignCosPhi(float v)
friend F32vec4 log(const F32vec4 &a)
float GetCosPhiPositive() const
float GetErr2DzDs() const
bool TransportToXWithMaterial(float x, float Bz, float maxSinPhi=.999)
void RotateXY(float alpha, float &x, float &y, float &sin) const
void InitDirection(float r0, float r1, float r2)
float GetDistXZ2(const PndCATrackParam &t) const
float GetS(float x, float y, float Bz) const
const float * GetPar() const
const float * Cov() const
float GetDist2(const PndCATrackParam &t) const
void SetPar(int i, float v)
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)
int hit(Int_t nEvents=0, TString inFile="sim.root", TString parFile="par.root", TString inDigi="digi.root", TString outFile="hit.root", Int_t timeBased=0)
bool TransportToX(float x, float Bz, float maxSinPhi=.999)
const float_v & Par(int i) const
static float ApproximateBetheBloch(float beta2)
const float * GetCov() const
static float BetheBlochGas(float bg)
void CalculateFitParameters(PndCATrackFitParam &par, float mass=0.13957)
PndCATrackParam(const TrackParamVector &v, int i)
float Err2QMomentum() const