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