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