HLT/TPCLib/AliHLTTPCFitter.cxx

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