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   AliTracker::FillResiduals(this,p,cov,c->GetVolumeId());
 239
 240   Int_t n=GetNumberOfClusters();
 241   fIndex[n]=index;
 242   SetNumberOfClusters(n+1);
 243   SetChi2(GetChi2()+chi2);
 244
 245   return kTRUE;
 246 }
 247
 248 Bool_t AliITStrackV2::Invariant() const {
 249   //------------------------------------------------------------------
 250   // This function is for debugging purpose only
 251   //------------------------------------------------------------------
 252   Int_t n=GetNumberOfClusters();
 253
 254   Double_t sP2=GetParameter()[2];
 255   if (TMath::Abs(sP2) >= kAlmost1){
 256      if (n>fgkWARN) Warning("Invariant","fP2=%f\n",sP2);
 257      return kFALSE;
 258   }
 259   Double_t sC00=GetCovariance()[0];
 260   if (sC00<=0 || sC00>9.) {
 261      if (n>fgkWARN) Warning("Invariant","fC00=%f\n",sC00);
 262      return kFALSE;
 263   }
 264   Double_t sC11=GetCovariance()[2];
 265   if (sC11<=0 || sC11>9.) {
 266      if (n>fgkWARN) Warning("Invariant","fC11=%f\n",sC11);
 267      return kFALSE;
 268   }
 269   Double_t sC22=GetCovariance()[5];
 270   if (sC22<=0 || sC22>1.) {
 271      if (n>fgkWARN) Warning("Invariant","fC22=%f\n",sC22);
 272      return kFALSE;
 273   }
 274   Double_t sC33=GetCovariance()[9];
 275   if (sC33<=0 || sC33>1.) {
 276      if (n>fgkWARN) Warning("Invariant","fC33=%f\n",sC33);
 277      return kFALSE;
 278   }
 279   Double_t sC44=GetCovariance()[14];
 280   if (sC44<=0 /*|| sC44>6e-5*/) {
 281      if (n>fgkWARN) Warning("Invariant","fC44=%f\n",sC44);
 282      return kFALSE;
 283   }
 284
 285   return kTRUE;
 286 }
 287
 288 //____________________________________________________________________________
 289 Bool_t AliITStrackV2::Propagate(Double_t alp,Double_t xk) {
 290   //------------------------------------------------------------------
 291   //This function propagates a track
 292   //------------------------------------------------------------------
 293   Double_t bz=GetBz();
 294   if (!AliExternalTrackParam::Propagate(alp,xk,bz)) return kFALSE;
 295
 296   if (!Invariant()) {
 297      AliWarning("Wrong invariant !");
 298      return kFALSE;
 299   }
 300
 301   return kTRUE;
 302 }
 303
 304 Bool_t AliITStrackV2::MeanBudgetToPrimVertex(Double_t xyz[3], Double_t step, Double_t &d) const {
 305
 306   //-------------------------------------------------------------------
 307   //  Get the mean material budget between the actual point and the
 308   //  primary vertex. (L.Gaudichet)
 309   //-------------------------------------------------------------------
 310
 311   Double_t cs=TMath::Cos(GetAlpha()), sn=TMath::Sin(GetAlpha());
 312   Double_t vertexX = xyz[0]*cs + xyz[1]*sn;
 313
 314   Int_t nstep = Int_t((GetX()-vertexX)/step);
 315   if (nstep<1) nstep = 1;
 316   step = (GetX()-vertexX)/nstep;
 317
 318   //  Double_t mparam[7], densMean=0, radLength=0, length=0;
 319   Double_t mparam[7];
 320   Double_t p1[3], p2[3], x = GetX(), bz = GetBz();
 321   GetXYZ(p1);
 322
 323   d=0.;
 324
 325   for (Int_t i=0; i<nstep; i++) {
 326     x  += step;
 327     if (!GetXYZAt(x, bz, p2)) return kFALSE;
 328     AliTracker::MeanMaterialBudget(p1, p2, mparam);
 329     if (mparam[1]>900000) return kFALSE;
 330     d  += mparam[1];
 331
 332     p1[0] = p2[0];
 333     p1[1] = p2[1];
 334     p1[2] = p2[2];
 335   }
 336
 337   return kTRUE;
 338 }
 339
 340 Bool_t AliITStrackV2::Improve(Double_t x0,Double_t xyz[3],Double_t ers[3]) {
 341   //------------------------------------------------------------------
 342   //This function improves angular track parameters
 343   //------------------------------------------------------------------
 344   Double_t cs=TMath::Cos(GetAlpha()), sn=TMath::Sin(GetAlpha());
 345 //Double_t xv = xyz[0]*cs + xyz[1]*sn; // vertex
 346   Double_t yv =-xyz[0]*sn + xyz[1]*cs; // in the
 347   Double_t zv = xyz[2];                // local frame
 348
 349   Double_t dy = Par(0) - yv, dz = Par(1) - zv;
 350   Double_t r2=GetX()*GetX() + dy*dy;
 351   Double_t p2=(1.+ GetTgl()*GetTgl())/(GetSigned1Pt()*GetSigned1Pt());
 352   Double_t beta2=p2/(p2 + GetMass()*GetMass());
 353   x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
 354   Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
 355   //Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*x0*9.36*2.33;
 356
 357   Double_t cnv=GetBz()*kB2C;
 358   {
 359     Double_t dummy = 4/r2 - GetC()*GetC();
 360     if (dummy < 0) return kFALSE;
 361     Double_t parp = 0.5*(GetC()*GetX() + dy*TMath::Sqrt(dummy));
 362     Double_t sigma2p = theta2*(1.- GetSnp()*GetSnp())*(1. + GetTgl()*GetTgl());
 363     sigma2p += Cov(0)/r2*(1.- dy*dy/r2)*(1.- dy*dy/r2);
 364     sigma2p += ers[1]*ers[1]/r2;
 365     sigma2p += 0.25*Cov(14)*cnv*cnv*GetX()*GetX();
 366     Double_t eps2p=sigma2p/(Cov(5) + sigma2p);
 367     Par(0) += Cov(3)/(Cov(5) + sigma2p)*(parp - GetSnp());
 368     Par(2)  = eps2p*GetSnp() + (1 - eps2p)*parp;
 369     Cov(5) *= eps2p;
 370     Cov(3) *= eps2p;
 371   }
 372   {
 373     Double_t parl=0.5*GetC()*dz/TMath::ASin(0.5*GetC()*TMath::Sqrt(r2));
 374     Double_t sigma2l=theta2;
 375     sigma2l += Cov(2)/r2 + Cov(0)*dy*dy*dz*dz/(r2*r2*r2);
 376     sigma2l += ers[2]*ers[2]/r2;
 377     Double_t eps2l = sigma2l/(Cov(9) + sigma2l);
 378     Par(1) += Cov(7 )/(Cov(9) + sigma2l)*(parl - Par(3));
 379     Par(4) += Cov(13)/(Cov(9) + sigma2l)*(parl - Par(3));
 380     Par(3)  = eps2l*Par(3) + (1-eps2l)*parl;
 381     Cov(9) *= eps2l;
 382     Cov(13)*= eps2l;
 383     Cov(7) *= eps2l;
 384   }
 385   if (!Invariant()) return kFALSE;
 386
 387   return kTRUE;
 388 }
 389
 390 void AliITStrackV2::CookdEdx(Double_t low, Double_t up) {
 391   //-----------------------------------------------------------------
 392   // This function calculates dE/dX within the "low" and "up" cuts.
 393   // Origin: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
 394   //-----------------------------------------------------------------
 395   // The clusters order is: SSD-2, SSD-1, SDD-2, SDD-1, SPD-2, SPD-1
 396
 397   Int_t i;
 398   Int_t nc=0;
 399   for (i=0; i<GetNumberOfClusters(); i++) {
 400     Int_t idx=GetClusterIndex(i);
 401     idx=(idx&0xf0000000)>>28;
 402     if (idx>1) nc++; // Take only SSD and SDD
 403   }
 404
 405   Int_t swap;//stupid sorting
 406   do {
 407     swap=0;
 408     for (i=0; i<nc-1; i++) {
 409       if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
 410       Float_t tmp=fdEdxSample[i];
 411       fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
 412       swap++;
 413     }
 414   } while (swap);
 415
 416   Int_t nl=Int_t(low*nc), nu=Int_t(up*nc); //b.b. to take two lowest dEdX
 417                                            // values from four ones choose
 418                                            // nu=2
 419   Float_t dedx=0;
 420   for (i=nl; i<nu; i++) dedx += fdEdxSample[i];
 421   if (nu-nl>0) dedx /= (nu-nl);
 422
 423   SetdEdx(dedx);
 424 }
 425
 426 Double_t AliITStrackV2::GetBz() const {
 427   //
 428   // returns Bz component of the magnetic field (kG)
 429   //
 430   if (AliTracker::UniformField()) return AliTracker::GetBz();
 431   Double_t r[3]; GetXYZ(r);
 432   return AliTracker::GetBz(r);
 433 }
 434
 435 //____________________________________________________________________________
 436 Bool_t AliITStrackV2::
 437 GetPhiZat(Double_t r, Double_t &phi, Double_t &z) const {
 438   //------------------------------------------------------------------
 439   // This function returns the global cylindrical (phi,z) of the track
 440   // position estimated at the radius r.
 441   // The track curvature is neglected.
 442   //------------------------------------------------------------------
 443   Double_t d=GetD(0.,0.);
 444   if (TMath::Abs(d) > r) return kFALSE;
 445
 446   Double_t rcurr=TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
 447   if (TMath::Abs(d) > rcurr) return kFALSE;
 448   Double_t phicurr=GetAlpha()+TMath::ASin(GetSnp());
 449
 450   phi=phicurr+TMath::ASin(d/r)-TMath::ASin(d/rcurr);
 451   z=GetZ()+GetTgl()*(TMath::Sqrt((r-d)*(r+d))-TMath::Sqrt((rcurr-d)*(rcurr+d)));
 452   return kTRUE;
 453 }
 454 //____________________________________________________________________________
 455 Bool_t AliITStrackV2::
 456 GetLocalXat(Double_t r,Double_t &xloc) const {
 457   //------------------------------------------------------------------
 458   // This function returns the local x of the track
 459   // position estimated at the radius r.
 460   // The track curvature is neglected.
 461   //------------------------------------------------------------------
 462   Double_t d=GetD(0.,0.);
 463   if (TMath::Abs(d) > r) return kFALSE;
 464
 465   Double_t rcurr=TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
 466   Double_t phicurr=GetAlpha()+TMath::ASin(GetSnp());
 467   Double_t phi=phicurr+TMath::ASin(d/r)-TMath::ASin(d/rcurr);
 468
 469   xloc=r*(TMath::Cos(phi)*TMath::Cos(GetAlpha())
 470          +TMath::Sin(phi)*TMath::Sin(GetAlpha()));
 471
 472   return kTRUE;
 473 }
 474