]>
Commit | Line | Data |
---|---|---|
a6c02c85 | 1 | // @(#) $Id$ |
4aa41877 | 2 | // Original: AliHLTConfMapFit.cxx,v 1.14 2005/06/14 10:55:21 cvetan |
a6c02c85 | 3 | |
2a083ac4 | 4 | /************************************************************************** |
9be2600f | 5 | * This file is property of and copyright by the ALICE HLT Project * |
6 | * ALICE Experiment at CERN, All rights reserved. * | |
2a083ac4 | 7 | * * |
9be2600f | 8 | * Primary Authors: Anders Vestbo, maintained by * |
9 | * Matthias Richter <Matthias.Richter@ift.uib.no> * | |
10 | * for The ALICE HLT Project. * | |
2a083ac4 | 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 | */ | |
a6c02c85 | 26 | |
e67b0680 | 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 | ||
a6c02c85 | 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 | ||
a6c02c85 | 41 | #if __GNUC__ >= 3 |
42 | using namespace std; | |
43 | #endif | |
44 | ||
e67b0680 | 45 | ClassImp(AliHLTTPCConfMapFit); |
a6c02c85 | 46 | |
2a083ac4 | 47 | AliHLTTPCConfMapFit::AliHLTTPCConfMapFit() |
48 | : | |
49 | fTrack(NULL), | |
50 | fVertex(NULL) | |
51 | { | |
52 | //constructor | |
53 | } | |
54 | ||
a6c02c85 | 55 | AliHLTTPCConfMapFit::AliHLTTPCConfMapFit(AliHLTTPCConfMapTrack *track,AliHLTTPCVertex *vertex) |
2a083ac4 | 56 | : |
57 | fTrack(track), | |
58 | fVertex(vertex) | |
59 | ||
a6c02c85 | 60 | { |
61 | //constructor | |
2a083ac4 | 62 | } |
63 | ||
2a083ac4 | 64 | AliHLTTPCConfMapFit::~AliHLTTPCConfMapFit() |
65 | { | |
66 | // destructor | |
a6c02c85 | 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 | ||
db16520a | 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 | ||
a6c02c85 | 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 | ||
e67b0680 | 118 | Int_t numOfHits = fTrack->GetNumberOfPoints(); |
a6c02c85 | 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 | } | |
e67b0680 | 132 | if(co!=numOfHits) |
a6c02c85 | 133 | LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapFit::FitCircle","TrackFit")<<AliHLTTPCLog::kDec<< |
e67b0680 | 134 | "Mismatch of hits. Counter: "<<co<<" nHits: "<<numOfHits<<ENDLOG; |
a6c02c85 | 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 | ||
e67b0680 | 235 | //for (hit_counter=0; hit_counter<numOfHits; hit_counter++) |
a6c02c85 | 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 | ||
e67b0680 | 300 | const Int_t ntry = 5 ; |
a6c02c85 | 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 | ||
1ee9dca4 | 417 | fTrack->SetY0err(0.03); |
31f11c63 | 418 | //set error for pT and Y. psi, Z and Tgl are set. |
419 | ||
a6c02c85 | 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(); | |
e67b0680 | 449 | //Int_t numOfHits = fTrack->GetNumberOfPoints(); |
a6c02c85 | 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 ; | |
e67b0680 | 465 | Double_t totalS ; |
a6c02c85 | 466 | |
467 | if ( fabs(localPsi) < 1. ) | |
468 | { | |
e67b0680 | 469 | totalS = 2.0 * radius * asin ( localPsi ) ; |
a6c02c85 | 470 | } |
471 | else | |
472 | { | |
e67b0680 | 473 | totalS = 2.0 * radius * AliHLTTPCTransform::Pi() ; |
a6c02c85 | 474 | } |
475 | ||
476 | AliHLTTPCConfMapPoint *previousHit = NULL; | |
477 | ||
478 | // FtfBaseHit *previousHit = 0 ; | |
479 | ||
480 | //for ( startLoop() ; done() ; nextHit() ) { | |
481 | Double_t dpsi,s; | |
482 | ||
e67b0680 | 483 | // for(hit_counter=0; hit_counter<numOfHits; hit_counter++) |
a6c02c85 | 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() ) | |
8fb3f728 | 489 | if(previousHit)// hits->First()) |
a6c02c85 | 490 | { |
8fb3f728 | 491 | // this means cHit != fTrack->GetFirstHit() |
a6c02c85 | 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 | |
e67b0680 | 500 | cHit->SetS(totalS); |
501 | // cHit->s = totalS ; | |
a6c02c85 | 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 | ||
e67b0680 | 533 | //for(hit_counter=0; hit_counter<numOfHits; hit_counter++) |
a6c02c85 | 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); | |
1ee9dca4 | 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 | ||
ec75ebfe | 576 | Double_t KTot = TMath::Sqrt((Kres * Kres) + (Kms * Kms)); |
577 | ||
578 | Double_t Pterr = (1/(0.3*AliHLTTPCTransform::GetBField()))*KTot; | |
1ee9dca4 | 579 | |
580 | fTrack->SetPterr(Pterr); | |
581 | ||
a6c02c85 | 582 | return 0 ; |
583 | } | |
db16520a | 584 | |
585 | ||
db16520a | 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 | ||
e67b0680 | 597 | Double_t s = 0.; |
598 | Double_t sx = 0.; | |
db16520a | 599 | |
e67b0680 | 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.; | |
db16520a | 605 | |
606 | Double_t chi2 = 0.; | |
607 | ||
e67b0680 | 608 | Int_t numOfHits = fTrack->GetNumberOfPoints(); |
db16520a | 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 | // ** | |
e67b0680 | 619 | s += currentHit->GetXYWeight(); |
620 | sx += currentHit->GetXYWeight() * currentHit->GetX(); | |
db16520a | 621 | } |
622 | ||
e67b0680 | 623 | if(co!=numOfHits) |
624 | LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapFit::FitLineXY","TrackFit") << "Mismatch of hits. Counter: "<<co<<" nHits: "<<numOfHits<<ENDLOG; | |
db16520a | 625 | |
e67b0680 | 626 | Double_t xav = (Double_t)sx / s; |
db16520a | 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; | |
e67b0680 | 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(); | |
db16520a | 638 | } |
639 | ||
e67b0680 | 640 | Double_t det = sPrime*sxxPrime + sxPrime*sxPrime; |
db16520a | 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 | ||
e67b0680 | 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 | |
db16520a | 651 | |
e67b0680 | 652 | Double_t sigma2b = (Double_t)1. / sxxPrime; |
653 | //-- Double_t sigma2aprime = (Double_t)1. /sPrime; | |
db16520a | 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 | ||
e67b0680 | 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.; | |
db16520a | 712 | |
713 | Double_t chi2 = 0.; | |
e67b0680 | 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.; | |
db16520a | 719 | |
8fb3f728 | 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 | |
db16520a | 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(); | |
8fb3f728 | 731 | if(previousHit) { |
732 | // this means currentHit != fTrack->GetFirstHit() | |
db16520a | 733 | Double_t dx = currentHit->GetX() - previousHit->GetX() ; |
734 | Double_t dy = currentHit->GetY() - previousHit->GetY() ; | |
e67b0680 | 735 | slength = previousHit->GetS() - (Double_t)sqrt ( dx*dx + dy*dy ); |
db16520a | 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(); | |
e67b0680 | 740 | slength = (Double_t)sqrt ( dx*dx + dy*dy ); |
db16520a | 741 | } |
742 | ||
e67b0680 | 743 | currentHit->SetS(slength); |
db16520a | 744 | |
745 | S += currentHit->GetZWeight(); | |
746 | Ss += currentHit->GetZWeight() * currentHit->GetS(); | |
8fb3f728 | 747 | |
748 | // Matthias 23.02.2011 | |
749 | // adding missing assignment, otherwise previousHit stays NULL | |
750 | previousHit=currentHit; | |
db16520a | 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 | ||
f7561f8d | 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; | |
e67b0680 | 766 | szPrime += currentHit->GetZWeight() * currentHit->GetZ(); |
f7561f8d | 767 | sszPrime += currentHit->GetZWeight() * lsPrime * currentHit->GetZ(); |
db16520a | 768 | } |
769 | ||
e67b0680 | 770 | Double_t det = sPrime*sssPrime + ssPrime*ssPrime; |
db16520a | 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 | ||
e67b0680 | 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 | |
db16520a | 781 | |
782 | Double_t a = aPrime - b*sav; | |
783 | ||
e67b0680 | 784 | Double_t sigma2b = (Double_t) 1. / sssPrime; |
785 | Double_t sigma2aprime = (Double_t) 1. /sPrime; | |
db16520a | 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 |