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