HLT/src/AliL3Fitter.cxx

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