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 "AliTracker.h"
  28 #include "AliESDtrack.h"
  29 #include "AliITStrackV2.h"
  30
  31 const Int_t AliITStrackV2::fgkWARN = 5;
  32
  33 ClassImp(AliITStrackV2)
  34
  35
  36 //____________________________________________________________________________
  37 AliITStrackV2::AliITStrackV2() : AliKalmanTrack(),
  38   fdEdx(0),
  39   fESDtrack(0)
  40 {
  41   for(Int_t i=0; i<2*AliITSgeomTGeo::kNLayers; i++) {fIndex[i]=-1; fModule[i]=-1;}
  42   for(Int_t i=0; i<4; i++) fdEdxSample[i]=0;
  43 }
  44
  45
  46 //____________________________________________________________________________
  47 AliITStrackV2::AliITStrackV2(AliESDtrack& t,Bool_t c) throw (const Char_t *) :
  48   AliKalmanTrack(),
  49   fdEdx(t.GetITSsignal()),
  50   fESDtrack(&t)
  51 {
  52   //------------------------------------------------------------------
  53   // Conversion ESD track -> ITS track.
  54   // If c==kTRUE, create the ITS track out of the constrained params.
  55   //------------------------------------------------------------------
  56   const AliExternalTrackParam *par=&t;
  57   if (c) {
  58     par=t.GetConstrainedParam();
  59     if (!par) throw "AliITStrackV2: conversion failed !\n";
  60   }
  61   Set(par->GetX(),par->GetAlpha(),par->GetParameter(),par->GetCovariance());
  62
  63   //if (!Invariant()) throw "AliITStrackV2: conversion failed !\n";
  64
  65   SetLabel(t.GetLabel());
  66   SetMass(t.GetMass());
  67   SetNumberOfClusters(t.GetITSclusters(fIndex));
  68
  69   if (t.GetStatus()&AliESDtrack::kTIME) {
  70     StartTimeIntegral();
  71     Double_t times[10]; t.GetIntegratedTimes(times); SetIntegratedTimes(times);
  72     SetIntegratedLength(t.GetIntegratedLength());
  73   }
  74
  75   for(Int_t i=0; i<4; i++) fdEdxSample[i]=0;
  76 }
  77
  78 void AliITStrackV2::UpdateESDtrack(ULong_t flags) const {
  79   fESDtrack->UpdateTrackParams(this,flags);
  80   // copy the module indices
  81   for(Int_t i=0;i<12;i++) {
  82     //   printf("     %d\n",GetModuleIndex(i));
  83     fESDtrack->SetITSModuleIndex(i,GetModuleIndex(i));
  84   }
  85 }
  86
  87 //____________________________________________________________________________
  88 AliITStrackV2::AliITStrackV2(const AliITStrackV2& t) :
  89   AliKalmanTrack(t),
  90   fdEdx(t.fdEdx),
  91   fESDtrack(t.fESDtrack)
  92 {
  93   //------------------------------------------------------------------
  94   //Copy constructor
  95   //------------------------------------------------------------------
  96   Int_t i;
  97   for (i=0; i<4; i++) fdEdxSample[i]=t.fdEdxSample[i];
  98   for (i=0; i<2*AliITSgeomTGeo::GetNLayers(); i++) {
  99     fIndex[i]=t.fIndex[i];
 100     fModule[i]=t.fModule[i];
 101   }
 102 }
 103
 104 //_____________________________________________________________________________
 105 Int_t AliITStrackV2::Compare(const TObject *o) const {
 106   //-----------------------------------------------------------------
 107   // This function compares tracks according to the their curvature
 108   //-----------------------------------------------------------------
 109   AliITStrackV2 *t=(AliITStrackV2*)o;
 110   //Double_t co=OneOverPt();
 111   //Double_t c =OneOverPt();
 112   Double_t co=t->GetSigmaY2()*t->GetSigmaZ2();
 113   Double_t c =GetSigmaY2()*GetSigmaZ2();
 114   if (c>co) return 1;
 115   else if (c<co) return -1;
 116   return 0;
 117 }
 118
 119 //____________________________________________________________________________
 120 Bool_t
 121 AliITStrackV2::PropagateToVertex(const AliESDVertex *v,Double_t d,Double_t x0)
 122 {
 123   //------------------------------------------------------------------
 124   //This function propagates a track to the minimal distance from the origin
 125   //------------------------------------------------------------------
 126   Double_t bz=GetBz();
 127   if (PropagateToDCA(v,bz,kVeryBig)) {
 128     Double_t xOverX0,xTimesRho;
 129     xOverX0 = d; xTimesRho = d*x0;
 130     if (CorrectForMeanMaterial(xOverX0,xTimesRho,kTRUE)) return kTRUE;
 131   }
 132   return kFALSE;
 133 }
 134
 135 //____________________________________________________________________________
 136 Bool_t AliITStrackV2::
 137 GetGlobalXYZat(Double_t xloc, Double_t &x, Double_t &y, Double_t &z) const {
 138   //------------------------------------------------------------------
 139   //This function returns a track position in the global system
 140   //------------------------------------------------------------------
 141   Double_t r[3];
 142   Bool_t rc=GetXYZAt(xloc, AliTracker::GetBz(), r);
 143   x=r[0]; y=r[1]; z=r[2];
 144   return rc;
 145 }
 146
 147 //_____________________________________________________________________________
 148 Double_t AliITStrackV2::GetPredictedChi2(const AliCluster *c) const {
 149   //-----------------------------------------------------------------
 150   // This function calculates a predicted chi2 increment.
 151   //-----------------------------------------------------------------
 152   Double_t p[2]={c->GetY(), c->GetZ()};
 153   Double_t cov[3]={c->GetSigmaY2(), 0., c->GetSigmaZ2()};
 154   return AliExternalTrackParam::GetPredictedChi2(p,cov);
 155 }
 156
 157 //____________________________________________________________________________
 158 Bool_t AliITStrackV2::PropagateTo(Double_t xk, Double_t d, Double_t x0) {
 159   //------------------------------------------------------------------
 160   //This function propagates a track
 161   //------------------------------------------------------------------
 162
 163   Double_t oldX=GetX(), oldY=GetY(), oldZ=GetZ();
 164
 165   Double_t bz=GetBz();
 166   if (!AliExternalTrackParam::PropagateTo(xk,bz)) return kFALSE;
 167   Double_t xOverX0,xTimesRho;
 168   xOverX0 = d; xTimesRho = d*x0;
 169   if (!CorrectForMeanMaterial(xOverX0,xTimesRho,kTRUE)) return kFALSE;
 170
 171   Double_t x=GetX(), y=GetY(), z=GetZ();
 172   if (IsStartedTimeIntegral() && x>oldX) {
 173     Double_t l2 = (x-oldX)*(x-oldX) + (y-oldY)*(y-oldY) + (z-oldZ)*(z-oldZ);
 174     AddTimeStep(TMath::Sqrt(l2));
 175   }
 176
 177   return kTRUE;
 178 }
 179
 180 //____________________________________________________________________________
 181 Bool_t AliITStrackV2::PropagateToTGeo(Double_t xToGo, Int_t nstep, Double_t &xOverX0, Double_t &xTimesRho) {
 182   //-------------------------------------------------------------------
 183   //  Propagates the track to a reference plane x=xToGo in n steps.
 184   //  These n steps are only used to take into account the curvature.
 185   //  The material is calculated with TGeo. (L.Gaudichet)
 186   //-------------------------------------------------------------------
 187
 188   Double_t startx = GetX(), starty = GetY(), startz = GetZ();
 189   Double_t sign = (startx<xToGo) ? -1.:1.;
 190   Double_t step = (xToGo-startx)/TMath::Abs(nstep);
 191
 192   Double_t start[3], end[3], mparam[7], bz = GetBz();
 193   Double_t x = startx;
 194
 195   for (Int_t i=0; i<nstep; i++) {
 196
 197     GetXYZ(start);   //starting global position
 198     x += step;
 199     if (!GetXYZAt(x, bz, end)) return kFALSE;
 200     if (!AliExternalTrackParam::PropagateTo(x, bz)) return kFALSE;
 201     AliTracker::MeanMaterialBudget(start, end, mparam);
 202     if (mparam[1]<900000) {
 203       xTimesRho = sign*mparam[4]*mparam[0];
 204       xOverX0   = mparam[1];
 205       if (!AliExternalTrackParam::CorrectForMeanMaterial(xOverX0,
 206                                       xTimesRho,GetMass())) return kFALSE;
 207     }
 208   }
 209
 210   if (IsStartedTimeIntegral() && GetX()>startx) {
 211     Double_t l2 = ( (GetX()-startx)*(GetX()-startx) +
 212                     (GetY()-starty)*(GetY()-starty) +
 213                     (GetZ()-startz)*(GetZ()-startz) );
 214     AddTimeStep(TMath::Sqrt(l2));
 215   }
 216
 217   return kTRUE;
 218 }
 219
 220 //____________________________________________________________________________
 221 Bool_t AliITStrackV2::Update(const AliCluster* c, Double_t chi2, Int_t index)
 222 {
 223   //------------------------------------------------------------------
 224   //This function updates track parameters
 225   //------------------------------------------------------------------
 226   Double_t p[2]={c->GetY(), c->GetZ()};
 227   Double_t cov[3]={c->GetSigmaY2(), 0., c->GetSigmaZ2()};
 228
 229   if (!AliExternalTrackParam::Update(p,cov)) return kFALSE;
 230
 231   if (!Invariant()) {
 232      AliWarning("Wrong invariant !");
 233      return kFALSE;
 234   }
 235
 236   if (chi2<0) return kTRUE;
 237
 238   Int_t n=GetNumberOfClusters();
 239   fIndex[n]=index;
 240   SetNumberOfClusters(n+1);
 241   SetChi2(GetChi2()+chi2);
 242
 243   return kTRUE;
 244 }
 245
 246 Bool_t AliITStrackV2::Invariant() const {
 247   //------------------------------------------------------------------
 248   // This function is for debugging purpose only
 249   //------------------------------------------------------------------
 250   Int_t n=GetNumberOfClusters();
 251
 252   Double_t sP2=GetParameter()[2];
 253   if (TMath::Abs(sP2) >= kAlmost1){
 254      if (n>fgkWARN) Warning("Invariant","fP2=%f\n",sP2);
 255      return kFALSE;
 256   }
 257   Double_t sC00=GetCovariance()[0];
 258   if (sC00<=0 || sC00>9.) {
 259      if (n>fgkWARN) Warning("Invariant","fC00=%f\n",sC00);
 260      return kFALSE;
 261   }
 262   Double_t sC11=GetCovariance()[2];
 263   if (sC11<=0 || sC11>9.) {
 264      if (n>fgkWARN) Warning("Invariant","fC11=%f\n",sC11);
 265      return kFALSE;
 266   }
 267   Double_t sC22=GetCovariance()[5];
 268   if (sC22<=0 || sC22>1.) {
 269      if (n>fgkWARN) Warning("Invariant","fC22=%f\n",sC22);
 270      return kFALSE;
 271   }
 272   Double_t sC33=GetCovariance()[9];
 273   if (sC33<=0 || sC33>1.) {
 274      if (n>fgkWARN) Warning("Invariant","fC33=%f\n",sC33);
 275      return kFALSE;
 276   }
 277   Double_t sC44=GetCovariance()[14];
 278   if (sC44<=0 /*|| sC44>6e-5*/) {
 279      if (n>fgkWARN) Warning("Invariant","fC44=%f\n",sC44);
 280      return kFALSE;
 281   }
 282
 283   return kTRUE;
 284 }
 285
 286 //____________________________________________________________________________
 287 Bool_t AliITStrackV2::Propagate(Double_t alp,Double_t xk) {
 288   //------------------------------------------------------------------
 289   //This function propagates a track
 290   //------------------------------------------------------------------
 291   Double_t bz=GetBz();
 292   if (!AliExternalTrackParam::Propagate(alp,xk,bz)) return kFALSE;
 293
 294   if (!Invariant()) {
 295      AliWarning("Wrong invariant !");
 296      return kFALSE;
 297   }
 298
 299   return kTRUE;
 300 }
 301
 302 Bool_t AliITStrackV2::MeanBudgetToPrimVertex(Double_t xyz[3], Double_t step, Double_t &d) const {
 303
 304   //-------------------------------------------------------------------
 305   //  Get the mean material budget between the actual point and the
 306   //  primary vertex. (L.Gaudichet)
 307   //-------------------------------------------------------------------
 308
 309   Double_t cs=TMath::Cos(GetAlpha()), sn=TMath::Sin(GetAlpha());
 310   Double_t vertexX = xyz[0]*cs + xyz[1]*sn;
 311
 312   Int_t nstep = Int_t((GetX()-vertexX)/step);
 313   if (nstep<1) nstep = 1;
 314   step = (GetX()-vertexX)/nstep;
 315
 316   //  Double_t mparam[7], densMean=0, radLength=0, length=0;
 317   Double_t mparam[7];
 318   Double_t p1[3], p2[3], x = GetX(), bz = GetBz();
 319   GetXYZ(p1);
 320
 321   d=0.;
 322
 323   for (Int_t i=0; i<nstep; i++) {
 324     x  += step;
 325     if (!GetXYZAt(x, bz, p2)) return kFALSE;
 326     AliTracker::MeanMaterialBudget(p1, p2, mparam);
 327     if (mparam[1]>900000) return kFALSE;
 328     d  += mparam[1];
 329
 330     p1[0] = p2[0];
 331     p1[1] = p2[1];
 332     p1[2] = p2[2];
 333   }
 334
 335   return kTRUE;
 336 }
 337
 338 Bool_t AliITStrackV2::Improve(Double_t x0,Double_t xyz[3],Double_t ers[3]) {
 339   //------------------------------------------------------------------
 340   //This function improves angular track parameters
 341   //------------------------------------------------------------------
 342   Double_t cs=TMath::Cos(GetAlpha()), sn=TMath::Sin(GetAlpha());
 343 //Double_t xv = xyz[0]*cs + xyz[1]*sn; // vertex
 344   Double_t yv =-xyz[0]*sn + xyz[1]*cs; // in the
 345   Double_t zv = xyz[2];                // local frame
 346
 347   Double_t dy = Par(0) - yv, dz = Par(1) - zv;
 348   Double_t r2=GetX()*GetX() + dy*dy;
 349   Double_t p2=(1.+ GetTgl()*GetTgl())/(GetSigned1Pt()*GetSigned1Pt());
 350   Double_t beta2=p2/(p2 + GetMass()*GetMass());
 351   x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
 352   Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
 353   //Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*x0*9.36*2.33;
 354
 355   Double_t cnv=GetBz()*kB2C;
 356   {
 357     Double_t dummy = 4/r2 - GetC()*GetC();
 358     if (dummy < 0) return kFALSE;
 359     Double_t parp = 0.5*(GetC()*GetX() + dy*TMath::Sqrt(dummy));
 360     Double_t sigma2p = theta2*(1.- GetSnp()*GetSnp())*(1. + GetTgl()*GetTgl());
 361     sigma2p += Cov(0)/r2*(1.- dy*dy/r2)*(1.- dy*dy/r2);
 362     sigma2p += ers[1]*ers[1]/r2;
 363     sigma2p += 0.25*Cov(14)*cnv*cnv*GetX()*GetX();
 364     Double_t eps2p=sigma2p/(Cov(5) + sigma2p);
 365     Par(0) += Cov(3)/(Cov(5) + sigma2p)*(parp - GetSnp());
 366     Par(2)  = eps2p*GetSnp() + (1 - eps2p)*parp;
 367     Cov(5) *= eps2p;
 368     Cov(3) *= eps2p;
 369   }
 370   {
 371     Double_t parl=0.5*GetC()*dz/TMath::ASin(0.5*GetC()*TMath::Sqrt(r2));
 372     Double_t sigma2l=theta2;
 373     sigma2l += Cov(2)/r2 + Cov(0)*dy*dy*dz*dz/(r2*r2*r2);
 374     sigma2l += ers[2]*ers[2]/r2;
 375     Double_t eps2l = sigma2l/(Cov(9) + sigma2l);
 376     Par(1) += Cov(7 )/(Cov(9) + sigma2l)*(parl - Par(3));
 377     Par(4) += Cov(13)/(Cov(9) + sigma2l)*(parl - Par(3));
 378     Par(3)  = eps2l*Par(3) + (1-eps2l)*parl;
 379     Cov(9) *= eps2l;
 380     Cov(13)*= eps2l;
 381     Cov(7) *= eps2l;
 382   }
 383   if (!Invariant()) return kFALSE;
 384
 385   return kTRUE;
 386 }
 387
 388 void AliITStrackV2::CookdEdx(Double_t low, Double_t up) {
 389   //-----------------------------------------------------------------
 390   // This function calculates dE/dX within the "low" and "up" cuts.
 391   // Origin: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
 392   //-----------------------------------------------------------------
 393   // The clusters order is: SSD-2, SSD-1, SDD-2, SDD-1, SPD-2, SPD-1
 394
 395   Int_t i;
 396   Int_t nc=0;
 397   for (i=0; i<GetNumberOfClusters(); i++) {
 398     Int_t idx=GetClusterIndex(i);
 399     idx=(idx&0xf0000000)>>28;
 400     if (idx>1) nc++; // Take only SSD and SDD
 401   }
 402
 403   Int_t swap;//stupid sorting
 404   do {
 405     swap=0;
 406     for (i=0; i<nc-1; i++) {
 407       if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
 408       Float_t tmp=fdEdxSample[i];
 409       fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
 410       swap++;
 411     }
 412   } while (swap);
 413
 414   Int_t nl=Int_t(low*nc), nu=Int_t(up*nc); //b.b. to take two lowest dEdX
 415                                            // values from four ones choose
 416                                            // nu=2
 417   Float_t dedx=0;
 418   for (i=nl; i<nu; i++) dedx += fdEdxSample[i];
 419   if (nu-nl>0) dedx /= (nu-nl);
 420
 421   SetdEdx(dedx);
 422 }
 423
 424 Double_t AliITStrackV2::GetBz() const {
 425   //
 426   // returns Bz component of the magnetic field (kG)
 427   //
 428   if (AliTracker::UniformField()) return AliTracker::GetBz();
 429   Double_t r[3]; GetXYZ(r);
 430   return AliTracker::GetBz(r);
 431 }
 432
 433 //____________________________________________________________________________
 434 Bool_t AliITStrackV2::
 435 GetPhiZat(Double_t r, Double_t &phi, Double_t &z) const {
 436   //------------------------------------------------------------------
 437   // This function returns the global cylindrical (phi,z) of the track
 438   // position estimated at the radius r.
 439   // The track curvature is neglected.
 440   //------------------------------------------------------------------
 441   Double_t d=GetD(0.,0.);
 442   if (TMath::Abs(d) > r) return kFALSE;
 443
 444   Double_t rcurr=TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
 445   if (TMath::Abs(d) > rcurr) return kFALSE;
 446   Double_t phicurr=GetAlpha()+TMath::ASin(GetSnp());
 447
 448   phi=phicurr+TMath::ASin(d/r)-TMath::ASin(d/rcurr);
 449   z=GetZ()+GetTgl()*(TMath::Sqrt(r*r-d*d) - TMath::Sqrt(rcurr*rcurr-d*d));
 450
 451   return kTRUE;
 452 }
 453 //____________________________________________________________________________
 454 Bool_t AliITStrackV2::
 455 GetLocalXat(Double_t r,Double_t &xloc) const {
 456   //------------------------------------------------------------------
 457   // This function returns the local x of the track
 458   // position estimated at the radius r.
 459   // The track curvature is neglected.
 460   //------------------------------------------------------------------
 461   Double_t d=GetD(0.,0.);
 462   if (TMath::Abs(d) > r) return kFALSE;
 463
 464   Double_t rcurr=TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
 465   Double_t phicurr=GetAlpha()+TMath::ASin(GetSnp());
 466   Double_t phi=phicurr+TMath::ASin(d/r)-TMath::ASin(d/rcurr);
 467
 468   xloc=r*(TMath::Cos(phi)*TMath::Cos(GetAlpha())
 469          +TMath::Sin(phi)*TMath::Sin(GetAlpha()));
 470
 471   return kTRUE;
 472 }
 473