HLT/TPCLib/tracking-ca/AliHLTTPCCATrackParam.cxx

   1 // $Id$
   2 //***************************************************************************
   3 // This file is property of and copyright by the ALICE HLT Project          *
   4 // ALICE Experiment at CERN, All rights reserved.                           *
   5 //                                                                          *
   6 // Primary Authors: Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de> *
   7 //                  Ivan Kisel <kisel@kip.uni-heidelberg.de>                *
   8 //                  for The ALICE HLT Project.                              *
   9 //                                                                          *
  10 // Permission to use, copy, modify and distribute this software and its     *
  11 // documentation strictly for non-commercial purposes is hereby granted     *
  12 // without fee, provided that the above copyright notice appears in all     *
  13 // copies and that both the copyright notice and this permission notice     *
  14 // appear in the supporting documentation. The authors make no claims       *
  15 // about the suitability of this software for any purpose. It is            *
  16 // provided "as is" without express or implied warranty.                    *
  17 //***************************************************************************
  18
  19 #include "AliHLTTPCCATrackParam.h"
  20 #include "TMath.h"
  21 #include "Riostream.h"
  22
  23 //ClassImp(AliHLTTPCCATrackParam)
  24
  25 //
  26 // Circle in XY:
  27 //
  28 // R  = 1/TMath::Abs(Kappa);
  29 // Xc = X - sin(Phi)/Kappa;
  30 // Yc = Y + cos(Phi)/Kappa;
  31 //
  32
  33
  34
  35 void AliHLTTPCCATrackParam::ConstructXY3( const Float_t x[3], const Float_t y[3],
  36                                           const Float_t sigmaY2[3], Float_t CosPhi0 )
  37 {
  38   //* Construct the track in XY plane by 3 points
  39
  40   Float_t x0 = x[0];
  41   Float_t y0 = y[0];
  42   Float_t x1 = x[1] - x0;
  43   Float_t y1 = y[1] - y0;
  44   Float_t x2 = x[2] - x0;
  45   Float_t y2 = y[2] - y0;
  46
  47   Float_t a1 = x1*x1 + y1*y1;
  48   Float_t a2 = x2*x2 + y2*y2;
  49   Float_t a = 2*(x1*y2 - y1*x2);
  50   Float_t lx =  a1*y2 - a2*y1;
  51   Float_t ly = -a1*x2 + a2*x1;
  52   Float_t l = TMath::Sqrt(lx*lx + ly*ly);
  53
  54   Float_t li = 1./l;
  55   Float_t li2 = li*li;
  56   Float_t li3 = li2*li;
  57   Float_t cosPhi = ly*li;
  58
  59   Float_t sinPhi = -lx*li;
  60   Float_t kappa = a/l;
  61
  62   Float_t dlx = a2 - a1;     //  D lx / D y0
  63   Float_t dly = -a;          //  D ly / D y0
  64   Float_t dA  = 2*(x2 - x1); // D a  / D y0
  65   Float_t dl = (lx*dlx + ly*dly)*li;
  66
  67   // D sinPhi,kappa / D y0
  68
  69   Float_t d0[2] = { -(dlx*ly-lx*dly)*ly*li3, (dA*l-a*dl)*li2 };
  70
  71   // D sinPhi,kappa / D y1
  72
  73   dlx = -a2 + 2*y1*y2;
  74   dly = -2*x2*y1;
  75   dA  = -2*x2;
  76   dl = (lx*dlx + ly*dly)*li;
  77
  78   Float_t d1[2] = { -(dlx*ly-lx*dly)*ly*li3, (dA*l-a*dl)*li2 };
  79
  80   // D sinPhi,kappa / D y2
  81
  82   dlx = a1 - 2*y1*y2;
  83   dly = -2*x1*y2;
  84   dA  = 2*x1;
  85   dl = (lx*dlx + ly*dly)*li;
  86
  87   Float_t d2[2] = { -(dlx*ly-lx*dly)*ly*li3, (dA*l-a*dl)*li2 };
  88
  89   if( CosPhi0*cosPhi <0 ){
  90     cosPhi = -cosPhi;
  91     sinPhi = -sinPhi;
  92     kappa = -kappa;
  93     d0[0] = -d0[0];
  94     d0[1] = -d0[1];
  95     d1[0] = -d1[0];
  96     d1[1] = -d1[1];
  97     d2[0] = -d2[0];
  98     d2[1] = -d2[1];
  99   }
 100
 101   X() = x0;
 102   Y() = y0;
 103   SinPhi() = sinPhi;
 104   Kappa() = kappa;
 105   CosPhi() = cosPhi;
 106
 107   Float_t s0 = sigmaY2[0];
 108   Float_t s1 = sigmaY2[1];
 109   Float_t s2 = sigmaY2[2];
 110
 111   fC[0] = s0;
 112   fC[1] = 0;
 113   fC[2] = 100.;
 114
 115   fC[3] = d0[0]*s0;
 116   fC[4] = 0;
 117   fC[5] = d0[0]*d0[0]*s0 + d1[0]*d1[0]*s1 + d2[0]*d2[0]*s2;
 118
 119   fC[6] = 0;
 120   fC[7] = 0;
 121   fC[8] = 0;
 122   fC[9] = 100.;
 123
 124   fC[10] = d0[1]*s0;
 125   fC[11] = 0;
 126   fC[12] = d0[0]*d0[1]*s0 + d1[0]*d1[1]*s1 + d2[0]*d2[1]*s2;
 127   fC[13] = 0;
 128   fC[14] = d0[1]*d0[1]*s0 + d1[1]*d1[1]*s1 + d2[1]*d2[1]*s2;
 129 }
 130
 131
 132 Float_t  AliHLTTPCCATrackParam::GetS( Float_t x, Float_t y ) const
 133 {
 134   //* Get XY path length to the given point
 135
 136   Float_t k  = GetKappa();
 137   Float_t ex = GetCosPhi();
 138   Float_t ey = GetSinPhi();
 139   x-= GetX();
 140   y-= GetY();
 141   Float_t dS = x*ex + y*ey;
 142   if( TMath::Abs(k)>1.e-4 ) dS = TMath::ATan2( k*dS, 1+k*(x*ey-y*ex) )/k;
 143   return dS;
 144 }
 145
 146 void  AliHLTTPCCATrackParam::GetDCAPoint( Float_t x, Float_t y, Float_t z,
 147                                            Float_t &xp, Float_t &yp, Float_t &zp ) const
 148 {
 149   //* Get the track point closest to the (x,y,z)
 150
 151   Float_t x0 = GetX();
 152   Float_t y0 = GetY();
 153   Float_t k  = GetKappa();
 154   Float_t ex = GetCosPhi();
 155   Float_t ey = GetSinPhi();
 156   Float_t dx = x - x0;
 157   Float_t dy = y - y0;
 158   Float_t ax = dx*k+ey;
 159   Float_t ay = dy*k-ex;
 160   Float_t a = sqrt( ax*ax+ay*ay );
 161   xp = x0 + (dx - ey*( (dx*dx+dy*dy)*k - 2*(-dx*ey+dy*ex) )/(a+1) )/a;
 162   yp = y0 + (dy + ex*( (dx*dx+dy*dy)*k - 2*(-dx*ey+dy*ex) )/(a+1) )/a;
 163   Float_t s = GetS(x,y);
 164   zp = GetZ() + GetDzDs()*s;
 165   if( TMath::Abs(k)>1.e-2 ){
 166     Float_t dZ = TMath::Abs( GetDzDs()*TMath::TwoPi()/k );
 167     if( dZ>.1 ){
 168       zp+= TMath::Nint((z-zp)/dZ)*dZ;
 169     }
 170   }
 171 }
 172
 173 void AliHLTTPCCATrackParam::ConstructXYZ3( const Float_t p0[5], const Float_t p1[5],
 174                                            const Float_t p2[5],
 175                                            Float_t CosPhi0, Float_t t0[] )
 176 {
 177   //* Construct the track in XYZ by 3 points
 178
 179   Float_t px[3]   = { p0[0], p1[0], p2[0] };
 180   Float_t py[3]   = { p0[1], p1[1], p2[1] };
 181   Float_t pz[3]   = { p0[2], p1[2], p2[2] };
 182   Float_t ps2y[3] = { p0[3]*p0[3], p1[3]*p1[3], p2[3]*p2[3] };
 183   Float_t ps2z[3] = { p0[4]*p0[4], p1[4]*p1[4], p2[4]*p2[4] };
 184
 185   Float_t kold = t0 ?t0[4] :0;
 186   ConstructXY3( px, py, ps2y, CosPhi0 );
 187
 188   Float_t pS[3] = { GetS(px[0],py[0]), GetS(px[1],py[1]), GetS(px[2],py[2]) };
 189   Float_t k = Kappa();
 190   if( TMath::Abs(k)>1.e-2 ){
 191     Float_t dS = TMath::Abs( TMath::TwoPi()/k );
 192     pS[1]+= TMath::Nint( (pS[0]-pS[1])/dS )*dS; // not more than half turn
 193     pS[2]+= TMath::Nint( (pS[1]-pS[2])/dS )*dS;
 194     if( t0 ){
 195       Float_t dZ = TMath::Abs(t0[3]*dS);
 196       if( TMath::Abs(dZ)>1. ){
 197         Float_t dsDz = 1./t0[3];
 198         if( kold*k<0 ) dsDz = -dsDz;
 199         Float_t s0 = (pz[0]-t0[1])*dsDz;
 200         Float_t s1 = (pz[1]-t0[1])*dsDz;
 201         Float_t s2 = (pz[2]-t0[1])*dsDz;
 202         pS[0]+= TMath::Nint( (s0-pS[0])/dS )*dS ;
 203         pS[1]+= TMath::Nint( (s1-pS[1])/dS )*dS ;
 204         pS[2]+= TMath::Nint( (s2-pS[2])/dS )*dS ;
 205       }
 206     }
 207   }
 208
 209   Float_t s = pS[0] + pS[1] + pS[2];
 210   Float_t z = pz[0] + pz[1] + pz[2];
 211   Float_t sz = pS[0]*pz[0] + pS[1]*pz[1] + pS[2]*pz[2];
 212   Float_t ss = pS[0]*pS[0] + pS[1]*pS[1] + pS[2]*pS[2];
 213
 214   Float_t a = 3*ss-s*s;
 215   Z() = (z*ss-sz*s)/a; // z0
 216   DzDs() = (3*sz-z*s)/a; // t = dz/ds
 217
 218   Float_t dz0[3] = {ss - pS[0]*s,ss - pS[1]*s,ss - pS[2]*s };
 219   Float_t dt [3] = {3*pS[0] - s, 3*pS[1] - s, 3*pS[2] - s };
 220
 221   fC[2] = (dz0[0]*dz0[0]*ps2z[0] + dz0[1]*dz0[1]*ps2z[1] + dz0[2]*dz0[2]*ps2z[2])/a/a;
 222   fC[7]= (dz0[0]*dt [0]*ps2z[0] + dz0[1]*dt [1]*ps2z[1] + dz0[2]*dt [2]*ps2z[2])/a/a;
 223   fC[9]= (dt [0]*dt [0]*ps2z[0] + dt [1]*dt [1]*ps2z[1] + dt [2]*dt [2]*ps2z[2])/a/a;
 224 }
 225
 226
 227 Bool_t  AliHLTTPCCATrackParam::TransportToX( Float_t x )
 228 {
 229   //* Transport the track parameters to X=x
 230
 231   Bool_t ret = 1;
 232
 233   Float_t x0  = X();
 234   //Float_t y0  = Y();
 235   Float_t k   = Kappa();
 236   Float_t ex = CosPhi();
 237   Float_t ey = SinPhi();
 238   Float_t dx = x - x0;
 239
 240   Float_t ey1 = k*dx + ey;
 241   Float_t ex1;
 242   if( TMath::Abs(ey1)>1 ){ // no intersection -> check the border
 243     ey1 = ( ey1>0 ) ?1 :-1;
 244     ex1 = 0;
 245     dx = ( TMath::Abs(k)>1.e-4) ? ( (ey1-ey)/k ) :0;
 246
 247     Float_t ddx = TMath::Abs(x0+dx - x)*k*k;
 248     Float_t hx[] = {0, -k, 1+ey };
 249     Float_t sx2 = hx[1]*hx[1]*fC[ 3] + hx[2]*hx[2]*fC[ 5];
 250     if( ddx*ddx>3.5*3.5*sx2 ) ret = 0; // x not withing the error
 251     ret = 0; // any case
 252     return ret;
 253   }else{
 254     ex1 = TMath::Sqrt(1 - ey1*ey1);
 255     if( ex<0 ) ex1 = -ex1;
 256   }
 257
 258   Float_t dx2 = dx*dx;
 259   CosPhi() = ex1;
 260   Float_t ss = ey+ey1;
 261   Float_t cc = ex+ex1;
 262   Float_t tg = 0;
 263   if( TMath::Abs(cc)>1.e-4 ) tg = ss/cc; // tan((phi1+phi)/2)
 264   else ret = 0;
 265   Float_t dy = dx*tg;
 266   Float_t dl = dx*TMath::Sqrt(1+tg*tg);
 267
 268   if( cc<0 ) dl = -dl;
 269   Float_t dSin = dl*k/2;
 270   if( dSin > 1 ) dSin = 1;
 271   if( dSin <-1 ) dSin = -1;
 272   Float_t dS = ( TMath::Abs(k)>1.e-4)  ? (2*TMath::ASin(dSin)/k) :dl;
 273   Float_t dz = dS*DzDs();
 274
 275   Float_t cci = 0, exi = 0, ex1i = 0;
 276   if( TMath::Abs(cc)>1.e-4 ) cci = 1./cc;
 277   else ret = 0;
 278   if( TMath::Abs(ex)>1.e-4 ) exi = 1./ex;
 279   else ret = 0;
 280   if( TMath::Abs(ex1)>1.e-4 ) ex1i = 1./ex1;
 281   else ret = 0;
 282
 283   if( !ret ) return ret;
 284
 285   X() += dx;
 286   fP[0]+= dy;
 287   fP[1]+= dz;
 288   fP[2] = ey1;
 289   fP[3] = fP[3];
 290   fP[4] = fP[4];
 291
 292   Float_t h2 = dx*(1+ ex*ex1 + ey*ey1 )*cci*exi*ex1i;
 293   Float_t h4 = dx2*(cc + ss*ey1*ex1i )*cci*cci;
 294
 295   Float_t c00 = fC[0];
 296   Float_t c10 = fC[1];
 297   Float_t c11 = fC[2];
 298   Float_t c20 = fC[3];
 299   Float_t c21 = fC[4];
 300   Float_t c22 = fC[5];
 301   Float_t c30 = fC[6];
 302   Float_t c31 = fC[7];
 303   Float_t c32 = fC[8];
 304   Float_t c33 = fC[9];
 305   Float_t c40 = fC[10];
 306   Float_t c41 = fC[11];
 307   Float_t c42 = fC[12];
 308   Float_t c43 = fC[13];
 309   Float_t c44 = fC[14];
 310
 311   //Float_t H0[5] = { 1,0, h2,  0, h4 };
 312   //Float_t H1[5] = { 0, 1, 0, dS,  0 };
 313   //Float_t H2[5] = { 0, 0, 1,  0, dx };
 314   //Float_t H3[5] = { 0, 0, 0,  1,  0 };
 315   //Float_t H4[5] = { 0, 0, 0,  0,  1 };
 316
 317
 318   fC[0]=( c00  + h2*h2*c22 + h4*h4*c44
 319           + 2*( h2*c20 + h4*c40 + h2*h4*c42 )  );
 320
 321   fC[1]= c10 + h2*c21 + h4*c41 + dS*(c30 + h2*c32 + h4*c43);
 322   fC[2]= c11 + 2*dS*c31 + dS*dS*c33;
 323
 324   fC[3]= c20 + h2*c22 + h4*c42 + dx*( c40 + h2*c42 + h4*c44);
 325   fC[4]= c21 + dS*c32 + dx*(c41 + dS*c43);
 326   fC[5]= c22 +2*dx*c42 + dx2*c44;
 327
 328   fC[6]= c30 + h2*c32 + h4*c43;
 329   fC[7]= c31 + dS*c33;
 330   fC[8]= c32 + dx*c43;
 331   fC[9]= c33;
 332
 333   fC[10]= c40 + h2*c42 + h4*c44;
 334   fC[11]= c41 + dS*c43;
 335   fC[12]= c42 + dx*c44;
 336   fC[13]= c43;
 337   fC[14]= c44;
 338
 339   return ret;
 340 }
 341
 342 Bool_t AliHLTTPCCATrackParam::TransportToXWithMaterial( Float_t x, Float_t Bz )
 343 {
 344   AliHLTTPCCATrackFitParam par;
 345   CalculateFitParameters( par, Bz );
 346   return TransportToXWithMaterial(x, par );
 347 }
 348
 349
 350 Bool_t  AliHLTTPCCATrackParam::TransportToXWithMaterial( Float_t x, AliHLTTPCCATrackFitParam &par )
 351 {
 352   //* Transport the track parameters to X=x
 353
 354   Bool_t ret = 1;
 355
 356   Float_t oldX=GetX();
 357
 358   Float_t x0  = X();
 359   //Float_t y0  = Y();
 360   Float_t k   = Kappa();
 361   Float_t ex = CosPhi();
 362   Float_t ey = SinPhi();
 363   Float_t dx = x - x0;
 364
 365   Float_t ey1 = k*dx + ey;
 366   Float_t ex1;
 367   if( TMath::Abs(ey1)>.99 ){ // no intersection -> check the border
 368     ey1 = ( ey1>0 ) ?1 :-1;
 369     ex1 = 0;
 370     dx = ( TMath::Abs(k)>1.e-4) ? ( (ey1-ey)/k ) :0;
 371
 372     Float_t ddx = TMath::Abs(x0+dx - x)*k*k;
 373     Float_t hx[] = {0, -k, 1+ey };
 374     Float_t sx2 = hx[1]*hx[1]*fC[ 3] + hx[2]*hx[2]*fC[ 5];
 375     if( ddx*ddx>3.5*3.5*sx2 ) ret = 0; // x not withing the error
 376     ret = 0; // any case
 377     return ret;
 378   }else{
 379     ex1 = TMath::Sqrt(1 - ey1*ey1);
 380     if( ex<0 ) ex1 = -ex1;
 381   }
 382
 383   Float_t dx2 = dx*dx;
 384   CosPhi() = ex1;
 385   Float_t ss = ey+ey1;
 386   Float_t cc = ex+ex1;
 387   Float_t tg = 0;
 388   if( TMath::Abs(cc)>1.e-4 ) tg = ss/cc; // tan((phi1+phi)/2)
 389   else ret = 0;
 390   Float_t dy = dx*tg;
 391   Float_t dl = dx*TMath::Sqrt(1+tg*tg);
 392
 393   if( cc<0 ) dl = -dl;
 394   Float_t dSin = dl*k/2;
 395   if( dSin > 1 ) dSin = 1;
 396   if( dSin <-1 ) dSin = -1;
 397   Float_t dS = ( TMath::Abs(k)>1.e-4)  ? (2*TMath::ASin(dSin)/k) :dl;
 398   Float_t dz = dS*DzDs();
 399
 400   Float_t cci = 0, exi = 0, ex1i = 0;
 401   if( TMath::Abs(cc)>1.e-4 ) cci = 1./cc;
 402   else ret = 0;
 403   if( TMath::Abs(ex)>1.e-4 ) exi = 1./ex;
 404   else ret = 0;
 405   if( TMath::Abs(ex1)>1.e-4 ) ex1i = 1./ex1;
 406   else ret = 0;
 407
 408   if( !ret ) return ret;
 409
 410   X() += dx;
 411   fP[0]+= dy;
 412   fP[1]+= dz;
 413   fP[2] = ey1;
 414   fP[3] = fP[3];
 415   fP[4] = fP[4];
 416
 417   Float_t h2 = dx*(1+ ex*ex1 + ey*ey1 )*cci*exi*ex1i;
 418   Float_t h4 = dx2*(cc + ss*ey1*ex1i )*cci*cci;
 419
 420   Float_t c00 = fC[0];
 421   Float_t c10 = fC[1];
 422   Float_t c11 = fC[2];
 423   Float_t c20 = fC[3];
 424   Float_t c21 = fC[4];
 425   Float_t c22 = fC[5];
 426   Float_t c30 = fC[6];
 427   Float_t c31 = fC[7];
 428   Float_t c32 = fC[8];
 429   Float_t c33 = fC[9];
 430   Float_t c40 = fC[10];
 431   Float_t c41 = fC[11];
 432   Float_t c42 = fC[12];
 433   Float_t c43 = fC[13];
 434   Float_t c44 = fC[14];
 435
 436   //Float_t H0[5] = { 1,0, h2,  0, h4 };
 437   //Float_t H1[5] = { 0, 1, 0, dS,  0 };
 438   //Float_t H2[5] = { 0, 0, 1,  0, dx };
 439   //Float_t H3[5] = { 0, 0, 0,  1,  0 };
 440   //Float_t H4[5] = { 0, 0, 0,  0,  1 };
 441
 442
 443   fC[0]=( c00  + h2*h2*c22 + h4*h4*c44
 444           + 2*( h2*c20 + h4*c40 + h2*h4*c42 )  );
 445
 446   fC[1]= c10 + h2*c21 + h4*c41 + dS*(c30 + h2*c32 + h4*c43);
 447   fC[2]= c11 + 2*dS*c31 + dS*dS*c33;
 448
 449   fC[3]= c20 + h2*c22 + h4*c42 + dx*( c40 + h2*c42 + h4*c44);
 450   fC[4]= c21 + dS*c32 + dx*(c41 + dS*c43);
 451   fC[5]= c22 +2*dx*c42 + dx2*c44;
 452
 453   fC[6]= c30 + h2*c32 + h4*c43;
 454   fC[7]= c31 + dS*c33;
 455   fC[8]= c32 + dx*c43;
 456   fC[9]= c33;
 457
 458   fC[10]= c40 + h2*c42 + h4*c44;
 459   fC[11]= c41 + dS*c43;
 460   fC[12]= c42 + dx*c44;
 461   fC[13]= c43;
 462   fC[14]= c44;
 463
 464   Float_t d = TMath::Sqrt(dS*dS + dz*dz );
 465
 466   if (oldX > GetX() ) d = -d;
 467   {
 468     Float_t rho=0.9e-3;
 469     Float_t radLen=28.94;
 470     CorrectForMeanMaterial(d*rho/radLen,d*rho,par);
 471   }
 472
 473   return ret;
 474 }
 475
 476
 477
 478 Float_t AliHLTTPCCATrackParam::ApproximateBetheBloch( Float_t beta2 )
 479 {
 480   //------------------------------------------------------------------
 481   // This is an approximation of the Bethe-Bloch formula with
 482   // the density effect taken into account at beta*gamma > 3.5
 483   // (the approximation is reasonable only for solid materials)
 484   //------------------------------------------------------------------
 485   if (beta2 >= 1) return 0;
 486
 487   if (beta2/(1-beta2)>3.5*3.5)
 488     return 0.153e-3/beta2*( log(3.5*5940)+0.5*log(beta2/(1-beta2)) - beta2);
 489   return 0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2);
 490 }
 491
 492
 493 void AliHLTTPCCATrackParam::CalculateFitParameters( AliHLTTPCCATrackFitParam &par, Float_t Bz, Float_t mass )
 494 {
 495   //*!
 496
 497   const Float_t kCLight = 0.000299792458;
 498   Float_t c = Bz*kCLight;
 499   Float_t p2 = (1.+ fP[3]*fP[3])*c*c;
 500   Float_t k2 = fP[4]*fP[4];
 501   Float_t beta2= p2 / (p2 + mass*mass*k2);
 502   Float_t Bethe = ApproximateBetheBloch(beta2);
 503
 504   Float_t P2 = (k2>1.e-8) ?p2/k2 :10000; // impuls 2
 505   par.fBethe = Bethe;
 506   par.fE = TMath::Sqrt( P2 + mass*mass);
 507   par.fTheta2 = 14.1*14.1/(beta2*p2*1e6)*k2;
 508   par.fEP2 = par.fE/p2*k2;
 509
 510   // Approximate energy loss fluctuation (M.Ivanov)
 511
 512   const Float_t knst=0.07; // To be tuned.
 513   par.fSigmadE2 = knst*par.fEP2*fP[4];
 514   par.fSigmadE2 = par.fSigmadE2 * par.fSigmadE2;
 515
 516   par.fK22 = (1. + fP[3]*fP[3]);
 517   par.fK33 = par.fK22*par.fK22;
 518   par.fK43 = fP[3]*fP[4]*par.fK22;
 519   par.fK44 = fP[3]*fP[3]*fP[4]*fP[4];
 520 }
 521
 522
 523 Bool_t AliHLTTPCCATrackParam::CorrectForMeanMaterial( Float_t xOverX0,  Float_t xTimesRho, AliHLTTPCCATrackFitParam &par )
 524 {
 525   //------------------------------------------------------------------
 526   // This function corrects the track parameters for the crossed material.
 527   // "xOverX0"   - X/X0, the thickness in units of the radiation length.
 528   // "xTimesRho" - is the product length*density (g/cm^2).
 529   //------------------------------------------------------------------
 530
 531   Float_t &fC22=fC[5];
 532   Float_t &fC33=fC[9];
 533   Float_t &fC40=fC[10];
 534   Float_t &fC41=fC[11];
 535   Float_t &fC42=fC[12];
 536   Float_t &fC43=fC[13];
 537   Float_t &fC44=fC[14];
 538
 539   //Energy losses************************
 540
 541   Float_t dE = par.fBethe*xTimesRho;
 542   if ( TMath::Abs(dE) > 0.3*par.fE ) return 0; //30% energy loss is too much!
 543   Float_t corr = (1.- par.fEP2*dE);
 544   if( corr<0.3 ) return 0;
 545   fP[4]*= corr;
 546   fC40*= corr;
 547   fC41*= corr;
 548   fC42*= corr;
 549   fC43*= corr;
 550   fC44*= corr*corr;
 551   fC44+= par.fSigmadE2*TMath::Abs(dE);
 552
 553
 554   //Multiple scattering******************
 555
 556   Float_t theta2 = par.fTheta2*TMath::Abs(xOverX0);
 557   fC22 += theta2*par.fK22*(1.- fP[2]*fP[2]);
 558   fC33 += theta2*par.fK33;
 559   fC43 += theta2*par.fK43;
 560   fC44 += theta2*par.fK44;
 561
 562   return 1;
 563 }
 564
 565
 566
 567 #include "Riostream.h"
 568
 569 Bool_t AliHLTTPCCATrackParam::Rotate( Float_t alpha )
 570 {
 571   //* Rotate the coordinate system in XY on the angle alpha
 572
 573   Float_t cA = TMath::Cos( alpha );
 574   Float_t sA = TMath::Sin( alpha );
 575   Float_t x = X(), y= Y(), sP= SinPhi(), cP= CosPhi();
 576   Float_t cosPhi = cP*cA + sP*sA;
 577   Float_t sinPhi =-cP*sA + sP*cA;
 578
 579   if( TMath::Abs(sinPhi)>.99 || TMath::Abs(cosPhi)<1.e-2 || TMath::Abs(cP)<1.e-2  ) return 0;
 580
 581   Float_t j0 = cP/cosPhi;
 582   Float_t j2 = cosPhi/cP;
 583
 584   X()      =   x*cA +  y*sA;
 585   Y()      =  -x*sA +  y*cA;
 586   CosPhi() =  cosPhi;
 587   SinPhi() =  sinPhi;
 588
 589
 590   //Float_t J[5][5] = { { j0, 0, 0,  0,  0 }, // Y
 591   //                      {  0, 1, 0,  0,  0 }, // Z
 592   //                      {  0, 0, j2, 0,  0 }, // SinPhi
 593   //                      {  0, 0, 0,  1,  0 }, // DzDs
 594   //                      {  0, 0, 0,  0,  1 } }; // Kappa
 595   //cout<<"alpha="<<alpha<<" "<<x<<" "<<y<<" "<<sP<<" "<<cP<<" "<<j0<<" "<<j2<<endl;
 596   //cout<<"      "<<fC[0]<<" "<<fC[1]<<" "<<fC[6]<<" "<<fC[10]<<" "<<fC[4]<<" "<<fC[5]<<" "<<fC[8]<<" "<<fC[12]<<endl;
 597   fC[0]*= j0*j0;
 598   fC[1]*= j0;
 599   //fC[3]*= j0;
 600   fC[6]*= j0;
 601   fC[10]*= j0;
 602
 603   //fC[3]*= j2;
 604   fC[4]*= j2;
 605   fC[5]*= j2*j2;
 606   fC[8]*= j2;
 607   fC[12]*= j2;
 608   //cout<<"      "<<fC[0]<<" "<<fC[1]<<" "<<fC[6]<<" "<<fC[10]<<" "<<fC[4]<<" "<<fC[5]<<" "<<fC[8]<<" "<<fC[12]<<endl;
 609   return 1;
 610 }
 611
 612
 613 Bool_t AliHLTTPCCATrackParam::Filter2( Float_t y, Float_t z, Float_t err2Y, Float_t err2Z )
 614 {
 615   //* Add the y,z measurement with the Kalman filter
 616
 617   Float_t
 618     c00 = fC[ 0],
 619     c10 = fC[ 1], c11 = fC[ 2],
 620     c20 = fC[ 3], c21 = fC[ 4],
 621     c30 = fC[ 6], c31 = fC[ 7],
 622     c40 = fC[10], c41 = fC[11];
 623
 624   Float_t
 625     z0 = y-fP[0],
 626     z1 = z-fP[1];
 627
 628   Float_t v[3] = {err2Y, 0, err2Z};
 629
 630   Float_t mS[3] = { c00+v[0], c10+v[1], c11+v[2] };
 631
 632   Float_t mSi[3];
 633   Float_t det = (mS[0]*mS[2] - mS[1]*mS[1]);
 634
 635   if( det < 1.e-8 ) return 0;
 636   det = 1./det;
 637   mSi[0] = mS[2]*det;
 638   mSi[1] = -mS[1]*det;
 639   mSi[2] = mS[0]*det;
 640
 641   // K = CHtS
 642
 643   Float_t k00, k01 , k10, k11, k20, k21, k30, k31, k40, k41;
 644
 645   k00 = c00*mSi[0] + c10*mSi[1]; k01 = c00*mSi[1] + c10*mSi[2];
 646   k10 = c10*mSi[0] + c11*mSi[1]; k11 = c10*mSi[1] + c11*mSi[2];
 647   k20 = c20*mSi[0] + c21*mSi[1]; k21 = c20*mSi[1] + c21*mSi[2];
 648   k30 = c30*mSi[0] + c31*mSi[1]; k31 = c30*mSi[1] + c31*mSi[2] ;
 649   k40 = c40*mSi[0] + c41*mSi[1]; k41 = c40*mSi[1] + c41*mSi[2] ;
 650
 651   Float_t sinPhi = fP[2] + k20*z0  + k21*z1  ;
 652   if( TMath::Abs(sinPhi)>=0.99 ) return 0;
 653
 654   fNDF  += 2;
 655   fChi2 += ( +(mSi[0]*z0 + mSi[1]*z1 )*z0
 656              +(mSi[1]*z0 + mSi[2]*z1 )*z1 );
 657
 658   fP[ 0]+= k00*z0  + k01*z1 ;
 659   fP[ 1]+= k10*z0  + k11*z1  ;
 660   fP[ 2] = sinPhi;
 661   fP[ 3]+= k30*z0  + k31*z1  ;
 662   fP[ 4]+= k40*z0  + k41*z1  ;
 663
 664
 665   fC[ 0]-= k00*c00 + k01*c10 ;
 666
 667   fC[ 1]-= k10*c00 + k11*c10 ;
 668   fC[ 2]-= k10*c10 + k11*c11 ;
 669
 670   fC[ 3]-= k20*c00 + k21*c10 ;
 671   fC[ 4]-= k20*c10 + k21*c11 ;
 672   fC[ 5]-= k20*c20 + k21*c21 ;
 673
 674   fC[ 6]-= k30*c00 + k31*c10 ;
 675   fC[ 7]-= k30*c10 + k31*c11 ;
 676   fC[ 8]-= k30*c20 + k31*c21 ;
 677   fC[ 9]-= k30*c30 + k31*c31 ;
 678
 679   fC[10]-= k40*c00 + k41*c10 ;
 680   fC[11]-= k40*c10 + k41*c11 ;
 681   fC[12]-= k40*c20 + k41*c21 ;
 682   fC[13]-= k40*c30 + k41*c31 ;
 683   fC[14]-= k40*c40 + k41*c41 ;
 684
 685   if( CosPhi()>=0 ){
 686     CosPhi() = TMath::Sqrt(1-SinPhi()*SinPhi());
 687   }else{
 688     CosPhi() = -TMath::Sqrt(1-SinPhi()*SinPhi());
 689   }
 690   return 1;
 691 }
 692
 693 void AliHLTTPCCATrackParam::FilterY( Float_t y, Float_t erry )
 694 {
 695   //* Add the y measurement with the Kalman filter
 696
 697   Float_t
 698     c00 = fC[ 0],
 699     c10 = fC[ 1],
 700     c20 = fC[ 3],
 701     c30 = fC[ 6],
 702     c40 = fC[10];
 703
 704   Float_t
 705     z0 = y-fP[0];
 706
 707   Float_t s = { c00+erry*erry };
 708   if( TMath::Abs(s)<1.e-4 ) return;
 709
 710   Float_t si = 1/s;
 711
 712   fNDF  += 1;
 713   fChi2 += si*z0*z0;
 714
 715   // K = CHtS
 716
 717   Float_t k0, k1 , k2, k3, k4;
 718
 719   k0 = c00*si;
 720   k1 = c10*si;
 721   k2 = c20*si;
 722   k3 = c30*si;
 723   k4 = c40*si;
 724
 725   Float_t sinPhi = fP[2] + k2*z0 ;
 726   if( TMath::Abs(sinPhi)>=0.99 ) return;
 727
 728   fP[ 0]+= k0*z0 ;
 729   fP[ 1]+= k1*z0 ;
 730   fP[ 2] = sinPhi;
 731   fP[ 3]+= k3*z0 ;
 732   fP[ 4]+= k4*z0 ;
 733
 734   fC[ 0]-= k0*c00;
 735
 736   fC[ 1]-= k1*c00;
 737   fC[ 2]-= k1*c10;
 738
 739   fC[ 3]-= k2*c00;
 740   fC[ 4]-= k2*c10;
 741   fC[ 5]-= k2*c20;
 742
 743   fC[ 6]-= k3*c00;
 744   fC[ 7]-= k3*c10;
 745   fC[ 8]-= k3*c20;
 746   fC[ 9]-= k3*c30;
 747
 748   fC[10]-= k4*c00;
 749   fC[11]-= k4*c10;
 750   fC[12]-= k4*c20;
 751   fC[13]-= k4*c30;
 752   fC[14]-= k4*c40;
 753
 754   if( CosPhi()>=0 ){
 755     CosPhi() = TMath::Sqrt(1-SinPhi()*SinPhi());
 756   }else{
 757     CosPhi() = -TMath::Sqrt(1-SinPhi()*SinPhi());
 758   }
 759
 760 }
 761
 762 void AliHLTTPCCATrackParam::FilterZ( Float_t z, Float_t errz )
 763 {
 764   //* Add the z measurement with the Kalman filter
 765
 766   Float_t
 767     c01 = fC[ 1],
 768     c11 = fC[ 2],
 769     c21 = fC[ 4],
 770     c31 = fC[ 7],
 771     c41 = fC[11];
 772
 773   Float_t
 774     z1 = z-fP[1];
 775
 776   Float_t s = c11 + errz*errz;
 777   if( TMath::Abs(s)<1.e-4 ) return;
 778
 779   Float_t si = 1./s;
 780
 781   fNDF  += 1;
 782   fChi2 += si*z1*z1;
 783
 784   // K = CHtS
 785
 786   Float_t k0, k1 , k2, k3, k4;
 787
 788   k0 = 0;//c01*si;
 789   k1 = c11*si;
 790   k2 = 0;//c21*si;
 791   k3 = c31*si;
 792   k4 = 0;//c41*si;
 793
 794   Float_t sinPhi = fP[2] + k2*z1  ;
 795   if( TMath::Abs(sinPhi)>=0.99 ) return;
 796
 797   fP[ 0]+= k0*z1 ;
 798   fP[ 1]+= k1*z1 ;
 799   fP[ 2] = sinPhi;
 800   fP[ 3]+= k3*z1 ;
 801   fP[ 4]+= k4*z1 ;
 802
 803
 804   fC[ 0]-= k0*c01 ;
 805
 806   fC[ 1]-= k1*c01 ;
 807   fC[ 2]-= k1*c11 ;
 808
 809   fC[ 3]-= k2*c01 ;
 810   fC[ 4]-= k2*c11 ;
 811   fC[ 5]-= k2*c21 ;
 812
 813   fC[ 6]-= k3*c01 ;
 814   fC[ 7]-= k3*c11 ;
 815   fC[ 8]-= k3*c21 ;
 816   fC[ 9]-= k3*c31 ;
 817
 818   fC[10]-= k4*c01 ;
 819   fC[11]-= k4*c11 ;
 820   fC[12]-= k4*c21 ;
 821   fC[13]-= k4*c31 ;
 822   fC[14]-= k4*c41 ;
 823
 824   if( CosPhi()>=0 ){
 825     CosPhi() = TMath::Sqrt(1-SinPhi()*SinPhi());
 826   }else{
 827     CosPhi() = -TMath::Sqrt(1-SinPhi()*SinPhi());
 828   }
 829
 830 }
 831
 832 void AliHLTTPCCATrackParam::Print() const
 833 {
 834   cout<<"track: "<<GetX()<<" "<<GetY()<<" "<<GetZ()<<" "<<GetSinPhi()<<" "<<GetDzDs()<<" "<<GetKappa()<<endl;
 835   cout<<"errs2: "<<GetErr2Y()<<" "<<GetErr2Z()<<" "<<GetErr2SinPhi()<<" "<<GetErr2DzDs()<<" "<<GetErr2Kappa()<<endl;
 836 }