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