HMPID/AliHMPIDtrack.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 #include "AliESDtrack.h"
  17 #include "AliTracker.h"
  18 #include "AliHMPIDtrack.h"
  19
  20 ClassImp(AliHMPIDtrack)
  21
  22 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  23 AliHMPIDtrack::AliHMPIDtrack():AliKalmanTrack()
  24 {
  25   //
  26   // def. ctor
  27   //
  28 }
  29 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  30 AliHMPIDtrack::AliHMPIDtrack(const AliHMPIDtrack& t):AliKalmanTrack(t)
  31 {
  32   //
  33   // cctor.
  34   //
  35  }
  36 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  37 AliHMPIDtrack::AliHMPIDtrack(const AliESDtrack& t):AliKalmanTrack()
  38 {
  39   //
  40   // Constructor from AliESDtrack
  41   //
  42   SetLabel(t.GetLabel());
  43   SetChi2(0.);
  44   SetMass(t.GetMass());
  45
  46   Set(t.GetX(),t.GetAlpha(),t.GetParameter(),t.GetCovariance());
  47
  48   if ((t.GetStatus()&AliESDtrack::kTIME) == 0) return;
  49   StartTimeIntegral();
  50   Double_t times[10]; t.GetIntegratedTimes(times); SetIntegratedTimes(times);
  51   SetIntegratedLength(t.GetIntegratedLength());
  52 }
  53 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  54 AliHMPIDtrack& AliHMPIDtrack::operator=(const AliHMPIDtrack &/*source*/)
  55 {
  56   // ass. op.
  57
  58   return *this;
  59 }
  60 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  61 Bool_t AliHMPIDtrack::PropagateTo(Double_t xk, Double_t xx0, Double_t xrho)
  62 {
  63   //
  64   // Propagates this track to a reference plane defined by "xk" [cm]
  65   // correcting for the mean crossed material.
  66   // Arguments:
  67   // "xx0"  - thickness/rad.length [units of the radiation length]
  68   // "xrho" - thickness*density    [g/cm^2]
  69   //  Returns: kTRUE if the track propagates to plane, else kFALSE
  70
  71   if (xk == GetX()) {
  72     return kTRUE;
  73   }
  74   Double_t bz   = GetBz();
  75   if (!AliExternalTrackParam::PropagateTo(xk,bz)) {
  76     return kFALSE;
  77   }
  78   if (!AliExternalTrackParam::CorrectForMeanMaterial(xx0,xrho,GetMass())) {
  79     return kFALSE;
  80   }
  81   return kTRUE;
  82 }//PropagateTo()
  83 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  84 Bool_t AliHMPIDtrack::Rotate(Double_t alpha, Bool_t absolute)
  85 {
  86   //
  87   // Rotates track parameters in R*phi plane
  88   // if absolute rotation alpha is in global system
  89   // otherwise alpha rotation is relative to the current rotation angle
  90   //
  91
  92   if (absolute) {
  93     alpha -= GetAlpha();
  94   }
  95
  96
  97   return AliExternalTrackParam::Rotate(GetAlpha()+alpha);
  98
  99 }
 100 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 101 Int_t AliHMPIDtrack::GetProlongation(Double_t xk, Double_t &y, Double_t &z)
 102 {
 103   //
 104   // Find a prolongation at given x
 105   // Return 0 if it does not exist
 106   //
 107
 108   Double_t bz = GetBz();
 109
 110   if (!AliExternalTrackParam::GetYAt(xk,bz,y)) {
 111     return 0;
 112   }
 113   if (!AliExternalTrackParam::GetZAt(xk,bz,z)) {
 114     return 0;
 115   }
 116
 117   return 1;
 118
 119 }
 120
 121 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 122 Bool_t   AliHMPIDtrack::PropagateToR(Double_t r,Double_t step)
 123 {
 124   //
 125   // Propagate track to the radial position
 126   // Rotation always connected to the last track position
 127   //
 128
 129   Double_t xyz0[3];
 130   Double_t xyz1[3];
 131   Double_t y;
 132   Double_t z;
 133
 134   Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
 135   // Direction +-
 136   Double_t dir    = (radius > r) ? -1.0 : 1.0;
 137
 138   for (Double_t x = radius+dir*step; dir*x < dir*r; x += dir*step) {
 139
 140     GetXYZ(xyz0);
 141     Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
 142     Rotate(alpha,kTRUE);
 143     GetXYZ(xyz0);
 144     if (!GetProlongation(x,y,z)) return kFALSE;
 145     xyz1[0] = x * TMath::Cos(alpha) + y * TMath::Sin(alpha);
 146     xyz1[1] = x * TMath::Sin(alpha) - y * TMath::Cos(alpha);
 147     xyz1[2] = z;
 148     Double_t param[7];
 149     AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
 150     if (param[1] <= 0) {
 151       param[1] = 100000000;
 152     }
 153     PropagateTo(x,param[1],param[0]*param[4]);
 154
 155   }
 156
 157   GetXYZ(xyz0);
 158   Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
 159   Rotate(alpha,kTRUE);
 160   GetXYZ(xyz0);
 161   if (!GetProlongation(r,y,z)) return kFALSE;
 162   xyz1[0] = r * TMath::Cos(alpha) + y * TMath::Sin(alpha);
 163   xyz1[1] = r * TMath::Sin(alpha) - y * TMath::Cos(alpha);
 164   xyz1[2] = z;
 165   Double_t param[7];
 166   AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
 167
 168   if (param[1] <= 0) {
 169     param[1] = 100000000;
 170   }
 171
 172   return PropagateTo(r,param[1],param[0]*param[4]);
 173 }
 174 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 175 Double_t AliHMPIDtrack::GetPredictedChi2(const AliCluster3D *c) const {
 176   //
 177   // Arguments: AliCluster3D
 178   // Returns:   Chi2 of track for the cluster
 179   Double_t      p[3]={c->GetX(),       c->GetY(),       c->GetZ()};
 180   Double_t  covyz[3]={c->GetSigmaY2(), c->GetSigmaYZ(), c->GetSigmaZ2()};
 181   Double_t covxyz[3]={c->GetSigmaX2(), c->GetSigmaXY(), c->GetSigmaXZ()};
 182   return AliExternalTrackParam::GetPredictedChi2(p, covyz, covxyz);
 183 }//GetPredictedChi2()
 184 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 185 Bool_t AliHMPIDtrack::PropagateTo(const AliCluster3D *c) {
 186   //
 187   // Arguments: AliCluster3D
 188   // Returns: kTRUE if the track propagates to the plane of the cluster
 189   Double_t      oldX=GetX(),   oldY=GetY(), oldZ=GetZ();
 190   Double_t      p[3]={c->GetX(), c->GetY(), c->GetZ()};
 191   Double_t  covyz[3]={c->GetSigmaY2(), c->GetSigmaYZ(), c->GetSigmaZ2()};
 192   Double_t covxyz[3]={c->GetSigmaX2(), c->GetSigmaXY(), c->GetSigmaXZ()};
 193   Double_t bz=GetBz();
 194
 195   if(!AliExternalTrackParam::PropagateTo(p, covyz, covxyz, bz)) return kFALSE;
 196   if(IsStartedTimeIntegral())
 197     {
 198       Double_t d = TMath::Sqrt((GetX()-oldX)*(GetX()-oldX) + (GetY()-oldY)*(GetY()-oldY) + (GetZ()-oldZ)*(GetZ()-oldZ));
 199       if (GetX()<oldX) d=-d;
 200       AddTimeStep(d);
 201   }
 202   return kTRUE;
 203 }//PropagateTo()
 204 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 205 Double_t AliHMPIDtrack::GetBz() const {
 206   //
 207   // Arguments: none
 208   // Returns: Bz component of the magnetic field (kG)
 209   //
 210   if (AliTracker::UniformField()) return AliTracker::GetBz();
 211   Double_t r[3]; GetXYZ(r);
 212   return AliTracker::GetBz(r);
 213 }//GetBz()
 214 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 215 Bool_t AliHMPIDtrack::Intersect(Double_t pnt[3], Double_t norm[3], Double_t bz) const {
 216   //+++++++++++++++++++++++++++++++++++++++++
 217   // Origin: K. Shileev (Kirill.Shileev@cern.ch)
 218   // Finds point of intersection (if exists) of the helix with the plane.
 219   // Stores result in fX and fP.
 220   // Arguments: planePoint,planeNorm - the plane defined by any plane's point
 221   // and vector, normal to the plane
 222   // Returns: kTrue if helix intersects the plane, kFALSE otherwise.
 223   //+++++++++++++++++++++++++++++++++++++++++
 224   Double_t x0[3]; GetXYZ(x0); //get track position in MARS
 225
 226   //estimates initial helix length up to plane
 227   Double_t s=(pnt[0]-x0[0])*norm[0] + (pnt[1]-x0[1])*norm[1] + (pnt[2]-x0[2])*norm[2];
 228   Double_t dist=99999,distPrev=dist;
 229   Double_t x[3],p[3];
 230   while(TMath::Abs(dist)>0.00001){
 231     //calculates helix at the distance s from x0 ALONG the helix
 232     Propagate(s,x,p,bz);
 233     //distance between current helix position and plane
 234     dist=(x[0]-pnt[0])*norm[0]+(x[1]-pnt[1])*norm[1]+(x[2]-pnt[2])*norm[2];
 235     if(TMath::Abs(dist) >= TMath::Abs(distPrev)) {return kFALSE;}
 236     distPrev=dist;
 237     s-=dist;
 238   }
 239   //on exit pnt is intersection point,norm is track vector at that point,
 240   //all in MARS
 241   for (Int_t i=0; i<3; i++) {pnt[i]=x[i]; norm[i]=p[i];}
 242   return kTRUE;
 243 }//Intersect()
 244 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 245 Bool_t AliHMPIDtrack::Intersect(AliHMPIDtrack *pTrk,Double_t pnt[3], Double_t norm[3]) {
 246   // Finds point of intersection (if exists) of the helix with the plane.
 247   // Stores result in fX and fP.
 248   // Arguments: planePoint,planeNorm - the plane defined by any plane's point
 249   // and vector, normal to the plane
 250   // Returns: kTrue if helix intersects the plane, kFALSE otherwise.
 251
 252   Double_t x0[3]; pTrk->GetXYZ(x0); //get track position in MARS
 253   Double_t dist=99999,distPrev=dist;
 254   Double_t x[3],p[3],
 255   pntrad= TMath::Sqrt(pnt[0]*pnt[0]+pnt[1]*pnt[1]);
 256   while(TMath::Abs(dist)> 0.000001){//0.00001){
 257     //calculates helix at the distance s from x0 ALONG the helix
 258     pTrk->PropagateTo(pntrad);pTrk->GetXYZ(x);pTrk->GetPxPyPz(p);
 259     //distance between current helix position and plane
 260     dist=(x[0]-pnt[0])*norm[0]+(x[1]-pnt[1])*norm[1]+(x[2]-pnt[2])*norm[2];
 261     pntrad=pntrad-dist*0.7;
 262     //Printf("--- 111 --- dist %lf",dist);
 263     if(TMath::Abs(2.0*dist) >= TMath::Abs(distPrev)) {return kFALSE;}
 264     distPrev=dist;
 265   }
 266   //on exit pnt is intersection point,norm is track vector at that point,
 267   //Printf("--- 222 --- dist %lf",dist);
 268   Printf("");
 269   for (Int_t i=0; i<3; i++) {pnt[i]=x[i]; norm[i]=p[i];}
 270   return kTRUE;
 271 }//Intersect()
 272 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 273 void AliHMPIDtrack::Propagate(Double_t len, Double_t x[3],Double_t p[3], Double_t bz) const {
 274   //+++++++++++++++++++++++++++++++++++++++++
 275   // Origin: K. Shileev (Kirill.Shileev@cern.ch)
 276   // Extrapolate track along simple helix in magnetic field
 277   // Arguments: len -distance alogn helix, [cm]
 278   //            bz  - mag field, [kGaus]
 279   // Returns: x and p contain extrapolated positon and momentum
 280   // The momentum returned for straight-line tracks is meaningless !
 281   //+++++++++++++++++++++++++++++++++++++++++
 282   GetXYZ(x);
 283   if (OneOverPt() < kAlmost0 || TMath::Abs(bz) < kAlmost0Field ){ //straight-line tracks
 284      Double_t unit[3]; GetDirection(unit);
 285      x[0]+=unit[0]*len;
 286      x[1]+=unit[1]*len;
 287      x[2]+=unit[2]*len;
 288
 289      p[0]=unit[0]/kAlmost0;
 290      p[1]=unit[1]/kAlmost0;
 291      p[2]=unit[2]/kAlmost0;
 292   } else {
 293      GetPxPyPz(p);
 294      Double_t pp=GetP();
 295      Double_t a = -kB2C*bz*GetSign(); ////////// what is kB2C
 296      Double_t rho = a/pp;
 297      x[0] += p[0]*TMath::Sin(rho*len)/a - p[1]*(1-TMath::Cos(rho*len))/a;
 298      x[1] += p[1]*TMath::Sin(rho*len)/a + p[0]*(1-TMath::Cos(rho*len))/a;
 299      x[2] += p[2]*len/pp;
 300      Double_t p0=p[0];
 301      p[0] = p0  *TMath::Cos(rho*len) - p[1]*TMath::Sin(rho*len);
 302      p[1] = p[1]*TMath::Cos(rho*len) + p0  *TMath::Sin(rho*len);
 303   }
 304 }//Propagate()
 305 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 306 Bool_t AliHMPIDtrack::Update(const AliHMPIDCluster *pClu, Double_t /*chisq*/, Int_t /*index*/)
 307 {
 308   //
 309   // Arguments: AliCluster3D, chi sq, and clu index
 310   // Returns: kTRUE if the track parameters are successfully updated
 311   Double_t p[2]={pClu->GetY(), pClu->GetZ()};
 312   Double_t cov[3]={pClu->GetSigmaY2(), 0., pClu->GetSigmaZ2()};
 313   if (!AliExternalTrackParam::Update(p,cov)) return kFALSE;
 314
 315   /*
 316   AliTracker::FillResiduals(this,p,cov,pClu->GetVolumeId());
 317
 318   Int_t n=GetNumberOfClusters();
 319   fIndex[n]=index;
 320   SetNumberOfClusters(n+1);
 321   SetChi2(GetChi2()+chisq);
 322 */
 323   return kTRUE;
 324
 325 }//Update()
 326 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++