ITS/AliITStrackV2.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////
  19 //                Implementation of the ITS track class
  20 //
  21 //          Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
  22 //     dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
  23 ///////////////////////////////////////////////////////////////////////////
  24 #include <TMath.h>
  25
  26 #include "AliCluster.h"
  27 #include "AliESDVertex.h"
  28 #include "AliITSReconstructor.h"
  29 #include "AliITStrackV2.h"
  30 #include "AliTracker.h"
  31 #include "AliLog.h"
  32 #include "AliPID.h"
  33
  34 const Int_t AliITStrackV2::fgkWARN = 5;
  35
  36 ClassImp(AliITStrackV2)
  37
  38
  39 //____________________________________________________________________________
  40 AliITStrackV2::AliITStrackV2() : AliKalmanTrack(),
  41   fCheckInvariant(kTRUE),
  42   fdEdx(0),
  43   fESDtrack(0)
  44 {
  45   for(Int_t i=0; i<2*AliITSgeomTGeo::kNLayers; i++) {fIndex[i]=-1; fModule[i]=-1;}
  46   for(Int_t i=0; i<AliITSgeomTGeo::kNLayers; i++) {fSharedWeight[i]=0;}
  47   for(Int_t i=0; i<4; i++) fdEdxSample[i]=0;
  48 }
  49
  50
  51 //____________________________________________________________________________
  52 AliITStrackV2::AliITStrackV2(AliESDtrack& t,Bool_t c):
  53   AliKalmanTrack(),
  54   fCheckInvariant(kTRUE),
  55   fdEdx(t.GetITSsignal()),
  56   fESDtrack(&t)
  57 {
  58   //------------------------------------------------------------------
  59   // Conversion ESD track -> ITS track.
  60   // If c==kTRUE, create the ITS track out of the constrained params.
  61   //------------------------------------------------------------------
  62   const AliExternalTrackParam *par=&t;
  63   if (c) {
  64     par=t.GetConstrainedParam();
  65     if (!par) AliError("AliITStrackV2: conversion failed !\n");
  66   }
  67   Set(par->GetX(),par->GetAlpha(),par->GetParameter(),par->GetCovariance());
  68
  69   SetLabel(t.GetITSLabel());
  70   SetMass(t.GetMassForTracking());
  71   SetNumberOfClusters(t.GetITSclusters(fIndex));
  72
  73   if (t.GetStatus()&AliESDtrack::kTIME) {
  74     StartTimeIntegral();
  75     Double_t times[AliPID::kSPECIESC];
  76     t.GetIntegratedTimes(times,AliPID::kSPECIESC);
  77     SetIntegratedTimes(times);
  78     SetIntegratedLength(t.GetIntegratedLength());
  79   }
  80
  81   for(Int_t i=0; i<AliITSgeomTGeo::kNLayers; i++) {fSharedWeight[i]=0;}
  82   for(Int_t i=0; i<4; i++) fdEdxSample[i]=0;
  83 }
  84
  85 //____________________________________________________________________________
  86 void AliITStrackV2::ResetClusters() {
  87   //------------------------------------------------------------------
  88   // Reset the array of attached clusters.
  89   //------------------------------------------------------------------
  90   for (Int_t i=0; i<2*AliITSgeomTGeo::kNLayers; i++) fIndex[i]=-1;
  91   for (Int_t i=0; i<AliITSgeomTGeo::kNLayers; i++) {fSharedWeight[i]=0;}
  92   SetChi2(0.);
  93   SetNumberOfClusters(0);
  94 }
  95
  96 //____________________________________________________________________________
  97 void AliITStrackV2::UpdateESDtrack(ULong_t flags) const {
  98   // Update track params
  99   fESDtrack->UpdateTrackParams(this,flags);
 100   //
 101   // set correctly the global label
 102   if (fESDtrack->IsOn(AliESDtrack::kTPCin)) {
 103     // for global track the GetLabel should be negative if
 104     // 1) GetTPCLabel<0
 105     // 2) this->GetLabel()<0
 106     // 3) GetTPCLabel() != this->GetLabel()
 107     int label = fESDtrack->GetTPCLabel();
 108     int itsLabel = GetLabel();
 109     if (label<0 || itsLabel<0 || itsLabel!=label) label = -TMath::Abs(label);
 110     fESDtrack->SetLabel(label);
 111   }
 112   //
 113   // copy the module indices
 114   Int_t i;
 115   for(i=0;i<2*AliITSgeomTGeo::kNLayers;i++) {
 116     //   printf("     %d\n",GetModuleIndex(i));
 117     fESDtrack->SetITSModuleIndex(i,GetModuleIndex(i));
 118   }
 119   // copy the map of shared clusters
 120   if(flags==AliESDtrack::kITSin) {
 121     UChar_t itsSharedMap=0;
 122     for(i=0;i<AliITSgeomTGeo::kNLayers;i++) {
 123       if(fSharedWeight[i]>0) SETBIT(itsSharedMap,i);
 124
 125     }
 126     fESDtrack->SetITSSharedMap(itsSharedMap);
 127   }
 128
 129   // copy the 4 dedx samples
 130   Double_t sdedx[4]={0.,0.,0.,0.};
 131   for(i=0; i<4; i++) sdedx[i]=fdEdxSample[i];
 132   fESDtrack->SetITSdEdxSamples(sdedx);
 133 }
 134
 135 //____________________________________________________________________________
 136 AliITStrackV2::AliITStrackV2(const AliITStrackV2& t) :
 137   AliKalmanTrack(t),
 138   fCheckInvariant(t.fCheckInvariant),
 139   fdEdx(t.fdEdx),
 140   fESDtrack(t.fESDtrack)
 141 {
 142   //------------------------------------------------------------------
 143   //Copy constructor
 144   //------------------------------------------------------------------
 145   Int_t i;
 146   for (i=0; i<4; i++) fdEdxSample[i]=t.fdEdxSample[i];
 147   for (i=0; i<2*AliITSgeomTGeo::GetNLayers(); i++) {
 148     fIndex[i]=t.fIndex[i];
 149     fModule[i]=t.fModule[i];
 150   }
 151   for (i=0; i<AliITSgeomTGeo::kNLayers; i++) {fSharedWeight[i]=t.fSharedWeight[i];}
 152 }
 153
 154 //_____________________________________________________________________________
 155 Int_t AliITStrackV2::Compare(const TObject *o) const {
 156   //-----------------------------------------------------------------
 157   // This function compares tracks according to the their curvature
 158   //-----------------------------------------------------------------
 159   AliITStrackV2 *t=(AliITStrackV2*)o;
 160   //Double_t co=OneOverPt();
 161   //Double_t c =OneOverPt();
 162   Double_t co=t->GetSigmaY2()*t->GetSigmaZ2();
 163   Double_t c =GetSigmaY2()*GetSigmaZ2();
 164   if (c>co) return 1;
 165   else if (c<co) return -1;
 166   return 0;
 167 }
 168
 169 //____________________________________________________________________________
 170 Bool_t
 171 AliITStrackV2::PropagateToVertex(const AliESDVertex *v,Double_t d,Double_t x0)
 172 {
 173   //------------------------------------------------------------------
 174   //This function propagates a track to the minimal distance from the origin
 175   //------------------------------------------------------------------
 176   Double_t bz=GetBz();
 177   if (PropagateToDCA(v,bz,kVeryBig)) {
 178     Double_t xOverX0,xTimesRho;
 179     xOverX0 = d; xTimesRho = d*x0;
 180     if (CorrectForMeanMaterial(xOverX0,xTimesRho,kTRUE)) return kTRUE;
 181   }
 182   return kFALSE;
 183 }
 184
 185 //____________________________________________________________________________
 186 Bool_t AliITStrackV2::
 187 GetGlobalXYZat(Double_t xloc, Double_t &x, Double_t &y, Double_t &z) const {
 188   //------------------------------------------------------------------
 189   //This function returns a track position in the global system
 190   //------------------------------------------------------------------
 191   Double_t r[3];
 192   Bool_t rc=GetXYZAt(xloc, GetBz(), r);
 193   x=r[0]; y=r[1]; z=r[2];
 194   return rc;
 195 }
 196
 197 //_____________________________________________________________________________
 198 Double_t AliITStrackV2::GetPredictedChi2(const AliCluster *c) const {
 199   //-----------------------------------------------------------------
 200   // This function calculates a predicted chi2 increment.
 201   //-----------------------------------------------------------------
 202   Double_t p[2]={c->GetY(), c->GetZ()};
 203   Double_t cov[3]={c->GetSigmaY2(), 0., c->GetSigmaZ2()};
 204   return AliExternalTrackParam::GetPredictedChi2(p,cov);
 205 }
 206
 207 //____________________________________________________________________________
 208 Bool_t AliITStrackV2::PropagateTo(Double_t xk, Double_t d, Double_t x0) {
 209   //------------------------------------------------------------------
 210   //This function propagates a track
 211   //------------------------------------------------------------------
 212
 213   Double_t oldX=GetX(), oldY=GetY(), oldZ=GetZ();
 214
 215   //Double_t bz=GetBz();
 216   //if (!AliExternalTrackParam::PropagateTo(xk,bz)) return kFALSE;
 217   Double_t b[3]; GetBxByBz(b);
 218   if (!AliExternalTrackParam::PropagateToBxByBz(xk,b)) return kFALSE;
 219   Double_t xOverX0,xTimesRho;
 220   xOverX0 = d; xTimesRho = d*x0;
 221   if (!CorrectForMeanMaterial(xOverX0,xTimesRho,kTRUE)) return kFALSE;
 222
 223   Double_t x=GetX(), y=GetY(), z=GetZ();
 224   if (IsStartedTimeIntegral() && x>oldX) {
 225     Double_t l2 = (x-oldX)*(x-oldX) + (y-oldY)*(y-oldY) + (z-oldZ)*(z-oldZ);
 226     AddTimeStep(TMath::Sqrt(l2));
 227   }
 228
 229   return kTRUE;
 230 }
 231
 232 //____________________________________________________________________________
 233 Bool_t AliITStrackV2::PropagateToTGeo(Double_t xToGo, Int_t nstep, Double_t &xOverX0, Double_t &xTimesRho, Bool_t addTime) {
 234   //-------------------------------------------------------------------
 235   //  Propagates the track to a reference plane x=xToGo in n steps.
 236   //  These n steps are only used to take into account the curvature.
 237   //  The material is calculated with TGeo. (L.Gaudichet)
 238   //-------------------------------------------------------------------
 239
 240   Double_t startx = GetX(), starty = GetY(), startz = GetZ();
 241   Double_t sign = (startx<xToGo) ? -1.:1.;
 242   Double_t step = (xToGo-startx)/TMath::Abs(nstep);
 243
 244   Double_t start[3], end[3], mparam[7];
 245   //Double_t bz = GetBz();
 246   Double_t b[3]; GetBxByBz(b);
 247   Double_t bz = b[2];
 248
 249   Double_t x = startx;
 250
 251   for (Int_t i=0; i<nstep; i++) {
 252
 253     GetXYZ(start);   //starting global position
 254     x += step;
 255     if (!GetXYZAt(x, bz, end)) return kFALSE;
 256     //if (!AliExternalTrackParam::PropagateTo(x, bz)) return kFALSE;
 257     if (!AliExternalTrackParam::PropagateToBxByBz(x, b)) return kFALSE;
 258     AliTracker::MeanMaterialBudget(start, end, mparam);
 259     xTimesRho = sign*mparam[4]*mparam[0];
 260     xOverX0   = mparam[1];
 261     if (mparam[1]<900000) {
 262       if (!AliExternalTrackParam::CorrectForMeanMaterial(xOverX0,
 263                            xTimesRho,GetMass())) return kFALSE;
 264     } else { // this happens when MeanMaterialBudget cannot cross a boundary
 265       return kFALSE;
 266     }
 267   }
 268
 269   if (addTime && IsStartedTimeIntegral() && GetX()>startx) {
 270     Double_t l2 = ( (GetX()-startx)*(GetX()-startx) +
 271                     (GetY()-starty)*(GetY()-starty) +
 272                     (GetZ()-startz)*(GetZ()-startz) );
 273     AddTimeStep(TMath::Sqrt(l2));
 274   }
 275
 276   return kTRUE;
 277 }
 278
 279 //____________________________________________________________________________
 280 Bool_t AliITStrackV2::Update(const AliCluster* c, Double_t chi2, Int_t index)
 281 {
 282   //------------------------------------------------------------------
 283   //This function updates track parameters
 284   //------------------------------------------------------------------
 285   Double_t p[2]={c->GetY(), c->GetZ()};
 286   Double_t cov[3]={c->GetSigmaY2(), c->GetSigmaYZ(), c->GetSigmaZ2()};
 287
 288   if (!AliExternalTrackParam::Update(p,cov)) return kFALSE;
 289
 290   Int_t n=GetNumberOfClusters();
 291   if (!Invariant()) {
 292     if (n>fgkWARN) AliDebug(1,"Wrong invariant !");
 293      return kFALSE;
 294   }
 295
 296   if (chi2<0) return kTRUE;
 297
 298   // fill residuals for ITS+TPC tracks
 299   if (fESDtrack) {
 300     if (fESDtrack->GetStatus()&AliESDtrack::kTPCin) {
 301       AliTracker::FillResiduals(this,p,cov,c->GetVolumeId());
 302     }
 303   }
 304
 305   fIndex[n]=index;
 306   SetNumberOfClusters(n+1);
 307   SetChi2(GetChi2()+chi2);
 308
 309   return kTRUE;
 310 }
 311
 312 Bool_t AliITStrackV2::Invariant() const {
 313   //------------------------------------------------------------------
 314   // This function is for debugging purpose only
 315   //------------------------------------------------------------------
 316   if(!fCheckInvariant) return kTRUE;
 317
 318   Int_t n=GetNumberOfClusters();
 319   static Float_t bz = GetBz();
 320   // take into account the misalignment error
 321   Float_t maxMisalErrY2=0,maxMisalErrZ2=0;
 322   //RS
 323   const AliITSRecoParam* recopar = AliITSReconstructor::GetRecoParam();
 324   if (!recopar) recopar = AliITSRecoParam::GetHighFluxParam();
 325
 326   for (Int_t lay=0; lay<AliITSgeomTGeo::kNLayers; lay++) {
 327     maxMisalErrY2 = TMath::Max(maxMisalErrY2,recopar->GetClusterMisalErrorY(lay,bz));
 328     maxMisalErrZ2 = TMath::Max(maxMisalErrZ2,recopar->GetClusterMisalErrorZ(lay,bz));
 329   }
 330   maxMisalErrY2 *= maxMisalErrY2;
 331   maxMisalErrZ2 *= maxMisalErrZ2;
 332   // this is because when we reset before refitting, we multiply the
 333   // matrix by 10
 334   maxMisalErrY2 *= 10.;
 335   maxMisalErrZ2 *= 10.;
 336
 337   Double_t sP2=GetParameter()[2];
 338   if (TMath::Abs(sP2) >= kAlmost1){
 339     if (n>fgkWARN) AliDebug(1,Form("fP2=%f\n",sP2));
 340      return kFALSE;
 341   }
 342   Double_t sC00=GetCovariance()[0];
 343   if (sC00<=0 || sC00>(9.+maxMisalErrY2)) {
 344     if (n>fgkWARN) AliDebug(1,Form("fC00=%f\n",sC00));
 345      return kFALSE;
 346   }
 347   Double_t sC11=GetCovariance()[2];
 348   if (sC11<=0 || sC11>(9.+maxMisalErrZ2)) {
 349     if (n>fgkWARN) AliDebug(1,Form("fC11=%f\n",sC11));
 350      return kFALSE;
 351   }
 352   Double_t sC22=GetCovariance()[5];
 353   if (sC22<=0 || sC22>1.) {
 354     if (n>fgkWARN) AliDebug(1,Form("fC22=%f\n",sC22));
 355      return kFALSE;
 356   }
 357   Double_t sC33=GetCovariance()[9];
 358   if (sC33<=0 || sC33>1.) {
 359     if (n>fgkWARN) AliDebug(1,Form("fC33=%f\n",sC33));
 360      return kFALSE;
 361   }
 362   Double_t sC44=GetCovariance()[14];
 363   if (sC44<=0 /*|| sC44>6e-5*/) {
 364     if (n>fgkWARN) AliDebug(1,Form("fC44=%f\n",sC44));
 365      return kFALSE;
 366   }
 367
 368   return kTRUE;
 369 }
 370
 371 //____________________________________________________________________________
 372 Bool_t AliITStrackV2::Propagate(Double_t alp,Double_t xk) {
 373   //------------------------------------------------------------------
 374   //This function propagates a track
 375   //------------------------------------------------------------------
 376   //Double_t bz=GetBz();
 377   //if (!AliExternalTrackParam::Propagate(alp,xk,bz)) return kFALSE;
 378   Double_t b[3]; GetBxByBz(b);
 379   if (!AliExternalTrackParam::PropagateBxByBz(alp,xk,b)) return kFALSE;
 380
 381   if (!Invariant()) {
 382     Int_t n=GetNumberOfClusters();
 383     if (n>fgkWARN) AliDebug(1,"Wrong invariant !");
 384     return kFALSE;
 385   }
 386
 387   return kTRUE;
 388 }
 389
 390 Bool_t AliITStrackV2::MeanBudgetToPrimVertex(Double_t xyz[3], Double_t step, Double_t &d) const {
 391
 392   //-------------------------------------------------------------------
 393   //  Get the mean material budget between the actual point and the
 394   //  primary vertex. (L.Gaudichet)
 395   //-------------------------------------------------------------------
 396
 397   Double_t cs=TMath::Cos(GetAlpha()), sn=TMath::Sin(GetAlpha());
 398   Double_t vertexX = xyz[0]*cs + xyz[1]*sn;
 399
 400   Int_t nstep = Int_t((GetX()-vertexX)/step);
 401   if (nstep<1) nstep = 1;
 402   step = (GetX()-vertexX)/nstep;
 403
 404   //  Double_t mparam[7], densMean=0, radLength=0, length=0;
 405   Double_t mparam[7];
 406   Double_t p1[3], p2[3], x = GetX(), bz = GetBz();
 407   GetXYZ(p1);
 408
 409   d=0.;
 410
 411   for (Int_t i=0; i<nstep; i++) {
 412     x  += step;
 413     if (!GetXYZAt(x, bz, p2)) return kFALSE;
 414     AliTracker::MeanMaterialBudget(p1, p2, mparam);
 415     if (mparam[1]>900000) return kFALSE;
 416     d  += mparam[1];
 417
 418     p1[0] = p2[0];
 419     p1[1] = p2[1];
 420     p1[2] = p2[2];
 421   }
 422
 423   return kTRUE;
 424 }
 425
 426 Bool_t AliITStrackV2::Improve(Double_t x0,Double_t xyz[3],Double_t ers[3]) {
 427   //------------------------------------------------------------------
 428   //This function improves angular track parameters
 429   //------------------------------------------------------------------
 430   //Store the initail track parameters
 431
 432   Double_t x = GetX();
 433   Double_t alpha = GetAlpha();
 434   Double_t par[] = {GetY(),GetZ(),GetSnp(),GetTgl(),GetSigned1Pt()};
 435   Double_t cov[] = {
 436     GetSigmaY2(),
 437     GetSigmaZY(),
 438     GetSigmaZ2(),
 439     GetSigmaSnpY(),
 440     GetSigmaSnpZ(),
 441     GetSigmaSnp2(),
 442     GetSigmaTglY(),
 443     GetSigmaTglZ(),
 444     GetSigmaTglSnp(),
 445     GetSigmaTgl2(),
 446     GetSigma1PtY(),
 447     GetSigma1PtZ(),
 448     GetSigma1PtSnp(),
 449     GetSigma1PtTgl(),
 450     GetSigma1Pt2()
 451   };
 452
 453
 454   Double_t cs=TMath::Cos(GetAlpha()), sn=TMath::Sin(GetAlpha());
 455   Double_t xv = xyz[0]*cs + xyz[1]*sn; // vertex
 456   Double_t yv =-xyz[0]*sn + xyz[1]*cs; // in the
 457   Double_t zv = xyz[2];                // local frame
 458
 459   Double_t dx = x - xv, dy = par[0] - yv, dz = par[1] - zv;
 460   Double_t r2=dx*dx + dy*dy;
 461   Double_t p2=(1.+ GetTgl()*GetTgl())/(GetSigned1Pt()*GetSigned1Pt());
 462   if (GetMass()<0) p2 *= 4; // q=2
 463   Double_t beta2=p2/(p2 + GetMass()*GetMass());
 464   x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
 465   Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
 466   //Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*x0*9.36*2.33;
 467
 468   Double_t bz=GetBz();
 469   Double_t cnv=bz*kB2C;
 470   Double_t curv=GetC(bz);
 471   {
 472     Double_t dummy = 4/r2 - curv*curv;
 473     if (dummy < 0) return kFALSE;
 474     Double_t parp = 0.5*(curv*dx + dy*TMath::Sqrt(dummy));
 475     Double_t sigma2p = theta2*(1.-GetSnp())*(1.+GetSnp())*(1. + GetTgl()*GetTgl());
 476     Double_t ovSqr2 = 1./TMath::Sqrt(r2);
 477     Double_t tfact = ovSqr2*(1.-dy*ovSqr2)*(1.+dy*ovSqr2);
 478     sigma2p += cov[0]*tfact*tfact;
 479     sigma2p += ers[1]*ers[1]/r2;
 480     sigma2p += 0.25*cov[14]*cnv*cnv*dx*dx;
 481     Double_t eps2p=sigma2p/(cov[5] + sigma2p);
 482     par[0] += cov[3]/(cov[5] + sigma2p)*(parp - GetSnp());
 483     par[2]  = eps2p*GetSnp() + (1 - eps2p)*parp;
 484     cov[5] *= eps2p;
 485     cov[3] *= eps2p;
 486   }
 487   {
 488     Double_t parl=0.5*curv*dz/TMath::ASin(0.5*curv*TMath::Sqrt(r2));
 489     Double_t sigma2l=theta2;
 490     sigma2l += cov[2]/r2 + cov[0]*dy*dy*dz*dz/(r2*r2*r2);
 491     sigma2l += ers[2]*ers[2]/r2;
 492     Double_t eps2l = sigma2l/(cov[9] + sigma2l);
 493     par[1] += cov[7 ]/(cov[9] + sigma2l)*(parl - par[3]);
 494     par[4] += cov[13]/(cov[9] + sigma2l)*(parl - par[3]);
 495     par[3]  = eps2l*par[3] + (1-eps2l)*parl;
 496     cov[9] *= eps2l;
 497     cov[13]*= eps2l;
 498     cov[7] *= eps2l;
 499   }
 500
 501   Set(x,alpha,par,cov);
 502
 503   if (!Invariant()) return kFALSE;
 504
 505   return kTRUE;
 506 }
 507
 508 void AliITStrackV2::CookdEdx(Double_t /*low*/, Double_t /*up*/) {
 509   //-----------------------------------------------------------------
 510   // This function calculates dE/dX within the "low" and "up" cuts.
 511   // Origin: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
 512   // Updated: F. Prino 8-June-2009
 513   //-----------------------------------------------------------------
 514   // The cluster order is: SDD-1, SDD-2, SSD-1, SSD-2
 515
 516   Int_t nc=0;
 517   Float_t dedx[4];
 518   for (Int_t il=0; il<4; il++) { // count good (>0) dE/dx values
 519     if(fdEdxSample[il]>0.){
 520       dedx[nc]= fdEdxSample[il];
 521       nc++;
 522     }
 523   }
 524   if(nc<1){
 525     SetdEdx(0.);
 526     return;
 527   }
 528
 529   Int_t swap; // sort in ascending order
 530   do {
 531     swap=0;
 532     for (Int_t i=0; i<nc-1; i++) {
 533       if (dedx[i]<=dedx[i+1]) continue;
 534       Float_t tmp=dedx[i];
 535       dedx[i]=dedx[i+1];
 536       dedx[i+1]=tmp;
 537       swap++;
 538     }
 539   } while (swap);
 540
 541
 542   Double_t sumamp=0,sumweight=0;
 543   Double_t weight[4]={1.,1.,0.,0.};
 544   if(nc==3) weight[1]=0.5;
 545   else if(nc<3) weight[1]=0.;
 546   for (Int_t i=0; i<nc; i++) {
 547     sumamp+= dedx[i]*weight[i];
 548     sumweight+=weight[i];
 549   }
 550   SetdEdx(sumamp/sumweight);
 551 }
 552
 553 //____________________________________________________________________________
 554 Bool_t AliITStrackV2::
 555 GetPhiZat(Double_t r, Double_t &phi, Double_t &z) const {
 556   //------------------------------------------------------------------
 557   // This function returns the global cylindrical (phi,z) of the track
 558   // position estimated at the radius r.
 559   // The track curvature is neglected.
 560   //------------------------------------------------------------------
 561   Double_t d=GetD(0.,0.);
 562   if (TMath::Abs(d) > r) {
 563     if (r>1e-1) return kFALSE;
 564     r = TMath::Abs(d);
 565   }
 566
 567   Double_t rcurr=TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
 568   if (TMath::Abs(d) > rcurr) return kFALSE;
 569   Double_t globXYZcurr[3]; GetXYZ(globXYZcurr);
 570   Double_t phicurr=TMath::ATan2(globXYZcurr[1],globXYZcurr[0]);
 571
 572   if (GetX()>=0.) {
 573     phi=phicurr+TMath::ASin(d/r)-TMath::ASin(d/rcurr);
 574   } else {
 575     phi=phicurr+TMath::ASin(d/r)+TMath::ASin(d/rcurr)-TMath::Pi();
 576   }
 577
 578   // return a phi in [0,2pi[
 579   if (phi<0.) phi+=2.*TMath::Pi();
 580   else if (phi>=2.*TMath::Pi()) phi-=2.*TMath::Pi();
 581   z=GetZ()+GetTgl()*(TMath::Sqrt((r-d)*(r+d))-TMath::Sqrt((rcurr-d)*(rcurr+d)));
 582   return kTRUE;
 583 }
 584 //____________________________________________________________________________
 585 Bool_t AliITStrackV2::
 586 GetLocalXat(Double_t r,Double_t &xloc) const {
 587   //------------------------------------------------------------------
 588   // This function returns the local x of the track
 589   // position estimated at the radius r.
 590   // The track curvature is neglected.
 591   //------------------------------------------------------------------
 592   Double_t d=GetD(0.,0.);
 593   if (TMath::Abs(d) > r) {
 594     if (r>1e-1) return kFALSE;
 595     r = TMath::Abs(d);
 596   }
 597
 598   Double_t rcurr=TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
 599   Double_t globXYZcurr[3]; GetXYZ(globXYZcurr);
 600   Double_t phicurr=TMath::ATan2(globXYZcurr[1],globXYZcurr[0]);
 601   Double_t phi;
 602   if (GetX()>=0.) {
 603     phi=phicurr+TMath::ASin(d/r)-TMath::ASin(d/rcurr);
 604   } else {
 605     phi=phicurr+TMath::ASin(d/r)+TMath::ASin(d/rcurr)-TMath::Pi();
 606   }
 607
 608   xloc=r*(TMath::Cos(phi)*TMath::Cos(GetAlpha())
 609          +TMath::Sin(phi)*TMath::Sin(GetAlpha()));
 610
 611   return kTRUE;
 612 }
 613
 614 //____________________________________________________________________________
 615 Bool_t AliITStrackV2::ImproveKalman(Double_t xyz[3],Double_t ers[3], const Double_t* xlMS, const Double_t* x2X0MS, Int_t nMS)
 616 {
 617   // Substitute the state of the track (p_{k|k},C_{k|k}) at the k-th measumerent by its
 618   // smoothed value from the k-th measurement + measurement at the vertex.
 619   // Account for the MS on nMS layers at x-postions xlMS with x/x0 = x2X0MS
 620   // p_{k|kv} = p_{k|k} + C_{k|k}*D^Tr_{k+1} B^{-1}_{k+1} ( vtx - D_{k+1}*p_{k|k})
 621   // C_{k|kv} = C_{k|k}*( I - D^Tr_{k+1} B^{-1}_{k+1} D_{k+1} C^Tr_{k|k})
 622   //
 623   // where D_{k} = H_{k} F_{k} with H being the matrix converting the tracks parameters
 624   // to measurements m_{k} = H_{k} p_{k} and F_{k} the matrix propagating the track between the
 625   // the point k-1 and k:  p_{k|k-1} = F_{k} p_{k-1|k-1}
 626   //
 627   // B_{k+1} = V_{k+1} + H_{k+1} C_{k+1|k} H^Tr_{k+1} with V_{k+1} being the error of the measurment
 628   // at point k+1 (i.e. vertex), and C_{k+1|k} - error matrix extrapolated from k-th measurement to
 629   // k+1 (vtx) and accounting for the MS inbetween
 630   //
 631   // H = {{1,0,0,0,0},{0,1,0,0,0}}
 632   //
 633   double covc[15], *cori = (double*) GetCovariance(),par[5] = {GetY(),GetZ(),GetSnp(),GetTgl(),GetSigned1Pt()},
 634     &c00=cori[0],
 635     &c01=cori[1],&c11=cori[2],
 636     &c02=cori[3],&c12=cori[4],&c22=cori[5],
 637     &c03=cori[6],&c13=cori[7],&c23=cori[8],&c33=cori[9],
 638     &c04=cori[10],&c14=cori[11],&c24=cori[12],&c34=cori[13],&c44=cori[14],
 639     // for smoothed cov matrix
 640     &cov00=covc[0],
 641     &cov01=covc[1],&cov11=covc[2],
 642     &cov02=covc[3],&cov12=covc[4],&cov22=covc[5],
 643     &cov03=covc[6],&cov13=covc[7],&cov23=covc[8],&cov33=covc[9],
 644     &cov04=covc[10],&cov14=covc[11],&cov24=covc[12],&cov34=covc[13],&cov44=covc[14];
 645   //
 646   double x = GetX(), alpha = GetAlpha();
 647   // vertex in the track frame
 648   double cs=TMath::Cos(alpha), sn=TMath::Sin(alpha);
 649   double xv = xyz[0]*cs + xyz[1]*sn, yv =-xyz[0]*sn + xyz[1]*cs, zv = xyz[2];
 650   double dx = xv - GetX();
 651   if (TMath::Abs(dx)<=kAlmost0)  return kTRUE;
 652   //
 653   double cnv=GetBz()*kB2C, x2r=cnv*par[4]*dx, f1=par[2], f2=f1+x2r;
 654   if (TMath::Abs(f1) >= kAlmost1 || TMath::Abs(f2) >= kAlmost1) {
 655     AliInfo(Form("Fail: %+e %+e",f1,f2));
 656     return kFALSE;
 657   }
 658   double r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2)), dx2r=dx/(r1+r2);
 659   // elements of matrix F_{k+1} (1s on diagonal)
 660   double f02 = 2*dx2r, f04 = cnv*dx*dx2r, f13/*, f24 = cnv*dx*/;
 661   if (TMath::Abs(x2r)<0.05) f13 = dx*r2+f2*(f1+f2)*dx2r; // see AliExternalTrackParam::PropagateTo
 662   else {
 663     double dy2dx = (f1+f2)/(r1+r2);
 664     f13 = 2*TMath::ASin(0.5*TMath::Sqrt(1+dy2dx*dy2dx)*x2r)/(cnv*par[4]);
 665   }
 666   // elements of matrix D_{k+1} = H_{k+1} * F_{k+1}
 667   // double d00 = 1., d11 = 1.;
 668   double &d02 = f02,  &d04 = f04, &d13 = f13;
 669   //
 670   // elements of matrix DC = D_{k+1}*C_{kk}^T
 671   double dc00 = c00+c02*d02+c04*d04,   dc10 = c01+c03*d13;
 672   double dc01 = c01+c12*d02+c14*d04,   dc11 = c11+c13*d13;
 673   double dc02 = c02+c22*d02+c24*d04,   dc12 = c12+c23*d13;
 674   double dc03 = c03+c23*d02+c34*d04,   dc13 = c13+c33*d13;
 675   double dc04 = c04+c24*d02+c44*d04,   dc14 = c14+c34*d13;
 676   //
 677   // difference between the vertex and the the track extrapolated to vertex
 678   yv -= par[0] + par[2]*d02 + par[4]*d04;
 679   zv -= par[1] + par[3]*d13;
 680   //
 681   // y,z part of the cov.matrix extrapolated to vtx (w/o MS contribution)
 682   // C_{k+1,k} = H F_{k+1} C_{k,k} F^Tr_{k+1} H^Tr = D C D^Tr
 683   double cv00 = dc00+dc02*d02+dc04*d04, cv01 = dc01+dc03*d13, cv11 = dc11+dc13*d13;
 684   //
 685   // add MS contribution layer by layer
 686   double xCurr = x;
 687   double p2Curr = par[2];
 688   //
 689   // precalculated factors of MS contribution matrix:
 690   double ms22t = (1. + par[3]*par[3]);
 691   double ms33t = ms22t*ms22t;
 692   double p34 = par[3]*par[4];
 693   double ms34t = p34*ms22t;
 694   double ms44t = p34*p34;
 695   //
 696   double p2=(1.+ par[3]*par[3])/(par[4]*par[4]);
 697   if (GetMass()<0) p2 *= 4; // q=2
 698   double beta2 = p2/(p2+GetMass()*GetMass());
 699   double theta2t = 14.1*14.1/(beta2*p2*1e6) * (1. + par[3]*par[3]);
 700   //
 701   // account for the MS in the layers between the last measurement and the vertex
 702   for (int il=0;il<nMS;il++) {
 703     double dfx = xlMS[il] - xCurr;
 704     xCurr = xlMS[il];
 705     p2Curr += dfx*cnv*par[4];   // p2 at the scattering layer
 706     double dxL=xv - xCurr;    // distance from scatering layer to vtx
 707     double x2rL=cnv*par[4]*dxL, f1L=p2Curr, f2L=f1L+x2rL;
 708     if (TMath::Abs(f1L) >= kAlmost1 || TMath::Abs(f2L) >= kAlmost1) {
 709       AliInfo(Form("FailMS at step %d of %d: dfx:%e dxL:%e %e %e",il,nMS,dfx,dxL,f1L,f2L));
 710       return kFALSE;
 711     }
 712     double r1L=TMath::Sqrt((1.-f1L)*(1.+f1L)), r2L=TMath::Sqrt((1.-f2L)*(1.+f2L)), dx2rL=dxL/(r1L+r2L);
 713     // elements of matrix for propagation from scatering layer to vertex
 714     double f02L = 2*dx2rL, f04L = cnv*dxL*dx2rL, f13L/*, f24L = cnv*dxL*/;
 715     if (TMath::Abs(x2rL)<0.05) f13L = dxL*r2L+f2L*(f1L+f2L)*dx2rL; // see AliExternalTrackParam::PropagateTo
 716     else {
 717       double dy2dxL = (f1L+f2L)/(r1L+r2L);
 718       f13L = 2*TMath::ASin(0.5*TMath::Sqrt(1+dy2dxL*dy2dxL)*x2rL)/(cnv*par[4]);
 719     }
 720     // MS contribution matrix:
 721     double theta2 = theta2t*TMath::Abs(x2X0MS[il]);
 722     double ms22 = theta2*(1.-p2Curr)*(1.+p2Curr)*ms22t;
 723     double ms33 = theta2*ms33t;
 724     double ms34 = theta2*ms34t;
 725     double ms44 = theta2*ms44t;
 726     //
 727     // add  H F MS F^Tr H^Tr to cv
 728     cv00 += f02L*f02L*ms22 + f04L*f04L*ms44;
 729     cv01 += f04L*f13L*ms34;
 730     cv11 += f13L*f13L*ms33;
 731   }
 732   //
 733   // inverse of matrix B
 734   double b11 = ers[1]*ers[1] + cv00;
 735   double b00 = ers[2]*ers[2] + cv11;
 736   double det = b11*b00 - cv01*cv01;
 737   if (TMath::Abs(det)<kAlmost0) {
 738     AliInfo(Form("Fail on det %e: %e %e %e",det,cv00,cv11,cv01));
 739     return kFALSE;
 740   }
 741   det = 1./det;
 742   b00 *= det; b11 *= det;
 743   double b01 = -cv01*det;
 744   //
 745   // elements of matrix DC^Tr * B^-1
 746   double dcb00 = b00*dc00+b01*dc10, dcb01 = b01*dc00+b11*dc10;
 747   double dcb10 = b00*dc01+b01*dc11, dcb11 = b01*dc01+b11*dc11;
 748   double dcb20 = b00*dc02+b01*dc12, dcb21 = b01*dc02+b11*dc12;
 749   double dcb30 = b00*dc03+b01*dc13, dcb31 = b01*dc03+b11*dc13;
 750   double dcb40 = b00*dc04+b01*dc14, dcb41 = b01*dc04+b11*dc14;
 751   //
 752   // p_{k|k+1} = p_{k|k} +  C_{k|k}*D^Tr_{k+1} B^{-1}_{k+1} ( vtx - D_{k+1}*p_{k|k})
 753   par[0] += dcb00*yv + dcb01*zv;
 754   par[1] += dcb10*yv + dcb11*zv;
 755   par[2] += dcb20*yv + dcb21*zv;
 756   par[3] += dcb30*yv + dcb31*zv;
 757   par[4] += dcb40*yv + dcb41*zv;
 758   //
 759   // C_{k|kv} = C_{k|k} - C_{k|k} D^Tr_{k+1} B^{-1}_{k+1} D_{k+1} C^Tr_{k|k})
 760   //
 761   cov00 = c00 - (dc00*dcb00 + dc10*dcb01);
 762   cov01 = c01 - (dc01*dcb00 + dc11*dcb01);
 763   cov02 = c02 - (dc02*dcb00 + dc12*dcb01);
 764   cov03 = c03 - (dc03*dcb00 + dc13*dcb01);
 765   cov04 = c04 - (dc04*dcb00 + dc14*dcb01);
 766   //
 767   cov11 = c11 - (dc01*dcb10 + dc11*dcb11);
 768   cov12 = c12 - (dc02*dcb10 + dc12*dcb11);
 769   cov13 = c13 - (dc03*dcb10 + dc13*dcb11);
 770   cov14 = c14 - (dc04*dcb10 + dc14*dcb11);
 771   //
 772   cov22 = c22 - (dc02*dcb20 + dc12*dcb21);
 773   cov23 = c23 - (dc03*dcb20 + dc13*dcb21);
 774   cov24 = c24 - (dc04*dcb20 + dc14*dcb21);
 775   //
 776   cov33 = c33 - (dc03*dcb30 + dc13*dcb31);
 777   cov34 = c34 - (dc04*dcb30 + dc14*dcb31);
 778   //
 779   cov44 = c44 - (dc04*dcb40 + dc14*dcb41);
 780   //
 781   Set(x,alpha,par,covc);
 782   if (!Invariant()) {
 783     AliInfo(Form("Fail on Invariant, X=%e",GetX()));
 784     return kFALSE;
 785   }
 786   return kTRUE;
 787   //
 788 }