HLT/TPCLib/AliHLTTPCFitter.cxx

   1 // @(#) $Id$
   2 // Original: AliHLTFitter.cxx,v 1.14 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, 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   AliHLTTPCFitter.cxx
  22     @author Anders Vestbo, maintained by Matthias Richter
  23     @date
  24     @brief  Fit class HLT for helix
  25 */
  26
  27 #include <math.h>
  28 #include "AliHLTTPCLogging.h"
  29 #include "AliHLTTPCFitter.h"
  30 #include "AliHLTTPCVertex.h"
  31 #include "AliHLTTPCTrack.h"
  32 #include "AliHLTTPCSpacePointData.h"
  33 #include "AliHLTTPCMemHandler.h"
  34 #include "AliHLTTPCTransform.h"
  35 //#include "AliHLTTPC.h"
  36
  37 #if __GNUC__ >= 3
  38 using namespace std;
  39 #endif
  40
  41 ClassImp(AliHLTTPCFitter)
  42
  43
  44 AliHLTTPCFitter::AliHLTTPCFitter()
  45   :
  46   fTrack(NULL),
  47   fVertex(NULL),
  48   fVertexConstraint(0)
  49 {
  50   //constructor
  51   // see header file for class documentation
  52   // or
  53   // refer to README to build package
  54   // or
  55   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  56   memset(fClusters,0,36*6*sizeof(AliHLTTPCSpacePointData*));
  57   memset(fNcl,0,36*6*sizeof(UInt_t));
  58 }
  59
  60 AliHLTTPCFitter::AliHLTTPCFitter(AliHLTTPCVertex *vertex,Bool_t vertexconstraint)
  61   :
  62   fTrack(NULL),
  63   fVertex(vertex),
  64   fVertexConstraint(vertexconstraint)
  65 {
  66   //constructor
  67   fTrack=0;
  68   memset(fClusters,0,36*6*sizeof(AliHLTTPCSpacePointData*));
  69 }
  70
  71 AliHLTTPCFitter::~AliHLTTPCFitter()
  72 {
  73   //destructor
  74   for(Int_t i=0; i<36; i++)
  75     {
  76       for(Int_t j=0; j<6; j++)
  77         {
  78           if(fClusters[i][j])
  79             delete [] fClusters[i][j];
  80         }
  81     }
  82 }
  83
  84 void AliHLTTPCFitter::LoadClusters(Char_t *path,Int_t event,Bool_t sp)
  85 {
  86   //load clusters
  87   Char_t fname[256];
  88   AliHLTTPCMemHandler *clusterfile[36][6];
  89   for(Int_t s=0; s<=35; s++)
  90     {
  91       for(Int_t p=0; p<6; p++)
  92         {
  93           Int_t patch;
  94           if(sp==kTRUE)
  95             patch=-1;
  96           else
  97             patch=p;
  98           if(fClusters[s][p])
  99             delete fClusters[s][p];
 100           fClusters[s][p] = 0;
 101           clusterfile[s][p] = new AliHLTTPCMemHandler();
 102           sprintf(fname,"%s/points_%d_%d_%d.raw",path,event,s,patch);
 103           if(!clusterfile[s][p]->SetBinaryInput(fname))
 104             {
 105               delete clusterfile[s][p];
 106               clusterfile[s][p] = 0;
 107               continue;
 108             }
 109           fClusters[s][p] = (AliHLTTPCSpacePointData*)clusterfile[s][p]->Allocate();
 110           clusterfile[s][p]->Binary2Memory(fNcl[s][p],fClusters[s][p]);
 111           clusterfile[s][p]->CloseBinaryInput();
 112           if(sp==kTRUE)
 113             break;
 114         }
 115     }
 116 }
 117
 118 void AliHLTTPCFitter::SortTrackClusters(AliHLTTPCTrack *track) const
 119 {
 120   //Sort the internal cluster list in each track with respect to row numbering.
 121   //This may be necessary when no conventional track follower has been
 122   //applied, in which the cluster list has been maintained in a more
 123   //arbitrary fashion.
 124
 125   Int_t nhits = track->GetNHits();
 126   Int_t *ids = (Int_t*)track->GetHitNumbers();
 127   Int_t *origids = new Int_t[nhits];
 128   Int_t *mk = new Int_t[nhits];
 129   Int_t k;
 130
 131   for(k=0; k<nhits; k++) {origids[k] = ids[k]; mk[k] = -1;}
 132
 133   Int_t slice,patch,id,padrow,maxrow,maxk;
 134   UInt_t pos;
 135   for(Int_t j=0; j<nhits; j++)
 136     {
 137       maxrow=-1;
 138       maxk=200;
 139       for(k=0; k<nhits; k++)
 140         {
 141           id=ids[k];
 142           if(id < 0) continue;
 143           slice = AliHLTTPCSpacePointData::GetSlice(id);
 144           patch = AliHLTTPCSpacePointData::GetPatch(id);
 145           pos = AliHLTTPCSpacePointData::GetNumber(id);
 146           AliHLTTPCSpacePointData *points = fClusters[slice][patch];
 147           padrow = points[pos].fPadRow;
 148           if(padrow > maxrow)
 149             {
 150               maxrow = padrow;
 151               maxk=k;
 152             }
 153         }
 154       mk[j]=maxk;
 155       ids[maxk]=-1;
 156     }
 157
 158   for(k=0; k<nhits; k++)
 159     ids[k] = origids[mk[k]];
 160   delete [] origids;
 161   delete [] mk;
 162 }
 163
 164 Int_t AliHLTTPCFitter::FitHelix(AliHLTTPCTrack *track)
 165 {
 166   //fit helix parameters
 167   fTrack = track;
 168   if(FitCircle())
 169     {
 170       LOG(AliHLTTPCLog::kError,"AliHLTTPCFitter::FitHelix","TrackFit")<<AliHLTTPCLog::kDec<<
 171         "Problems during circle fit"<<ENDLOG;
 172       return 1;
 173     }
 174   if(FitLine())
 175     {
 176       LOG(AliHLTTPCLog::kError,"AliHLTTPCFitter::FitHelix","TrackFit")<<AliHLTTPCLog::kDec<<
 177         "Problems during line fit"<<ENDLOG;
 178       return 1;
 179     }
 180   return 0;
 181 }
 182
 183 Int_t AliHLTTPCFitter::FitCircle()
 184 {
 185   //-----------------------------------------------------------------
 186   //Fits circle parameters using algorithm
 187   //described by ChErnov and Oskov in Computer Physics
 188   //Communications.
 189   //
 190   //Written in FORTRAN by Jawluen Tang, Physics department , UT-Austin
 191   //Moved to C by Pablo Yepes
 192   //Moved to AliROOT by ASV.
 193   //------------------------------------------------------------------
 194
 195   Double_t wsum  = 0.0 ;
 196   Double_t xav   = 0.0 ;
 197   Double_t yav   = 0.0 ;
 198
 199   //
 200   //     Loop over hits calculating average
 201   Double_t * fXYWeight = new Double_t[(fTrack->GetNHits())];
 202   UInt_t *hitnum = fTrack->GetHitNumbers();
 203   for(Int_t i=0; i<fTrack->GetNHits(); i++)
 204     {
 205       UInt_t id = hitnum[i];
 206       Int_t slice = AliHLTTPCSpacePointData::GetSlice(id);
 207       Int_t patch = AliHLTTPCSpacePointData::GetPatch(id);
 208       UInt_t pos = AliHLTTPCSpacePointData::GetNumber(id);
 209       AliHLTTPCSpacePointData *points = fClusters[slice][patch];
 210       fXYWeight[i] = 1./ (Double_t)(points[pos].fSigmaY2 + points[pos].fSigmaY2);
 211       wsum += fXYWeight[i];
 212       xav += fXYWeight[i]*points[pos].fX;
 213       yav += fXYWeight[i]*points[pos].fY;
 214     }
 215   if (fVertexConstraint == kTRUE)
 216     {
 217       wsum += fVertex->GetXYWeight() ;
 218       xav  += fVertex->GetX() ;
 219       yav  += fVertex->GetY() ;
 220     }
 221
 222   xav = xav / wsum ;
 223   yav = yav / wsum ;
 224 //
 225 //  CALCULATE <X**2>, <XY>, AND <Y**2> WITH <X> = 0, & <Y> = 0
 226 //
 227   Double_t xxav  = 0.0 ;
 228   Double_t xyav  = 0.0 ;
 229   Double_t yyav  = 0.0 ;
 230   Double_t xi, yi ;
 231
 232   for(Int_t i=0; i<fTrack->GetNHits(); i++)
 233     {
 234       UInt_t id = hitnum[i];
 235       Int_t slice = AliHLTTPCSpacePointData::GetSlice(id);
 236       Int_t patch = AliHLTTPCSpacePointData::GetPatch(id);
 237       UInt_t pos = AliHLTTPCSpacePointData::GetNumber(id);
 238       AliHLTTPCSpacePointData *points = fClusters[slice][patch];
 239
 240       xi = points[pos].fX -xav;
 241       yi        = points[pos].fY - yav ;
 242       xxav     += xi * xi * fXYWeight[i];
 243       xyav     += xi * yi * fXYWeight[i];
 244       yyav     += yi * yi * fXYWeight[i];
 245     }
 246
 247   if (fVertexConstraint == kTRUE)
 248     {
 249       xi        = fVertex->GetX() - xav ;
 250       yi        = fVertex->GetY() - yav ;
 251       xxav     += xi * xi * fVertex->GetXYWeight() ;
 252       xyav     += xi * yi * fVertex->GetXYWeight() ;
 253       yyav     += yi * yi * fVertex->GetXYWeight() ;
 254     }
 255   xxav = xxav / wsum ;
 256   xyav = xyav / wsum ;
 257   yyav = yyav / wsum ;
 258 //
 259 //-->  ROTATE COORDINATES SO THAT <XY> = 0
 260 //
 261 //-->  SIGN(C**2 - S**2) = SIGN(XXAV - YYAV) >
 262 //-->  &                                     > ==> NEW : (XXAV-YYAV) > 0
 263 //-->  SIGN(S) = SIGN(XYAV)                  >
 264
 265   Double_t a = fabs( xxav - yyav ) ;
 266   Double_t b = 4.0 * xyav * xyav ;
 267
 268   Double_t asqpb  = a * a + b  ;
 269   Double_t rasqpb = sqrt ( asqpb) ;
 270
 271   Double_t splus  = 1.0 + a / rasqpb ;
 272   Double_t sminus = b / (asqpb * splus) ;
 273
 274   splus  = sqrt (0.5 * splus ) ;
 275   sminus = sqrt (0.5 * sminus) ;
 276 //
 277 //->  FIRST REQUIRE : SIGN(C**2 - S**2) = SIGN(XXAV - YYAV)
 278 //
 279   Double_t sinrot, cosrot ;
 280   if ( xxav <= yyav ) {
 281          cosrot = sminus ;
 282          sinrot = splus  ;
 283   }
 284   else {
 285           cosrot = splus ;
 286           sinrot = sminus ;
 287   }
 288 //
 289 //->  REQUIRE : SIGN(S) = SIGN(XYAV) * SIGN(C) (ASSUMING SIGN(C) > 0)
 290 //
 291   if ( xyav < 0.0 ) sinrot = - sinrot ;
 292 //
 293 //-->  WE NOW HAVE THE SMALLEST ANGLE THAT GUARANTEES <X**2> > <Y**2>
 294 //-->  TO GET THE SIGN OF THE CHARGE RIGHT, THE NEW X-AXIS MUST POINT
 295 //-->  OUTWARD FROM THE ORGIN.  WE ARE FREE TO CHANGE SIGNS OF BOTH
 296 //-->  COSROT AND SINROT SIMULTANEOUSLY TO ACCOMPLISH THIS.
 297 //
 298 //-->  CHOOSE SIGN OF C WISELY TO BE ABLE TO GET THE SIGN OF THE CHARGE
 299 //
 300   if ( cosrot*xav+sinrot*yav < 0.0 ) {
 301           cosrot = -cosrot ;
 302           sinrot = -sinrot ;
 303   }
 304 //
 305 //->  NOW GET <R**2> AND RSCALE= SQRT(<R**2>)
 306 //
 307   Double_t rrav   = xxav + yyav ;
 308   Double_t rscale = sqrt(rrav) ;
 309
 310   xxav   = 0.0 ;
 311   yyav   = 0.0 ;
 312   xyav   = 0.0 ;
 313   Double_t xrrav = 0.0 ;
 314   Double_t yrrav = 0.0 ;
 315   Double_t rrrrav  = 0.0 ;
 316
 317   Double_t xixi, yiyi, riri, wiriri, xold, yold ;
 318
 319   for(Int_t i=0; i<fTrack->GetNHits(); i++)
 320     {
 321       UInt_t id = hitnum[i];
 322       Int_t slice = AliHLTTPCSpacePointData::GetSlice(id);
 323       Int_t patch = AliHLTTPCSpacePointData::GetPatch(id);
 324       UInt_t pos = AliHLTTPCSpacePointData::GetNumber(id);
 325       AliHLTTPCSpacePointData *points = fClusters[slice][patch];
 326
 327       xold = points[pos].fX - xav ;
 328       yold = points[pos].fY - yav ;
 329       //
 330       //-->  ROTATE SO THAT <XY> = 0 & DIVIDE BY RSCALE SO THAT <R**2> = 1
 331       //
 332       xi = (  cosrot * xold + sinrot * yold ) / rscale ;
 333       yi = ( -sinrot * xold + cosrot * yold ) / rscale ;
 334
 335       xixi   = xi * xi ;
 336       yiyi   = yi * yi ;
 337       riri   = xixi + yiyi ;
 338       wiriri = fXYWeight[i] * riri ;
 339
 340       xyav   += fXYWeight[i] * xi * yi ;
 341       xxav   += fXYWeight[i] * xixi ;
 342       yyav   += fXYWeight[i] * yiyi ;
 343
 344       xrrav  += wiriri * xi ;
 345       yrrav  += wiriri * yi ;
 346       rrrrav += wiriri * riri ;
 347     }
 348 //
 349 //   Include vertex if required
 350 //
 351   if (fVertexConstraint == kTRUE)
 352     {
 353       xold = fVertex->GetX() - xav ;
 354         yold = fVertex->GetY() - yav ;
 355         //
 356         //-->  ROTATE SO THAT <XY> = 0 & DIVIDE BY RSCALE SO THAT <R**2> = 1
 357         //
 358         xi = (  cosrot * xold + sinrot * yold ) / rscale ;
 359         yi = ( -sinrot * xold + cosrot * yold ) / rscale ;
 360
 361         xixi   = xi * xi ;
 362         yiyi   = yi * yi ;
 363         riri   = xixi + yiyi ;
 364         wiriri = fVertex->GetXYWeight() * riri ;
 365
 366         xyav   += fVertex->GetXYWeight() * xi * yi ;
 367         xxav   += fVertex->GetXYWeight() * xixi ;
 368         yyav   += fVertex->GetXYWeight() * yiyi ;
 369
 370         xrrav  += wiriri * xi ;
 371         yrrav  += wiriri * yi ;
 372         rrrrav += wiriri * riri ;
 373   }
 374   //
 375   //
 376   //
 377   //-->  DIVIDE BY WSUM TO MAKE AVERAGES
 378   //
 379   xxav    = xxav   / wsum ;
 380   yyav    = yyav   / wsum ;
 381   xrrav   = xrrav  / wsum ;
 382   yrrav   = yrrav  / wsum ;
 383   rrrrav  = rrrrav / wsum ;
 384   xyav    = xyav   / wsum ;
 385
 386   Int_t const kntry = 5 ;
 387 //
 388 //-->  USE THESE TO GET THE COEFFICIENTS OF THE 4-TH ORDER POLYNIMIAL
 389 //-->  DON'T PANIC - THE THIRD ORDER TERM IS ZERO !
 390 //
 391   Double_t xrrxrr = xrrav * xrrav ;
 392   Double_t yrryrr = yrrav * yrrav ;
 393   Double_t rrrrm1 = rrrrav - 1.0  ;
 394   Double_t xxyy   = xxav  * yyav  ;
 395
 396   Double_t c0  =          rrrrm1*xxyy - xrrxrr*yyav - yrryrr*xxav ;
 397   Double_t c1  =        - rrrrm1      + xrrxrr      + yrryrr   - 4.0*xxyy ;
 398   Double_t c2  =   4.0  + rrrrm1                               - 4.0*xxyy ;
 399   Double_t c4  = - 4.0  ;
 400 //
 401 //-->  COEFFICIENTS OF THE DERIVATIVE - USED IN NEWTON-RAPHSON ITERATIONS
 402 //
 403   Double_t c2d =   2.0 * c2 ;
 404   Double_t c4d =   4.0 * c4 ;
 405 //
 406 //-->  0'TH VALUE OF LAMDA - LINEAR INTERPOLATION BETWEEN P(0) & P(YYAV)
 407 //
 408 //   LAMDA = YYAV * C0 / (C0 + YRRSQ * (XXAV-YYAV))
 409   Double_t lamda  = 0.0 ;
 410   Double_t dlamda = 0.0 ;
 411 //
 412   Double_t chiscl = wsum * rscale * rscale ;
 413   Double_t dlamax = 0.001 / chiscl ;
 414
 415   Double_t p, pd ;
 416   for ( int itry = 1 ; itry <= kntry ; itry++ ) {
 417      p      = c0 + lamda * (c1 + lamda * (c2 + lamda * lamda * c4 )) ;
 418      pd     = (c1 + lamda * (c2d + lamda * lamda * c4d)) ;
 419      dlamda = -p / pd ;
 420      lamda  = lamda + dlamda ;
 421      if (fabs(dlamda)<   dlamax) break ;
 422   }
 423
 424   //Double_t chi2 = (Double_t)(chiscl * lamda) ;
 425   //fTrack->SetChiSq1(chi2);
 426   // Double_t dchisq = chiscl * dlamda ;
 427   //
 428   //-->  NOW CALCULATE THE MATRIX ELEMENTS FOR ALPHA, BETA & KAPPA
 429   //
 430   Double_t h11   = xxav  -     lamda ;
 431   Double_t h14   = xrrav ;
 432   Double_t h22   = yyav  -     lamda ;
 433   Double_t h24   = yrrav ;
 434   Double_t h34   = 1.0   + 2.0*lamda ;
 435   if ( h11 == 0.0 || h22 == 0.0 ){
 436     LOG(AliHLTTPCLog::kError,"AliHLTTPCFitter::FitCircle","TrackFit")<<AliHLTTPCLog::kDec<<
 437       "Problems fitting circle"<<ENDLOG;
 438     return 1 ;
 439   }
 440   Double_t rootsq = (h14*h14)/(h11*h11) + 4.0*h34 ;
 441
 442   Double_t ratio, kappa, beta ;
 443   if ( fabs(h22) > fabs(h24) ) {
 444      ratio  = h24 / h22 ;
 445      rootsq = ratio * ratio + rootsq ;
 446      kappa = 1.0 / sqrt(rootsq) ;
 447      beta  = - ratio * kappa ;
 448   }
 449   else {
 450      ratio  = h22 / h24 ;
 451      rootsq = 1.0 + ratio * ratio * rootsq ;
 452      beta  = 1.0 / sqrt(rootsq) ;
 453      if ( h24 > 0 ) beta = - beta ;
 454      kappa = -ratio * beta ;
 455   }
 456   Double_t alpha = - (h14/h11) * kappa ;
 457 //
 458 //-->  transform these into the lab coordinate system
 459 //-->  first get kappa and back to real dimensions
 460 //
 461   Double_t kappa1 = kappa / rscale ;
 462   Double_t dbro   = 0.5   / kappa1 ;
 463 //
 464 //-->  next rotate alpha and beta and scale
 465 //
 466   Double_t alphar = (cosrot * alpha - sinrot * beta)* dbro ;
 467   Double_t betar  = (sinrot * alpha + cosrot * beta)* dbro ;
 468 //
 469 //-->  then translate by (xav,yav)
 470 //
 471   Double_t acent  = (double)(xav - alphar) ;
 472   Double_t bcent  = (double)(yav - betar ) ;
 473   Double_t radius = (double)dbro ;
 474 //
 475 //   Get charge
 476 //
 477   Int_t q = ( ( yrrav < 0 ) ? 1 : -1 ) ;
 478   fTrack->SetCharge(q);
 479
 480   //Set the first point on the track to the space point coordinates of the innermost track
 481   //This will be updated to lie on the fit later on (AliHLTTPCTrack::UpdateToFirstPoint).
 482   Double_t x0,y0,psi,pt ;
 483   Int_t lastid=fTrack->GetNHits()-1;
 484   UInt_t id = hitnum[lastid];
 485   Int_t slice = AliHLTTPCSpacePointData::GetSlice(id);
 486   Int_t patch = AliHLTTPCSpacePointData::GetPatch(id);
 487   UInt_t pos = AliHLTTPCSpacePointData::GetNumber(id);
 488   AliHLTTPCSpacePointData *points = fClusters[slice][patch];
 489   x0   =  points[pos].fX;
 490   y0   =  points[pos].fY;
 491   fTrack->SetFirstPoint(x0,y0,0); //Z-value is set in FitLine
 492
 493   //Set the remaining fit parameters
 494   psi  = (Double_t)atan2(bcent-y0,acent-x0) ;
 495   psi  = psi + q * 0.5F * AliHLTTPCTransform::Pi() ;
 496   if ( psi < 0 ) psi = psi + 2*AliHLTTPCTransform::Pi();
 497
 498   pt   = (Double_t)(AliHLTTPCTransform::GetBFact() * AliHLTTPCTransform::GetBField() * radius ) ;
 499   fTrack->SetPsi(psi);
 500   fTrack->SetPt(pt);
 501   fTrack->SetRadius(radius);
 502   fTrack->SetCenterX(acent);
 503   fTrack->SetCenterY(bcent);
 504 //
 505 //    Get errors from fast fit
 506 //
 507   //if ( getPara()->getErrors ) getErrorsCircleFit ( acent, bcent, radius ) ;
 508 //
 509   delete [] fXYWeight;
 510   return 0 ;
 511 }
 512
 513 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 514 //    Fit Line in s-z plane
 515 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 516 Int_t AliHLTTPCFitter::FitLine ( )
 517 {
 518   //
 519   //Initialization
 520   //
 521   Double_t sum = 0.F ;
 522   Double_t ss  = 0.F ;
 523   Double_t sz  = 0.F ;
 524   Double_t sss = 0.F ;
 525   Double_t ssz = 0.F ;
 526   //
 527   //find sum , sums ,sumz, sumss
 528   //
 529   Double_t dx, dy ;
 530   Double_t radius = (Double_t)(fTrack->GetPt() / ( AliHLTTPCTransform::GetBFact() * AliHLTTPCTransform::GetBField() ) ) ;
 531
 532   Double_t * fS = new Double_t[(fTrack->GetNHits())];
 533   Double_t *fZWeight = new Double_t[fTrack->GetNHits()];
 534   UInt_t *hitnum = fTrack->GetHitNumbers();
 535   if (0)//fVertexConstraint==kTRUE)
 536     {
 537       UInt_t id = hitnum[0];
 538       Int_t slice = AliHLTTPCSpacePointData::GetSlice(id);
 539       Int_t patch = AliHLTTPCSpacePointData::GetPatch(id);
 540       UInt_t pos = AliHLTTPCSpacePointData::GetNumber(id);
 541       AliHLTTPCSpacePointData *points = fClusters[slice][patch];
 542
 543       dx = points[pos].fX - fVertex->GetX();
 544       dy = points[pos].fY - fVertex->GetY();
 545     }
 546   else
 547     {
 548       UInt_t id = hitnum[0];
 549       Int_t slice = AliHLTTPCSpacePointData::GetSlice(id);
 550       Int_t patch = AliHLTTPCSpacePointData::GetPatch(id);
 551       UInt_t posf = AliHLTTPCSpacePointData::GetNumber(id);
 552       AliHLTTPCSpacePointData *pointsf = fClusters[slice][patch];
 553       id = hitnum[(fTrack->GetNHits()-1)];
 554       slice = AliHLTTPCSpacePointData::GetSlice(id);
 555       patch = AliHLTTPCSpacePointData::GetPatch(id);
 556       UInt_t posl = AliHLTTPCSpacePointData::GetNumber(id);
 557       AliHLTTPCSpacePointData *pointsl = fClusters[slice][patch];
 558       dx = pointsf[posf].fX - pointsl[posl].fX;
 559       dy = pointsf[posf].fY - pointsl[posl].fY;
 560     }
 561
 562   Double_t localPsi = 0.5F * sqrt ( dx*dx + dy*dy ) / radius ;
 563   Double_t totals ;
 564
 565   if ( fabs(localPsi) < 1. )
 566     {
 567       totals = 2.0 * radius * asin ( localPsi ) ;
 568     }
 569   else
 570     {
 571       totals = 2.0 * radius * AliHLTTPCTransform::Pi() ;
 572     }
 573
 574   Double_t dpsi,s;
 575
 576   for(Int_t i=0; i<fTrack->GetNHits(); i++)
 577     {
 578       UInt_t id = hitnum[i];
 579       Int_t slice = AliHLTTPCSpacePointData::GetSlice(id);
 580       Int_t patch = AliHLTTPCSpacePointData::GetPatch(id);
 581       UInt_t pos = AliHLTTPCSpacePointData::GetNumber(id);
 582       AliHLTTPCSpacePointData *points = fClusters[slice][patch];
 583
 584       fZWeight[i] = 1./(Double_t)(points[pos].fSigmaZ2);
 585       if(i>0)
 586         {
 587           id = hitnum[i-1];
 588           slice = AliHLTTPCSpacePointData::GetSlice(id);
 589           patch = AliHLTTPCSpacePointData::GetPatch(id);
 590           UInt_t lastpos = AliHLTTPCSpacePointData::GetNumber(id);
 591           AliHLTTPCSpacePointData *lastpoints = fClusters[slice][patch];
 592           dx = points[pos].fX -lastpoints[lastpos].fX;
 593           dy = points[pos].fY -lastpoints[lastpos].fY;
 594           dpsi = 0.5 * (Double_t)sqrt ( dx*dx + dy*dy ) / radius ;
 595           if(fabs(dpsi) > 1)
 596             return 1;
 597           fTrack->SetPsierr(dpsi);
 598           s = fS[i-1] - 2.0 * radius * (Double_t)asin ( dpsi ) ;
 599           fS[i]=s;
 600         }
 601       else
 602         fS[i]=totals;
 603
 604       sum += fZWeight[i];
 605       ss  += fZWeight[i] * fS[i];
 606       sz  += fZWeight[i] * points[pos].fZ;
 607       sss += fZWeight[i] * fS[i] * fS[i];
 608       ssz += fZWeight[i] * fS[i] * points[pos].fZ;
 609
 610     }
 611
 612
 613   Double_t chi2,det = sum * sss - ss * ss;
 614   if ( fabs(det) < 1e-20)
 615     {
 616       chi2 = 99999.F ;
 617       //fTrack->SetChiSq2(chi2);
 618       return 0 ;
 619     }
 620
 621   //Compute the best fitted parameters A,B
 622   Double_t tanl,z0,dtanl,dz0;
 623
 624   tanl = (Double_t)((sum * ssz - ss * sz ) / det );
 625   z0   = (Double_t)((sz * sss - ssz * ss ) / det );
 626
 627   fTrack->SetTgl(tanl);
 628   fTrack->SetZ0(z0);
 629
 630   //calculate chi-square
 631   chi2 = 0.;
 632   Double_t r1 ;
 633
 634   for(Int_t i=0; i<fTrack->GetNHits(); i++)
 635     {
 636       UInt_t id = hitnum[i];
 637       Int_t slice = AliHLTTPCSpacePointData::GetSlice(id);
 638       Int_t patch = AliHLTTPCSpacePointData::GetPatch(id);
 639       UInt_t pos = AliHLTTPCSpacePointData::GetNumber(id);
 640       AliHLTTPCSpacePointData *points = fClusters[slice][patch];
 641       r1   = points[pos].fZ - tanl * fS[i] - z0 ;
 642       chi2 += (Double_t) ( (Double_t)(fZWeight[i]) * (r1 * r1) );
 643     }
 644
 645   //fTrack->SetChiSq2(chi2);
 646   //
 647   //calculate estimated variance
 648   //varsq=chi/(double(n)-2.)
 649   //calculate covariance matrix
 650   //siga=sqrt(varsq*sxx/det)
 651   //sigb=sqrt(varsq*sum/det)
 652   //
 653   dtanl = (Double_t) ( sum / det );
 654   dz0   = (Double_t) ( sss / det );
 655
 656   fTrack->SetTglerr(dtanl);
 657   fTrack->SetZ0err(dz0);
 658   delete [] fZWeight;
 659   delete [] fS;
 660   return 0 ;
 661 }