]>
Commit | Line | Data |
---|---|---|
a6c02c85 | 1 | // @(#) $Id$ |
4aa41877 | 2 | // Original: AliHLTFitter.cxx,v 1.14 2005/06/14 10:55:21 cvetan |
a6c02c85 | 3 | |
6235cd38 | 4 | /************************************************************************** |
9be2600f | 5 | * This file is property of and copyright by the ALICE HLT Project * |
6 | * ALICE Experiment at CERN, All rights reserved. * | |
6235cd38 | 7 | * * |
9be2600f | 8 | * Primary Authors: Anders Vestbo, maintained by * |
9 | * Matthias Richter <Matthias.Richter@ift.uib.no> * | |
10 | * for The ALICE HLT Project. * | |
6235cd38 | 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 AliHLTTPCFitter.cxx | |
22 | @author Anders Vestbo, maintained by Matthias Richter | |
23 | @date | |
24 | @brief Fit class HLT for helix | |
a6c02c85 | 25 | */ |
26 | ||
27 | #include <math.h> | |
a6c02c85 | 28 | #include "AliHLTTPCLogging.h" |
29 | #include "AliHLTTPCFitter.h" | |
30 | #include "AliHLTTPCVertex.h" | |
31 | #include "AliHLTTPCTrack.h" | |
32 | #include "AliHLTTPCSpacePointData.h" | |
a6c02c85 | 33 | #include "AliHLTTPCTransform.h" |
a7c32962 | 34 | //#include "AliHLTTPC.h" |
a6c02c85 | 35 | |
36 | #if __GNUC__ >= 3 | |
37 | using namespace std; | |
38 | #endif | |
39 | ||
40 | ClassImp(AliHLTTPCFitter) | |
41 | ||
42 | ||
43 | AliHLTTPCFitter::AliHLTTPCFitter() | |
6235cd38 | 44 | : |
45 | fTrack(NULL), | |
46 | fVertex(NULL), | |
47 | fVertexConstraint(0) | |
a6c02c85 | 48 | { |
49 | //constructor | |
e419b223 | 50 | // see header file for class documentation |
51 | // or | |
52 | // refer to README to build package | |
53 | // or | |
54 | // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt | |
a6c02c85 | 55 | memset(fClusters,0,36*6*sizeof(AliHLTTPCSpacePointData*)); |
6235cd38 | 56 | memset(fNcl,0,36*6*sizeof(UInt_t)); |
57 | } | |
58 | ||
a6c02c85 | 59 | AliHLTTPCFitter::AliHLTTPCFitter(AliHLTTPCVertex *vertex,Bool_t vertexconstraint) |
e419b223 | 60 | : |
61 | fTrack(NULL), | |
62 | fVertex(vertex), | |
63 | fVertexConstraint(vertexconstraint) | |
a6c02c85 | 64 | { |
65 | //constructor | |
66 | fTrack=0; | |
a6c02c85 | 67 | memset(fClusters,0,36*6*sizeof(AliHLTTPCSpacePointData*)); |
68 | } | |
69 | ||
70 | AliHLTTPCFitter::~AliHLTTPCFitter() | |
71 | { | |
72 | //destructor | |
73 | for(Int_t i=0; i<36; i++) | |
74 | { | |
75 | for(Int_t j=0; j<6; j++) | |
76 | { | |
77 | if(fClusters[i][j]) | |
78 | delete [] fClusters[i][j]; | |
79 | } | |
80 | } | |
81 | } | |
82 | ||
a6c02c85 | 83 | void AliHLTTPCFitter::SortTrackClusters(AliHLTTPCTrack *track) const |
84 | { | |
85 | //Sort the internal cluster list in each track with respect to row numbering. | |
86 | //This may be necessary when no conventional track follower has been | |
87 | //applied, in which the cluster list has been maintained in a more | |
88 | //arbitrary fashion. | |
89 | ||
90 | Int_t nhits = track->GetNHits(); | |
91 | Int_t *ids = (Int_t*)track->GetHitNumbers(); | |
92 | Int_t *origids = new Int_t[nhits]; | |
93 | Int_t *mk = new Int_t[nhits]; | |
94 | Int_t k; | |
95 | ||
96 | for(k=0; k<nhits; k++) {origids[k] = ids[k]; mk[k] = -1;} | |
97 | ||
98 | Int_t slice,patch,id,padrow,maxrow,maxk; | |
99 | UInt_t pos; | |
100 | for(Int_t j=0; j<nhits; j++) | |
101 | { | |
102 | maxrow=-1; | |
103 | maxk=200; | |
104 | for(k=0; k<nhits; k++) | |
105 | { | |
106 | id=ids[k]; | |
107 | if(id < 0) continue; | |
a371a266 | 108 | slice = AliHLTTPCSpacePointData::GetSlice(id); |
109 | patch = AliHLTTPCSpacePointData::GetPatch(id); | |
110 | pos = AliHLTTPCSpacePointData::GetNumber(id); | |
a6c02c85 | 111 | AliHLTTPCSpacePointData *points = fClusters[slice][patch]; |
112 | padrow = points[pos].fPadRow; | |
113 | if(padrow > maxrow) | |
114 | { | |
115 | maxrow = padrow; | |
116 | maxk=k; | |
117 | } | |
118 | } | |
119 | mk[j]=maxk; | |
120 | ids[maxk]=-1; | |
121 | } | |
122 | ||
123 | for(k=0; k<nhits; k++) | |
124 | ids[k] = origids[mk[k]]; | |
125 | delete [] origids; | |
126 | delete [] mk; | |
127 | } | |
128 | ||
129 | Int_t AliHLTTPCFitter::FitHelix(AliHLTTPCTrack *track) | |
130 | { | |
131 | //fit helix parameters | |
132 | fTrack = track; | |
133 | if(FitCircle()) | |
134 | { | |
135 | LOG(AliHLTTPCLog::kError,"AliHLTTPCFitter::FitHelix","TrackFit")<<AliHLTTPCLog::kDec<< | |
136 | "Problems during circle fit"<<ENDLOG; | |
137 | return 1; | |
138 | } | |
139 | if(FitLine()) | |
140 | { | |
141 | LOG(AliHLTTPCLog::kError,"AliHLTTPCFitter::FitHelix","TrackFit")<<AliHLTTPCLog::kDec<< | |
142 | "Problems during line fit"<<ENDLOG; | |
143 | return 1; | |
144 | } | |
145 | return 0; | |
146 | } | |
147 | ||
148 | Int_t AliHLTTPCFitter::FitCircle() | |
149 | { | |
150 | //----------------------------------------------------------------- | |
151 | //Fits circle parameters using algorithm | |
152 | //described by ChErnov and Oskov in Computer Physics | |
153 | //Communications. | |
154 | // | |
155 | //Written in FORTRAN by Jawluen Tang, Physics department , UT-Austin | |
156 | //Moved to C by Pablo Yepes | |
157 | //Moved to AliROOT by ASV. | |
158 | //------------------------------------------------------------------ | |
379ff4f0 | 159 | |
160 | if (!fTrack) return -1; | |
a6c02c85 | 161 | |
162 | Double_t wsum = 0.0 ; | |
163 | Double_t xav = 0.0 ; | |
164 | Double_t yav = 0.0 ; | |
165 | ||
166 | // | |
167 | // Loop over hits calculating average | |
168 | Double_t * fXYWeight = new Double_t[(fTrack->GetNHits())]; | |
169 | UInt_t *hitnum = fTrack->GetHitNumbers(); | |
379ff4f0 | 170 | if (!fXYWeight || !hitnum) { |
171 | if (fXYWeight) delete [] fXYWeight; | |
172 | return -1; | |
173 | } | |
174 | memset(fXYWeight, 0, fTrack->GetNHits()*sizeof(Double_t)); | |
a6c02c85 | 175 | for(Int_t i=0; i<fTrack->GetNHits(); i++) |
176 | { | |
177 | UInt_t id = hitnum[i]; | |
a371a266 | 178 | Int_t slice = AliHLTTPCSpacePointData::GetSlice(id); |
179 | Int_t patch = AliHLTTPCSpacePointData::GetPatch(id); | |
180 | UInt_t pos = AliHLTTPCSpacePointData::GetNumber(id); | |
a6c02c85 | 181 | AliHLTTPCSpacePointData *points = fClusters[slice][patch]; |
182 | fXYWeight[i] = 1./ (Double_t)(points[pos].fSigmaY2 + points[pos].fSigmaY2); | |
183 | wsum += fXYWeight[i]; | |
184 | xav += fXYWeight[i]*points[pos].fX; | |
185 | yav += fXYWeight[i]*points[pos].fY; | |
186 | } | |
187 | if (fVertexConstraint == kTRUE) | |
188 | { | |
189 | wsum += fVertex->GetXYWeight() ; | |
190 | xav += fVertex->GetX() ; | |
191 | yav += fVertex->GetY() ; | |
192 | } | |
193 | ||
194 | xav = xav / wsum ; | |
195 | yav = yav / wsum ; | |
196 | // | |
197 | // CALCULATE <X**2>, <XY>, AND <Y**2> WITH <X> = 0, & <Y> = 0 | |
198 | // | |
199 | Double_t xxav = 0.0 ; | |
200 | Double_t xyav = 0.0 ; | |
201 | Double_t yyav = 0.0 ; | |
202 | Double_t xi, yi ; | |
203 | ||
204 | for(Int_t i=0; i<fTrack->GetNHits(); i++) | |
205 | { | |
206 | UInt_t id = hitnum[i]; | |
a371a266 | 207 | Int_t slice = AliHLTTPCSpacePointData::GetSlice(id); |
208 | Int_t patch = AliHLTTPCSpacePointData::GetPatch(id); | |
209 | UInt_t pos = AliHLTTPCSpacePointData::GetNumber(id); | |
a6c02c85 | 210 | AliHLTTPCSpacePointData *points = fClusters[slice][patch]; |
211 | ||
212 | xi = points[pos].fX -xav; | |
213 | yi = points[pos].fY - yav ; | |
214 | xxav += xi * xi * fXYWeight[i]; | |
215 | xyav += xi * yi * fXYWeight[i]; | |
216 | yyav += yi * yi * fXYWeight[i]; | |
217 | } | |
218 | ||
219 | if (fVertexConstraint == kTRUE) | |
220 | { | |
221 | xi = fVertex->GetX() - xav ; | |
222 | yi = fVertex->GetY() - yav ; | |
223 | xxav += xi * xi * fVertex->GetXYWeight() ; | |
224 | xyav += xi * yi * fVertex->GetXYWeight() ; | |
225 | yyav += yi * yi * fVertex->GetXYWeight() ; | |
226 | } | |
227 | xxav = xxav / wsum ; | |
228 | xyav = xyav / wsum ; | |
229 | yyav = yyav / wsum ; | |
230 | // | |
231 | //--> ROTATE COORDINATES SO THAT <XY> = 0 | |
232 | // | |
233 | //--> SIGN(C**2 - S**2) = SIGN(XXAV - YYAV) > | |
234 | //--> & > ==> NEW : (XXAV-YYAV) > 0 | |
235 | //--> SIGN(S) = SIGN(XYAV) > | |
236 | ||
237 | Double_t a = fabs( xxav - yyav ) ; | |
238 | Double_t b = 4.0 * xyav * xyav ; | |
239 | ||
240 | Double_t asqpb = a * a + b ; | |
241 | Double_t rasqpb = sqrt ( asqpb) ; | |
242 | ||
243 | Double_t splus = 1.0 + a / rasqpb ; | |
244 | Double_t sminus = b / (asqpb * splus) ; | |
245 | ||
246 | splus = sqrt (0.5 * splus ) ; | |
247 | sminus = sqrt (0.5 * sminus) ; | |
248 | // | |
249 | //-> FIRST REQUIRE : SIGN(C**2 - S**2) = SIGN(XXAV - YYAV) | |
250 | // | |
251 | Double_t sinrot, cosrot ; | |
252 | if ( xxav <= yyav ) { | |
253 | cosrot = sminus ; | |
254 | sinrot = splus ; | |
255 | } | |
256 | else { | |
257 | cosrot = splus ; | |
258 | sinrot = sminus ; | |
259 | } | |
260 | // | |
261 | //-> REQUIRE : SIGN(S) = SIGN(XYAV) * SIGN(C) (ASSUMING SIGN(C) > 0) | |
262 | // | |
263 | if ( xyav < 0.0 ) sinrot = - sinrot ; | |
264 | // | |
265 | //--> WE NOW HAVE THE SMALLEST ANGLE THAT GUARANTEES <X**2> > <Y**2> | |
266 | //--> TO GET THE SIGN OF THE CHARGE RIGHT, THE NEW X-AXIS MUST POINT | |
267 | //--> OUTWARD FROM THE ORGIN. WE ARE FREE TO CHANGE SIGNS OF BOTH | |
268 | //--> COSROT AND SINROT SIMULTANEOUSLY TO ACCOMPLISH THIS. | |
269 | // | |
270 | //--> CHOOSE SIGN OF C WISELY TO BE ABLE TO GET THE SIGN OF THE CHARGE | |
271 | // | |
272 | if ( cosrot*xav+sinrot*yav < 0.0 ) { | |
273 | cosrot = -cosrot ; | |
274 | sinrot = -sinrot ; | |
275 | } | |
276 | // | |
277 | //-> NOW GET <R**2> AND RSCALE= SQRT(<R**2>) | |
278 | // | |
279 | Double_t rrav = xxav + yyav ; | |
280 | Double_t rscale = sqrt(rrav) ; | |
281 | ||
282 | xxav = 0.0 ; | |
283 | yyav = 0.0 ; | |
284 | xyav = 0.0 ; | |
285 | Double_t xrrav = 0.0 ; | |
286 | Double_t yrrav = 0.0 ; | |
287 | Double_t rrrrav = 0.0 ; | |
288 | ||
289 | Double_t xixi, yiyi, riri, wiriri, xold, yold ; | |
290 | ||
291 | for(Int_t i=0; i<fTrack->GetNHits(); i++) | |
292 | { | |
293 | UInt_t id = hitnum[i]; | |
a371a266 | 294 | Int_t slice = AliHLTTPCSpacePointData::GetSlice(id); |
295 | Int_t patch = AliHLTTPCSpacePointData::GetPatch(id); | |
296 | UInt_t pos = AliHLTTPCSpacePointData::GetNumber(id); | |
a6c02c85 | 297 | AliHLTTPCSpacePointData *points = fClusters[slice][patch]; |
298 | ||
299 | xold = points[pos].fX - xav ; | |
300 | yold = points[pos].fY - yav ; | |
301 | // | |
302 | //--> ROTATE SO THAT <XY> = 0 & DIVIDE BY RSCALE SO THAT <R**2> = 1 | |
303 | // | |
304 | xi = ( cosrot * xold + sinrot * yold ) / rscale ; | |
305 | yi = ( -sinrot * xold + cosrot * yold ) / rscale ; | |
306 | ||
307 | xixi = xi * xi ; | |
308 | yiyi = yi * yi ; | |
309 | riri = xixi + yiyi ; | |
310 | wiriri = fXYWeight[i] * riri ; | |
311 | ||
312 | xyav += fXYWeight[i] * xi * yi ; | |
313 | xxav += fXYWeight[i] * xixi ; | |
314 | yyav += fXYWeight[i] * yiyi ; | |
315 | ||
316 | xrrav += wiriri * xi ; | |
317 | yrrav += wiriri * yi ; | |
318 | rrrrav += wiriri * riri ; | |
319 | } | |
320 | // | |
321 | // Include vertex if required | |
322 | // | |
323 | if (fVertexConstraint == kTRUE) | |
324 | { | |
325 | xold = fVertex->GetX() - xav ; | |
326 | yold = fVertex->GetY() - yav ; | |
327 | // | |
328 | //--> ROTATE SO THAT <XY> = 0 & DIVIDE BY RSCALE SO THAT <R**2> = 1 | |
329 | // | |
330 | xi = ( cosrot * xold + sinrot * yold ) / rscale ; | |
331 | yi = ( -sinrot * xold + cosrot * yold ) / rscale ; | |
332 | ||
333 | xixi = xi * xi ; | |
334 | yiyi = yi * yi ; | |
335 | riri = xixi + yiyi ; | |
336 | wiriri = fVertex->GetXYWeight() * riri ; | |
337 | ||
338 | xyav += fVertex->GetXYWeight() * xi * yi ; | |
339 | xxav += fVertex->GetXYWeight() * xixi ; | |
340 | yyav += fVertex->GetXYWeight() * yiyi ; | |
341 | ||
342 | xrrav += wiriri * xi ; | |
343 | yrrav += wiriri * yi ; | |
344 | rrrrav += wiriri * riri ; | |
345 | } | |
346 | // | |
347 | // | |
348 | // | |
349 | //--> DIVIDE BY WSUM TO MAKE AVERAGES | |
350 | // | |
351 | xxav = xxav / wsum ; | |
352 | yyav = yyav / wsum ; | |
353 | xrrav = xrrav / wsum ; | |
354 | yrrav = yrrav / wsum ; | |
355 | rrrrav = rrrrav / wsum ; | |
356 | xyav = xyav / wsum ; | |
357 | ||
358 | Int_t const kntry = 5 ; | |
359 | // | |
360 | //--> USE THESE TO GET THE COEFFICIENTS OF THE 4-TH ORDER POLYNIMIAL | |
361 | //--> DON'T PANIC - THE THIRD ORDER TERM IS ZERO ! | |
362 | // | |
363 | Double_t xrrxrr = xrrav * xrrav ; | |
364 | Double_t yrryrr = yrrav * yrrav ; | |
365 | Double_t rrrrm1 = rrrrav - 1.0 ; | |
366 | Double_t xxyy = xxav * yyav ; | |
367 | ||
368 | Double_t c0 = rrrrm1*xxyy - xrrxrr*yyav - yrryrr*xxav ; | |
369 | Double_t c1 = - rrrrm1 + xrrxrr + yrryrr - 4.0*xxyy ; | |
370 | Double_t c2 = 4.0 + rrrrm1 - 4.0*xxyy ; | |
371 | Double_t c4 = - 4.0 ; | |
372 | // | |
373 | //--> COEFFICIENTS OF THE DERIVATIVE - USED IN NEWTON-RAPHSON ITERATIONS | |
374 | // | |
375 | Double_t c2d = 2.0 * c2 ; | |
376 | Double_t c4d = 4.0 * c4 ; | |
377 | // | |
378 | //--> 0'TH VALUE OF LAMDA - LINEAR INTERPOLATION BETWEEN P(0) & P(YYAV) | |
379 | // | |
380 | // LAMDA = YYAV * C0 / (C0 + YRRSQ * (XXAV-YYAV)) | |
381 | Double_t lamda = 0.0 ; | |
382 | Double_t dlamda = 0.0 ; | |
383 | // | |
384 | Double_t chiscl = wsum * rscale * rscale ; | |
385 | Double_t dlamax = 0.001 / chiscl ; | |
386 | ||
387 | Double_t p, pd ; | |
388 | for ( int itry = 1 ; itry <= kntry ; itry++ ) { | |
389 | p = c0 + lamda * (c1 + lamda * (c2 + lamda * lamda * c4 )) ; | |
390 | pd = (c1 + lamda * (c2d + lamda * lamda * c4d)) ; | |
391 | dlamda = -p / pd ; | |
392 | lamda = lamda + dlamda ; | |
393 | if (fabs(dlamda)< dlamax) break ; | |
394 | } | |
395 | ||
396 | //Double_t chi2 = (Double_t)(chiscl * lamda) ; | |
397 | //fTrack->SetChiSq1(chi2); | |
398 | // Double_t dchisq = chiscl * dlamda ; | |
399 | // | |
400 | //--> NOW CALCULATE THE MATRIX ELEMENTS FOR ALPHA, BETA & KAPPA | |
401 | // | |
402 | Double_t h11 = xxav - lamda ; | |
403 | Double_t h14 = xrrav ; | |
404 | Double_t h22 = yyav - lamda ; | |
405 | Double_t h24 = yrrav ; | |
406 | Double_t h34 = 1.0 + 2.0*lamda ; | |
407 | if ( h11 == 0.0 || h22 == 0.0 ){ | |
408 | LOG(AliHLTTPCLog::kError,"AliHLTTPCFitter::FitCircle","TrackFit")<<AliHLTTPCLog::kDec<< | |
409 | "Problems fitting circle"<<ENDLOG; | |
379ff4f0 | 410 | delete [] fXYWeight; |
a6c02c85 | 411 | return 1 ; |
412 | } | |
413 | Double_t rootsq = (h14*h14)/(h11*h11) + 4.0*h34 ; | |
414 | ||
415 | Double_t ratio, kappa, beta ; | |
416 | if ( fabs(h22) > fabs(h24) ) { | |
417 | ratio = h24 / h22 ; | |
418 | rootsq = ratio * ratio + rootsq ; | |
419 | kappa = 1.0 / sqrt(rootsq) ; | |
420 | beta = - ratio * kappa ; | |
421 | } | |
422 | else { | |
423 | ratio = h22 / h24 ; | |
424 | rootsq = 1.0 + ratio * ratio * rootsq ; | |
425 | beta = 1.0 / sqrt(rootsq) ; | |
426 | if ( h24 > 0 ) beta = - beta ; | |
427 | kappa = -ratio * beta ; | |
428 | } | |
429 | Double_t alpha = - (h14/h11) * kappa ; | |
430 | // | |
431 | //--> transform these into the lab coordinate system | |
432 | //--> first get kappa and back to real dimensions | |
433 | // | |
434 | Double_t kappa1 = kappa / rscale ; | |
435 | Double_t dbro = 0.5 / kappa1 ; | |
436 | // | |
437 | //--> next rotate alpha and beta and scale | |
438 | // | |
439 | Double_t alphar = (cosrot * alpha - sinrot * beta)* dbro ; | |
440 | Double_t betar = (sinrot * alpha + cosrot * beta)* dbro ; | |
441 | // | |
442 | //--> then translate by (xav,yav) | |
443 | // | |
444 | Double_t acent = (double)(xav - alphar) ; | |
445 | Double_t bcent = (double)(yav - betar ) ; | |
446 | Double_t radius = (double)dbro ; | |
447 | // | |
448 | // Get charge | |
449 | // | |
450 | Int_t q = ( ( yrrav < 0 ) ? 1 : -1 ) ; | |
451 | fTrack->SetCharge(q); | |
452 | ||
453 | //Set the first point on the track to the space point coordinates of the innermost track | |
454 | //This will be updated to lie on the fit later on (AliHLTTPCTrack::UpdateToFirstPoint). | |
455 | Double_t x0,y0,psi,pt ; | |
456 | Int_t lastid=fTrack->GetNHits()-1; | |
457 | UInt_t id = hitnum[lastid]; | |
a371a266 | 458 | Int_t slice = AliHLTTPCSpacePointData::GetSlice(id); |
459 | Int_t patch = AliHLTTPCSpacePointData::GetPatch(id); | |
460 | UInt_t pos = AliHLTTPCSpacePointData::GetNumber(id); | |
a6c02c85 | 461 | AliHLTTPCSpacePointData *points = fClusters[slice][patch]; |
462 | x0 = points[pos].fX; | |
463 | y0 = points[pos].fY; | |
464 | fTrack->SetFirstPoint(x0,y0,0); //Z-value is set in FitLine | |
465 | ||
466 | //Set the remaining fit parameters | |
467 | psi = (Double_t)atan2(bcent-y0,acent-x0) ; | |
468 | psi = psi + q * 0.5F * AliHLTTPCTransform::Pi() ; | |
469 | if ( psi < 0 ) psi = psi + 2*AliHLTTPCTransform::Pi(); | |
470 | ||
471 | pt = (Double_t)(AliHLTTPCTransform::GetBFact() * AliHLTTPCTransform::GetBField() * radius ) ; | |
472 | fTrack->SetPsi(psi); | |
473 | fTrack->SetPt(pt); | |
474 | fTrack->SetRadius(radius); | |
475 | fTrack->SetCenterX(acent); | |
476 | fTrack->SetCenterY(bcent); | |
477 | // | |
478 | // Get errors from fast fit | |
479 | // | |
480 | //if ( getPara()->getErrors ) getErrorsCircleFit ( acent, bcent, radius ) ; | |
481 | // | |
482 | delete [] fXYWeight; | |
483 | return 0 ; | |
484 | } | |
485 | ||
486 | //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
487 | // Fit Line in s-z plane | |
488 | //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
489 | Int_t AliHLTTPCFitter::FitLine ( ) | |
490 | { | |
379ff4f0 | 491 | // Fit Line in s-z plane |
492 | if (!fTrack) return -1; | |
493 | ||
a6c02c85 | 494 | // |
495 | //Initialization | |
496 | // | |
497 | Double_t sum = 0.F ; | |
498 | Double_t ss = 0.F ; | |
499 | Double_t sz = 0.F ; | |
500 | Double_t sss = 0.F ; | |
501 | Double_t ssz = 0.F ; | |
502 | // | |
503 | //find sum , sums ,sumz, sumss | |
504 | // | |
505 | Double_t dx, dy ; | |
506 | Double_t radius = (Double_t)(fTrack->GetPt() / ( AliHLTTPCTransform::GetBFact() * AliHLTTPCTransform::GetBField() ) ) ; | |
507 | ||
508 | Double_t * fS = new Double_t[(fTrack->GetNHits())]; | |
509 | Double_t *fZWeight = new Double_t[fTrack->GetNHits()]; | |
510 | UInt_t *hitnum = fTrack->GetHitNumbers(); | |
379ff4f0 | 511 | if (!fS || !fZWeight || !hitnum) { |
512 | if (fS) delete [] fS; | |
513 | if (fZWeight) delete [] fZWeight; | |
514 | return -1; | |
515 | } | |
516 | memset(fS, 0, fTrack->GetNHits()*sizeof(Double_t)); | |
517 | memset(fZWeight, 0, fTrack->GetNHits()*sizeof(Double_t)); | |
a6c02c85 | 518 | if (0)//fVertexConstraint==kTRUE) |
519 | { | |
520 | UInt_t id = hitnum[0]; | |
a371a266 | 521 | Int_t slice = AliHLTTPCSpacePointData::GetSlice(id); |
522 | Int_t patch = AliHLTTPCSpacePointData::GetPatch(id); | |
523 | UInt_t pos = AliHLTTPCSpacePointData::GetNumber(id); | |
a6c02c85 | 524 | AliHLTTPCSpacePointData *points = fClusters[slice][patch]; |
525 | ||
526 | dx = points[pos].fX - fVertex->GetX(); | |
527 | dy = points[pos].fY - fVertex->GetY(); | |
528 | } | |
529 | else | |
530 | { | |
531 | UInt_t id = hitnum[0]; | |
a371a266 | 532 | Int_t slice = AliHLTTPCSpacePointData::GetSlice(id); |
533 | Int_t patch = AliHLTTPCSpacePointData::GetPatch(id); | |
534 | UInt_t posf = AliHLTTPCSpacePointData::GetNumber(id); | |
a6c02c85 | 535 | AliHLTTPCSpacePointData *pointsf = fClusters[slice][patch]; |
536 | id = hitnum[(fTrack->GetNHits()-1)]; | |
a371a266 | 537 | slice = AliHLTTPCSpacePointData::GetSlice(id); |
538 | patch = AliHLTTPCSpacePointData::GetPatch(id); | |
539 | UInt_t posl = AliHLTTPCSpacePointData::GetNumber(id); | |
a6c02c85 | 540 | AliHLTTPCSpacePointData *pointsl = fClusters[slice][patch]; |
541 | dx = pointsf[posf].fX - pointsl[posl].fX; | |
542 | dy = pointsf[posf].fY - pointsl[posl].fY; | |
543 | } | |
544 | ||
545 | Double_t localPsi = 0.5F * sqrt ( dx*dx + dy*dy ) / radius ; | |
546 | Double_t totals ; | |
547 | ||
548 | if ( fabs(localPsi) < 1. ) | |
549 | { | |
550 | totals = 2.0 * radius * asin ( localPsi ) ; | |
551 | } | |
552 | else | |
553 | { | |
554 | totals = 2.0 * radius * AliHLTTPCTransform::Pi() ; | |
555 | } | |
556 | ||
557 | Double_t dpsi,s; | |
558 | ||
559 | for(Int_t i=0; i<fTrack->GetNHits(); i++) | |
560 | { | |
561 | UInt_t id = hitnum[i]; | |
a371a266 | 562 | Int_t slice = AliHLTTPCSpacePointData::GetSlice(id); |
563 | Int_t patch = AliHLTTPCSpacePointData::GetPatch(id); | |
564 | UInt_t pos = AliHLTTPCSpacePointData::GetNumber(id); | |
a6c02c85 | 565 | AliHLTTPCSpacePointData *points = fClusters[slice][patch]; |
566 | ||
567 | fZWeight[i] = 1./(Double_t)(points[pos].fSigmaZ2); | |
568 | if(i>0) | |
569 | { | |
570 | id = hitnum[i-1]; | |
a371a266 | 571 | slice = AliHLTTPCSpacePointData::GetSlice(id); |
572 | patch = AliHLTTPCSpacePointData::GetPatch(id); | |
573 | UInt_t lastpos = AliHLTTPCSpacePointData::GetNumber(id); | |
a6c02c85 | 574 | AliHLTTPCSpacePointData *lastpoints = fClusters[slice][patch]; |
575 | dx = points[pos].fX -lastpoints[lastpos].fX; | |
576 | dy = points[pos].fY -lastpoints[lastpos].fY; | |
577 | dpsi = 0.5 * (Double_t)sqrt ( dx*dx + dy*dy ) / radius ; | |
379ff4f0 | 578 | if(fabs(dpsi) > 1) { |
579 | delete [] fS; | |
580 | delete [] fZWeight; | |
a6c02c85 | 581 | return 1; |
379ff4f0 | 582 | } |
a6c02c85 | 583 | fTrack->SetPsierr(dpsi); |
584 | s = fS[i-1] - 2.0 * radius * (Double_t)asin ( dpsi ) ; | |
585 | fS[i]=s; | |
586 | } | |
587 | else | |
588 | fS[i]=totals; | |
589 | ||
590 | sum += fZWeight[i]; | |
591 | ss += fZWeight[i] * fS[i]; | |
592 | sz += fZWeight[i] * points[pos].fZ; | |
593 | sss += fZWeight[i] * fS[i] * fS[i]; | |
594 | ssz += fZWeight[i] * fS[i] * points[pos].fZ; | |
595 | ||
596 | } | |
597 | ||
598 | ||
599 | Double_t chi2,det = sum * sss - ss * ss; | |
600 | if ( fabs(det) < 1e-20) | |
601 | { | |
602 | chi2 = 99999.F ; | |
603 | //fTrack->SetChiSq2(chi2); | |
379ff4f0 | 604 | delete [] fS; |
605 | delete [] fZWeight; | |
a6c02c85 | 606 | return 0 ; |
607 | } | |
608 | ||
609 | //Compute the best fitted parameters A,B | |
610 | Double_t tanl,z0,dtanl,dz0; | |
611 | ||
612 | tanl = (Double_t)((sum * ssz - ss * sz ) / det ); | |
613 | z0 = (Double_t)((sz * sss - ssz * ss ) / det ); | |
614 | ||
615 | fTrack->SetTgl(tanl); | |
616 | fTrack->SetZ0(z0); | |
617 | ||
618 | //calculate chi-square | |
619 | chi2 = 0.; | |
620 | Double_t r1 ; | |
621 | ||
622 | for(Int_t i=0; i<fTrack->GetNHits(); i++) | |
623 | { | |
624 | UInt_t id = hitnum[i]; | |
a371a266 | 625 | Int_t slice = AliHLTTPCSpacePointData::GetSlice(id); |
626 | Int_t patch = AliHLTTPCSpacePointData::GetPatch(id); | |
627 | UInt_t pos = AliHLTTPCSpacePointData::GetNumber(id); | |
a6c02c85 | 628 | AliHLTTPCSpacePointData *points = fClusters[slice][patch]; |
629 | r1 = points[pos].fZ - tanl * fS[i] - z0 ; | |
630 | chi2 += (Double_t) ( (Double_t)(fZWeight[i]) * (r1 * r1) ); | |
631 | } | |
632 | ||
633 | //fTrack->SetChiSq2(chi2); | |
634 | // | |
635 | //calculate estimated variance | |
636 | //varsq=chi/(double(n)-2.) | |
637 | //calculate covariance matrix | |
638 | //siga=sqrt(varsq*sxx/det) | |
639 | //sigb=sqrt(varsq*sum/det) | |
640 | // | |
641 | dtanl = (Double_t) ( sum / det ); | |
642 | dz0 = (Double_t) ( sss / det ); | |
643 | ||
644 | fTrack->SetTglerr(dtanl); | |
645 | fTrack->SetZ0err(dz0); | |
646 | delete [] fZWeight; | |
647 | delete [] fS; | |
648 | return 0 ; | |
649 | } |