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