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