HLT/TPCLib/AliHLTTPCTrack.cxx

   1 // @(#) $Id$
   2 // Original: AliHLTTrack.cxx,v 1.32 2005/06/14 10:55:21 cvetan
   3
   4 /**************************************************************************
   5  * This file is property of and copyright by the ALICE HLT Project        *
   6  * ALICE Experiment at CERN, All rights reserved.                         *
   7  *                                                                        *
   8  * Primary Authors: Anders Vestbo, Uli Frankenfeld, maintained by         *
   9  *                  Matthias Richter <Matthias.Richter@ift.uib.no>        *
  10  *                  for The ALICE HLT Project.                            *
  11  *                                                                        *
  12  * Permission to use, copy, modify and distribute this software and its   *
  13  * documentation strictly for non-commercial purposes is hereby granted   *
  14  * without fee, provided that the above copyright notice appears in all   *
  15  * copies and that both the copyright notice and this permission notice   *
  16  * appear in the supporting documentation. The authors make no claims     *
  17  * about the suitability of this software for any purpose. It is          *
  18  * provided "as is" without express or implied warranty.                  *
  19  **************************************************************************/
  20
  21 /** @file   AliHLTTPCTrack.cxx
  22     @author Anders Vestbo, Uli Frankenfeld, maintained by Matthias Richter
  23     @date
  24     @brief  HLT TPC track implementation (conformal mapping) */
  25
  26
  27 #include "AliHLTTPCLogging.h"
  28 #include "AliHLTTPCTrack.h"
  29 #include "AliHLTTPCTransform.h"
  30 #include "AliHLTTPCVertex.h"
  31 #include "AliHLTTPCSpacePointData.h"
  32
  33 #if __GNUC__ >= 3
  34 using namespace std;
  35 #endif
  36
  37 ClassImp(AliHLTTPCTrack)
  38
  39
  40 AliHLTTPCTrack::AliHLTTPCTrack()
  41   :
  42   fNHits(0),
  43   fMCid(-1),
  44   fKappa(0),
  45   fRadius(0),
  46   fCenterX(0),
  47   fCenterY(0),
  48   fFromMainVertex(0),
  49   fSector(0),
  50   fQ(0),
  51
  52   fTanl(0),
  53   fPsi(0),
  54   fPt(0),
  55   fLength(0),
  56
  57   fPterr(0),
  58   fPsierr(0),
  59   fZ0err(0),
  60   fY0err(0),
  61   fTanlerr(0),
  62
  63   fPhi0(0),
  64   fR0(0),
  65   fZ0(0),
  66
  67   //  fPoint({0,0,0}),
  68   fPointPsi(0),
  69
  70   fIsPoint(0),
  71   fIsLocal(true),
  72   //  fRowRange({0,0}),
  73
  74   fPID(0),
  75   fId(-1)
  76 {
  77   //Constructor
  78   fRowRange[0]=0;
  79   fRowRange[1]=0;
  80   fPoint[0]=0;
  81   fPoint[1]=0;
  82   fPoint[2]=0;
  83
  84   SetFirstPoint(0,0,0);
  85   SetLastPoint(0,0,0);
  86   memset(fHitNumbers,0,fgkHitArraySize*sizeof(UInt_t));
  87 }
  88
  89 void AliHLTTPCTrack::Copy(AliHLTTPCTrack *tpt)
  90 {
  91   //setter
  92   SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
  93   SetPhi0(tpt->GetPhi0());
  94   SetKappa(tpt->GetKappa());
  95   SetNHits(tpt->GetNHits());
  96   SetFirstPoint(tpt->GetFirstPointX(),tpt->GetFirstPointY(),tpt->GetFirstPointZ());
  97   SetLastPoint(tpt->GetLastPointX(),tpt->GetLastPointY(),tpt->GetLastPointZ());
  98   SetPt(tpt->GetPt());
  99   SetPsi(tpt->GetPsi());
 100   SetTgl(tpt->GetTgl());
 101   SetPterr(tpt->GetPterr());
 102   SetPsierr(tpt->GetPsierr());
 103   SetTglerr(tpt->GetTglerr());
 104   SetZ0err(tpt->GetZ0err());
 105   SetY0err(tpt->GetY0err());
 106   SetCharge(tpt->GetCharge());
 107   SetHits(tpt->GetNHits(),(UInt_t *)tpt->GetHitNumbers());
 108   SetMCid(tpt->GetMCid());
 109   SetPID(tpt->GetPID());
 110   SetSector(tpt->GetSector());
 111   SetId( tpt->GetId());
 112 }
 113
 114 Int_t AliHLTTPCTrack::Compare(const AliHLTTPCTrack *track) const
 115 {
 116   // compare tracks
 117   if(track->GetNHits() < GetNHits()) return 1;
 118   if(track->GetNHits() > GetNHits()) return -1;
 119   return 0;
 120 }
 121
 122 AliHLTTPCTrack::~AliHLTTPCTrack()
 123 {
 124   //Nothing to do
 125 }
 126
 127 Double_t AliHLTTPCTrack::GetP() const
 128 {
 129   // Returns total momentum.
 130   return fabs(GetPt())*sqrt(1. + GetTgl()*GetTgl());
 131 }
 132
 133 Double_t AliHLTTPCTrack::GetPseudoRapidity() const
 134 { //get pseudo rap
 135   return 0.5 * log((GetP() + GetPz()) / (GetP() - GetPz()));
 136 }
 137
 138 Double_t AliHLTTPCTrack::GetRapidity() const
 139 {
 140   //get rap
 141   const Double_t kmpi = 0.13957;
 142   return 0.5 * log((kmpi + GetPz()) / (kmpi - GetPz()));
 143 }
 144
 145 void AliHLTTPCTrack::Rotate(Int_t slice,Bool_t tolocal)
 146 {
 147   //Rotate track to global parameters
 148   //If flag tolocal is set, the track is rotated
 149   //to local coordinates.
 150
 151   Float_t psi[1] = {GetPsi()};
 152   if(!tolocal)
 153     AliHLTTPCTransform::Local2GlobalAngle(psi,slice);
 154   else
 155     AliHLTTPCTransform::Global2LocalAngle(psi,slice);
 156   SetPsi(psi[0]);
 157   Float_t first[3];
 158   first[0] = GetFirstPointX();
 159   first[1] = GetFirstPointY();
 160   first[2] = GetFirstPointZ();
 161   if(!tolocal)
 162     AliHLTTPCTransform::Local2Global(first,slice);
 163   else
 164     AliHLTTPCTransform::Global2LocHLT(first,slice);
 165   //AliHLTTPCTransform::Global2Local(first,slice,kTRUE);
 166
 167   SetFirstPoint(first[0],first[1],first[2]);
 168   Float_t last[3];
 169   last[0] = GetLastPointX();
 170   last[1] = GetLastPointY();
 171   last[2] = GetLastPointZ();
 172   if(!tolocal)
 173     AliHLTTPCTransform::Local2Global(last,slice);
 174   else
 175     AliHLTTPCTransform::Global2LocHLT(last,slice);
 176   //AliHLTTPCTransform::Global2Local(last,slice,kTRUE);
 177   SetLastPoint(last[0],last[1],last[2]);
 178
 179   Float_t center[3] = {GetCenterX(),GetCenterY(),0};
 180   if(!tolocal)
 181     AliHLTTPCTransform::Local2Global(center,slice);
 182   else
 183     AliHLTTPCTransform::Global2LocHLT(center,slice);
 184   //AliHLTTPCTransform::Global2Local(center,slice,kTRUE);
 185   SetCenterX(center[0]);
 186   SetCenterY(center[1]);
 187
 188   SetPhi0(atan2(fFirstPoint[1],fFirstPoint[0]));
 189   SetR0(sqrt(fFirstPoint[0]*fFirstPoint[0]+fFirstPoint[1]*fFirstPoint[1]));
 190
 191   if(!tolocal)
 192     fIsLocal=kFALSE;
 193   else
 194     fIsLocal=kTRUE;
 195 }
 196
 197 void AliHLTTPCTrack::CalculateHelix()
 198 {
 199   // fit assigned clusters to helix
 200   // for straight line fit
 201   if (AliHLTTPCTransform::GetBFieldValue() == 0.0 ){
 202     fRadius = 999999;  //just zero
 203
 204     SetPhi0(atan2(fFirstPoint[1],fFirstPoint[0]));
 205     SetR0(sqrt(fFirstPoint[0]*fFirstPoint[0]+fFirstPoint[1]*fFirstPoint[1]));
 206   }
 207   // for helix fit
 208   else {
 209     //Calculate Radius, CenterX and CenterY from Psi, X0, Y0
 210     fRadius = fPt / (AliHLTTPCTransform::GetBFieldValue());
 211     if(fRadius) fKappa = -fQ*1./fRadius;
 212     else fRadius = 999999;  //just zero
 213     Double_t trackPhi0 = fPsi + fQ * AliHLTTPCTransform::PiHalf();
 214
 215     fCenterX = fFirstPoint[0] - fRadius *  cos(trackPhi0);
 216     fCenterY = fFirstPoint[1] - fRadius *  sin(trackPhi0);
 217
 218     SetPhi0(atan2(fFirstPoint[1],fFirstPoint[0]));
 219     SetR0(sqrt(fFirstPoint[0]*fFirstPoint[0]+fFirstPoint[1]*fFirstPoint[1]));
 220   }
 221 }
 222
 223 Double_t AliHLTTPCTrack::GetCrossingAngle(Int_t padrow,Int_t slice)
 224 {
 225   //Calculate the crossing angle between track and given padrow.
 226   //Take the dot product of the tangent vector of the track, and
 227   //vector perpendicular to the padrow.
 228   //In order to do this, we need the tangent vector to the track at the
 229   //point. This is done by rotating the radius vector by 90 degrees;
 230   //rotation matrix: (  0  1 )
 231   //                 ( -1  0 )
 232
 233   Float_t angle=0;//Angle perpendicular to the padrow in local coordinates
 234   if(slice>=0)//Global coordinates
 235     {
 236       AliHLTTPCTransform::Local2GlobalAngle(&angle,slice);
 237       if(!CalculateReferencePoint(angle,AliHLTTPCTransform::Row2X(padrow)))
 238         cerr<<"AliHLTTPCTrack::GetCrossingAngle : Track does not cross line in slice "<<slice<<" row "<<padrow<<endl;
 239     }
 240   else //should be in local coordinates
 241     {
 242       Float_t xyz[3];
 243       GetCrossingPoint(padrow,xyz);
 244       fPoint[0] = xyz[0];
 245       fPoint[1] = xyz[1];
 246       fPoint[2] = xyz[2];
 247     }
 248
 249   Double_t tangent[2];
 250
 251   tangent[0] = (fPoint[1] - GetCenterY())/GetRadius();
 252   tangent[1] = -1.*(fPoint[0] - GetCenterX())/GetRadius();
 253
 254   Double_t perppadrow[2] = {cos(angle),sin(angle)};
 255   Double_t cosbeta = fabs(tangent[0]*perppadrow[0] + tangent[1]*perppadrow[1]);
 256   if(cosbeta > 1) cosbeta=1;
 257   return acos(cosbeta);
 258 }
 259
 260 Bool_t AliHLTTPCTrack::GetCrossingPoint(Int_t padrow,Float_t *xyz)
 261 {
 262   //Assumes the track is given in local coordinates
 263   if(!IsLocal())
 264     {
 265       cerr<<"GetCrossingPoint: Track is given on global coordinates"<<endl;
 266       return false;
 267     }
 268
 269   Double_t xHit = AliHLTTPCTransform::Row2X(padrow);
 270
 271 //if (xHit < xyz[0]){
 272 //    LOG(AliHLTTPCLog::kError,"AliHLTTPCTRACK::GetCrossingPoint","")<< "Track doesn't cross padrow "
 273 //                              << padrow <<"(x=" << xHit << "). Smallest x=" << xyz[0] << ENDLOG;
 274 //      return false;
 275 //}
 276
 277   // for straight line fit
 278   if (AliHLTTPCTransform::GetBFieldValue() == 0.0 ){
 279
 280     Double_t yHit = GetFirstPointY() + (Double_t) tan( GetPsi() ) * (xHit - GetFirstPointX());
 281
 282     Double_t s = (xHit - GetFirstPointX())*(xHit - GetFirstPointX()) + (yHit - GetFirstPointY())*(yHit - GetFirstPointY());
 283
 284     Double_t zHit = GetFirstPointZ() + s * GetTgl();
 285
 286     xyz[0] = xHit;
 287     xyz[1] = yHit;
 288     xyz[2] = zHit;
 289   }
 290   // for helix fit
 291     else {
 292       xyz[0] = xHit;
 293       Double_t aa = (xHit - GetCenterX())*(xHit - GetCenterX());
 294       Double_t r2 = GetRadius()*GetRadius();
 295       if(aa > r2)
 296         return false;
 297
 298       Double_t aa2 = sqrt(r2 - aa);
 299       Double_t y1 = GetCenterY() + aa2;
 300       Double_t y2 = GetCenterY() - aa2;
 301       xyz[1] = y1;
 302       if(fabs(y2) < fabs(y1)) xyz[1] = y2;
 303
 304       Double_t yHit = xyz[1];
 305       Double_t angle1 = atan2((yHit - GetCenterY()),(xHit - GetCenterX()));
 306       if(angle1 < 0) angle1 += 2.*AliHLTTPCTransform::Pi();
 307       Double_t angle2 = atan2((GetFirstPointY() - GetCenterY()),(GetFirstPointX() - GetCenterX()));
 308       if(angle2 < 0) angle2 += AliHLTTPCTransform::TwoPi();
 309
 310       Double_t diffangle = angle1 - angle2;
 311       diffangle = fmod(diffangle,AliHLTTPCTransform::TwoPi());
 312       if((GetCharge()*diffangle) > 0) diffangle = diffangle - GetCharge()*AliHLTTPCTransform::TwoPi();
 313
 314       Double_t stot = fabs(diffangle)*GetRadius();
 315
 316       Double_t zHit = GetFirstPointZ() + stot*GetTgl();
 317
 318       xyz[2] = zHit;
 319     }
 320
 321   return true;
 322 }
 323
 324 Bool_t AliHLTTPCTrack::CalculateReferencePoint(Double_t angle,Double_t radius)
 325 {
 326   // Global coordinate: crossing point with y = ax+ b;
 327   // a=tan(angle-AliHLTTPCTransform::PiHalf());
 328   //
 329   const Double_t krr=radius; //position of reference plane
 330   const Double_t kxr = cos(angle) * krr;
 331   const Double_t kyr = sin(angle) * krr;
 332
 333   Double_t a = tan(angle-AliHLTTPCTransform::PiHalf());
 334   Double_t b = kyr - a * kxr;
 335
 336   Double_t pp=(fCenterX+a*fCenterY-a*b)/(1+pow(a,2));
 337   Double_t qq=(pow(fCenterX,2)+pow(fCenterY,2)-2*fCenterY*b+pow(b,2)-pow(fRadius,2))/(1+pow(a,2));
 338
 339   Double_t racine = pp*pp-qq;
 340   if(racine<0) return IsPoint(kFALSE);      //no Point
 341
 342   Double_t rootRacine = sqrt(racine);
 343   Double_t x0 = pp+rootRacine;
 344   Double_t x1 = pp-rootRacine;
 345   Double_t y0 = a*x0 + b;
 346   Double_t y1 = a*x1 + b;
 347
 348   Double_t diff0 = sqrt(pow(x0-kxr,2)+pow(y0-kyr,2));
 349   Double_t diff1 = sqrt(pow(x1-kxr,2)+pow(y1-kyr,2));
 350
 351   if(diff0<diff1){
 352     fPoint[0]=x0;
 353     fPoint[1]=y0;
 354   }
 355   else{
 356     fPoint[0]=x1;
 357     fPoint[1]=y1;
 358   }
 359
 360   Double_t pointPhi0  = atan2(fPoint[1]-fCenterY,fPoint[0]-fCenterX);
 361   Double_t trackPhi0  = atan2(fFirstPoint[1]-fCenterY,fFirstPoint[0]-fCenterX);
 362   if(fabs(trackPhi0-pointPhi0)>AliHLTTPCTransform::Pi()){
 363     if(trackPhi0<pointPhi0) trackPhi0 += AliHLTTPCTransform::TwoPi();
 364     else                    pointPhi0 += AliHLTTPCTransform::TwoPi();
 365   }
 366   Double_t stot = -fQ * (pointPhi0-trackPhi0) * fRadius ;
 367   fPoint[2]   = fFirstPoint[2] + stot * fTanl;
 368
 369   fPointPsi = pointPhi0 - fQ * AliHLTTPCTransform::PiHalf();
 370   if(fPointPsi<0.)  fPointPsi+= AliHLTTPCTransform::TwoPi();
 371   fPointPsi = fmod(fPointPsi, AliHLTTPCTransform::TwoPi());
 372
 373   return IsPoint(kTRUE);
 374 }
 375
 376 Bool_t AliHLTTPCTrack::CalculateEdgePoint(Double_t angle)
 377 {
 378   // Global coordinate: crossing point with y = ax; a=tan(angle);
 379   //
 380   Double_t rmin=AliHLTTPCTransform::Row2X(AliHLTTPCTransform::GetFirstRow(-1));  //min Radius of TPC
 381   Double_t rmax=AliHLTTPCTransform::Row2X(AliHLTTPCTransform::GetLastRow(-1)); //max Radius of TPC
 382
 383   Double_t a = tan(angle);
 384   Double_t pp=(fCenterX+a*fCenterY)/(1+pow(a,2));
 385   Double_t qq=(pow(fCenterX,2)+pow(fCenterY,2)-pow(fRadius,2))/(1+pow(a,2));
 386   Double_t racine = pp*pp-qq;
 387   if(racine<0) return IsPoint(kFALSE);     //no Point
 388   Double_t rootRacine = sqrt(racine);
 389   Double_t x0 = pp+rootRacine;
 390   Double_t x1 = pp-rootRacine;
 391   Double_t y0 = a*x0;
 392   Double_t y1 = a*x1;
 393
 394   Double_t r0 = sqrt(pow(x0,2)+pow(y0,2));
 395   Double_t r1 = sqrt(pow(x1,2)+pow(y1,2));
 396   //find the right crossing point:
 397   //inside the TPC modules
 398   Bool_t ok0 = kFALSE;
 399   Bool_t ok1 = kFALSE;
 400
 401   if(r0>rmin&&r0<rmax){
 402     Double_t da=atan2(y0,x0);
 403     if(da<0) da+=AliHLTTPCTransform::TwoPi();
 404     if(fabs(da-angle)<0.5)
 405       ok0 = kTRUE;
 406   }
 407   if(r1>rmin&&r1<rmax){
 408     Double_t da=atan2(y1,x1);
 409     if(da<0) da+=AliHLTTPCTransform::TwoPi();
 410     if(fabs(da-angle)<0.5)
 411       ok1 = kTRUE;
 412   }
 413   if(!(ok0||ok1)) return IsPoint(kFALSE);   //no Point
 414
 415   if(ok0&&ok1){
 416     Double_t diff0 = sqrt(pow(fFirstPoint[0]-x0,2)+pow(fFirstPoint[1]-y0,2));
 417     Double_t diff1 = sqrt(pow(fFirstPoint[0]-x1,2)+pow(fFirstPoint[1]-y1,2));
 418     if(diff0<diff1) ok1 = kFALSE; //use ok0
 419     else ok0 = kFALSE;            //use ok1
 420   }
 421   if(ok0){fPoint[0]=x0; fPoint[1]=y0;}
 422   else   {fPoint[0]=x1; fPoint[1]=y1;}
 423
 424   Double_t pointPhi0  = atan2(fPoint[1]-fCenterY,fPoint[0]-fCenterX);
 425   Double_t trackPhi0  = atan2(fFirstPoint[1]-fCenterY,fFirstPoint[0]-fCenterX);
 426   if(fabs(trackPhi0-pointPhi0)>AliHLTTPCTransform::Pi()){
 427     if(trackPhi0<pointPhi0) trackPhi0 += AliHLTTPCTransform::TwoPi();
 428     else                    pointPhi0 += AliHLTTPCTransform::TwoPi();
 429   }
 430   Double_t stot = -fQ * (pointPhi0-trackPhi0) * fRadius ;
 431   fPoint[2]   = fFirstPoint[2] + stot * fTanl;
 432
 433   fPointPsi = pointPhi0 - fQ * AliHLTTPCTransform::PiHalf();
 434   if(fPointPsi<0.)  fPointPsi+= AliHLTTPCTransform::TwoPi();
 435   fPointPsi = fmod(fPointPsi, AliHLTTPCTransform::TwoPi());
 436
 437   return IsPoint(kTRUE);
 438 }
 439
 440 Bool_t AliHLTTPCTrack::CalculatePoint(Double_t xplane)
 441 {
 442   // Local coordinate: crossing point with x plane
 443   //
 444   Double_t racine = pow(fRadius,2)-pow(xplane-fCenterX,2);
 445   if(racine<0) return IsPoint(kFALSE);
 446   Double_t rootRacine = sqrt(racine);
 447
 448   Double_t y0 = fCenterY + rootRacine;
 449   Double_t y1 = fCenterY - rootRacine;
 450   //Double_t diff0 = sqrt(pow(fFirstPoint[0]-xplane)+pow(fFirstPoint[1]-y0));
 451   //Double_t diff1 = sqrt(pow(fFirstPoint[0]-xplane)+pow(fFirstPoint[1]-y1));
 452   Double_t diff0 = fabs(y0-fFirstPoint[1]);
 453   Double_t diff1 = fabs(y1-fFirstPoint[1]);
 454
 455   fPoint[0]=xplane;
 456   if(diff0<diff1) fPoint[1]=y0;
 457   else            fPoint[1]=y1;
 458
 459   Double_t pointPhi0  = atan2(fPoint[1]-fCenterY,fPoint[0]-fCenterX);
 460   Double_t trackPhi0  = atan2(fFirstPoint[1]-fCenterY,fFirstPoint[0]-fCenterX);
 461   if(fabs(trackPhi0-pointPhi0)>AliHLTTPCTransform::Pi()){
 462     if(trackPhi0<pointPhi0) trackPhi0 += AliHLTTPCTransform::TwoPi();
 463     else                    pointPhi0 += AliHLTTPCTransform::TwoPi();
 464   }
 465   Double_t stot = -fQ * (pointPhi0-trackPhi0) * fRadius ;
 466   fPoint[2]   = fFirstPoint[2] + stot * fTanl;
 467
 468   fPointPsi = pointPhi0 - fQ * AliHLTTPCTransform::PiHalf();
 469   if(fPointPsi<0.)  fPointPsi+= AliHLTTPCTransform::TwoPi();
 470   fPointPsi = fmod(fPointPsi, AliHLTTPCTransform::TwoPi());
 471
 472   return IsPoint(kTRUE);
 473 }
 474
 475 void AliHLTTPCTrack::UpdateToFirstPoint()
 476 {
 477   //Update track parameters to the innermost point on the track.
 478   //This means that the parameters of the track will be given in the point
 479   //of closest approach to the first innermost point, i.e. the point
 480   //lying on the track fit (and not the coordinates of the innermost point itself).
 481   //This function assumes that fFirstPoint is already set to the coordinates of the innermost
 482   //assigned cluster.
 483   //
 484   //During the helix-fit, the first point on the track is set to the coordinates
 485   //of the innermost assigned cluster. This may be ok, if you just want a fast
 486   //estimate of the "global" track parameters; such as the momentum etc.
 487   //However, if you later on want to do more precise local calculations, such
 488   //as impact parameter, residuals etc, you need to give the track parameters
 489   //according to the actual fit.
 490   // for straight line fit
 491   if (AliHLTTPCTransform::GetBFieldValue() == 0.0 ){
 492     Double_t xc = GetCenterX() - GetFirstPointX();
 493     Double_t yc = GetCenterY() - GetFirstPointY();
 494
 495     Double_t xn = (Double_t) sin( GetPsi() );
 496     Double_t yn = -1. * (Double_t) cos( GetPsi() );
 497
 498     Double_t d = xc*xn + yc*yn;
 499
 500     Double_t distx = d * xn;
 501     Double_t disty = d * yn;
 502
 503     Double_t point[2];
 504
 505     point[0] = distx + GetFirstPointX();
 506     point[1] = disty + GetFirstPointY();
 507
 508     //Update the track parameters
 509     SetR0(sqrt(point[0]*point[0]+point[1]*point[1]));
 510     SetPhi0(atan2(point[1],point[0]));
 511     SetFirstPoint(point[0],point[1],GetZ0());
 512   }
 513   // for helix fit
 514   else {
 515     Double_t xc = GetCenterX() - GetFirstPointX();
 516     Double_t yc = GetCenterY() - GetFirstPointY();
 517
 518     Double_t distx1 = xc*(1 + GetRadius()/sqrt(xc*xc + yc*yc));
 519     Double_t disty1 = yc*(1 + GetRadius()/sqrt(xc*xc + yc*yc));
 520     Double_t distance1 = sqrt(distx1*distx1 + disty1*disty1);
 521
 522     Double_t distx2 = xc*(1 - GetRadius()/sqrt(xc*xc + yc*yc));
 523     Double_t disty2 = yc*(1 - GetRadius()/sqrt(xc*xc + yc*yc));
 524     Double_t distance2 = sqrt(distx2*distx2 + disty2*disty2);
 525
 526     //Choose the closest:
 527     Double_t point[2];
 528     if(distance1 < distance2)
 529       {
 530         point[0] = distx1 + GetFirstPointX();
 531         point[1] = disty1 + GetFirstPointY();
 532       }
 533     else
 534       {
 535         point[0] = distx2 + GetFirstPointX();
 536         point[1] = disty2 + GetFirstPointY();
 537       }
 538
 539     Double_t pointpsi = atan2(point[1]-GetCenterY(),point[0]-GetCenterX());
 540     pointpsi -= GetCharge()*AliHLTTPCTransform::PiHalf();
 541     if(pointpsi < 0) pointpsi += AliHLTTPCTransform::TwoPi();
 542
 543     //Update the track parameters
 544     SetR0(sqrt(point[0]*point[0]+point[1]*point[1]));
 545     SetPhi0(atan2(point[1],point[0]));
 546     SetFirstPoint(point[0],point[1],GetZ0());
 547     SetPsi(pointpsi);
 548   }
 549 }
 550
 551 void AliHLTTPCTrack::GetClosestPoint(AliHLTTPCVertex *vertex,Double_t &closestX,Double_t &closestY,Double_t &closestZ)
 552 {
 553   //Calculate the point of closest approach to the vertex
 554   //This function calculates the minimum distance from the helix to the vertex, and choose
 555   //the corresponding point lying on the helix as the point of closest approach.
 556
 557   Double_t xc = GetCenterX() - vertex->GetX();
 558   Double_t yc = GetCenterY() - vertex->GetY();
 559
 560   Double_t distx1 = xc*(1 + GetRadius()/sqrt(xc*xc + yc*yc));
 561   Double_t disty1 = yc*(1 + GetRadius()/sqrt(xc*xc + yc*yc));
 562   Double_t distance1 = sqrt(distx1*distx1 + disty1*disty1);
 563
 564   Double_t distx2 = xc*(1 - GetRadius()/sqrt(xc*xc + yc*yc));
 565   Double_t disty2 = yc*(1 - GetRadius()/sqrt(xc*xc + yc*yc));
 566   Double_t distance2 = sqrt(distx2*distx2 + disty2*disty2);
 567
 568   //Choose the closest:
 569   if(distance1 < distance2)
 570     {
 571       closestX = distx1 + vertex->GetX();
 572       closestY = disty1 + vertex->GetY();
 573     }
 574   else
 575     {
 576       closestX = distx2 + vertex->GetX();
 577       closestY = disty2 + vertex->GetY();
 578     }
 579
 580   //Get the z coordinate:
 581   Double_t angle1 = atan2((closestY-GetCenterY()),(closestX-GetCenterX()));
 582   if(angle1 < 0) angle1 = angle1 + AliHLTTPCTransform::TwoPi();
 583
 584   Double_t angle2 = atan2((GetFirstPointY()-GetCenterY()),(GetFirstPointX()-GetCenterX()));
 585   if(angle2 < 0) angle2 = angle2 + AliHLTTPCTransform::TwoPi();
 586
 587   Double_t diffAngle = angle1 - angle2;
 588   diffAngle = fmod(diffAngle,AliHLTTPCTransform::TwoPi());
 589
 590   if((GetCharge()*diffAngle) < 0) diffAngle = diffAngle + GetCharge()*AliHLTTPCTransform::TwoPi();
 591   Double_t stot = fabs(diffAngle)*GetRadius();
 592   closestZ = GetFirstPointZ() - stot*GetTgl();
 593 }
 594
 595 void AliHLTTPCTrack::Print(Option_t* /*option*/) const
 596 {
 597 //print out parameters of track
 598
 599  LOG(AliHLTTPCLog::kInformational,"AliHLTTPCTrack::Print","Print values")
 600     <<"NH="<<fNHits<<" "<<fMCid<<" K="<<fKappa<<" R="<<fRadius<<" Cx="<<fCenterX<<" Cy="<<fCenterY<<" MVT="
 601     <<fFromMainVertex<<" Row0="<<fRowRange[0]<<" Row1="<<fRowRange[1]<<" Sector="<<fSector<<" Q="<<fQ<<" TgLam="
 602     <<fTanl<<" psi="<<fPsi<<" pt="<<fPt<<" L="<<fLength<<" "<<fPterr<<" "<<fPsierr<<" "<<fZ0err<<" "
 603     <<fTanlerr<<" phi0="<<fPhi0<<" R0="<<fR0<<" Z0="<<fZ0<<" X0="<<fFirstPoint[0]<<" Y0="<<fFirstPoint[1]<<" Z0="
 604     <<fFirstPoint[2]<<" XL="<<fLastPoint[0]<<" YL="<<fLastPoint[1]<<" ZL="<<fLastPoint[2]<<" "
 605     <<fPoint[0]<<" "<<fPoint[1]<<" "<<fPoint[2]<<" "<<fPointPsi<<" "<<fIsPoint<<" local="
 606     <<fIsLocal<<" "<<fPID<<ENDLOG;
 607
 608 }
 609
 610 int AliHLTTPCTrack::Convert2AliKalmanTrack()
 611 {
 612   // The method has been copied from AliHLTHoughKalmanTrack and adapted
 613   // to the TPC conformal mapping track parametrization
 614   int iResult=0;
 615
 616   // sector A00 starts at 3 o'clock, sectors are counted counterclockwise
 617   // median of sector 00 is at 10 degrees, median of sector A04 at 90
 618   //
 619
 620   Double_t charge=(double) GetCharge();
 621   Double_t param[5];
 622   param[1] = GetFirstPointZ();
 623   param[3] = GetTgl();
 624   param[4] = charge*(1.0/GetPt());
 625
 626   Double_t alpha, phi, xl, yl;
 627
 628   // rotate to local coordinates if necessary
 629
 630   if(GetSector() == -1){ // track in global coordinates
 631
 632     alpha = TMath::ATan2(GetFirstPointY(),GetFirstPointX());
 633     double sinAlpha = TMath::Sin(alpha);
 634     double cosAlpha = TMath::Cos(alpha);
 635
 636     phi = GetPsi() - alpha;
 637     xl =  GetFirstPointX()*cosAlpha + GetFirstPointY()*sinAlpha;
 638     yl = -GetFirstPointX()*sinAlpha + GetFirstPointY()*cosAlpha;
 639
 640   } else{ // track in local coordinates
 641
 642     alpha = (GetSector()+0.5)*(TMath::TwoPi()/18);
 643     phi = GetPsi();
 644     xl = GetFirstPointX();
 645     yl = GetFirstPointY();
 646
 647     // normalize alpha to [-Pi,+Pi]
 648
 649     alpha = alpha - TMath::TwoPi() * TMath::Floor( alpha /TMath::TwoPi()+.5);
 650   }
 651
 652   // extra rotation to keep phi in the range (-Pi/2,+Pi/2)
 653   {
 654     const Double_t kMaxPhi = TMath::PiOver2() - 10./180.*TMath::Pi();
 655
 656     // normalize phi to [-Pi,+Pi]
 657
 658     phi = phi - TMath::TwoPi() * TMath::Floor( phi /TMath::TwoPi()+.5);
 659
 660     if( phi >= kMaxPhi )
 661       {
 662         alpha += TMath::PiOver2();
 663         phi   -= TMath::PiOver2();
 664         Double_t xtmp = xl;
 665         xl =  yl;
 666         yl = -xtmp;
 667       }
 668     else if( phi <= -kMaxPhi )
 669       {
 670         alpha += -TMath::PiOver2();
 671         phi   -= -TMath::PiOver2();
 672         Double_t xtmp = xl;
 673         xl = -yl;
 674         yl =  xtmp;
 675       }
 676   }
 677
 678   param[0] = yl;
 679   param[2] = TMath::Sin(phi);
 680
 681   //covariance matrix
 682   Double_t cov[15]={
 683     GetY0err(),                         //Error in Y (Y and X are the same)
 684     0.,  GetZ0err(),                    //Error in Z
 685     0.,  0.,  GetPsierr(),              //Error for Psi
 686     0.,  0.,  0.,  GetTglerr(),         //Error for Tgl
 687     0.,  0.,  0.,  0.,  GetPterr()      //Error for Pt
 688   };
 689
 690   Int_t nCluster = GetNHits();
 691   fdEdx=0;
 692
 693   // the Set function was not available in earlier versions, check done
 694   // during configure; for the AliRoot build, by default ON
 695 #ifdef EXTERNALTRACKPARAM_V1
 696 #warning track conversion to ESD format needs AliRoot version > v4-05-04
 697   //TODO (Feb 07): make this a real warning when logging system is adapted
 698   //HLTWarning("track conversion to ESD format needs AliRoot version > v4-05-04");
 699 #else
 700   Set(xl,alpha,param,cov);
 701   SetNumberOfClusters(nCluster);
 702   SetChi2(0.);
 703   SetFakeRatio(0.);
 704   SetMass(0.13957);
 705 #endif
 706
 707   return iResult;
 708 }
 709
 710 void AliHLTTPCTrack::SetHits(Int_t nhits,UInt_t *hits)
 711 {
 712   // set hit array
 713   if (!hits) return;
 714   if (nhits>fgkHitArraySize) {
 715     LOG(AliHLTTPCLog::kWarning,"AliHLTTPCTrack::SetHits","too many hits")
 716       << "too many hits (" << nhits << ") for hit array of size " << fgkHitArraySize << ENDLOG;
 717     SetNHits(fgkHitArraySize);
 718   } else {
 719     SetNHits(nhits);
 720   }
 721   memcpy(fHitNumbers,hits,fNHits*sizeof(UInt_t));
 722 }
 723
 724 Double_t AliHLTTPCTrack::GetLengthXY() const
 725 {
 726   //calculates the length of the arc in XY-plane. This is the length of the track in XY-plane.
 727   //Using a^2 = b^2 + c^2 - 2bc * cosA for finding the angle between first and last point.
 728   //Length of arc is arc = r*A. Where A is the angle between first and last point.
 729
 730   Double_t dx = GetLastPointX()-GetFirstPointX();
 731   Double_t dy = GetLastPointY()-GetFirstPointY();
 732   Double_t a = TMath::Sqrt((dx*dx)+(dy*dy));
 733   Double_t r = GetRadius();
 734   Double_t r2 = r*r;
 735
 736   Double_t A = TMath::ACos((r2+r2-(a*a))/(2*r2));
 737
 738   return r*A;
 739 }
 740
 741 Double_t AliHLTTPCTrack::GetLengthTot() const
 742 {
 743   //Calculates the length of the track in 3D
 744
 745
 746
 747
 748
 749   return 100.0;
 750
 751 }
 752
 753 int AliHLTTPCTrack::CheckConsistency()
 754 {
 755   // Check consistency of all members
 756   int iResult=0;
 757   if (CheckDoubleMember(&fPterr,   0., "fPterr")<0) iResult=-EDOM;
 758   if (CheckDoubleMember(&fPsierr,  0., "fPsierr")<0) iResult=-EDOM;
 759   if (CheckDoubleMember(&fZ0err,   0., "fZ0err")<0) iResult=-EDOM;
 760   if (CheckDoubleMember(&fY0err,   0., "fY0err")<0) iResult=-EDOM;
 761   if (CheckDoubleMember(&fTanlerr, 0., "fTanlerr")<0) iResult=-EDOM;
 762   return iResult;
 763 }
 764
 765 int AliHLTTPCTrack::CheckDoubleMember(double* pMember, double def, const char* name) const
 766 {
 767   // Check consistency of a Double member
 768   if (!pMember) return -EINVAL;
 769   if (TMath::Abs(*pMember)>kVeryBig) {
 770     LOG(AliHLTTPCLog::kWarning,"AliHLTTPCTrack","member consistency")
 771       << "invalid Double number %f" << *pMember << " in member " << name << ENDLOG;
 772     *pMember=def;
 773     return -EDOM;
 774   }
 775   return 0;
 776 }