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