L3 becomes HLT
[u/mrichter/AliRoot.git] / HLT / TPCLib / AliHLTTPCConfMapFit.cxx
1 // @(#) $Id$
2 // Original: AliHLTConfMapFit.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 #include "AliHLTTPCRootTypes.h"
8 #include "AliHLTTPCLogging.h"
9 #include "AliHLTTPCVertex.h"
10 #include "AliHLTTPCConfMapTrack.h"
11 #include "AliHLTTPCConfMapPoint.h"
12 #include "AliHLTTPCTransform.h"
13 #include "AliHLTTPCConfMapFit.h"
14
15 /** \class AliHLTTPCConfMapFit
16 <pre>
17 //_____________________________________________________________
18 // AliHLTTPCConfMapFit
19 //
20 // Fit class for conformal mapping tracking
21 </pre>
22 */
23
24 #if __GNUC__ >= 3
25 using namespace std;
26 #endif
27
28 ClassImp(AliHLTTPCConfMapFit)
29
30
31 AliHLTTPCConfMapFit::AliHLTTPCConfMapFit(AliHLTTPCConfMapTrack *track,AliHLTTPCVertex *vertex)
32 {
33   //constructor
34   fTrack = track;
35   fVertex = vertex;
36 }
37
38 Int_t AliHLTTPCConfMapFit::FitHelix()
39 {
40   //fit the helix
41   if(FitCircle())
42     {
43       LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapFit::FitHelix","TrackFit")<<AliHLTTPCLog::kDec<<
44         "Problems during circle fit"<<ENDLOG;
45       return 1;
46     }
47   if(FitLine())
48     {
49       LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapFit::FitHelix","TrackFit")<<AliHLTTPCLog::kDec<<
50         "Problems during line fit"<<ENDLOG;
51       return 1;
52     }
53   return 0;
54 }
55
56 // #### -B0-CHANGE-START == JMT
57
58 Int_t AliHLTTPCConfMapFit::FitStraightLine() {
59     //fit the straight line 
60     if(FitLineXY()) {
61         LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapFit::FitStraightLine","TrackFit")<<AliHLTTPCLog::kDec<<
62             "Problems during stright line fit in XY plane"<<ENDLOG;
63         return 1;
64     }
65     if(FitLineSZ()){
66         LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapFit::FitStraightLine","TrackFit")<<AliHLTTPCLog::kDec<<
67             "Problems during stright line fit in SZ plane"<<ENDLOG;
68         return 1;
69     }
70     return 0;
71 }
72
73 // #### -B0-CHANGE-END == JMT
74
75 Int_t AliHLTTPCConfMapFit::FitCircle()
76 {
77   //-----------------------------------------------------------------
78   //Fits circle parameters using algorithm
79   //described by ChErnov and Oskov in Computer Physics
80   //Communications.
81   // 
82   //Written in FORTRAN by Jawluen Tang, Physics department , UT-Austin 
83   //Moved to C by Pablo Yepes
84   //Moved to AliROOT by ASV.
85   //------------------------------------------------------------------
86   
87   Double_t wsum  = 0.0 ;
88   Double_t xav   = 0.0 ;
89   Double_t yav   = 0.0 ;
90   
91   Int_t num_of_hits = fTrack->GetNumberOfPoints();
92   //
93   //     Loop over hits calculating average
94   Int_t co=0;
95   
96   for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())
97     {
98       co++;
99       AliHLTTPCConfMapPoint *cHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
100       cHit->SetXYWeight( 1./ (Double_t)(cHit->GetXerr()*cHit->GetXerr() + cHit->GetYerr()*cHit->GetYerr()) );
101       wsum      += cHit->GetXYWeight() ;
102       xav       += cHit->GetXYWeight() * cHit->GetX() ;
103       yav       += cHit->GetXYWeight() * cHit->GetY() ;
104     }
105   if(co!=num_of_hits) 
106     LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapFit::FitCircle","TrackFit")<<AliHLTTPCLog::kDec<<
107       "Mismatch of hits. Counter: "<<co<<" nHits: "<<num_of_hits<<ENDLOG;
108   if (fTrack->ComesFromMainVertex() == true)
109     {    
110       wsum += fVertex->GetXYWeight() ;
111       xav  += fVertex->GetX() ;
112       yav  += fVertex->GetY() ;
113     }
114   
115   xav = xav / wsum ;
116   yav = yav / wsum ;
117 //
118 //  CALCULATE <X**2>, <XY>, AND <Y**2> WITH <X> = 0, & <Y> = 0
119 //
120   Double_t xxav  = 0.0 ;
121   Double_t xyav  = 0.0 ; 
122   Double_t yyav  = 0.0 ;
123   Double_t xi, yi ;
124
125   for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())
126     { 
127       //AliHLTTPCConfMapPoint *cHit = (AliHLTTPCConfMapPoint *)hits->At(hit_counter);
128       AliHLTTPCConfMapPoint *cHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
129       xi        = cHit->GetX() - xav ;
130       yi        = cHit->GetY() - yav ;
131       xxav     += xi * xi * cHit->GetXYWeight() ;
132       xyav     += xi * yi * cHit->GetXYWeight() ;
133       yyav     += yi * yi * cHit->GetXYWeight() ;
134     }
135   
136   if (fTrack->ComesFromMainVertex() == true)
137     {
138       xi        = fVertex->GetX() - xav ;
139       yi        = fVertex->GetY() - yav ;
140       xxav     += xi * xi * fVertex->GetXYWeight() ;
141       xyav     += xi * yi * fVertex->GetXYWeight() ;
142       yyav     += yi * yi * fVertex->GetXYWeight() ; 
143     }
144   xxav = xxav / wsum ;
145   xyav = xyav / wsum ;
146   yyav = yyav / wsum ;
147 //
148 //-->  ROTATE COORDINATES SO THAT <XY> = 0
149 //
150 //-->  SIGN(C**2 - S**2) = SIGN(XXAV - YYAV) >
151 //-->  &                                     > ==> NEW : (XXAV-YYAV) > 0
152 //-->  SIGN(S) = SIGN(XYAV)                  >
153
154   Double_t a = fabs( xxav - yyav ) ;
155   Double_t b = 4.0 * xyav * xyav ;
156
157   Double_t asqpb  = a * a + b  ;
158   Double_t rasqpb = sqrt ( asqpb) ;
159
160   Double_t splus  = 1.0 + a / rasqpb ;
161   Double_t sminus = b / (asqpb * splus) ;
162
163   splus  = sqrt (0.5 * splus ) ;
164   sminus = sqrt (0.5 * sminus) ;
165 //
166 //->  FIRST REQUIRE : SIGN(C**2 - S**2) = SIGN(XXAV - YYAV)
167 //
168   Double_t sinrot, cosrot ;
169   if ( xxav <= yyav ) {
170          cosrot = sminus ;
171          sinrot = splus  ;
172   }
173   else {
174           cosrot = splus ;
175           sinrot = sminus ;
176   }
177 //
178 //->  REQUIRE : SIGN(S) = SIGN(XYAV) * SIGN(C) (ASSUMING SIGN(C) > 0)
179 //
180   if ( xyav < 0.0 ) sinrot = - sinrot ;
181 //
182 //-->  WE NOW HAVE THE SMALLEST ANGLE THAT GUARANTEES <X**2> > <Y**2>
183 //-->  TO GET THE SIGN OF THE CHARGE RIGHT, THE NEW X-AXIS MUST POINT
184 //-->  OUTWARD FROM THE ORGIN.  WE ARE FREE TO CHANGE SIGNS OF BOTH
185 //-->  COSROT AND SINROT SIMULTANEOUSLY TO ACCOMPLISH THIS.
186 //
187 //-->  CHOOSE SIGN OF C WISELY TO BE ABLE TO GET THE SIGN OF THE CHARGE
188 //
189   if ( cosrot*xav+sinrot*yav < 0.0 ) {
190           cosrot = -cosrot ;
191           sinrot = -sinrot ;
192   }
193 //
194 //->  NOW GET <R**2> AND RSCALE= SQRT(<R**2>)
195 //
196   Double_t rrav   = xxav + yyav ;
197   Double_t rscale = sqrt(rrav) ;
198
199   xxav   = 0.0 ;
200   yyav   = 0.0 ;
201   xyav   = 0.0 ;
202   Double_t xrrav = 0.0 ;
203   Double_t yrrav = 0.0 ;
204   Double_t rrrrav  = 0.0 ;
205
206   Double_t xixi, yiyi, riri, wiriri, xold, yold ;
207   
208   //for (hit_counter=0; hit_counter<num_of_hits; hit_counter++) 
209   for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())  
210     { 
211       //AliHLTTPCConfMapPoint *cHit = (AliHLTTPCConfMapPoint*)hits->At(hit_counter);  
212       AliHLTTPCConfMapPoint* cHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
213
214       xold = cHit->GetX() - xav ;
215       yold = cHit->GetY() - yav ;
216       //
217       //-->  ROTATE SO THAT <XY> = 0 & DIVIDE BY RSCALE SO THAT <R**2> = 1
218       //
219       xi = (  cosrot * xold + sinrot * yold ) / rscale ;
220       yi = ( -sinrot * xold + cosrot * yold ) / rscale ;
221       
222       xixi   = xi * xi ;
223       yiyi   = yi * yi ;
224       riri   = xixi + yiyi ;
225       wiriri = cHit->GetXYWeight() * riri ;
226       
227       xyav   += cHit->GetXYWeight() * xi * yi ;
228       xxav   += cHit->GetXYWeight() * xixi ;
229       yyav   += cHit->GetXYWeight() * yiyi ;
230       
231       xrrav  += wiriri * xi ;
232       yrrav  += wiriri * yi ;
233       rrrrav += wiriri * riri ;
234     }
235   //
236 //   Include vertex if required
237 //
238   if (fTrack->ComesFromMainVertex() == true)
239     {
240         xold = fVertex->GetX() - xav ;
241         yold = fVertex->GetY() - yav ;
242         //
243         //-->  ROTATE SO THAT <XY> = 0 & DIVIDE BY RSCALE SO THAT <R**2> = 1
244         //
245         xi = (  cosrot * xold + sinrot * yold ) / rscale ;
246         yi = ( -sinrot * xold + cosrot * yold ) / rscale ;
247         
248         xixi   = xi * xi ;
249         yiyi   = yi * yi ;
250         riri   = xixi + yiyi ;
251         wiriri = fVertex->GetXYWeight() * riri ;
252
253         xyav   += fVertex->GetXYWeight() * xi * yi ;
254         xxav   += fVertex->GetXYWeight() * xixi ;
255         yyav   += fVertex->GetXYWeight() * yiyi ;
256
257         xrrav  += wiriri * xi ;
258         yrrav  += wiriri * yi ;
259         rrrrav += wiriri * riri ;
260   }
261   //
262   //    
263   //
264   //-->  DIVIDE BY WSUM TO MAKE AVERAGES
265   //
266   xxav    = xxav   / wsum ;
267   yyav    = yyav   / wsum ;
268   xrrav   = xrrav  / wsum ;
269   yrrav   = yrrav  / wsum ;
270   rrrrav  = rrrrav / wsum ;
271   xyav    = xyav   / wsum ;
272
273   Int_t const ntry = 5 ;
274 //
275 //-->  USE THESE TO GET THE COEFFICIENTS OF THE 4-TH ORDER POLYNIMIAL
276 //-->  DON'T PANIC - THE THIRD ORDER TERM IS ZERO !
277 //
278   Double_t xrrxrr = xrrav * xrrav ;
279   Double_t yrryrr = yrrav * yrrav ;
280   Double_t rrrrm1 = rrrrav - 1.0  ;
281   Double_t xxyy   = xxav  * yyav  ;        
282
283   Double_t c0  =          rrrrm1*xxyy - xrrxrr*yyav - yrryrr*xxav ;
284   Double_t c1  =        - rrrrm1      + xrrxrr      + yrryrr   - 4.0*xxyy ;        
285   Double_t c2  =   4.0  + rrrrm1                               - 4.0*xxyy ;           
286   Double_t c4  = - 4.0  ;                
287 //
288 //-->  COEFFICIENTS OF THE DERIVATIVE - USED IN NEWTON-RAPHSON ITERATIONS
289 //
290   Double_t c2d =   2.0 * c2 ;
291   Double_t c4d =   4.0 * c4 ;
292 //
293 //-->  0'TH VALUE OF LAMDA - LINEAR INTERPOLATION BETWEEN P(0) & P(YYAV)
294 //
295 //   LAMDA = YYAV * C0 / (C0 + YRRSQ * (XXAV-YYAV))
296   Double_t lamda  = 0.0 ;
297   Double_t dlamda = 0.0 ;
298 //
299   Double_t chiscl = wsum * rscale * rscale ;
300   Double_t dlamax = 0.001 / chiscl ;   
301    
302   Double_t p, pd ;
303   for ( int itry = 1 ; itry <= ntry ; itry++ ) {
304      p      = c0 + lamda * (c1 + lamda * (c2 + lamda * lamda * c4 )) ;
305      pd     = (c1 + lamda * (c2d + lamda * lamda * c4d)) ;
306      dlamda = -p / pd ;
307      lamda  = lamda + dlamda ;
308      if (fabs(dlamda)<   dlamax) break ;
309   }
310
311   Double_t chi2 = (Double_t)(chiscl * lamda) ;
312  
313   fTrack->SetChiSq1(chi2);
314   // Double_t dchisq = chiscl * dlamda ;             
315 //
316 //-->  NOW CALCULATE THE MATRIX ELEMENTS FOR ALPHA, BETA & KAPPA
317 //
318   Double_t h11   = xxav  -     lamda ;
319   Double_t h14   = xrrav ;
320   Double_t h22   = yyav  -     lamda ; 
321   Double_t h24   = yrrav ;
322   Double_t h34   = 1.0   + 2.0*lamda ;
323   if ( h11 == 0.0 || h22 == 0.0 ){
324     LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapFit::FitCircle","TrackFit")<<AliHLTTPCLog::kDec<<
325       "Problems fitting circle"<<ENDLOG;
326     return 1 ;
327   }
328   Double_t rootsq = (h14*h14)/(h11*h11) + 4.0*h34 ;
329
330   Double_t ratio, kappa, beta ;
331   if ( fabs(h22) > fabs(h24) ) {
332      ratio  = h24 / h22 ;
333      rootsq = ratio * ratio + rootsq ;
334      kappa = 1.0 / sqrt(rootsq) ;
335      beta  = - ratio * kappa ;
336   }
337   else {
338      ratio  = h22 / h24 ;
339      rootsq = 1.0 + ratio * ratio * rootsq ;
340      beta  = 1.0 / sqrt(rootsq) ;
341      if ( h24 > 0 ) beta = - beta ;
342      kappa = -ratio * beta ;
343   }            
344   Double_t alpha = - (h14/h11) * kappa ;
345 //
346 //-->  transform these into the lab coordinate system
347 //-->  first get kappa and back to real dimensions
348 //
349   Double_t kappa1 = kappa / rscale ;
350   Double_t dbro   = 0.5   / kappa1 ;
351 //
352 //-->  next rotate alpha and beta and scale
353 //
354   Double_t alphar = (cosrot * alpha - sinrot * beta)* dbro ;
355   Double_t betar  = (sinrot * alpha + cosrot * beta)* dbro ;
356 //
357 //-->  then translate by (xav,yav)
358 //
359   Double_t acent  = (double)(xav - alphar) ;
360   Double_t bcent  = (double)(yav - betar ) ;
361   Double_t radius = (double)dbro ;
362 //
363 //   Get charge
364 //
365   Int_t q = ( ( yrrav < 0 ) ? 1 : -1 ) ;
366
367   fTrack->SetCharge(q);
368   
369   
370   //Set the first point on the track to the space point coordinates of the innermost track
371   //This will be updated to lie on the fit later on (AliHLTTPCTrack::UpdateToFirstPoint).
372   Double_t x0,y0,psi,pt ;
373   AliHLTTPCConfMapPoint *lHit = (AliHLTTPCConfMapPoint*)fTrack->GetLastHit();
374   x0 = lHit->GetX();
375   y0 = lHit->GetY();
376   fTrack->SetFirstPoint(x0,y0,0); //Z-value is set in FitLine
377
378   psi  = (Double_t)atan2(bcent-y0,acent-x0) ;
379   psi  = psi + q * AliHLTTPCTransform::PiHalf();
380   if ( psi < 0 ) psi = psi + AliHLTTPCTransform::TwoPi();
381   pt   = (Double_t)(AliHLTTPCTransform::GetBFieldValue() * radius ) ;
382   
383   //Update the track parameters with the parameters from this fit:
384   fTrack->SetPsi(psi);
385   fTrack->SetPt(pt);
386   fTrack->SetRadius(radius);
387   fTrack->SetCenterX(acent);
388   fTrack->SetCenterY(bcent);
389
390   //
391 //    Get errors from fast fit
392 //
393   //if ( getPara()->getErrors ) getErrorsCircleFit ( acent, bcent, radius ) ;
394 //
395   return 0 ;
396   
397 }
398
399 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
400 //    Fit Line in s-z plane
401 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
402 Int_t AliHLTTPCConfMapFit::FitLine ( )
403 {
404   //
405   //Initialization 
406   //
407   Double_t sum = 0.F ;
408   Double_t ss  = 0.F ;
409   Double_t sz  = 0.F ;
410   Double_t sss = 0.F ;
411   Double_t ssz = 0.F ;
412   //
413   //find sum , sums ,sumz, sumss 
414   // 
415   Double_t dx, dy ;
416   Double_t radius = (Double_t)(fTrack->GetPt() / AliHLTTPCTransform::GetBFieldValue() ) ;
417
418   //TObjArray *hits = fTrack->GetHits();
419   //Int_t num_of_hits = fTrack->GetNumberOfPoints();
420
421   if (0)// fTrack->ComesFromMainVertex() == true ) 
422     {
423       dx = ((AliHLTTPCConfMapPoint*)fTrack->GetFirstHit())->GetX() - fVertex->GetX();
424       dy = ((AliHLTTPCConfMapPoint*)fTrack->GetFirstHit())->GetY() - fVertex->GetY() ;
425     }
426   else 
427     {
428       dx = ((AliHLTTPCConfMapPoint *)fTrack->GetFirstHit())->GetX() - ((AliHLTTPCConfMapPoint *)fTrack->GetLastHit())->GetX() ;
429       dy = ((AliHLTTPCConfMapPoint *)fTrack->GetFirstHit())->GetY() - ((AliHLTTPCConfMapPoint *)fTrack->GetLastHit())->GetY() ;
430       //dx = ((AliHLTTPCConfMapPoint *)hits->First())->GetX() - ((AliHLTTPCConfMapPoint *)hits->Last())->GetX() ;
431       //dy = ((AliHLTTPCConfMapPoint *)hits->First())->GetY() - ((AliHLTTPCConfMapPoint *)hits->Last())->GetY() ;
432     }
433   
434   Double_t localPsi = 0.5F * sqrt ( dx*dx + dy*dy ) / radius ;
435   Double_t total_s ;
436   
437   if ( fabs(localPsi) < 1. ) 
438     {
439       total_s = 2.0 * radius * asin ( localPsi ) ;
440     } 
441   else 
442     { 
443       total_s = 2.0 * radius * AliHLTTPCTransform::Pi() ;
444     } 
445   
446   AliHLTTPCConfMapPoint *previousHit = NULL;
447
448   // FtfBaseHit *previousHit = 0  ;
449   
450   //for ( startLoop() ; done() ; nextHit() ) {
451   Double_t dpsi,s;
452
453   //  for(hit_counter=0; hit_counter<num_of_hits; hit_counter++)
454   for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())  
455     {
456       // AliHLTTPCConfMapPoint *cHit = (AliHLTTPCConfMapPoint*)hits->At(hit_counter);
457       AliHLTTPCConfMapPoint *cHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
458       // if ( GetCurrentHit() != GetFirstHit() ) 
459       if(cHit != fTrack->GetFirstHit())//  hits->First())
460         {
461           dx   = cHit->GetX() - previousHit->GetX() ;
462           dy   = cHit->GetY() - previousHit->GetY() ;
463           dpsi = 0.5 * (Double_t)sqrt ( dx*dx + dy*dy ) / radius ;
464           fTrack->SetPsierr(dpsi);
465           s = previousHit->GetS() - 2.0 * radius * (Double_t)asin ( dpsi ) ;
466           cHit->SetS(s);
467         }
468       else
469         cHit->SetS(total_s);
470       //        cHit->s = total_s ;
471     
472       sum += cHit->GetZWeight() ;
473       ss  += cHit->GetZWeight() * cHit->GetS() ;
474       sz  += cHit->GetZWeight() * cHit->GetZ() ;
475       sss += cHit->GetZWeight() * cHit->GetS() * cHit->GetS() ;
476       ssz += cHit->GetZWeight() * cHit->GetS() * cHit->GetZ() ;
477       previousHit = cHit ;
478     }
479   
480   Double_t chi2,det = sum * sss - ss * ss;
481   if ( fabs(det) < 1e-20)
482     { 
483       chi2 = 99999.F ;
484       fTrack->SetChiSq2(chi2);
485       return 0 ;
486     }
487   
488   //Compute the best fitted parameters A,B
489   Double_t tanl,z0,dtanl,dz0;
490
491   tanl = (Double_t)((sum * ssz - ss * sz ) / det );
492   z0   = (Double_t)((sz * sss - ssz * ss ) / det );
493
494   fTrack->SetTgl(tanl);
495   fTrack->SetZ0(z0);
496   
497   //     calculate chi-square 
498   
499   chi2 = 0.;
500   Double_t r1 ;
501   
502   //for(hit_counter=0; hit_counter<num_of_hits; hit_counter++)
503   for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())  
504     {
505       //AliHLTTPCConfMapPoint *cHit = (AliHLTTPCConfMapPoint*)hits->At(hit_counter);
506       AliHLTTPCConfMapPoint *cHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
507       r1   = cHit->GetZ() - tanl * cHit->GetS() - z0 ;
508       chi2 += (Double_t) ( (Double_t)cHit->GetZWeight() * (r1 * r1) );
509     }
510   fTrack->SetChiSq2(chi2);
511   //
512   //     calculate estimated variance
513   //      varsq=chi/(double(n)-2.) 
514   //     calculate covariance matrix 
515   //      siga=sqrt(varsq*sxx/det) 
516   //      sigb=sqrt(varsq*sum/det) 
517   //
518   dtanl = (Double_t) ( sum / det );
519   dz0   = (Double_t) ( sss / det );
520   
521   fTrack->SetTglerr(dtanl);
522   fTrack->SetZ0err(dz0);
523   
524   return 0 ;
525
526
527
528 // #### -B0-CHANGE-START == JMT
529
530 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
531 //    Straight Line Fit  in x-y plane
532 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
533 Int_t AliHLTTPCConfMapFit::FitLineXY ( ){
534     // -----------------------------------------------------------------------------
535     // Implementation after Numerical Recipes in C, 2nd Edtion, chapter 15.2, p. 661
536     // with y = b*x + a
537     // and Data Analysis for Physical Science Students, Luis Lyons, chapter 2.4 p.51 
538     // with y = a' + bx' , x' = x - <x>
539     // -----------------------------------------------------------------------------
540
541     Double_t S = 0.;
542     Double_t Sx = 0.;
543
544     Double_t SPrime = 0.;
545     Double_t SxPrime = 0.;
546     Double_t SxxPrime = 0.;
547     Double_t SyPrime = 0.;
548     Double_t SxyPrime = 0.;
549
550     Double_t chi2 = 0.;
551
552     Int_t num_of_hits = fTrack->GetNumberOfPoints();
553
554     Int_t co=0;
555     
556     // - Loop over hits calculating average : xav
557     for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit()) {
558         co++;
559         AliHLTTPCConfMapPoint *currentHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
560         // ** maybe not necessary, already done in ConfMapPoint
561         currentHit->SetXYWeight( 1./ (Double_t)(currentHit->GetXerr()*currentHit->GetXerr() + currentHit->GetYerr()*currentHit->GetYerr()) );
562         // **
563         S   += currentHit->GetXYWeight();
564         Sx  += currentHit->GetXYWeight() * currentHit->GetX();
565     }   
566     
567     if(co!=num_of_hits) 
568         LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapFit::FitLineXY","TrackFit") << "Mismatch of hits. Counter: "<<co<<" nHits: "<<num_of_hits<<ENDLOG;
569
570     Double_t xav = (Double_t)Sx / S;
571     
572     // Calculate weighted means
573     for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit()) {
574         AliHLTTPCConfMapPoint *currentHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
575
576         Double_t xPrime =  currentHit->GetX() - xav;
577         SPrime   += currentHit->GetXYWeight();
578         SxPrime  += currentHit->GetXYWeight() * xPrime;
579         SxxPrime += currentHit->GetXYWeight() * xPrime * xPrime;
580         SyPrime  += currentHit->GetXYWeight() * currentHit->GetY();
581         SxyPrime += currentHit->GetXYWeight() * xPrime * currentHit->GetY();
582     }
583
584     Double_t det = SPrime*SxxPrime + SxPrime*SxPrime;
585
586     if (fabs(det) < 1e-20) { 
587         LOG(AliHLTTPCLog::kDebug,"AliHLTTPCConfMapFit::FitLineXY","TrackFit") << "Determinant == 0" << ENDLOG;  
588         chi2 = 99999.F ;
589         fTrack->SetChiSq1(chi2);
590         return -1 ;
591     }
592
593     Double_t b   = (Double_t)(SPrime*SxyPrime - SxPrime*SyPrime) / det;        // line parameter b
594     Double_t aPrime   = (Double_t)(SxxPrime*SyPrime - SxPrime*SxyPrime) / det; // line parameter a
595
596     Double_t sigma2b = (Double_t)1. / SxxPrime;
597     //-- Double_t sigma2aprime = (Double_t)1. /SPrime;
598
599     // Get gradient angle psi of line in xy plane
600     Double_t psi  = (Double_t) atan(b) ; 
601
602     // Calculate chi2
603     for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit()) {
604         AliHLTTPCConfMapPoint *currentHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
605         Double_t tempchi = currentHit->GetY() - aPrime - b*(currentHit->GetX() - xav);
606         chi2 += tempchi*tempchi*currentHit->GetXYWeight() ;
607     }
608
609     Double_t a = aPrime - b*xav;
610
611
612     // Set TrackParameter
613     fTrack->SetChiSq1(chi2);
614     fTrack->SetPsi(psi);
615     fTrack->SetPsierr(sigma2b);   
616     fTrack->SetCenterX(0.);    // Set to point on the track (for UpdateToFirstPoint)
617     fTrack->SetCenterY(a);     // Set to point on the track (for UpdateToFirstPoint)
618
619     //Set the first point on the track to the space point coordinates of the innermost track
620     //This will be updated to lie on the fit later on (AliHLTTPCTrack::UpdateToFirstPoint).
621     AliHLTTPCConfMapPoint *lastHit = (AliHLTTPCConfMapPoint*)fTrack->GetLastHit();
622     Double_t x0 = lastHit->GetX();
623     Double_t y0 = lastHit->GetY();
624     fTrack->SetFirstPoint(x0,y0,0); //Z-value is set in FitLineSZ
625     
626
627     //Set Defaults
628     fTrack->SetRadius(-1.);
629     fTrack->SetCharge(1);
630     fTrack->SetPt(-1.);
631   
632
633     return 0;
634 }
635
636
637 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
638 //    Straight Line Fit  in s-z plane
639 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
640 Int_t AliHLTTPCConfMapFit::FitLineSZ ( ){
641     // -----------------------------------------------------------------------------
642     // Implementation after Numerical Recipes in C, 2nd Edtion, chapter 15.2, p. 661
643     // with z = b*s + a
644     // and Data Analysis for Physical Science Students, Luis Lyons, chapter 2.4 p.51 
645     // with z = a' + bs' , s' = s - <s>
646     // -----------------------------------------------------------------------------
647
648     Double_t S = 0.;
649     Double_t Ss = 0.;
650
651     Double_t SPrime = 0.;
652     Double_t SsPrime = 0.;
653     Double_t SssPrime = 0.;
654     Double_t SzPrime = 0.;
655     Double_t SszPrime = 0.;
656
657     Double_t chi2 = 0.;
658     Double_t s = 0.;
659
660     AliHLTTPCConfMapPoint *previousHit = NULL;
661   
662     // - Loop over hits calculating length in xy-plane: s
663     // - Loop over hits calculating average : sav
664     for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit()) {
665         AliHLTTPCConfMapPoint *currentHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
666         if(currentHit != fTrack->GetFirstHit()) {
667             Double_t dx = currentHit->GetX() - previousHit->GetX() ;
668             Double_t dy = currentHit->GetY() - previousHit->GetY() ;
669             s = previousHit->GetS() - (Double_t)sqrt ( dx*dx + dy*dy );
670         }
671         else{
672             Double_t dx = ((AliHLTTPCConfMapPoint *)fTrack->GetFirstHit())->GetX() - ((AliHLTTPCConfMapPoint *)fTrack->GetLastHit())->GetX(); 
673             Double_t dy = ((AliHLTTPCConfMapPoint *)fTrack->GetFirstHit())->GetY() - ((AliHLTTPCConfMapPoint *)fTrack->GetLastHit())->GetY();
674             s = (Double_t)sqrt ( dx*dx + dy*dy );
675         }
676
677         currentHit->SetS(s);
678
679         S   += currentHit->GetZWeight();
680         Ss  += currentHit->GetZWeight() * currentHit->GetS();
681     }
682
683     Double_t sav = (Double_t)Ss / S;
684
685     // Calculate weighted means
686     for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit()) {
687         AliHLTTPCConfMapPoint *currentHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
688
689         Double_t sPrime =  currentHit->GetS() - sav;
690         SPrime   += currentHit->GetZWeight();
691         SsPrime  += currentHit->GetZWeight() * sPrime;
692         SssPrime += currentHit->GetZWeight() * sPrime * sPrime;
693         SzPrime  += currentHit->GetZWeight() * currentHit->GetZ();
694         SszPrime += currentHit->GetZWeight() * sPrime * currentHit->GetZ();
695     }
696
697     Double_t det = SPrime*SssPrime + SsPrime*SsPrime;
698
699     if (fabs(det) < 1e-20) { 
700         LOG(AliHLTTPCLog::kDebug,"AliHLTTPCConfMapFit::FitLineSZ","TrackFit") << "Determinant == 0" << ENDLOG;  
701         chi2 = 99999.F ;
702         fTrack->SetChiSq2(chi2);
703         return -1 ;
704     }
705
706     Double_t b   = (Double_t)(SPrime*SszPrime - SsPrime*SzPrime) / det;        // line parameter b
707     Double_t aPrime   = (Double_t)(SssPrime*SzPrime - SsPrime*SszPrime) / det; // line parameter a
708
709     Double_t a = aPrime - b*sav;
710
711     Double_t sigma2b = (Double_t) 1. / SssPrime;
712     Double_t sigma2aprime = (Double_t) 1. /SPrime;
713
714     Double_t sigma2a = sigma2aprime + sav*sav * sigma2b*sigma2b;
715
716     // Calculate chi2
717     for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit()) {
718         AliHLTTPCConfMapPoint *currentHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
719         Double_t tempchi = currentHit->GetZ() - aPrime - b*(currentHit->GetS() - sav);
720         chi2 += tempchi*tempchi*currentHit->GetZWeight() ;
721     }
722
723     // Set TrackParameter
724     fTrack->SetChiSq2(chi2);
725     fTrack->SetTgl(b);
726     fTrack->SetZ0(a);
727     fTrack->SetTglerr(sigma2b);
728 //  fTrack->SetZ0err(sigma2aprime);   // maybe subject to check
729     fTrack->SetZ0err(sigma2a);        // maybe subject to check
730     return 0;
731 }
732
733 // #### -B0-CHANGE-END == JMT