HLT/src/AliL3Fitter.cxx

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