ITS/UPGRADE/ITSUpgradeRec/AliITSUSeed.cxx

   1 #include <TString.h>
   2 #include <TMath.h>
   3 #include "AliITSUSeed.h"
   4 #include "AliLog.h"
   5 #include "AliESDtrack.h"
   6 using namespace TMath;
   7
   8 ClassImp(AliITSUSeed)
   9
  10 //_________________________________________________________________________
  11 AliITSUSeed::AliITSUSeed()
  12 : fHitsPattern(0)
  13   ,fNChildren(0)
  14   ,fClID(0)
  15   ,fChi2Glo(0)
  16   ,fChi2Cl(0)
  17   ,fChi2Penalty(0)
  18   ,fChi2Match(0)
  19   ,fChi2ITSSA(0)
  20   ,fParent(0)
  21 #ifdef _ITSU_TUNING_MODE_ // this is used only for tuning histo filling
  22   ,fOrdBranch(0)
  23   ,fOrdCand(0)
  24 #endif
  25 {
  26   // def c-tor
  27   ResetFMatrix();
  28   for (int i=kNKElem;i--;) fKMatrix[i] = 0.;
  29   for (int i=kNRElem;i--;) fRMatrix[i] = 0.;
  30   fResid[0]=0.;
  31   fResid[1]=0.;
  32   fCovIYZ[0]=0.;
  33   fCovIYZ[1]=0.;
  34   fCovIYZ[2]=0.;
  35 }
  36
  37 //_________________________________________________________________________
  38 AliITSUSeed::~AliITSUSeed()
  39 {
  40   // d-tor
  41 }
  42
  43 //_________________________________________________________________________
  44 AliITSUSeed::AliITSUSeed(const AliITSUSeed& src)
  45   :AliExternalTrackParam(src)
  46   ,fHitsPattern(src.fHitsPattern)
  47   ,fNChildren(src.fNChildren)
  48   ,fClID(src.fClID)
  49   ,fChi2Glo(src.fChi2Glo)
  50   ,fChi2Cl(src.fChi2Cl)
  51   ,fChi2Penalty(src.fChi2Penalty)
  52   ,fChi2Match(src.fChi2Match)
  53   ,fChi2ITSSA(src.fChi2ITSSA)
  54   ,fParent(src.fParent)
  55 #ifdef _ITSU_TUNING_MODE_ // this is used only for tuning histo filling
  56   ,fOrdBranch(src.fOrdBranch)
  57   ,fOrdCand(src.fOrdCand)
  58 #endif
  59 {
  60   // def c-tor
  61   for (int i=kNFElem;i--;) fFMatrix[i] = src.fFMatrix[i];
  62   for (int i=kNKElem;i--;) fKMatrix[i] = src.fKMatrix[i];
  63   for (int i=kNRElem;i--;) fRMatrix[i] = src.fRMatrix[i];
  64   fResid[0]=src.fResid[0];
  65   fResid[1]=src.fResid[1];
  66   fCovIYZ[0]=src.fCovIYZ[0];
  67   fCovIYZ[1]=src.fCovIYZ[1];
  68   fCovIYZ[2]=src.fCovIYZ[2];
  69   //
  70 }
  71
  72 //_________________________________________________________________________
  73 AliITSUSeed &AliITSUSeed::operator=(const AliITSUSeed& src)
  74 {
  75   // def c-tor
  76   if (this == &src) return *this;
  77   this->~AliITSUSeed();
  78   new(this) AliITSUSeed(src);
  79   return *this;
  80 }
  81
  82 //_________________________________________________________________________
  83 void AliITSUSeed::Print(Option_t* opt) const
  84 {
  85   // print seed info
  86   int lr,cl = GetLrCluster(lr);
  87   printf("%cLr%d Nchild: %3d Cl:%4d Chi2Glo:%7.2f(%7.2f) Chi2Cl:%7.2f Penalty: %7.2f Mtc:%6.3f Bwd:%6.3f",IsKilled() ? '-':' ',
  88          lr,GetNChildren(),cl,GetChi2Glo(),GetChi2GloNrm(),GetChi2Cl(), GetChi2Penalty(), GetChi2ITSTPC(), GetChi2ITSSA());
  89   printf(" |");
  90   int lrc=0;
  91   const AliITSUSeed *sdc = this;
  92   while(1) {
  93     if (lrc<lr) printf(".");
  94     else {
  95       sdc = sdc->GetParent(lrc);
  96       if (!sdc) break;
  97       printf("%c",sdc->GetClusterID()<0 ? '.': (sdc->IsFake() ? '-':'+'));
  98     }
  99     lrc++;
 100   }
 101   printf("|\n");
 102   TString opts = opt; opts.ToLower();
 103   if (opts.Contains("etp")) AliExternalTrackParam::Print();
 104   if (opts.Contains("parent") && GetParent()) GetParent()->Print(opt);
 105 }
 106
 107 //______________________________________________________________________________
 108 void AliITSUSeed::InitFromSeed(const AliExternalTrackParam* seed)
 109 {
 110   // init seed from ESD track
 111   TObject::Clear();
 112   AliExternalTrackParam::operator=(*seed);
 113   ResetFMatrix();
 114   fHitsPattern = 0;
 115   fClID = 0;
 116   fNChildren = 0;
 117   fChi2Glo = fChi2Cl = fChi2Penalty = 0;
 118   fParent = 0; //!!!
 119 }
 120
 121
 122 //______________________________________________________________________________
 123 Float_t AliITSUSeed::GetChi2GloNrm() const
 124 {
 125   int ndf = 2*GetNLayersHit() - 5;
 126   return (ndf>0 ? fChi2Glo/ndf : fChi2Glo) + fChi2Penalty;
 127 }
 128
 129
 130 //______________________________________________________________________________
 131 Int_t AliITSUSeed::Compare(const TObject* obj)  const
 132 {
 133   // compare clusters accodring to specific mode
 134   const AliITSUSeed* sd = (const AliITSUSeed*)obj;
 135   const Float_t kTol = 1e-5;
 136   if (!IsKilled() && sd->IsKilled()) return -1;
 137   if ( IsKilled() &&!sd->IsKilled()) return  1;
 138   //
 139   float chi2This  = GetChi2GloNrm();
 140   float chi2Other = sd->GetChi2GloNrm();
 141
 142   if      (chi2This+kTol<chi2Other) return -1;
 143   else if (chi2This-kTol>chi2Other) return  1;
 144   return 0;
 145 }
 146
 147 //______________________________________________________________________________
 148 Bool_t AliITSUSeed::IsEqual(const TObject* obj)  const
 149 {
 150   // compare clusters accodring to specific mode
 151   const AliITSUSeed* sd = (const AliITSUSeed*)obj;
 152   const Float_t kTol = 1e-5;
 153   if (IsKilled() != sd->IsKilled()) return kFALSE;
 154   return Abs(GetChi2GloNrm() - sd->GetChi2GloNrm())<kTol;
 155 }
 156
 157 //______________________________________________________________________________
 158 Bool_t AliITSUSeed::PropagateToX(Double_t xk, Double_t b)
 159 {
 160   // Propagate this track to the plane X=xk (cm) in the field "b" (kG)
 161   Double_t dx=xk-fX;
 162   if (TMath::Abs(dx)<=kAlmost0)  return kTRUE;
 163
 164   Double_t crv=GetC(b);
 165   if (TMath::Abs(b) < kAlmost0Field) crv=0.;
 166   Double_t x2r = crv*dx;
 167   Double_t f1=fP[2], f2=f1 + x2r;
 168   if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
 169   if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
 170   if (TMath::Abs(fP[4])< kAlmost0) return kFALSE;
 171
 172   Double_t &fP0=fP[0], &fP1=fP[1], &fP2=fP[2], &fP3=fP[3], &fP4=fP[4];
 173   Double_t
 174   &fC00=fC[0],
 175   &fC10=fC[1],   &fC11=fC[2],
 176   &fC20=fC[3],   &fC21=fC[4],   &fC22=fC[5],
 177   &fC30=fC[6],   &fC31=fC[7],   &fC32=fC[8],   &fC33=fC[9],
 178   &fC40=fC[10],  &fC41=fC[11],  &fC42=fC[12],  &fC43=fC[13], &fC44=fC[14];
 179
 180   Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
 181   if (TMath::Abs(r1)<kAlmost0)  return kFALSE;
 182   if (TMath::Abs(r2)<kAlmost0)  return kFALSE;
 183
 184   fX=xk;
 185   double dy2dx = (f1+f2)/(r1+r2);
 186   fP0 += dx*dy2dx;
 187   if (TMath::Abs(x2r)<0.05) {
 188     fP1 += dx*(r2 + f2*dy2dx)*fP3;  // Many thanks to P.Hristov !
 189     fP2 += x2r;
 190   }
 191   else {
 192     // for small dx/R the linear apporximation of the arc by the segment is OK,
 193     // but at large dx/R the error is very large and leads to incorrect Z propagation
 194     // angle traversed delta = 2*asin(dist_start_end / R / 2), hence the arc is: R*deltaPhi
 195     // The dist_start_end is obtained from sqrt(dx^2+dy^2) = x/(r1+r2)*sqrt(2+f1*f2+r1*r2)
 196     // Similarly, the rotation angle in linear in dx only for dx<<R
 197     double chord = dx*TMath::Sqrt(1+dy2dx*dy2dx);   // distance from old position to new one
 198     double rot = 2*TMath::ASin(0.5*chord*crv); // angular difference seen from the circle center
 199     fP1 += rot/crv*fP3;
 200     fP2  = TMath::Sin(rot + TMath::ASin(fP2));
 201   }
 202
 203   //f = F - 1
 204   double r1i = 1./r1;
 205   double r2i = 1./r2;
 206   double tg1 = f1*r1i;
 207   double tg2 = f2*r2i;
 208   double v0 = 1. + dy2dx*tg2;
 209   double v1 = (r1i+r2i)*(dy2dx*(tg1+tg2)+2);
 210   double v2 = (r1i+r2i)*v0;
 211   //
 212   double f24 = dx*crv/fP4;
 213   double f02 = dx*v1;
 214   double f04 = dx*v2*f24;
 215   double f12 = dx*fP3*    (f2*v1+dy2dx-tg2);
 216   double f13 = dx*r2*v0;
 217   double f14 = dx*f24*fP3*(f2*v2+dy2dx-tg2);
 218   //
 219   //b = C*ft
 220   Double_t b00=f02*fC20 + f04*fC40, b01=f12*fC20 + f14*fC40 + f13*fC30;
 221   Double_t b02=f24*fC40;
 222   Double_t b10=f02*fC21 + f04*fC41, b11=f12*fC21 + f14*fC41 + f13*fC31;
 223   Double_t b12=f24*fC41;
 224   Double_t b20=f02*fC22 + f04*fC42, b21=f12*fC22 + f14*fC42 + f13*fC32;
 225   Double_t b22=f24*fC42;
 226   Double_t b40=f02*fC42 + f04*fC44, b41=f12*fC42 + f14*fC44 + f13*fC43;
 227   Double_t b42=f24*fC44;
 228   Double_t b30=f02*fC32 + f04*fC43, b31=f12*fC32 + f14*fC43 + f13*fC33;
 229   Double_t b32=f24*fC43;
 230
 231   //a = f*b = f*C*ft
 232   Double_t a00=f02*b20+f04*b40,a01=f02*b21+f04*b41,a02=f02*b22+f04*b42;
 233   Double_t a11=f12*b21+f14*b41+f13*b31,a12=f12*b22+f14*b42+f13*b32;
 234   Double_t a22=f24*b42;
 235
 236   //F*C*Ft = C + (b + bt + a)
 237   fC00 += b00 + b00 + a00;
 238   fC10 += b10 + b01 + a01;
 239   fC20 += b20 + b02 + a02;
 240   fC30 += b30;
 241   fC40 += b40;
 242   fC11 += b11 + b11 + a11;
 243   fC21 += b21 + b12 + a12;
 244   fC31 += b31;
 245   fC41 += b41;
 246   fC22 += b22 + b22 + a22;
 247   fC32 += b32;
 248   fC42 += b42;
 249   //
 250   // update stored transformation matrix   F = Fnew*Fold
 251   fFMatrix[kF04] += f04 + f24*fFMatrix[kF02];
 252   fFMatrix[kF14] += f14 + f24*fFMatrix[kF12];
 253   fFMatrix[kF02] += f02;
 254   fFMatrix[kF12] += f12;
 255   fFMatrix[kF13] += f13;
 256   fFMatrix[kF24] += f24;
 257   //
 258   CheckCovariance();
 259
 260   return kTRUE;
 261 }
 262
 263 //__________________________________________________________________
 264 Int_t AliITSUSeed::GetClusterIndex(Int_t ind) const
 265 {
 266   // get ind-th cluster index
 267   int ncl = 0;
 268   const AliITSUSeed* seed = this;
 269   while(seed) {
 270     if ( seed->HasCluster() && (ncl++==ind) ) return seed->GetLrClusterID();//GetClusterID();
 271     seed = (AliITSUSeed*)seed->GetParent();
 272   }
 273   return -1;
 274   //
 275 }
 276
 277 //______________________________________________________________________________
 278 Bool_t AliITSUSeed::RotateToAlpha(Double_t alpha)
 279 {
 280   // Transform this track to the local coord. system rotated
 281   // by angle "alpha" (rad) with respect to the global coord. system.
 282   //
 283   if (TMath::Abs(fP[2]) >= kAlmost1) {
 284      AliError(Form("Precondition is not satisfied: |sin(phi)|>1 ! %f",fP[2]));
 285      return kFALSE;
 286   }
 287   //
 288   if      (alpha < -TMath::Pi()) alpha += 2*TMath::Pi();
 289   else if (alpha >= TMath::Pi()) alpha -= 2*TMath::Pi();
 290   //
 291   Double_t &fP0=fP[0];
 292   Double_t &fP2=fP[2];
 293   Double_t &fC00=fC[0];
 294   Double_t &fC10=fC[1];
 295   Double_t &fC20=fC[3];
 296   Double_t &fC21=fC[4];
 297   Double_t &fC22=fC[5];
 298   Double_t &fC30=fC[6];
 299   Double_t &fC32=fC[8];
 300   Double_t &fC40=fC[10];
 301   Double_t &fC42=fC[12];
 302   //
 303   Double_t x=fX;
 304   Double_t ca=TMath::Cos(alpha-fAlpha), sa=TMath::Sin(alpha-fAlpha);
 305   Double_t sf=fP2, cf=TMath::Sqrt((1.- fP2)*(1.+fP2)); // Improve precision
 306   // RS: check if rotation does no invalidate track model (cos(local_phi)>=0, i.e. particle
 307   // direction in local frame is along the X axis
 308   if ((cf*ca+sf*sa)<0) {
 309     AliDebug(1,Form("Rotation failed: local cos(phi) would become %.2f",cf*ca+sf*sa));
 310     return kFALSE;
 311   }
 312   //
 313   Double_t tmp=sf*ca - cf*sa;
 314
 315   if (TMath::Abs(tmp) >= kAlmost1) {
 316      if (TMath::Abs(tmp) > 1.+ Double_t(FLT_EPSILON))
 317         AliWarning(Form("Rotation failed ! %.10e",tmp));
 318      return kFALSE;
 319   }
 320   fAlpha = alpha;
 321   fX =  x*ca + fP0*sa;
 322   fP0= -x*sa + fP0*ca;
 323   fP2=  tmp;
 324
 325   if (TMath::Abs(cf)<kAlmost0) {
 326     AliError(Form("Too small cosine value %f",cf));
 327     cf = kAlmost0;
 328   }
 329
 330   Double_t rr=(ca+sf/cf*sa);
 331
 332   fC00 *= (ca*ca);
 333   fC10 *= ca;
 334   fC20 *= ca*rr;
 335   fC21 *= rr;
 336   fC22 *= rr*rr;
 337   fC30 *= ca;
 338   fC32 *= rr;
 339   fC40 *= ca;
 340   fC42 *= rr;
 341   //
 342   fRMatrix[kR00] = ca;
 343   fRMatrix[kR22] = rr;
 344   //
 345   CheckCovariance();
 346
 347   return kTRUE;
 348 }
 349
 350 //______________________________________________________________________________
 351 Bool_t AliITSUSeed::GetTrackingXAtXAlpha(double xOther, double alpOther, double bz, double &xdst)
 352 {
 353   // calculate X and Y in the tracking frame of the track, corresponding to other X,Alpha tracking
 354   double ca=TMath::Cos(alpOther-fAlpha), sa=TMath::Sin(alpOther-fAlpha);
 355   double &y=fP[0], &sf=fP[2], cf=Sqrt((1.-sf)*(1.+sf));
 356   double eta = xOther - fX*ca - y*sa;
 357   double xi  = sf*ca - cf*sa;
 358   if (Abs(xi)>= kAlmost1) return kFALSE;
 359   double nu  = xi + GetC(bz)*eta;
 360   if (Abs(nu)>= kAlmost1) return kFALSE;
 361   xdst = xOther*ca - sa*( y*ca-fX*sa + eta*(xi+nu)/(Sqrt((1.-xi)*(1.+xi)) + Sqrt((1.-nu)*(1.+nu))) );
 362   return kTRUE;
 363 }
 364
 365 //____________________________________________________________________
 366 Double_t AliITSUSeed::GetPredictedChi2(Double_t p[2],Double_t cov[3])
 367 {
 368   // Estimate the chi2 of the space point "p" with the cov. matrix "cov"
 369   // Store info needed for update and smoothing
 370   Double_t sdd = fC[0] + cov[0];
 371   Double_t sdz = fC[1] + cov[1];
 372   Double_t szz = fC[2] + cov[2];
 373   Double_t det = sdd*szz - sdz*sdz;
 374   if (TMath::Abs(det) < kAlmost0) return kVeryBig;
 375   det = 1./det;
 376   fCovIYZ[0] =  szz*det;
 377   fCovIYZ[1] = -sdz*det;
 378   fCovIYZ[2] =  sdd*det;
 379   double &dy = fResid[0] = p[0] - fP[0];
 380   double &dz = fResid[1] = p[1] - fP[1];
 381   //
 382   return dy*(dy*fCovIYZ[0]+dz*fCovIYZ[1]) + dz*(dy*fCovIYZ[1]+dz*(fCovIYZ[2]));
 383   //
 384 }
 385
 386 //____________________________________________________________________
 387 Bool_t AliITSUSeed::Update()
 388 {
 389   // Update the track parameters with the measurement stored during GetPredictedChi2
 390   //
 391   Double_t &fP0=fP[0], &fP1=fP[1], &fP2=fP[2], &fP3=fP[3], &fP4=fP[4],
 392     &fC00=fC[kS00],
 393     &fC10=fC[kS10],  &fC11=fC[kS11],
 394     &fC20=fC[kS20],  &fC21=fC[kS21],  &fC22=fC[kS22],
 395     &fC30=fC[kS30],  &fC31=fC[kS31],  &fC32=fC[kS32],  &fC33=fC[kS33],
 396     &fC40=fC[kS40],  &fC41=fC[kS41],  &fC42=fC[kS42],  &fC43=fC[kS43], &fC44=fC[kS44];
 397   //
 398   double &r00=fCovIYZ[0],&r01=fCovIYZ[1],&r11=fCovIYZ[2];
 399   double &dy=fResid[0], &dz=fResid[1];
 400   //
 401   // store info needed for smoothing in the fKMatrix
 402   double &k00 = fKMatrix[kK00] = fC00*r00+fC10*r01;
 403   double &k01 = fKMatrix[kK01] = fC00*r01+fC10*r11;
 404   double &k10 = fKMatrix[kK10] = fC10*r00+fC11*r01;
 405   double &k11 = fKMatrix[kK11] = fC10*r01+fC11*r11;
 406   double &k20 = fKMatrix[kK20] = fC20*r00+fC21*r01;
 407   double &k21 = fKMatrix[kK21] = fC20*r01+fC21*r11;
 408   double &k30 = fKMatrix[kK30] = fC30*r00+fC31*r01;
 409   double &k31 = fKMatrix[kK31] = fC30*r01+fC31*r11;
 410   double &k40 = fKMatrix[kK40] = fC40*r00+fC41*r01;
 411   double &k41 = fKMatrix[kK41] = fC40*r01+fC41*r11;
 412   //
 413   Double_t sf=fP2 + k20*dy + k21*dz;
 414   if (TMath::Abs(sf) > kAlmost1) return kFALSE;
 415
 416   fP0 += k00*dy + k01*dz;
 417   fP1 += k10*dy + k11*dz;
 418   fP2  = sf;
 419   fP3 += k30*dy + k31*dz;
 420   fP4 += k40*dy + k41*dz;
 421   //
 422   Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
 423   Double_t c12=fC21, c13=fC31, c14=fC41;
 424
 425   fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
 426   fC20-=k00*c02+k01*c12;   fC30-=k00*c03+k01*c13;
 427   fC40-=k00*c04+k01*c14;
 428
 429   fC11-=k10*c01+k11*fC11;
 430   fC21-=k10*c02+k11*c12;   fC31-=k10*c03+k11*c13;
 431   fC41-=k10*c04+k11*c14;
 432
 433   fC22-=k20*c02+k21*c12;   fC32-=k20*c03+k21*c13;
 434   fC42-=k20*c04+k21*c14;
 435
 436   fC33-=k30*c03+k31*c13;
 437   fC43-=k30*c04+k31*c14;
 438
 439   fC44-=k40*c04+k41*c14;
 440   //
 441   CheckCovariance();
 442   //
 443   return kTRUE;
 444 }
 445
 446
 447 //____________________________________________________________________
 448 Bool_t AliITSUSeed::Smooth(Double_t vecL[5],Double_t matL[15])
 449 {
 450   // Prepare MBF smoothing auxiliary params for smoothing at prev. point:
 451   // \hat{l_N} = 0
 452   // \hat{L_N} = 0
 453   // \tilde{l_j} = -H_j^T N_{j}^{-1} z_j + B_{j}^T \hat{l_j}
 454   // \tilde{L_j} =  H_j^T N_{j}^{-1} H_j + B_j^T \hat{L_j} B_j
 455   // \hat{l_j} = F_j^T \tilde{l_{j+1}}
 456   // \hat{L_j} = F_j^T \tilde{L_{j+1}} F_j
 457   //
 458   // P_{j/N} = P_{j/j} - P_{j/j} \hat{L_j} P_{j/j}
 459   // \hat{x_{j/N}} = \hat{x_{j/j}} - P_{j/j} \hat{l_j}
 460   //
 461   // N^-1 = fCovIYZ
 462   // z = fResid
 463   // B = I - K H
 464   // H = {{1,0,0,0,0},{0,1,0,0,0}}
 465   //
 466   // calc. \tilde{l_j}
 467   //
 468   if (GetClusterID()<0) return kTRUE;
 469   //
 470
 471   double
 472     &k00=fKMatrix[kK00],&k01=fKMatrix[kK01],
 473     &k10=fKMatrix[kK10],&k11=fKMatrix[kK11],
 474     &k20=fKMatrix[kK20],&k21=fKMatrix[kK21],
 475     &k30=fKMatrix[kK30],&k31=fKMatrix[kK31],
 476     &k40=fKMatrix[kK40],&k41=fKMatrix[kK41];
 477   double
 478     &l00=matL[kS00],
 479     &l10=matL[kS10], &l11=matL[kS11],
 480     &l20=matL[kS20], &l21=matL[kS21], &l22=matL[kS22],
 481     &l30=matL[kS30], &l31=matL[kS31], &l32=matL[kS32], &l33=matL[kS33],
 482     &l40=matL[kS40], &l41=matL[kS41], &l42=matL[kS42], &l43=matL[kS43], &l44=matL[kS44];
 483   //
 484   // calculate correction
 485   double corrVec[5]={0},corrMat[15]={0};
 486   corrVec[0] = fC[kS00]*vecL[0] + fC[kS10]*vecL[1] + fC[kS20]*vecL[2] + fC[kS30]*vecL[3] + fC[kS40]*vecL[4];
 487   corrVec[1] = fC[kS10]*vecL[0] + fC[kS11]*vecL[1] + fC[kS21]*vecL[2] + fC[kS31]*vecL[3] + fC[kS41]*vecL[4];
 488   corrVec[2] = fC[kS20]*vecL[0] + fC[kS21]*vecL[1] + fC[kS22]*vecL[2] + fC[kS32]*vecL[3] + fC[kS42]*vecL[4];
 489   corrVec[3] = fC[kS30]*vecL[0] + fC[kS31]*vecL[1] + fC[kS32]*vecL[2] + fC[kS33]*vecL[3] + fC[kS43]*vecL[4];
 490   corrVec[4] = fC[kS40]*vecL[0] + fC[kS41]*vecL[1] + fC[kS42]*vecL[2] + fC[kS43]*vecL[3] + fC[kS44]*vecL[4];
 491   //
 492   double *crm = ProdABA(fC,matL);
 493   for (int i=0;i<15;i++) corrMat[i] = crm[i];
 494
 495   double vcL0 = vecL[0], vcL1 = vecL[1];
 496   vecL[0] -= k00*vcL0+k10*vcL1+k20*vecL[2]+k30*vecL[3]+k40*vecL[4] + fCovIYZ[0]*fResid[0] + fCovIYZ[1]*fResid[1];
 497   vecL[1] -= k01*vcL0+k11*vcL1+k21*vecL[2]+k31*vecL[3]+k41*vecL[4] + fCovIYZ[1]*fResid[0] + fCovIYZ[2]*fResid[1];
 498
 499   /*
 500   double vcL0 = vecL[0], vcL1 = vecL[1];
 501   vecL[0] -= k00*vcL0+k10*vcL1+fKMatrix[kK20]*vecL[2]+k30*vecL[3]+k40*vecL[4] + fCovIYZ[0]*fResid[0] + fCovIYZ[1]*fResid[1];
 502   vecL[1] -= k01*vcL0+fKMatrix[kK11]*vcL1+k21*vecL[2]+k31*vecL[3]+k41*vecL[4] + fCovIYZ[1]*fResid[0] + fCovIYZ[2]*fResid[1];
 503   vecL[3] += fFMatrix[kF13]*vecL[1];
 504   vecL[4]  = fFMatrix[kF04]*vecL[0] + fFMatrix[kF14]*vecL[1] + fFMatrix[kF24]*vecL[2] + fFMatrix[kF44]*vecL[4];
 505   vecL[2] += fFMatrix[kF02]*vecL[0] + fFMatrix[kF12]*vecL[1];
 506   //
 507   */
 508   // and \hat{l_j} in one go
 509
 510   // L = H^T * sg * H + (I-KH)^T * L * (I - KH)
 511   double v00 =  k00*l00+k10*l10+k20*l20+k30*l30+k40*l40;
 512   double v10 =  k00*l10+k10*l11+k20*l21+k30*l31+k40*l41;
 513   double v20 =  k00*l20+k10*l21+k20*l22+k30*l32+k40*l42;
 514   double v30 =  k00*l30+k10*l31+k20*l32+k30*l33+k40*l43;
 515   double v40 =  k00*l40+k10*l41+k20*l42+k30*l43+k40*l44;
 516   //
 517   double v01 =  k01*l00+k11*l10+k21*l20+k31*l30+k41*l40;
 518   double v11 =  k01*l10+k11*l11+k21*l21+k31*l31+k41*l41;
 519   double v21 =  k01*l20+k11*l21+k21*l22+k31*l32+k41*l42;
 520   double v31 =  k01*l30+k11*l31+k21*l32+k31*l33+k41*l43;
 521   double v41 =  k01*l40+k11*l41+k21*l42+k31*l43+k41*l44;
 522   //
 523   // (H^T * K^T * L * K * H) - (L * K * H) - (H^T * K^T * L) + (H^T*N^-1*H)
 524   l00 += k00*v00 + k10*v10 + k20*v20 + k30*v30 + k40*v40 - v00 - v00 + fCovIYZ[0];
 525   l10 += k01*v00 + k11*v10 + k21*v20 + k31*v30 + k41*v40 - v01 - v10 + fCovIYZ[1];
 526   l11 += k01*v01 + k11*v11 + k21*v21 + k31*v31 + k41*v41 - v11 - v11 + fCovIYZ[2];
 527   //
 528   l20 -= v20;
 529   l21 -= v21;
 530   l30 -= v30;
 531   l31 -= v31;
 532   l40 -= v40;
 533   l41 -= v41;
 534   //
 535   printf("CorrMt:\n");
 536   printf("%+e\n%+e %+e\n%+e %+e %+e\n%+e %+e %+e %+e\n%+e %+e %+e %+e %+e\n",
 537          corrMat[kS00],corrMat[kS10],corrMat[kS11],corrMat[kS20],corrMat[kS21],corrMat[kS22],
 538          corrMat[kS30],corrMat[kS31],corrMat[kS32],corrMat[kS33],
 539          corrMat[kS40],corrMat[kS41],corrMat[kS42],corrMat[kS43],corrMat[kS44]);
 540
 541   printf("SMcorr: %+e %+e %+e %+e %+e\n",corrVec[0],corrVec[1],corrVec[2],corrVec[3],corrVec[4]);
 542
 543   printf("State : "); this->AliExternalTrackParam::Print("");
 544   //
 545   printf("\nBefore transport back (RotElems: %+e %+e)\n",fRMatrix[kR00],fRMatrix[kR22]);
 546   printf("Res: %+e %+e | Err: %+e %+e %+e\n",fResid[0],fResid[1],fCovIYZ[0],fCovIYZ[1],fCovIYZ[2]);
 547   printf("Lr%d VecL: ",GetLayerID()); for (int i=0;i<5;i++) printf("%+e ",vecL[i]); printf("\n");
 548   //
 549   printf("%+e\n%+e %+e\n%+e %+e %+e\n%+e %+e %+e %+e\n%+e %+e %+e %+e %+e\n",
 550          matL[kS00],matL[kS10],matL[kS11],matL[kS20],matL[kS21],matL[kS22],
 551          matL[kS30],matL[kS31],matL[kS32],matL[kS33],matL[kS40],matL[kS41],matL[kS42],matL[kS43],matL[kS44]);
 552   //
 553   printf("F: "); for (int i=0;i<kNFElem;i++) printf("%+e ",fFMatrix[i]); printf("\n");
 554   printf("K: "); for (int i=0;i<kNKElem;i++) printf("%+e ",fKMatrix[i]); printf("\n");
 555   //
 556   // apply rotation matrix (diagonal)
 557   vecL[0] *= fRMatrix[kR00];
 558   vecL[2] *= fRMatrix[kR22];
 559   //
 560   l00 *= fRMatrix[kR00]*fRMatrix[kR00];
 561   l10 *= fRMatrix[kR00];
 562   l20 *= fRMatrix[kR22]*fRMatrix[kR00];
 563   l21 *= fRMatrix[kR22];
 564   l22 *= fRMatrix[kR22]*fRMatrix[kR22];
 565   l30 *= fRMatrix[kR00];
 566   l32 *= fRMatrix[kR22];
 567   l40 *= fRMatrix[kR00];
 568   l42 *= fRMatrix[kR22];
 569   //
 570   // Apply translation matrix F^T. Note, that fFMatrix keeps non-trivial elems of F-1 = f, except the e-loss coeff f44
 571   // We need F^T*L* F = L + (L*f) + (L*f)^T + f^T * (L*f)
 572   //
 573   double
 574     &f02=fFMatrix[kF02],&f04=fFMatrix[kF04],
 575     &f12=fFMatrix[kF12],&f13=fFMatrix[kF13],&f14=fFMatrix[kF14],
 576     &f24=fFMatrix[kF24],
 577     f44 =fFMatrix[kF44];
 578   //
 579   vecL[4]  = f04*vecL[0]+f14*vecL[1]+f24*vecL[2]+f44*vecL[4];
 580   vecL[3] += f13*vecL[1];
 581   vecL[2] += f02*vecL[0]+f12*vecL[1];
 582   //
 583   f44 -= 1.0; // !!!!!
 584   //
 585   //b = L*f
 586   Double_t b02=l00*f02+l10*f12, b03=l10*f13, b04=l00*f04+l10*f14+l20*f24+l40*f44;
 587   Double_t b12=l10*f02+l11*f12, b13=l11*f13, b14=l10*f04+l11*f14+l21*f24+l41*f44;
 588   Double_t b22=l20*f02+l21*f12, b23=l21*f13, b24=l20*f04+l21*f14+l22*f24+l42*f44;
 589   Double_t b32=l30*f02+l31*f12, b33=l31*f13, b34=l30*f04+l31*f14+l32*f24+l43*f44;
 590   Double_t b42=l40*f02+l41*f12, b43=l41*f13, b44=l40*f04+l41*f14+l42*f24+l44*f44;
 591   //
 592   //a = f^T * b = f^T * L * f, profit from symmetry
 593   Double_t a22=f02*b02+f12*b12, a33=f13*b13, a44=f04*b04+f14*b14+f24*b24+f44*b44,
 594     a32=f13*b12, //= a23=f02*b03+f12*b13,
 595     a42=f02*b04+f12*b14, //f04*b02+f14*b12+f24*b22+f44*b42 = a24
 596     a43=f13*b14;         //f04*b03+f14*b13+f24*b23+f44*b43 = a34
 597   //
 598   // F^T*L* F = L + (b + b^T + a)
 599   l44 += b44 + b44 + a44;
 600   l43 += b43 + b34 + a43;
 601   l42 += b42 + b24 + a42;
 602   l41 += b14;
 603   l40 += b04;
 604   l33 += b33 + b33 + a33;
 605   l32 += b32 + b23 + a32;
 606   l31 += b13;
 607   l30 += b03;
 608   l22 += b22 + b23 + a22;
 609   l21 += b12;
 610   l20 += b02;
 611   //
 612   printf("After transport back\n");
 613   printf("Lr%d VecL: ",GetLayerID()); for (int i=0;i<5;i++) printf("%+e ",vecL[i]); printf("\n");
 614   //
 615   printf("%+e\n%+e %+e\n%+e %+e %+e\n%+e %+e %+e %+e\n%+e %+e %+e %+e %+e\n",
 616          matL[kS00],matL[kS10],matL[kS11],matL[kS20],matL[kS21],matL[kS22],
 617          matL[kS30],matL[kS31],matL[kS32],matL[kS33],matL[kS40],matL[kS41],matL[kS42],matL[kS43],matL[kS44]);
 618
 619   return kTRUE;
 620 }
 621
 622 //____________________________________________________________________
 623 Double_t* AliITSUSeed::ProdABA(const double a[15],const double b[15]) const
 624 {
 625   // product of symmetric matrices A*B*A
 626   //
 627   const Short_t knd[5][5] = {
 628     {kS00,kS10,kS20,kS30,kS40},
 629     {kS10,kS11,kS21,kS31,kS41},
 630     {kS20,kS21,kS22,kS32,kS42},
 631     {kS30,kS31,kS32,kS33,kS43},
 632     {kS40,kS41,kS42,kS43,kS44}
 633   };
 634   //
 635   static double aba[15];
 636   // 1) ba = B*A
 637   double ba[5][5];
 638   for (int i=5;i--;) for (int j=5;j--;) {
 639       ba[i][j] = 0;
 640       for (int k=5;k--;) ba[i][j] += b[knd[i][k]]*a[knd[k][j]];
 641     }
 642   //
 643   // 2) A * ba, lower triangle only
 644   for (int i=5;i--;) for (int j=i+1;j--;) {
 645       aba[knd[i][j]] = 0;
 646       for (int k=5;k--;) aba[knd[i][j]] += a[knd[i][k]]*ba[k][j];
 647     }
 648   //
 649   return &aba[0];
 650 }
 651
 652 //____________________________________________________________________
 653 Bool_t AliITSUSeed::ContainsFake() const
 654 {
 655   // check if the full branch containes a fake cluster
 656   const AliITSUSeed* seed = this;
 657   while(seed) {
 658     if ( seed->IsFake() ) return kTRUE;
 659     seed = (AliITSUSeed*)seed->GetParent();
 660   }
 661   return kFALSE;
 662 }
 663
 664 //__________________________________________________________________
 665 Int_t AliITSUSeed::FetchClusterInfo(Int_t *clIDarr) const
 666 {
 667   // fill cl.id's in the array. The clusters of layer L will be set at slots
 668   // clID[2L] (and clID[2L+1] if there is an extra cluster).
 669   Int_t lr,ncl=0;
 670   const AliITSUSeed* seed = this;
 671   do {
 672     int clID = seed->GetLrCluster(lr);
 673     if (clID>=0) {
 674       lr<<=1;
 675       clIDarr[ clIDarr[lr]<0 ? lr : lr+1 ] = clID;
 676       ncl++;
 677     }
 678   } while ((seed=(AliITSUSeed*)seed->GetParent()));
 679   return ncl;
 680 }
 681
 682 /*
 683 //____________________________________________________________________
 684 Bool_t AliITSUSeed::Smooth(Double_t vecL[5],Double_t matL[15])
 685 {
 686   // Prepare MBF smoothing auxiliary params for smoothing at prev. point:
 687   // \hat{l_N} = 0
 688   // \hat{L_N} = 0
 689   // \tilde{l_j} = -H_j^T N_{j}^{-1} z_j + B_{j}^T \hat{l_j}
 690   // \tilde{L_j} =  H_j^T N_{j}^{-1} H_j + B_j^T \hat{L_j} B_j
 691   // \hat{l_j} = F_j^T \tilde{l_{j+1}}
 692   // \hat{L_j} = F_j^T \tilde{L_{j+1}} F_j
 693   //
 694   // P_{j/N} = P_{j/j} - P_{j/j} \hat{L_j} P_{j/j}
 695   // \hat{x_{j/N}} = \hat{x_{j/j}} - P_{j/j} \hat{l_j}
 696   //
 697   // N^-1 = fCovIYZ
 698   // z = fResid
 699   // B = I - K H
 700   // H = {{1,0,0,0,0},{0,1,0,0,0}}
 701   //
 702   // calc. \tilde{l_j} and \hat{l_j} in one go
 703   //
 704   if (GetClusterID()<0) return kTRUE;
 705   //
 706   double
 707     &k00=fKMatrix[kK00],&k01=fKMatrix[kK01],
 708     &k10=fKMatrix[kK10],&k11=fKMatrix[kK11],
 709     &k20=fKMatrix[kK20],&k21=fKMatrix[kK21],
 710     &k30=fKMatrix[kK30],&k31=fKMatrix[kK31],
 711     &k40=fKMatrix[kK40],&k41=fKMatrix[kK41];
 712   double
 713     &matL00=matL[kS00],
 714     &matL10=matL[kS01],  &matL11=matL[kS11],
 715     &matL20=matL[kS20],  &matL21=matL[kS21],  &matL22=matL[kS22],
 716     &matL30=matL[kS30],  &matL31=matL[kS31],  &matL32=matL[kS32],  &matL33=matL[kS33],
 717     &matL40=matL[kS40],  &matL41=matL[kS41],  &matL42=matL[kS42],  &matL43=matL[kS43], &matL44=matL[kS44];
 718   //
 719   double vcL0 = vecL[0], vcL1 = vecL[1];
 720   vecL[0] -= k00*vcL0+k10*vcL1+fKMatrix[kK20]*vecL[2]+k30*vecL[3]+k40*vecL[4] + fCovIYZ[0]*fResid[0] + fCovIYZ[1]*fResid[1];
 721   vecL[1] -= k01*vcL0+fKMatrix[kK11]*vcL1+k21*vecL[2]+k31*vecL[3]+k41*vecL[4] + fCovIYZ[1]*fResid[0] + fCovIYZ[2]*fResid[1];
 722   vecL[3] += fFMatrix[kF13]*vecL[1];
 723   vecL[4] =  fFMatrix[kF04]*vecL[0] + fFMatrix[kF14]*vecL[1] + fFMatrix[kF24]*vecL[2] + fFMatrix[kF44]*vecL[4];
 724   vecL[2] += fFMatrix[kF02]*vecL[0] + fFMatrix[kF12]*vecL[1];
 725   //
 726
 727   // L = H^T * sg * H + (I-KH)^T * L * (I - KH)
 728   double v00 =  k00*matL00+k10*matL10+k20*matL20+k30*matL30+k40*matL40;
 729   double v10 =  k00*matL10+k10*matL11+k20*matL21+k30*matL31+k40*matL41;
 730   double v20 =  k00*matL20+k10*matL12+k20*matL22+k30*matL32+k40*matL42;
 731   double v30 =  k00*matL30+k10*matL13+k20*matL23+k30*matL33+k40*matL43;
 732   double v40 =  k00*matL40+k10*matL14+k20*matL24+k30*matL34+k40*matL44;
 733   //
 734   double v01 =  k01*matL00+k11*matL10+k21*matL20+k31*matL30+k41*matL40;
 735   double v11 =  k01*matL01+k11*matL11+k21*matL21+k31*matL31+k41*matL41;
 736   double v21 =  k01*matL02+k11*matL12+k21*matL22+k31*matL32+k41*matL42;
 737   double v31 =  k01*matL03+k11*matL13+k21*matL23+k31*matL33+k41*matL43;
 738   double v41 =  k01*matL04+k11*matL14+k21*matL24+k31*matL34+k41*matL44;
 739   //
 740   double t00 =  k00*matL00+k10*matL01+k20*matL02+k30*matL03+k40*matL04;
 741   double t10 =  k00*matL10+k10*matL11+k20*matL12+k30*matL13+k40*matL14;
 742   double t20 =  k00*matL20+k10*matL21+k20*matL22+k30*matL23+k40*matL24;
 743   double t30 =  k00*matL30+k10*matL31+k20*matL32+k30*matL33+k40*matL34;
 744   double t40 =  k00*matL40+k10*matL41+k20*matL42+k30*matL43+k40*matL44;
 745   //
 746   double t01 =  k01*matL00+k11*matL01+k21*matL02+k31*matL03+k41*matL04;
 747   double t11 =  k01*matL10+k11*matL11+k21*matL12+k31*matL13+k41*matL14;
 748   double t21 =  k01*matL20+k11*matL21+k21*matL22+k31*matL23+k41*matL24;
 749   double t31 =  k01*matL30+k11*matL31+k21*matL32+k31*matL33+k41*matL34;
 750   double t41 =  k01*matL40+k11*matL41+k21*matL42+k31*matL43+k41*matL44;
 751   //
 752   // (H^T * K^T * L * K * H) - (L * K * H) - (H^T * K^T * L) + (H^T*N^-1*H)
 753   matL00 += k00*v00+k10*v10+k20*v20*k30*v30+k40*v40 - t00 - v00 + fCovIYZ[0];
 754   matL01 += k01*v00+k11*v10+k21*v20*k31*v30+k41*v40 - t01 - v10 + fCovIYZ[1];
 755   matL10 += k00*v01+k10*v11+k20*v21*k30*v31+k40*v41 - t10 - v01 + fCovIYZ[1];
 756   matL11 += k01*v01+k11*v11+k21*v21*k31*v31+k41*v41 - t11 - v11 + fCovIYZ[2];
 757   //
 758   matL20 -= t20;
 759   matL21 -= t21;
 760   matL30 -= t30;
 761   matL31 -= t31;
 762   matL40 -= t40;
 763   matL41 -= t41;
 764   //
 765   matL02 -= v20;
 766   matL03 -= v30;
 767   matL04 -= v40;
 768   matL12 -= v21;
 769   matL13 -= v31;
 770   matL14 -= v41;
 771   //
 772   printf("Lr%d VecL: ",GetLayerID()); for (int i=0;i<5;i++) printf("%+e ",vecL[i]); printf("\n");
 773   printf("F: "); for (int i=0;i<kNFElem;i++) printf("%+e ",fFMatrix[i]); printf("\n");
 774   printf("K: "); for (int i=0;i<kNKElem;i++) printf("%+e ",fKMatrix[i]); printf("\n");
 775   //
 776   for (int j=0;j<5;j++) {
 777     for (int i=0;i<5;i++) printf("%+e ",matL[j][i]); printf("\n");
 778   }
 779   //
 780   return kTRUE;
 781 }
 782
 783  */
 784