]>
Commit | Line | Data |
---|---|---|
75bd7f81 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | /* $Id$ */ | |
17 | ||
18 | /////////////////////////////////////////////////////////////////////////////// | |
19 | // // | |
20 | // The seed of a local TRD track // | |
21 | // // | |
75bd7f81 | 22 | /////////////////////////////////////////////////////////////////////////////// |
23 | ||
ad670fba | 24 | #include "TMath.h" |
25 | #include "TLinearFitter.h" | |
26 | ||
75bd7f81 | 27 | #include "AliTRDseed.h" |
28 | #include "AliTRDcluster.h" | |
29 | #include "AliTRDtracker.h" | |
30 | ||
75bd7f81 | 31 | ClassImp(AliTRDseed) |
32 | ||
33 | //_____________________________________________________________________________ | |
ad670fba | 34 | AliTRDseed::AliTRDseed() |
35 | :TObject() | |
36 | ,fTilt(0) | |
37 | ,fPadLength(0) | |
38 | ,fX0(0) | |
39 | ,fSigmaY(0) | |
40 | ,fSigmaY2(0) | |
41 | ,fMeanz(0) | |
42 | ,fZProb(0) | |
43 | ,fN(0) | |
44 | ,fN2(0) | |
45 | ,fNUsed(0) | |
46 | ,fFreq(0) | |
47 | ,fNChange(0) | |
48 | ,fMPads(0) | |
49 | ,fC(0) | |
50 | ,fCC(0) | |
51 | ,fChi2(1.0e10) | |
52 | ,fChi2Z(1.0e10) | |
75bd7f81 | 53 | { |
54 | // | |
55 | // Default constructor | |
56 | // | |
57 | ||
58 | for (Int_t i = 0; i < 25; i++) { | |
ad670fba | 59 | fX[i] = 0; // x position |
60 | fY[i] = 0; // y position | |
61 | fZ[i] = 0; // z position | |
62 | fIndexes[i] = 0; // Indexes | |
63 | fClusters[i] = 0x0; // Clusters | |
64 | fUsable[i] = 0; // Indication - usable cluster | |
65 | } | |
66 | ||
67 | for (Int_t i = 0; i < 2; i++) { | |
68 | fYref[i] = 0; // Reference y | |
69 | fZref[i] = 0; // Reference z | |
70 | fYfit[i] = 0; // Y fit position +derivation | |
71 | fYfitR[i] = 0; // Y fit position +derivation | |
72 | fZfit[i] = 0; // Z fit position | |
73 | fZfitR[i] = 0; // Z fit position | |
74 | fLabels[i] = 0; // Labels | |
75 | } | |
76 | ||
77 | } | |
78 | ||
79 | //_____________________________________________________________________________ | |
80 | AliTRDseed::AliTRDseed(const AliTRDseed &s) | |
81 | :TObject(s) | |
82 | ,fTilt(0) | |
83 | ,fPadLength(0) | |
84 | ,fX0(0) | |
85 | ,fSigmaY(0) | |
86 | ,fSigmaY2(0) | |
87 | ,fMeanz(0) | |
88 | ,fZProb(0) | |
89 | ,fN(0) | |
90 | ,fN2(0) | |
91 | ,fNUsed(0) | |
92 | ,fFreq(0) | |
93 | ,fNChange(0) | |
94 | ,fMPads(0) | |
95 | ,fC(0) | |
96 | ,fCC(0) | |
97 | ,fChi2(0) | |
98 | ,fChi2Z(0) | |
99 | { | |
100 | // | |
101 | // Copy constructor | |
102 | // | |
103 | ||
104 | for (Int_t i = 0; i < 25; i++) { | |
105 | fX[i] = 0; // x position | |
106 | fY[i] = 0; // y position | |
107 | fZ[i] = 0; // z position | |
108 | fIndexes[i] = 0; // Indexes | |
109 | fClusters[i] = 0x0; // Clusters | |
110 | fUsable[i] = 0; // Indication - usable cluster | |
75bd7f81 | 111 | } |
ad670fba | 112 | |
75bd7f81 | 113 | for (Int_t i = 0; i < 2; i++) { |
ad670fba | 114 | fYref[i] = 0; // Reference y |
115 | fZref[i] = 0; // Reference z | |
116 | fYfit[i] = 0; // Y fit position +derivation | |
117 | fYfitR[i] = 0; // Y fit position +derivation | |
118 | fZfit[i] = 0; // Z fit position | |
119 | fZfitR[i] = 0; // Z fit position | |
120 | fLabels[i] = 0; // Labels | |
75bd7f81 | 121 | } |
122 | ||
123 | } | |
124 | ||
125 | //_____________________________________________________________________________ | |
126 | void AliTRDseed::Reset() | |
127 | { | |
128 | // | |
129 | // Reset seed | |
130 | // | |
131 | ||
132 | for (Int_t i = 0; i < 25; i++) { | |
ad670fba | 133 | fX[i] = 0; // X position |
134 | fY[i] = 0; // Y position | |
135 | fZ[i] = 0; // Z position | |
136 | fIndexes[i] = 0; // Indexes | |
137 | fClusters[i] = 0; // Clusters | |
75bd7f81 | 138 | fUsable[i] = kFALSE; |
139 | } | |
ad670fba | 140 | |
75bd7f81 | 141 | for (Int_t i = 0; i < 2; i++) { |
ad670fba | 142 | fYref[i] = 0; // Reference y |
143 | fZref[i] = 0; // Reference z | |
144 | fYfit[i] = 0; // Y fit position +derivation | |
145 | fYfitR[i] = 0; // Y fit position +derivation | |
146 | fZfit[i] = 0; // Z fit position | |
147 | fZfitR[i] = 0; // Z fit position | |
148 | fLabels[i] = -1; // Labels | |
75bd7f81 | 149 | } |
ad670fba | 150 | fSigmaY = 0; // "Robust" sigma in y |
151 | fSigmaY2 = 0; // "Robust" sigma in y | |
152 | fMeanz = 0; // Mean vaue of z | |
153 | fZProb = 0; // Max probbable z | |
75bd7f81 | 154 | fMPads = 0; |
ad670fba | 155 | fN = 0; // Number of associated clusters |
156 | fN2 = 0; // Number of not crossed | |
157 | fNUsed = 0; // Number of used clusters | |
158 | fNChange = 0; // Change z counter | |
75bd7f81 | 159 | |
160 | } | |
161 | ||
162 | //_____________________________________________________________________________ | |
163 | void AliTRDseed::CookLabels() | |
164 | { | |
165 | // | |
166 | // Cook 2 labels for seed | |
167 | // | |
168 | ||
169 | Int_t labels[200]; | |
170 | Int_t out[200]; | |
ad670fba | 171 | Int_t nlab = 0; |
172 | ||
173 | for (Int_t i = 0; i < 25; i++) { | |
75bd7f81 | 174 | if (!fClusters[i]) continue; |
ad670fba | 175 | for (Int_t ilab = 0; ilab < 3; ilab++) { |
176 | if (fClusters[i]->GetLabel(ilab) >= 0) { | |
75bd7f81 | 177 | labels[nlab] = fClusters[i]->GetLabel(ilab); |
178 | nlab++; | |
179 | } | |
180 | } | |
181 | } | |
ad670fba | 182 | |
75bd7f81 | 183 | Int_t nlab2 = AliTRDtracker::Freq(nlab,labels,out,kTRUE); |
184 | fLabels[0] = out[0]; | |
ad670fba | 185 | if ((nlab2 > 1) && |
186 | (out[3] > 1)) { | |
187 | fLabels[1] = out[2]; | |
188 | } | |
75bd7f81 | 189 | |
190 | } | |
191 | ||
192 | //_____________________________________________________________________________ | |
193 | void AliTRDseed::UseClusters() | |
194 | { | |
195 | // | |
196 | // Use clusters | |
197 | // | |
198 | ||
ad670fba | 199 | for (Int_t i = 0; i < 25; i++) { |
75bd7f81 | 200 | if (!fClusters[i]) continue; |
201 | if (!(fClusters[i]->IsUsed())) fClusters[i]->Use(); | |
202 | } | |
203 | ||
204 | } | |
205 | ||
206 | //_____________________________________________________________________________ | |
207 | void AliTRDseed::Update() | |
208 | { | |
209 | // | |
210 | // Update the seed. | |
211 | // | |
212 | ||
213 | const Float_t kRatio = 0.8; | |
214 | const Int_t kClmin = 6; | |
215 | const Float_t kmaxtan = 2; | |
216 | ||
ad670fba | 217 | if (TMath::Abs(fYref[1]) > kmaxtan) return; // Track inclined too much |
218 | ||
219 | Float_t sigmaexp = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction | |
220 | Float_t ycrosscor = fPadLength * fTilt * 0.5; // Y correction for crossing | |
221 | fNChange = 0; | |
222 | ||
223 | Double_t sumw; | |
224 | Double_t sumwx; | |
225 | Double_t sumwx2; | |
226 | Double_t sumwy; | |
227 | Double_t sumwxy; | |
228 | Double_t sumwz; | |
229 | Double_t sumwxz; | |
230 | ||
231 | Int_t zints[25]; // Histograming of the z coordinate | |
232 | // Get 1 and second max probable coodinates in z | |
233 | Int_t zouts[50]; | |
234 | Float_t allowedz[25]; // Allowed z for given time bin | |
235 | Float_t yres[25]; // Residuals from reference | |
236 | Float_t anglecor = fTilt * fZref[1]; // Correction to the angle | |
75bd7f81 | 237 | |
238 | ||
ad670fba | 239 | fN = 0; |
240 | fN2 = 0; | |
241 | for (Int_t i = 0; i < 25; i++) { | |
242 | yres[i] = 10000.0; | |
75bd7f81 | 243 | if (!fClusters[i]) continue; |
ad670fba | 244 | yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i]; // Residual y |
75bd7f81 | 245 | zints[fN] = Int_t(fZ[i]); |
246 | fN++; | |
247 | } | |
248 | ||
ad670fba | 249 | if (fN < kClmin) return; |
75bd7f81 | 250 | Int_t nz = AliTRDtracker::Freq(fN,zints,zouts,kFALSE); |
251 | fZProb = zouts[0]; | |
ad670fba | 252 | if (nz <= 1) zouts[3] = 0; |
253 | if (zouts[1] + zouts[3] < kClmin) return; | |
75bd7f81 | 254 | |
ad670fba | 255 | // Z distance bigger than pad - length |
256 | if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) { | |
257 | zouts[3]=0; | |
258 | } | |
75bd7f81 | 259 | |
260 | Int_t breaktime = -1; | |
261 | Bool_t mbefore = kFALSE; | |
262 | Int_t cumul[25][2]; | |
ad670fba | 263 | Int_t counts[2] = { 0, 0 }; |
75bd7f81 | 264 | |
ad670fba | 265 | if (zouts[3] >= 3) { |
266 | ||
75bd7f81 | 267 | // |
ad670fba | 268 | // Find the break time allowing one chage on pad-rows |
269 | // with maximal numebr of accepted clusters | |
75bd7f81 | 270 | // |
ad670fba | 271 | fNChange = 1; |
272 | for (Int_t i = 0; i < 25; i++) { | |
75bd7f81 | 273 | cumul[i][0] = counts[0]; |
274 | cumul[i][1] = counts[1]; | |
ad670fba | 275 | if (TMath::Abs(fZ[i]-zouts[0]) < 2) counts[0]++; |
276 | if (TMath::Abs(fZ[i]-zouts[2]) < 2) counts[1]++; | |
75bd7f81 | 277 | } |
ad670fba | 278 | Int_t maxcount = 0; |
279 | for (Int_t i = 0; i < 24; i++) { | |
280 | Int_t after = cumul[24][0] - cumul[i][0]; | |
75bd7f81 | 281 | Int_t before = cumul[i][1]; |
ad670fba | 282 | if (after + before > maxcount) { |
283 | maxcount = after + before; | |
284 | breaktime = i; | |
285 | mbefore = kFALSE; | |
75bd7f81 | 286 | } |
ad670fba | 287 | after = cumul[24][1] - cumul[i][1]; |
75bd7f81 | 288 | before = cumul[i][0]; |
ad670fba | 289 | if (after + before > maxcount) { |
290 | maxcount = after + before; | |
291 | breaktime = i; | |
292 | mbefore = kTRUE; | |
75bd7f81 | 293 | } |
294 | } | |
ad670fba | 295 | |
296 | breaktime -= 1; | |
297 | ||
75bd7f81 | 298 | } |
ad670fba | 299 | |
300 | for (Int_t i = 0; i < 25; i++) { | |
301 | if (i > breaktime) allowedz[i] = mbefore ? zouts[2] : zouts[0]; | |
302 | if (i <= breaktime) allowedz[i] = (!mbefore) ? zouts[2] : zouts[0]; | |
75bd7f81 | 303 | } |
ad670fba | 304 | |
305 | if (((allowedz[0] > allowedz[24]) && (fZref[1] < 0)) || | |
306 | ((allowedz[0] < allowedz[24]) && (fZref[1] > 0))) { | |
75bd7f81 | 307 | // |
ad670fba | 308 | // Tracklet z-direction not in correspondance with track z direction |
75bd7f81 | 309 | // |
ad670fba | 310 | fNChange = 0; |
311 | for (Int_t i = 0; i < 25; i++) { | |
312 | allowedz[i] = zouts[0]; // Only longest taken | |
75bd7f81 | 313 | } |
314 | } | |
315 | ||
ad670fba | 316 | if (fNChange > 0) { |
75bd7f81 | 317 | // |
ad670fba | 318 | // Cross pad -row tracklet - take the step change into account |
75bd7f81 | 319 | // |
ad670fba | 320 | for (Int_t i = 0; i < 25; i++) { |
75bd7f81 | 321 | if (!fClusters[i]) continue; |
ad670fba | 322 | if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue; |
323 | yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i]; // Residual y | |
324 | if (TMath::Abs(fZ[i] - fZProb) > 2) { | |
325 | if (fZ[i] > fZProb) yres[i] += fTilt * fPadLength; | |
326 | if (fZ[i] < fZProb) yres[i] -= fTilt * fPadLength; | |
75bd7f81 | 327 | } |
328 | } | |
329 | } | |
330 | ||
331 | Double_t yres2[25]; | |
ad670fba | 332 | Double_t mean; |
333 | Double_t sigma; | |
334 | for (Int_t i = 0; i < 25; i++) { | |
75bd7f81 | 335 | if (!fClusters[i]) continue; |
ad670fba | 336 | if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue; |
337 | yres2[fN2] = yres[i]; | |
75bd7f81 | 338 | fN2++; |
339 | } | |
ad670fba | 340 | if (fN2 < kClmin) { |
75bd7f81 | 341 | fN2 = 0; |
342 | return; | |
343 | } | |
344 | EvaluateUni(fN2,yres2,mean,sigma,Int_t(fN2*kRatio-2)); | |
ad670fba | 345 | if (sigma < sigmaexp * 0.8) { |
346 | sigma = sigmaexp; | |
347 | } | |
75bd7f81 | 348 | fSigmaY = sigma; |
ad670fba | 349 | |
350 | // Reset sums | |
351 | sumw = 0; | |
352 | sumwx = 0; | |
353 | sumwx2 = 0; | |
354 | sumwy = 0; | |
355 | sumwxy = 0; | |
356 | sumwz = 0; | |
357 | sumwxz = 0; | |
358 | ||
359 | fN2 = 0; | |
360 | fMeanz = 0; | |
361 | fMPads = 0; | |
362 | ||
363 | for (Int_t i = 0; i < 25; i++) { | |
364 | ||
365 | fUsable[i] = kFALSE; | |
75bd7f81 | 366 | if (!fClusters[i]) continue; |
ad670fba | 367 | if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue; |
368 | if (TMath::Abs(yres[i] - mean) > 4.0 * sigma) continue; | |
75bd7f81 | 369 | fUsable[i] = kTRUE; |
370 | fN2++; | |
ad670fba | 371 | fMPads += fClusters[i]->GetNPads(); |
372 | Float_t weight = 1.0; | |
373 | if (fClusters[i]->GetNPads() > 4) weight = 0.5; | |
374 | if (fClusters[i]->GetNPads() > 5) weight = 0.2; | |
375 | ||
75bd7f81 | 376 | Double_t x = fX[i]; |
ad670fba | 377 | sumw += weight; |
378 | sumwx += x * weight; | |
379 | sumwx2 += x*x * weight; | |
380 | sumwy += weight * yres[i]; | |
381 | sumwxy += weight * (yres[i]) * x; | |
382 | sumwz += weight * fZ[i]; | |
383 | sumwxz += weight * fZ[i] * x; | |
384 | ||
75bd7f81 | 385 | } |
ad670fba | 386 | |
387 | if (fN2 < kClmin){ | |
75bd7f81 | 388 | fN2 = 0; |
389 | return; | |
390 | } | |
ad670fba | 391 | fMeanz = sumwz / sumw; |
392 | Float_t correction = 0; | |
393 | if (fNChange > 0) { | |
394 | // Tracklet on boundary | |
395 | if (fMeanz < fZProb) correction = ycrosscor; | |
396 | if (fMeanz > fZProb) correction = -ycrosscor; | |
75bd7f81 | 397 | } |
ad670fba | 398 | |
399 | Double_t det = sumw * sumwx2 - sumwx * sumwx; | |
400 | fYfitR[0] = (sumwx2 * sumwy - sumwx * sumwxy) / det; | |
401 | fYfitR[1] = (sumw * sumwxy - sumwx * sumwy) / det; | |
75bd7f81 | 402 | |
ad670fba | 403 | fSigmaY2 = 0; |
404 | for (Int_t i = 0; i < 25; i++) { | |
75bd7f81 | 405 | if (!fUsable[i]) continue; |
ad670fba | 406 | Float_t delta = yres[i] - fYfitR[0] - fYfitR[1] * fX[i]; |
407 | fSigmaY2 += delta*delta; | |
75bd7f81 | 408 | } |
ad670fba | 409 | fSigmaY2 = TMath::Sqrt(fSigmaY2 / Float_t(fN2-2)); |
410 | ||
411 | fZfitR[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det; | |
412 | fZfitR[1] = (sumw * sumwxz - sumwx * sumwz) / det; | |
413 | fZfit[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det; | |
414 | fZfit[1] = (sumw * sumwxz - sumwx * sumwz) / det; | |
415 | fYfitR[0] += fYref[0] + correction; | |
416 | fYfitR[1] += fYref[1]; | |
417 | fYfit[0] = fYfitR[0]; | |
418 | fYfit[1] = fYfitR[1]; | |
75bd7f81 | 419 | |
420 | UpdateUsed(); | |
421 | ||
422 | } | |
423 | ||
424 | //_____________________________________________________________________________ | |
425 | void AliTRDseed::UpdateUsed() | |
426 | { | |
427 | // | |
428 | // Update used seed | |
429 | // | |
430 | ||
431 | fNUsed = 0; | |
ad670fba | 432 | for (Int_t i = 0; i < 25; i++) { |
433 | if (!fClusters[i]) continue; | |
434 | if ((fClusters[i]->IsUsed())) fNUsed++; | |
75bd7f81 | 435 | } |
436 | ||
437 | } | |
438 | ||
439 | //_____________________________________________________________________________ | |
440 | void AliTRDseed::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean | |
441 | , Double_t &sigma, Int_t hh) | |
442 | { | |
443 | // | |
444 | // Robust estimator in 1D case MI version | |
445 | // | |
446 | // For the univariate case | |
447 | // estimates of location and scatter are returned in mean and sigma parameters | |
448 | // the algorithm works on the same principle as in multivariate case - | |
449 | // it finds a subset of size hh with smallest sigma, and then returns mean and | |
450 | // sigma of this subset | |
451 | // | |
452 | ||
ad670fba | 453 | if (hh == 0) { |
454 | hh = (nvectors + 2) / 2; | |
455 | } | |
456 | ||
457 | Double_t faclts[] = { 2.6477, 2.5092, 2.3826, 2.2662, 2.1587 | |
458 | , 2.0589, 1.9660, 1.879, 1.7973, 1.7203 | |
459 | , 1.6473 }; | |
460 | Int_t *index = new Int_t[nvectors]; | |
75bd7f81 | 461 | TMath::Sort(nvectors, data, index, kFALSE); |
462 | ||
ad670fba | 463 | Int_t nquant = TMath::Min(Int_t(Double_t(((hh * 1.0 / nvectors) - 0.5) * 40)) + 1,11); |
75bd7f81 | 464 | Double_t factor = faclts[nquant-1]; |
465 | ||
ad670fba | 466 | Double_t sumx = 0.0; |
467 | Double_t sumx2 = 0.0; | |
75bd7f81 | 468 | Int_t bestindex = -1; |
ad670fba | 469 | Double_t bestmean = 0.0; |
470 | Double_t bestsigma = data[index[nvectors-1]] - data[index[0]]; // Maximal possible sigma | |
471 | for (Int_t i = 0; i < hh; i++) { | |
75bd7f81 | 472 | sumx += data[index[i]]; |
473 | sumx2 += data[index[i]]*data[index[i]]; | |
474 | } | |
475 | ||
ad670fba | 476 | Double_t norm = 1.0 / Double_t(hh); |
477 | Double_t norm2 = 1.0 / Double_t(hh - 1); | |
478 | for (Int_t i = hh; i < nvectors; i++) { | |
479 | ||
75bd7f81 | 480 | Double_t cmean = sumx*norm; |
ad670fba | 481 | Double_t csigma = (sumx2 - hh*cmean*cmean) * norm2; |
482 | if (csigma < bestsigma) { | |
75bd7f81 | 483 | bestmean = cmean; |
484 | bestsigma = csigma; | |
ad670fba | 485 | bestindex = i - hh; |
75bd7f81 | 486 | } |
487 | ||
ad670fba | 488 | sumx += data[index[i]] - data[index[i-hh]]; |
489 | sumx2 += data[index[i]]*data[index[i]] - data[index[i-hh]]*data[index[i-hh]]; | |
490 | ||
75bd7f81 | 491 | } |
492 | ||
ad670fba | 493 | Double_t bstd = factor * TMath::Sqrt(TMath::Abs(bestsigma)); |
75bd7f81 | 494 | mean = bestmean; |
495 | sigma = bstd; | |
ad670fba | 496 | |
75bd7f81 | 497 | delete [] index; |
498 | ||
499 | } | |
500 | ||
501 | //_____________________________________________________________________________ | |
502 | Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror) | |
503 | { | |
504 | // | |
505 | // Fit the Rieman tilt | |
506 | // | |
507 | ||
ad670fba | 508 | // Fitting with tilting pads - kz not fixed |
509 | TLinearFitter fitterT2(4,"hyp4"); | |
75bd7f81 | 510 | fitterT2.StoreData(kTRUE); |
ad670fba | 511 | Float_t xref2 = (cseed[2].fX0 + cseed[3].fX0) * 0.5; // Reference x0 for z |
512 | ||
513 | Int_t npointsT = 0; | |
75bd7f81 | 514 | fitterT2.ClearPoints(); |
ad670fba | 515 | |
516 | for (Int_t iLayer = 0; iLayer < 6; iLayer++) { | |
517 | ||
75bd7f81 | 518 | if (!cseed[iLayer].IsOK()) continue; |
519 | Double_t tilt = cseed[iLayer].fTilt; | |
520 | ||
ad670fba | 521 | for (Int_t itime = 0; itime < 25; itime++) { |
522 | ||
75bd7f81 | 523 | if (!cseed[iLayer].fUsable[itime]) continue; |
524 | // x relative to the midle chamber | |
ad670fba | 525 | Double_t x = cseed[iLayer].fX[itime] + cseed[iLayer].fX0 - xref2; |
526 | Double_t y = cseed[iLayer].fY[itime]; | |
527 | Double_t z = cseed[iLayer].fZ[itime]; | |
75bd7f81 | 528 | |
529 | // | |
ad670fba | 530 | // Tilted rieman |
75bd7f81 | 531 | // |
532 | Double_t uvt[6]; | |
ad670fba | 533 | Double_t x2 = cseed[iLayer].fX[itime] + cseed[iLayer].fX0; // Global x |
534 | Double_t t = 1.0 / (x2*x2 + y*y); | |
535 | uvt[1] = t; | |
536 | uvt[0] = 2.0 * x2 * uvt[1]; | |
537 | uvt[2] = 2.0 * tilt * uvt[1]; | |
538 | uvt[3] = 2.0 * tilt *uvt[1] * x; | |
539 | uvt[4] = 2.0 * (y + tilt * z) * uvt[1]; | |
540 | ||
541 | Double_t error = 2.0 * uvt[1]; | |
542 | if (terror) { | |
543 | error *= cseed[iLayer].fSigmaY; | |
544 | } | |
545 | else { | |
546 | error *= 0.2; //Default error | |
547 | } | |
75bd7f81 | 548 | fitterT2.AddPoint(uvt,uvt[4],error); |
549 | npointsT++; | |
ad670fba | 550 | |
75bd7f81 | 551 | } |
ad670fba | 552 | |
75bd7f81 | 553 | } |
ad670fba | 554 | |
75bd7f81 | 555 | fitterT2.Eval(); |
556 | Double_t rpolz0 = fitterT2.GetParameter(3); | |
557 | Double_t rpolz1 = fitterT2.GetParameter(4); | |
ad670fba | 558 | |
75bd7f81 | 559 | // |
ad670fba | 560 | // Linear fitter - not possible to make boundaries |
75bd7f81 | 561 | // non accept non possible z and dzdx combination |
562 | // | |
ad670fba | 563 | Bool_t acceptablez = kTRUE; |
564 | for (Int_t iLayer = 0; iLayer < 6; iLayer++) { | |
565 | if (cseed[iLayer].IsOK()) { | |
566 | Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].fX0 - xref2); | |
567 | if (TMath::Abs(cseed[iLayer].fZProb - zT2) > cseed[iLayer].fPadLength * 0.5 + 1.0) { | |
75bd7f81 | 568 | acceptablez = kFALSE; |
ad670fba | 569 | } |
75bd7f81 | 570 | } |
571 | } | |
ad670fba | 572 | if (!acceptablez) { |
573 | Double_t zmf = cseed[2].fZref[0] + cseed[2].fZref[1] * (xref2 - cseed[2].fX0); | |
574 | Double_t dzmf = (cseed[2].fZref[1] + cseed[3].fZref[1]) * 0.5; | |
75bd7f81 | 575 | fitterT2.FixParameter(3,zmf); |
576 | fitterT2.FixParameter(4,dzmf); | |
577 | fitterT2.Eval(); | |
578 | fitterT2.ReleaseParameter(3); | |
579 | fitterT2.ReleaseParameter(4); | |
580 | rpolz0 = fitterT2.GetParameter(3); | |
581 | rpolz1 = fitterT2.GetParameter(4); | |
582 | } | |
583 | ||
ad670fba | 584 | Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT); |
75bd7f81 | 585 | Double_t params[3]; |
ad670fba | 586 | params[0] = fitterT2.GetParameter(0); |
587 | params[1] = fitterT2.GetParameter(1); | |
588 | params[2] = fitterT2.GetParameter(2); | |
589 | Double_t curvature = 1.0 + params[1] * params[1] - params[2] * params[0]; | |
590 | ||
591 | for (Int_t iLayer = 0; iLayer < 6; iLayer++) { | |
592 | ||
593 | Double_t x = cseed[iLayer].fX0; | |
594 | Double_t y = 0; | |
595 | Double_t dy = 0; | |
596 | Double_t z = 0; | |
597 | Double_t dz = 0; | |
598 | ||
75bd7f81 | 599 | // y |
ad670fba | 600 | Double_t res2 = (x * params[0] + params[1]); |
601 | res2 *= res2; | |
602 | res2 = 1.0 - params[2]*params[0] + params[1]*params[1] - res2; | |
603 | if (res2 >= 0) { | |
75bd7f81 | 604 | res2 = TMath::Sqrt(res2); |
ad670fba | 605 | y = (1.0 - res2) / params[0]; |
75bd7f81 | 606 | } |
ad670fba | 607 | |
75bd7f81 | 608 | //dy |
ad670fba | 609 | Double_t x0 = -params[1] / params[0]; |
610 | if (-params[2]*params[0] + params[1]*params[1] + 1 > 0) { | |
611 | Double_t rm1 = params[0] / TMath::Sqrt(-params[2]*params[0] + params[1]*params[1] + 1); | |
612 | if (1.0/(rm1*rm1) - (x-x0) * (x-x0) > 0.0) { | |
613 | Double_t res = (x - x0) / TMath::Sqrt(1.0 / (rm1*rm1) - (x-x0)*(x-x0)); | |
614 | if (params[0] < 0) res *= -1.0; | |
75bd7f81 | 615 | dy = res; |
616 | } | |
617 | } | |
ad670fba | 618 | z = rpolz0 + rpolz1 * (x - xref2); |
75bd7f81 | 619 | dz = rpolz1; |
620 | cseed[iLayer].fYref[0] = y; | |
621 | cseed[iLayer].fYref[1] = dy; | |
622 | cseed[iLayer].fZref[0] = z; | |
623 | cseed[iLayer].fZref[1] = dz; | |
ad670fba | 624 | cseed[iLayer].fC = curvature; |
75bd7f81 | 625 | |
626 | } | |
627 | ||
628 | return chi2TR; | |
629 | ||
630 | } |