TPC/AliTPCtrack.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 TPC track class
  20 //        This class is used by the AliTPCtracker class
  21 //      Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
  22 //-----------------------------------------------------------------
  23
  24 #include <Riostream.h>
  25
  26 #include "AliTPCtrack.h"
  27 #include "AliCluster.h"
  28 #include "AliTracker.h"
  29 #include "AliESDtrack.h"
  30 #include "AliESDVertex.h"
  31 #include "TTreeStream.h"
  32
  33 ClassImp(AliTPCtrack)
  34
  35 //_________________________________________________________________________
  36 AliTPCtrack::AliTPCtrack():
  37   AliKalmanTrack(),
  38   fdEdx(0),
  39   fSdEdx(1e10),
  40   fNFoundable(0),
  41   fBConstrain(kFALSE),
  42   fLastPoint(-1),
  43   fFirstPoint(-1),
  44   fRemoval(0),
  45   fTrackType(0),
  46   fLab2(-1),
  47   fNShared(0),
  48   fReference()
  49 {
  50   //-------------------------------------------------
  51   // default constructor
  52   //-------------------------------------------------
  53   for (Int_t i=0; i<kMaxRow;i++) fIndex[i]=-2;
  54   for (Int_t i=0; i<4;i++) fPoints[i]=0.;
  55   for (Int_t i=0; i<12;i++) fKinkPoint[i]=0.;
  56   for (Int_t i=0; i<3;i++) fKinkIndexes[i]=0;
  57   for (Int_t i=0; i<3;i++) fV0Indexes[i]=0;
  58 }
  59
  60 //_________________________________________________________________________
  61
  62
  63
  64 AliTPCtrack::AliTPCtrack(Double_t x, Double_t alpha, const Double_t p[5],
  65                          const Double_t cov[15], Int_t index) :
  66   AliKalmanTrack(),
  67   fdEdx(0),
  68   fSdEdx(1e10),
  69   fNFoundable(0),
  70   fBConstrain(kFALSE),
  71   fLastPoint(0),
  72   fFirstPoint(0),
  73   fRemoval(0),
  74   fTrackType(0),
  75   fLab2(0),
  76   fNShared(0),
  77   fReference()
  78 {
  79   //-----------------------------------------------------------------
  80   // This is the main track constructor.
  81   //-----------------------------------------------------------------
  82   Double_t cnv=1./(AliTracker::GetBz()*kB2C);
  83
  84   Double_t pp[5]={
  85     p[0],
  86     p[1],
  87     x*p[4] - p[2],
  88     p[3],
  89     p[4]*cnv
  90   };
  91
  92   Double_t c22 = x*x*cov[14] - 2*x*cov[12] + cov[5];
  93   Double_t c32 = x*cov[13] - cov[8];
  94   Double_t c20 = x*cov[10] - cov[3],
  95            c21 = x*cov[11] - cov[4], c42 = x*cov[14] - cov[12];
  96
  97   Double_t cc[15]={
  98     cov[0 ],
  99     cov[1 ],     cov[2 ],
 100     c20,         c21,         c22,
 101     cov[6 ],     cov[7 ],     c32,     cov[9 ],
 102     cov[10]*cnv, cov[11]*cnv, c42*cnv, cov[13]*cnv, cov[14]*cnv*cnv
 103   };
 104
 105   Double_t mostProbablePt=AliExternalTrackParam::GetMostProbablePt();
 106   Double_t p0=TMath::Sign(1/mostProbablePt,pp[4]);
 107   Double_t w0=cc[14]/(cc[14] + p0*p0), w1=p0*p0/(cc[14] + p0*p0);
 108   pp[4] = w0*p0 + w1*pp[4];
 109   cc[10]*=w1; cc[11]*=w1; cc[12]*=w1; cc[13]*=w1; cc[14]*=w1;
 110
 111   Set(x,alpha,pp,cc);
 112
 113   SetNumberOfClusters(1);
 114
 115   fIndex[0]=index;
 116   for (Int_t i=1; i<kMaxRow;i++) fIndex[i]=-2;
 117   for (Int_t i=0; i<4;i++) fPoints[i]=0.;
 118   for (Int_t i=0; i<12;i++) fKinkPoint[i]=0.;
 119   for (Int_t i=0; i<3;i++) fKinkIndexes[i]=0;
 120   for (Int_t i=0; i<3;i++) fV0Indexes[i]=0;
 121 }
 122
 123 //_____________________________________________________________________________
 124 AliTPCtrack::AliTPCtrack(const AliESDtrack& t, TTreeSRedirector *pcstream) :
 125   AliKalmanTrack(),
 126   fdEdx(t.GetTPCsignal()),
 127   fSdEdx(1e10),
 128   fNFoundable(0),
 129   fBConstrain(kFALSE),
 130   fLastPoint(0),
 131   fFirstPoint(0),
 132   fRemoval(0),
 133   fTrackType(0),
 134   fLab2(0),
 135   fNShared(0),
 136   fReference()
 137 {
 138   //-----------------------------------------------------------------
 139   // Conversion AliESDtrack -> AliTPCtrack.
 140   //-----------------------------------------------------------------
 141   const Double_t kmaxC[4]={10,10,0.1,0.1};  // cuts on the rms /fP0,fP1,fP2,fP3
 142   SetNumberOfClusters(t.GetTPCclusters(fIndex));
 143   SetLabel(t.GetLabel());
 144   SetMass(t.GetMass());
 145   for (Int_t i=0; i<4;i++) fPoints[i]=0.;
 146   for (Int_t i=0; i<12;i++) fKinkPoint[i]=0.;
 147   for (Int_t i=0; i<3;i++) fKinkIndexes[i]=0;
 148   for (Int_t i=0; i<3;i++) fV0Indexes[i]=0;
 149   //
 150   // choose parameters to start
 151   //
 152   Int_t reject=0;
 153   AliExternalTrackParam param(t);
 154   const AliExternalTrackParam  *tpcout=(t.GetFriendTrack())? ((AliESDfriendTrack*)(t.GetFriendTrack()))->GetTPCOut():0;
 155   const AliExternalTrackParam  *tpcin = t.GetInnerParam();
 156   const AliExternalTrackParam  *tpc=(tpcout)?tpcout:tpcin;
 157   if (!tpc) tpc=&param;
 158   Bool_t isOK=kTRUE;
 159   if (param.GetCovariance()[0]>kmaxC[0]*kmaxC[0]) isOK=kFALSE;
 160   if (param.GetCovariance()[2]>kmaxC[1]*kmaxC[1]) isOK=kFALSE;
 161   if (param.GetCovariance()[5]>kmaxC[2]*kmaxC[2]) isOK=kFALSE;
 162   if (param.GetCovariance()[9]>kmaxC[3]*kmaxC[3]) isOK=kFALSE;
 163   if (!isOK){
 164     param=*tpc;
 165     isOK=kTRUE;
 166     reject=1;
 167   }
 168   param.Rotate(tpc->GetAlpha());
 169   isOK=AliTracker::PropagateTrackToBxByBz(&param,tpc->GetX(),t.GetMass(),2.,kFALSE);
 170   if (param.GetCovariance()[0]>kmaxC[0]*kmaxC[0]) isOK=kFALSE;
 171   if (param.GetCovariance()[2]>kmaxC[1]*kmaxC[1]) isOK=kFALSE;
 172   if (param.GetCovariance()[5]>kmaxC[2]*kmaxC[2]) isOK=kFALSE;
 173   if (param.GetCovariance()[9]>kmaxC[3]*kmaxC[3]) isOK=kFALSE;
 174   if (!isOK){
 175     param=*tpc;
 176     isOK=kTRUE;
 177     reject=2;
 178   }
 179   if (reject>0){
 180     param.ResetCovariance(4.);  // reset covariance if start from backup param
 181   }
 182   //
 183   //
 184   if (pcstream){
 185     AliExternalTrackParam dummy;
 186     AliExternalTrackParam *ptpc=(AliExternalTrackParam *)tpc;
 187     if (!ptpc) ptpc=&dummy;
 188     AliESDtrack *esd= (AliESDtrack *)&t;
 189     (*pcstream)<<"trackP"<<
 190       "reject="<<reject<<   // flag - rejection of current esd track parameters
 191       "esd.="<<esd<<        // original esd track
 192       "tr.="<<&param<<      // starting track parameters
 193       "out.="<<ptpc<<       // backup tpc parameters
 194       "\n";
 195   }
 196
 197   Set(param.GetX(),param.GetAlpha(),param.GetParameter(),param.GetCovariance());
 198
 199   if ((t.GetStatus()&AliESDtrack::kTIME) == 0) return;
 200   StartTimeIntegral();
 201   Double_t times[10]; t.GetIntegratedTimes(times); SetIntegratedTimes(times);
 202   SetIntegratedLength(t.GetIntegratedLength());
 203 }
 204
 205 //_____________________________________________________________________________
 206 AliTPCtrack::AliTPCtrack(const AliTPCtrack& t) :
 207   AliKalmanTrack(t),
 208   fdEdx(t.fdEdx),
 209   fSdEdx(t.fSdEdx),
 210   fNFoundable(t.fNFoundable),
 211   fBConstrain(t.fBConstrain),
 212   fLastPoint(t.fLastPoint),
 213   fFirstPoint(t.fFirstPoint),
 214   fRemoval(t.fRemoval),
 215   fTrackType(t.fTrackType),
 216   fLab2(t.fLab2),
 217   fNShared(t.fNShared),
 218   fReference(t.fReference)
 219
 220 {
 221   //-----------------------------------------------------------------
 222   // This is a track copy constructor.
 223   //-----------------------------------------------------------------
 224   Set(t.GetX(),t.GetAlpha(),t.GetParameter(),t.GetCovariance());
 225
 226   for (Int_t i=0; i<kMaxRow; i++) fIndex[i]=t.fIndex[i];
 227   for (Int_t i=0; i<4;i++) fPoints[i]=t.fPoints[i];
 228   for (Int_t i=0; i<12;i++) fKinkPoint[i]=t.fKinkPoint[i];
 229   for (Int_t i=0; i<3;i++) fKinkIndexes[i]=t.fKinkIndexes[i];
 230   for (Int_t i=0; i<3;i++) fV0Indexes[i]=t.fV0Indexes[i];
 231 }
 232
 233 AliTPCtrack& AliTPCtrack::operator=(const AliTPCtrack& o){
 234   if(this!=&o){
 235     AliKalmanTrack::operator=(o);
 236     fdEdx = o.fdEdx;
 237     for(Int_t i = 0;i<kMaxRow;++i)fIndex[i] = o.fIndex[i];
 238     for(Int_t i = 0;i<4;++i)fPoints[i] = o.fPoints[i];
 239     fSdEdx = o.fSdEdx;
 240     fNFoundable = o.fNFoundable;
 241     fBConstrain = o.fBConstrain;
 242     fLastPoint  = o.fLastPoint;
 243     fFirstPoint = o.fFirstPoint;
 244     fTrackType  = o.fTrackType;
 245     fLab2       = o.fLab2;
 246     fNShared    = o.fNShared;
 247     fReference  = o.fReference;
 248     for(Int_t i = 0;i<12;++i) fKinkPoint[i] = o.fKinkPoint[i];
 249
 250     for(Int_t i = 0;i<3;++i){
 251       fKinkIndexes[i] = o.fKinkIndexes[i];
 252       fV0Indexes[i] = o.fV0Indexes[i];
 253     }
 254   }
 255   return *this;
 256
 257 }
 258
 259
 260 //_____________________________________________________________________________
 261 Int_t AliTPCtrack::Compare(const TObject *o) const {
 262   //-----------------------------------------------------------------
 263   // This function compares tracks according to the their curvature
 264   //-----------------------------------------------------------------
 265   AliTPCtrack *t=(AliTPCtrack*)o;
 266   //Double_t co=t->OneOverPt();
 267   //Double_t c = OneOverPt();
 268   Double_t co=t->GetSigmaY2()*t->GetSigmaZ2();
 269   Double_t c =GetSigmaY2()*GetSigmaZ2();
 270   if (c>co) return 1;
 271   else if (c<co) return -1;
 272   return 0;
 273 }
 274
 275 Double_t AliTPCtrack::GetPredictedChi2(const AliCluster *c) const {
 276   //-----------------------------------------------------------------
 277   // This function calculates a predicted chi2 increment.
 278   //-----------------------------------------------------------------
 279   Double_t p[2]={c->GetY(), c->GetZ()};
 280   Double_t cov[3]={c->GetSigmaY2(), 0., c->GetSigmaZ2()};
 281   return AliExternalTrackParam::GetPredictedChi2(p,cov);
 282 }
 283
 284 //_____________________________________________________________________________
 285 Bool_t AliTPCtrack::PropagateTo(Double_t xk,Double_t rho,Double_t x0) {
 286   //-----------------------------------------------------------------
 287   //  This function propagates a track to a reference plane x=xk.
 288   //  rho - density of the crossed matrial (g/cm^3)
 289   //  x0  - radiation length of the crossed material (g/cm^2)
 290   //-----------------------------------------------------------------
 291   //
 292   Double_t bz=GetBz();
 293   Double_t zat=0;
 294   if (!GetZAt(xk, bz,zat)) return kFALSE;
 295   if (TMath::Abs(zat)>250.){
 296     // Don't propagate track outside of the fiducial volume - material budget not proper one
 297     //
 298     //AliWarning("Propagate outside of fiducial volume");
 299     return kFALSE;
 300   }
 301
 302   Double_t oldX=GetX(), oldY=GetY(), oldZ=GetZ();
 303   //if (!AliExternalTrackParam::PropagateTo(xk,bz)) return kFALSE;
 304   Double_t b[3]; GetBxByBz(b);
 305   if (!AliExternalTrackParam::PropagateToBxByBz(xk,b)) return kFALSE;
 306
 307   Double_t d = TMath::Sqrt((GetX()-oldX)*(GetX()-oldX) +
 308                            (GetY()-oldY)*(GetY()-oldY) +
 309                            (GetZ()-oldZ)*(GetZ()-oldZ));
 310   if (IsStartedTimeIntegral() && GetX()>oldX) AddTimeStep(d);
 311
 312   if (oldX < xk) d = -d;
 313   if (!AliExternalTrackParam::CorrectForMeanMaterial(d*rho/x0,d*rho,GetMass(),
 314       kFALSE,AliExternalTrackParam::BetheBlochGas)) return kFALSE;
 315
 316   return kTRUE;
 317 }
 318
 319 //_____________________________________________________________________________
 320 Bool_t
 321 AliTPCtrack::PropagateToVertex(const AliESDVertex *v,Double_t rho,Double_t x0)
 322 {
 323   //-----------------------------------------------------------------
 324   // This function propagates tracks to the vertex
 325   // rho - density of the crossed matrial (g/cm3)
 326   // x0  - radiation length of the crossed material (g/cm2)
 327   //-----------------------------------------------------------------
 328   Double_t oldX=GetX(), oldY=GetY(), oldZ=GetZ();
 329
 330   //Double_t bz=GetBz();
 331   //if (!PropagateToDCA(v,bz,kVeryBig)) return kFALSE;
 332   Double_t b[3]; GetBxByBz(b);
 333   if (!PropagateToDCABxByBz(v,b,kVeryBig)) return kFALSE;
 334
 335   Double_t d = TMath::Sqrt((GetX()-oldX)*(GetX()-oldX) +
 336                            (GetY()-oldY)*(GetY()-oldY) +
 337                            (GetZ()-oldZ)*(GetZ()-oldZ));
 338
 339   if (oldX < GetX()) d = -d;
 340   if (!AliExternalTrackParam::CorrectForMeanMaterial(d*rho/x0,d*rho,GetMass(),
 341       kFALSE,AliExternalTrackParam::BetheBlochGas)) return kFALSE;
 342
 343   return kTRUE;
 344 }
 345
 346 //_____________________________________________________________________________
 347 Bool_t AliTPCtrack::Update(const AliCluster *c, Double_t chisq, Int_t index) {
 348   //-----------------------------------------------------------------
 349   // This function associates a cluster with this track.
 350   //-----------------------------------------------------------------
 351   Double_t p[2]={c->GetY(), c->GetZ()};
 352   Double_t cov[3]={c->GetSigmaY2(), 0., c->GetSigmaZ2()};
 353
 354   if (!AliExternalTrackParam::Update(p,cov)) return kFALSE;
 355
 356   AliTracker::FillResiduals(this,p,cov,c->GetVolumeId());
 357
 358   Int_t n=GetNumberOfClusters();
 359   fIndex[n]=index;
 360   SetNumberOfClusters(n+1);
 361   SetChi2(GetChi2()+chisq);
 362
 363   return kTRUE;
 364 }
 365
 366 ////////////////////////////////////////////////////////////////////////
 367 // MI ADDITION
 368
 369 Float_t AliTPCtrack::Density(Int_t row0, Int_t row1)
 370 {
 371   //
 372   // calculate cluster density
 373   Int_t good  = 0;
 374   Int_t found = 0;
 375   //if (row0<fFirstPoint) row0 = fFirstPoint;
 376   if (row1>fLastPoint) row1 = fLastPoint;
 377
 378
 379   for (Int_t i=row0;i<=row1;i++){
 380     //    Int_t index = fClusterIndex[i];
 381     Int_t index = fIndex[i];
 382     if (index!=-1)  good++;
 383     if (index>0)    found++;
 384   }
 385   Float_t density=0;
 386   if (good>0) density = Float_t(found)/Float_t(good);
 387   return density;
 388 }
 389
 390
 391 Float_t AliTPCtrack::Density2(Int_t row0, Int_t row1)
 392 {
 393   //
 394   // calculate cluster density
 395   Int_t good  = 0;
 396   Int_t found = 0;
 397   //
 398   for (Int_t i=row0;i<=row1;i++){
 399     Int_t index = fIndex[i];
 400     if (index!=-1)  good++;
 401     if (index>0)    found++;
 402   }
 403   Float_t density=0;
 404   if (good>0) density = Float_t(found)/Float_t(good);
 405   return density;
 406 }
 407
 408 void  AliTPCtrack::UpdatePoints()
 409 {
 410   //--------------------------------------------------
 411   //calculates first ,amx dens and last points
 412   //--------------------------------------------------
 413   Float_t density[160];
 414   for (Int_t i=0;i<160;i++) density[i]=-1.;
 415   fPoints[0]= 160;
 416   fPoints[1] = -1;
 417   //
 418   Int_t ngood=0;
 419   Int_t undeff=0;
 420   Int_t nall =0;
 421   Int_t range=20;
 422   for (Int_t i=0;i<160;i++){
 423     Int_t last = i-range;
 424     if (nall<range) nall++;
 425     if (last>=0){
 426       if (fIndex[last]>0&& (fIndex[last]&0x8000)==0) ngood--;
 427       if (fIndex[last]==-1) undeff--;
 428     }
 429     if (fIndex[i]>0&& (fIndex[i]&0x8000)==0)   ngood++;
 430     if (fIndex[i]==-1) undeff++;
 431     if (nall==range &&undeff<range/2) density[i-range/2] = Float_t(ngood)/Float_t(nall-undeff);
 432   }
 433   Float_t maxdens=0;
 434   Int_t indexmax =0;
 435   for (Int_t i=0;i<160;i++){
 436     if (density[i]<0) continue;
 437     if (density[i]>maxdens){
 438       maxdens=density[i];
 439       indexmax=i;
 440     }
 441   }
 442   //
 443   //max dens point
 444   fPoints[3] = maxdens;
 445   fPoints[1] = indexmax;
 446   //
 447   // last point
 448   for (Int_t i=indexmax;i<160;i++){
 449     if (density[i]<0) continue;
 450     if (density[i]<maxdens/2.) {
 451       break;
 452     }
 453     fPoints[2]=i;
 454   }
 455   //
 456   // first point
 457   for (Int_t i=indexmax;i>0;i--){
 458     if (density[i]<0) continue;
 459     if (density[i]<maxdens/2.) {
 460       break;
 461     }
 462     fPoints[0]=i;
 463   }
 464   //
 465 }
 466