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