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