]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HLT/TPCLib/AliHLTTPCConfMapFit.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
L3 becomes HLT
[u/mrichter/AliRoot.git] / HLT / TPCLib / AliHLTTPCConfMapFit.cxx
CommitLineData
a6c02c85 1// @(#) $Id$
4aa41877 2// Original: AliHLTConfMapFit.cxx,v 1.14 2005/06/14 10:55:21 cvetan
a6c02c85 3
4// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
5//*-- Copyright &copy ALICE HLT Group
6
a6c02c85 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
25using namespace std;
26#endif
27
28ClassImp(AliHLTTPCConfMapFit)
29
30
31AliHLTTPCConfMapFit::AliHLTTPCConfMapFit(AliHLTTPCConfMapTrack *track,AliHLTTPCVertex *vertex)
32{
33 //constructor
34 fTrack = track;
35 fVertex = vertex;
36}
37
38Int_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
db16520a 56// #### -B0-CHANGE-START == JMT
57
58Int_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
a6c02c85 75Int_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//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
402Int_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}
db16520a 526
527
528// #### -B0-CHANGE-START == JMT
529
530//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
531// Straight Line Fit in x-y plane
532//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
533Int_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//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
640Int_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