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