L3 becomes HLT
[u/mrichter/AliRoot.git] / HLT / TPCLib / AliHLTTPCFitter.cxx
1 // @(#) $Id$
2 // Original: AliHLTFitter.cxx,v 1.14 2005/06/14 10:55:21 cvetan 
3
4 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
5 //*-- Copyright &copy ALICE HLT Group 
6
7 /** \class AliHLTTPCFitter
8 <pre>
9 //_____________________________________________________________
10 // AliHLTTPCFitter
11 //
12 // Fit class HLT for helix
13 </pre>
14 */
15
16 #include <math.h>
17 #include "AliHLTTPCLogging.h"
18 #include "AliHLTTPCFitter.h"
19 #include "AliHLTTPCVertex.h"
20 #include "AliHLTTPCTrack.h"
21 #include "AliHLTTPCSpacePointData.h"
22 #include "AliHLTTPCMemHandler.h"
23 #include "AliHLTTPCTransform.h"
24 #include "AliHLTTPC.h"
25
26 #if __GNUC__ >= 3
27 using namespace std;
28 #endif
29
30 ClassImp(AliHLTTPCFitter)
31
32
33 AliHLTTPCFitter::AliHLTTPCFitter()
34 {
35   //constructor
36   fTrack=0;
37   fVertex=0;
38   memset(fClusters,0,36*6*sizeof(AliHLTTPCSpacePointData*));
39 }
40
41 AliHLTTPCFitter::AliHLTTPCFitter(AliHLTTPCVertex *vertex,Bool_t vertexconstraint)
42 {
43   //constructor
44   fTrack=0;
45   fVertex = vertex;
46   fVertexConstraint=vertexconstraint;
47   memset(fClusters,0,36*6*sizeof(AliHLTTPCSpacePointData*));
48 }
49
50 AliHLTTPCFitter::~AliHLTTPCFitter()
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 AliHLTTPCFitter::LoadClusters(Char_t *path,Int_t event,Bool_t sp)
64 {
65   //load clusters
66   Char_t fname[256];
67   AliHLTTPCMemHandler *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 AliHLTTPCMemHandler();
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] = (AliHLTTPCSpacePointData*)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 AliHLTTPCFitter::SortTrackClusters(AliHLTTPCTrack *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           AliHLTTPCSpacePointData *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 AliHLTTPCFitter::FitHelix(AliHLTTPCTrack *track)
144 {
145   //fit helix parameters
146   fTrack = track;
147   if(FitCircle())
148     {
149       LOG(AliHLTTPCLog::kError,"AliHLTTPCFitter::FitHelix","TrackFit")<<AliHLTTPCLog::kDec<<
150         "Problems during circle fit"<<ENDLOG;
151       return 1;
152     }
153   if(FitLine())
154     {
155       LOG(AliHLTTPCLog::kError,"AliHLTTPCFitter::FitHelix","TrackFit")<<AliHLTTPCLog::kDec<<
156         "Problems during line fit"<<ENDLOG;
157       return 1;
158     }
159   return 0;
160 }
161
162 Int_t AliHLTTPCFitter::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       AliHLTTPCSpacePointData *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       AliHLTTPCSpacePointData *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       AliHLTTPCSpacePointData *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(AliHLTTPCLog::kError,"AliHLTTPCFitter::FitCircle","TrackFit")<<AliHLTTPCLog::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 (AliHLTTPCTrack::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   AliHLTTPCSpacePointData *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 * AliHLTTPCTransform::Pi() ;
475   if ( psi < 0 ) psi = psi + 2*AliHLTTPCTransform::Pi();
476   
477   pt   = (Double_t)(AliHLTTPCTransform::GetBFact() * AliHLTTPCTransform::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 AliHLTTPCFitter::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() / ( AliHLTTPCTransform::GetBFact() * AliHLTTPCTransform::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       AliHLTTPCSpacePointData *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       AliHLTTPCSpacePointData *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       AliHLTTPCSpacePointData *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 * AliHLTTPCTransform::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       AliHLTTPCSpacePointData *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           AliHLTTPCSpacePointData *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       AliHLTTPCSpacePointData *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