removing deprecated components from registration
[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(previousHit)//  hits->First())
490         {
491           // this means cHit != fTrack->GetFirstHit()
492           dx   = cHit->GetX() - previousHit->GetX() ;
493           dy   = cHit->GetY() - previousHit->GetY() ;
494           dpsi = 0.5 * (Double_t)sqrt ( dx*dx + dy*dy ) / radius ;
495           fTrack->SetPsierr(dpsi);
496           s = previousHit->GetS() - 2.0 * radius * (Double_t)asin ( dpsi ) ;
497           cHit->SetS(s);
498         }
499       else
500         cHit->SetS(totalS);
501       //        cHit->s = totalS ;
502     
503       sum += cHit->GetZWeight() ;
504       ss  += cHit->GetZWeight() * cHit->GetS() ;
505       sz  += cHit->GetZWeight() * cHit->GetZ() ;
506       sss += cHit->GetZWeight() * cHit->GetS() * cHit->GetS() ;
507       ssz += cHit->GetZWeight() * cHit->GetS() * cHit->GetZ() ;
508       previousHit = cHit ;
509     }
510   
511   Double_t chi2,det = sum * sss - ss * ss;
512   if ( fabs(det) < 1e-20)
513     { 
514       chi2 = 99999.F ;
515       fTrack->SetChiSq2(chi2);
516       return 0 ;
517     }
518   
519   //Compute the best fitted parameters A,B
520   Double_t tanl,z0,dtanl,dz0;
521
522   tanl = (Double_t)((sum * ssz - ss * sz ) / det );
523   z0   = (Double_t)((sz * sss - ssz * ss ) / det );
524
525   fTrack->SetTgl(tanl);
526   fTrack->SetZ0(z0);
527   
528   //     calculate chi-square 
529   
530   chi2 = 0.;
531   Double_t r1 ;
532   
533   //for(hit_counter=0; hit_counter<numOfHits; hit_counter++)
534   for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())  
535     {
536       //AliHLTTPCConfMapPoint *cHit = (AliHLTTPCConfMapPoint*)hits->At(hit_counter);
537       AliHLTTPCConfMapPoint *cHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
538       r1   = cHit->GetZ() - tanl * cHit->GetS() - z0 ;
539       chi2 += (Double_t) ( (Double_t)cHit->GetZWeight() * (r1 * r1) );
540     }
541   fTrack->SetChiSq2(chi2);
542   //
543   //     calculate estimated variance
544   //      varsq=chi/(double(n)-2.) 
545   //     calculate covariance matrix 
546   //      siga=sqrt(varsq*sxx/det) 
547   //      sigb=sqrt(varsq*sum/det) 
548   //
549   dtanl = (Double_t) ( sum / det );
550   dz0   = (Double_t) ( sss / det );
551   
552   fTrack->SetTglerr(dtanl);
553   fTrack->SetZ0err(dz0);
554
555     //The calculation of pT comes from "Perticle Physics Booklet": 
556   //24.11 Measurement of particle momenta in a uniform magnetic field.
557   //(dk)^2=(dk_res)^2 + (dk_ms)^2
558   //for now k_ms is 0. Need to find length of track in 3D. 
559
560   Double_t lengthXY = fTrack->GetLengthXY();
561   //Double_t lengthTot = fTrack->GetLengthTot();
562   //Double_t beta = fTrack->GetP()/TMath::Sqrt((fTrack->GetP()*fTrack->GetP())+(0.13957*0.13957));
563   //Double_t lambda = TMath::ATan(fTrack->GetTgl());
564   Double_t lengthXY2 = lengthXY*lengthXY;
565   Int_t nCluster4 = fTrack->GetNHits()+4;
566
567   Double_t Kres = 0.03/lengthXY2;
568   Kres = Kres * TMath::Sqrt(720/nCluster4);
569
570   //Double_t Kres = (0.03/(lengthXY*lengthXY))*TMath::Sqrt(720/(fTrack->GetNHits()+4));
571
572   //Double_t d = lengthTot*fTrack->GetP()*beta*TMath::Cos(lambda)*TMath::Cos(lambda);
573   //Double_t Kms = (0.016/d)*TMath::Sqrt(lengthTot/24.0);
574   Double_t Kms = 0.0;
575
576   Double_t KTot = TMath::Sqrt((Kres * Kres) + (Kms * Kms));
577   
578   Double_t Pterr = (1/(0.3*AliHLTTPCTransform::GetBField()))*KTot;
579
580   fTrack->SetPterr(Pterr);
581
582   return 0 ;
583
584
585
586 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
587 //    Straight Line Fit  in x-y plane
588 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
589 Int_t AliHLTTPCConfMapFit::FitLineXY ( ){
590     // -----------------------------------------------------------------------------
591     // Implementation after Numerical Recipes in C, 2nd Edtion, chapter 15.2, p. 661
592     // with y = b*x + a
593     // and Data Analysis for Physical Science Students, Luis Lyons, chapter 2.4 p.51 
594     // with y = a' + bx' , x' = x - <x>
595     // -----------------------------------------------------------------------------
596
597     Double_t s = 0.;
598     Double_t sx = 0.;
599
600     Double_t sPrime = 0.;
601     Double_t sxPrime = 0.;
602     Double_t sxxPrime = 0.;
603     Double_t syPrime = 0.;
604     Double_t sxyPrime = 0.;
605
606     Double_t chi2 = 0.;
607
608     Int_t numOfHits = fTrack->GetNumberOfPoints();
609
610     Int_t co=0;
611     
612     // - Loop over hits calculating average : xav
613     for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit()) {
614         co++;
615         AliHLTTPCConfMapPoint *currentHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
616         // ** maybe not necessary, already done in ConfMapPoint
617         currentHit->SetXYWeight( 1./ (Double_t)(currentHit->GetXerr()*currentHit->GetXerr() + currentHit->GetYerr()*currentHit->GetYerr()) );
618         // **
619         s   += currentHit->GetXYWeight();
620         sx  += currentHit->GetXYWeight() * currentHit->GetX();
621     }   
622     
623     if(co!=numOfHits) 
624         LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapFit::FitLineXY","TrackFit") << "Mismatch of hits. Counter: "<<co<<" nHits: "<<numOfHits<<ENDLOG;
625
626     Double_t xav = (Double_t)sx / s;
627     
628     // Calculate weighted means
629     for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit()) {
630         AliHLTTPCConfMapPoint *currentHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
631
632         Double_t xPrime =  currentHit->GetX() - xav;
633         sPrime   += currentHit->GetXYWeight();
634         sxPrime  += currentHit->GetXYWeight() * xPrime;
635         sxxPrime += currentHit->GetXYWeight() * xPrime * xPrime;
636         syPrime  += currentHit->GetXYWeight() * currentHit->GetY();
637         sxyPrime += currentHit->GetXYWeight() * xPrime * currentHit->GetY();
638     }
639
640     Double_t det = sPrime*sxxPrime + sxPrime*sxPrime;
641
642     if (fabs(det) < 1e-20) { 
643         LOG(AliHLTTPCLog::kDebug,"AliHLTTPCConfMapFit::FitLineXY","TrackFit") << "Determinant == 0" << ENDLOG;  
644         chi2 = 99999.F ;
645         fTrack->SetChiSq1(chi2);
646         return -1 ;
647     }
648
649     Double_t b   = (Double_t)(sPrime*sxyPrime - sxPrime*syPrime) / det;        // line parameter b
650     Double_t aPrime   = (Double_t)(sxxPrime*syPrime - sxPrime*sxyPrime) / det; // line parameter a
651
652     Double_t sigma2b = (Double_t)1. / sxxPrime;
653     //-- Double_t sigma2aprime = (Double_t)1. /sPrime;
654
655     // Get gradient angle psi of line in xy plane
656     Double_t psi  = (Double_t) atan(b) ; 
657
658     // Calculate chi2
659     for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit()) {
660         AliHLTTPCConfMapPoint *currentHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
661         Double_t tempchi = currentHit->GetY() - aPrime - b*(currentHit->GetX() - xav);
662         chi2 += tempchi*tempchi*currentHit->GetXYWeight() ;
663     }
664
665     Double_t a = aPrime - b*xav;
666
667
668     // Set TrackParameter
669     fTrack->SetChiSq1(chi2);
670     fTrack->SetPsi(psi);
671     fTrack->SetPsierr(sigma2b);   
672     fTrack->SetCenterX(0.);    // Set to point on the track (for UpdateToFirstPoint)
673     fTrack->SetCenterY(a);     // Set to point on the track (for UpdateToFirstPoint)
674
675     //Set the first point on the track to the space point coordinates of the innermost track
676     //This will be updated to lie on the fit later on (AliHLTTPCTrack::UpdateToFirstPoint).
677     AliHLTTPCConfMapPoint *lastHit = (AliHLTTPCConfMapPoint*)fTrack->GetLastHit();
678     Double_t x0 = lastHit->GetX();
679     Double_t y0 = lastHit->GetY();
680     fTrack->SetFirstPoint(x0,y0,0); //Z-value is set in FitLineSZ
681     
682
683     //Set Defaults
684     fTrack->SetRadius(-1.);
685     fTrack->SetCharge(1);
686     fTrack->SetPt(-1.);
687   
688
689     return 0;
690 }
691
692
693 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
694 //    Straight Line Fit  in s-z plane
695 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
696 Int_t AliHLTTPCConfMapFit::FitLineSZ ( ){
697     // -----------------------------------------------------------------------------
698     // Implementation after Numerical Recipes in C, 2nd Edtion, chapter 15.2, p. 661
699     // with z = b*s + a
700     // and Data Analysis for Physical Science Students, Luis Lyons, chapter 2.4 p.51 
701     // with z = a' + bs' , s' = s - <s>
702     // -----------------------------------------------------------------------------
703
704     Double_t S = 0.;
705     Double_t Ss = 0.;
706
707     Double_t sPrime = 0.;
708     Double_t ssPrime = 0.;
709     Double_t sssPrime = 0.;
710     Double_t szPrime = 0.;
711     Double_t sszPrime = 0.;
712
713     Double_t chi2 = 0.;
714
715     // Matthias 16.10.2007
716     // what's that!!! local variables 's' and 'S'
717     // change Double_t s = 0.; -> slength
718     Double_t slength = 0.;
719
720     // Matthias 23.02.2011
721     // this looks like a bug: previousHit is initialized, not changed
722     // in the loop, but dereferenced for all but the first hit
723     // changing the condition and adding an assignment at the end
724     // of the loop
725     AliHLTTPCConfMapPoint *previousHit = NULL;
726   
727     // - Loop over hits calculating length in xy-plane: s
728     // - Loop over hits calculating average : sav
729     for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit()) {
730         AliHLTTPCConfMapPoint *currentHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
731         if(previousHit) {
732             // this means currentHit != fTrack->GetFirstHit()
733             Double_t dx = currentHit->GetX() - previousHit->GetX() ;
734             Double_t dy = currentHit->GetY() - previousHit->GetY() ;
735             slength = previousHit->GetS() - (Double_t)sqrt ( dx*dx + dy*dy );
736         }
737         else{
738             Double_t dx = ((AliHLTTPCConfMapPoint *)fTrack->GetFirstHit())->GetX() - ((AliHLTTPCConfMapPoint *)fTrack->GetLastHit())->GetX(); 
739             Double_t dy = ((AliHLTTPCConfMapPoint *)fTrack->GetFirstHit())->GetY() - ((AliHLTTPCConfMapPoint *)fTrack->GetLastHit())->GetY();
740             slength = (Double_t)sqrt ( dx*dx + dy*dy );
741         }
742
743         currentHit->SetS(slength);
744
745         S   += currentHit->GetZWeight();
746         Ss  += currentHit->GetZWeight() * currentHit->GetS();
747
748         // Matthias 23.02.2011
749         // adding missing assignment, otherwise previousHit stays NULL
750         previousHit=currentHit;
751     }
752
753     Double_t sav = (Double_t)Ss / S;
754
755     // Calculate weighted means
756     for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit()) {
757         AliHLTTPCConfMapPoint *currentHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
758
759         // Matthias 20.05.2008
760         // here was a shadowed variable, sPrime is formerly defined
761         // renamed it to lsPrime ('local')
762         Double_t lsPrime =  currentHit->GetS() - sav;
763         lsPrime   += currentHit->GetZWeight();
764         ssPrime  += currentHit->GetZWeight() * lsPrime;
765         sssPrime += currentHit->GetZWeight() * lsPrime * lsPrime;
766         szPrime  += currentHit->GetZWeight() * currentHit->GetZ();
767         sszPrime += currentHit->GetZWeight() * lsPrime * currentHit->GetZ();
768     }
769
770     Double_t det = sPrime*sssPrime + ssPrime*ssPrime;
771
772     if (fabs(det) < 1e-20) { 
773         LOG(AliHLTTPCLog::kDebug,"AliHLTTPCConfMapFit::FitLineSZ","TrackFit") << "Determinant == 0" << ENDLOG;  
774         chi2 = 99999.F ;
775         fTrack->SetChiSq2(chi2);
776         return -1 ;
777     }
778
779     Double_t b   = (Double_t)(sPrime*sszPrime - ssPrime*szPrime) / det;        // line parameter b
780     Double_t aPrime   = (Double_t)(sssPrime*szPrime - ssPrime*sszPrime) / det; // line parameter a
781
782     Double_t a = aPrime - b*sav;
783
784     Double_t sigma2b = (Double_t) 1. / sssPrime;
785     Double_t sigma2aprime = (Double_t) 1. /sPrime;
786
787     Double_t sigma2a = sigma2aprime + sav*sav * sigma2b*sigma2b;
788
789     // Calculate chi2
790     for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit()) {
791         AliHLTTPCConfMapPoint *currentHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
792         Double_t tempchi = currentHit->GetZ() - aPrime - b*(currentHit->GetS() - sav);
793         chi2 += tempchi*tempchi*currentHit->GetZWeight() ;
794     }
795
796     // Set TrackParameter
797     fTrack->SetChiSq2(chi2);
798     fTrack->SetTgl(b);
799     fTrack->SetZ0(a);
800     fTrack->SetTglerr(sigma2b);
801 //  fTrack->SetZ0err(sigma2aprime);   // maybe subject to check
802     fTrack->SetZ0err(sigma2a);        // maybe subject to check
803     return 0;
804 }
805