]>
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 | ||
4fad09c9 | 27 | #include "AliMathBase.h" |
28 | ||
75bd7f81 | 29 | #include "AliTRDseed.h" |
30 | #include "AliTRDcluster.h" | |
31 | #include "AliTRDtracker.h" | |
2985ffcb | 32 | #include "AliTRDtrackerV1.h" |
75bd7f81 | 33 | |
75bd7f81 | 34 | ClassImp(AliTRDseed) |
fbb2ea06 | 35 | |
75bd7f81 | 36 | //_____________________________________________________________________________ |
ad670fba | 37 | AliTRDseed::AliTRDseed() |
38 | :TObject() | |
39 | ,fTilt(0) | |
40 | ,fPadLength(0) | |
41 | ,fX0(0) | |
42 | ,fSigmaY(0) | |
43 | ,fSigmaY2(0) | |
44 | ,fMeanz(0) | |
45 | ,fZProb(0) | |
46 | ,fN(0) | |
47 | ,fN2(0) | |
48 | ,fNUsed(0) | |
49 | ,fFreq(0) | |
50 | ,fNChange(0) | |
51 | ,fMPads(0) | |
52 | ,fC(0) | |
53 | ,fCC(0) | |
54 | ,fChi2(1.0e10) | |
55 | ,fChi2Z(1.0e10) | |
75bd7f81 | 56 | { |
57 | // | |
58 | // Default constructor | |
59 | // | |
60 | ||
e4f2f73d | 61 | for (Int_t i = 0; i < knTimebins; i++) { |
ad670fba | 62 | fX[i] = 0; // x position |
63 | fY[i] = 0; // y position | |
64 | fZ[i] = 0; // z position | |
65 | fIndexes[i] = 0; // Indexes | |
4d6aee34 | 66 | fClusters[i] = NULL; // Clusters |
ad670fba | 67 | fUsable[i] = 0; // Indication - usable cluster |
68 | } | |
69 | ||
70 | for (Int_t i = 0; i < 2; i++) { | |
71 | fYref[i] = 0; // Reference y | |
72 | fZref[i] = 0; // Reference z | |
73 | fYfit[i] = 0; // Y fit position +derivation | |
74 | fYfitR[i] = 0; // Y fit position +derivation | |
75 | fZfit[i] = 0; // Z fit position | |
76 | fZfitR[i] = 0; // Z fit position | |
77 | fLabels[i] = 0; // Labels | |
78 | } | |
79 | ||
80 | } | |
81 | ||
82 | //_____________________________________________________________________________ | |
83 | AliTRDseed::AliTRDseed(const AliTRDseed &s) | |
84 | :TObject(s) | |
4fad09c9 | 85 | ,fTilt(s.fTilt) |
86 | ,fPadLength(s.fPadLength) | |
87 | ,fX0(s.fX0) | |
88 | ,fSigmaY(s.fSigmaY) | |
89 | ,fSigmaY2(s.fSigmaY2) | |
90 | ,fMeanz(s.fMeanz) | |
91 | ,fZProb(s.fZProb) | |
92 | ,fN(s.fN) | |
93 | ,fN2(s.fN2) | |
94 | ,fNUsed(s.fNUsed) | |
95 | ,fFreq(s.fFreq) | |
96 | ,fNChange(s.fNChange) | |
97 | ,fMPads(s.fMPads) | |
98 | ,fC(s.fC) | |
99 | ,fCC(s.fCC) | |
100 | ,fChi2(s.fChi2) | |
101 | ,fChi2Z(s.fChi2Z) | |
ad670fba | 102 | { |
103 | // | |
104 | // Copy constructor | |
105 | // | |
106 | ||
e4f2f73d | 107 | for (Int_t i = 0; i < knTimebins; i++) { |
4fad09c9 | 108 | fX[i] = s.fX[i]; // x position |
109 | fY[i] = s.fY[i]; // y position | |
110 | fZ[i] = s.fZ[i]; // z position | |
111 | fIndexes[i] = s.fIndexes[i]; // Indexes | |
112 | fClusters[i] = s.fClusters[i]; // Clusters | |
113 | fUsable[i] = s.fUsable[i]; // Indication - usable cluster | |
75bd7f81 | 114 | } |
ad670fba | 115 | |
75bd7f81 | 116 | for (Int_t i = 0; i < 2; i++) { |
4fad09c9 | 117 | fYref[i] = s.fYref[i]; // Reference y |
118 | fZref[i] = s.fZref[i]; // Reference z | |
119 | fYfit[i] = s.fYfit[i]; // Y fit position +derivation | |
120 | fYfitR[i] = s.fYfitR[i]; // Y fit position +derivation | |
121 | fZfit[i] = s.fZfit[i]; // Z fit position | |
122 | fZfitR[i] = s.fZfitR[i]; // Z fit position | |
123 | fLabels[i] = s.fLabels[i]; // Labels | |
75bd7f81 | 124 | } |
125 | ||
126 | } | |
127 | ||
7e5909fe | 128 | //_____________________________________________________________________________ |
129 | AliTRDseed &AliTRDseed::operator=(const AliTRDseed &s) | |
130 | { | |
131 | // | |
132 | // Assignment operator | |
133 | // | |
134 | ||
135 | if (this != &s) { | |
136 | ((AliTRDseed &) s).Copy(*this); | |
137 | } | |
138 | ||
139 | return *this; | |
140 | ||
141 | } | |
142 | ||
e4f2f73d | 143 | //_____________________________________________________________________________ |
0906e73e | 144 | void AliTRDseed::Copy(TObject &o) const |
d76231c8 | 145 | { |
0906e73e | 146 | //printf("AliTRDseed::Copy()\n"); |
d76231c8 | 147 | |
0906e73e | 148 | AliTRDseed &seed = (AliTRDseed &)o; |
149 | ||
0906e73e | 150 | seed.fTilt = fTilt; |
e4f2f73d | 151 | seed.fPadLength = fPadLength; |
152 | seed.fX0 = fX0; | |
153 | seed.fSigmaY = fSigmaY; | |
154 | seed.fSigmaY2 = fSigmaY2; | |
155 | seed.fMeanz = fMeanz; | |
156 | seed.fZProb = fZProb; | |
157 | seed.fN = fN; | |
158 | seed.fN2 = fN2; | |
159 | seed.fNUsed = fNUsed; | |
160 | seed.fFreq = fFreq; | |
161 | seed.fNChange = fNChange; | |
162 | seed.fMPads = fMPads; | |
163 | seed.fC = fC; | |
164 | seed.fCC = fCC; | |
165 | seed.fChi2 = fChi2; | |
166 | seed.fChi2Z = fChi2Z; | |
0906e73e | 167 | for (Int_t i = 0; i < knTimebins; i++) { |
e4f2f73d | 168 | seed.fX[i] = fX[i]; |
169 | seed.fY[i] = fY[i]; | |
170 | seed.fZ[i] = fZ[i]; | |
171 | seed.fIndexes[i] = fIndexes[i]; | |
172 | seed.fClusters[i] = fClusters[i]; | |
173 | seed.fUsable[i] = fUsable[i]; | |
174 | } | |
175 | ||
176 | for (Int_t i = 0; i < 2; i++) { | |
177 | seed.fYref[i] = fYref[i]; | |
178 | seed.fZref[i] = fZref[i]; | |
179 | seed.fYfit[i] = fYfit[i]; | |
180 | seed.fYfitR[i] = fYfitR[i]; | |
181 | seed.fZfit[i] = fZfit[i]; | |
182 | seed.fZfitR[i] = fZfitR[i]; | |
183 | seed.fLabels[i] = fLabels[i]; | |
184 | } | |
185 | ||
186 | TObject::Copy(seed); | |
187 | ||
188 | } | |
189 | ||
75bd7f81 | 190 | //_____________________________________________________________________________ |
191 | void AliTRDseed::Reset() | |
192 | { | |
193 | // | |
194 | // Reset seed | |
195 | // | |
196 | ||
e4f2f73d | 197 | for (Int_t i = 0; i < knTimebins; i++) { |
ad670fba | 198 | fX[i] = 0; // X position |
199 | fY[i] = 0; // Y position | |
200 | fZ[i] = 0; // Z position | |
201 | fIndexes[i] = 0; // Indexes | |
4d6aee34 | 202 | fClusters[i] = NULL; // Clusters |
75bd7f81 | 203 | fUsable[i] = kFALSE; |
204 | } | |
ad670fba | 205 | |
75bd7f81 | 206 | for (Int_t i = 0; i < 2; i++) { |
ad670fba | 207 | fYref[i] = 0; // Reference y |
208 | fZref[i] = 0; // Reference z | |
209 | fYfit[i] = 0; // Y fit position +derivation | |
210 | fYfitR[i] = 0; // Y fit position +derivation | |
211 | fZfit[i] = 0; // Z fit position | |
212 | fZfitR[i] = 0; // Z fit position | |
213 | fLabels[i] = -1; // Labels | |
75bd7f81 | 214 | } |
ad670fba | 215 | fSigmaY = 0; // "Robust" sigma in y |
216 | fSigmaY2 = 0; // "Robust" sigma in y | |
217 | fMeanz = 0; // Mean vaue of z | |
218 | fZProb = 0; // Max probbable z | |
75bd7f81 | 219 | fMPads = 0; |
ad670fba | 220 | fN = 0; // Number of associated clusters |
221 | fN2 = 0; // Number of not crossed | |
222 | fNUsed = 0; // Number of used clusters | |
223 | fNChange = 0; // Change z counter | |
75bd7f81 | 224 | |
225 | } | |
226 | ||
227 | //_____________________________________________________________________________ | |
228 | void AliTRDseed::CookLabels() | |
229 | { | |
230 | // | |
231 | // Cook 2 labels for seed | |
232 | // | |
233 | ||
234 | Int_t labels[200]; | |
235 | Int_t out[200]; | |
ad670fba | 236 | Int_t nlab = 0; |
237 | ||
2985ffcb | 238 | for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) { |
75bd7f81 | 239 | if (!fClusters[i]) continue; |
ad670fba | 240 | for (Int_t ilab = 0; ilab < 3; ilab++) { |
241 | if (fClusters[i]->GetLabel(ilab) >= 0) { | |
75bd7f81 | 242 | labels[nlab] = fClusters[i]->GetLabel(ilab); |
243 | nlab++; | |
244 | } | |
245 | } | |
246 | } | |
ad670fba | 247 | |
75bd7f81 | 248 | Int_t nlab2 = AliTRDtracker::Freq(nlab,labels,out,kTRUE); |
249 | fLabels[0] = out[0]; | |
ad670fba | 250 | if ((nlab2 > 1) && |
251 | (out[3] > 1)) { | |
252 | fLabels[1] = out[2]; | |
253 | } | |
75bd7f81 | 254 | |
255 | } | |
256 | ||
257 | //_____________________________________________________________________________ | |
258 | void AliTRDseed::UseClusters() | |
259 | { | |
260 | // | |
261 | // Use clusters | |
262 | // | |
263 | ||
2985ffcb | 264 | for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) { |
75bd7f81 | 265 | if (!fClusters[i]) continue; |
266 | if (!(fClusters[i]->IsUsed())) fClusters[i]->Use(); | |
267 | } | |
268 | ||
269 | } | |
270 | ||
271 | //_____________________________________________________________________________ | |
272 | void AliTRDseed::Update() | |
273 | { | |
274 | // | |
275 | // Update the seed. | |
276 | // | |
277 | ||
0906e73e | 278 | |
bcb6fb78 | 279 | |
280 | // linear fit on the y direction with respect to the reference direction. | |
281 | // The residuals for each x (x = xc - x0) are deduced from: | |
282 | // dy = y - yt (1) | |
283 | // the tilting correction is written : | |
284 | // y = yc + h*(zc-zt) (2) | |
285 | // yt = y0+dy/dx*x (3) | |
286 | // zt = z0+dz/dx*x (4) | |
287 | // from (1),(2),(3) and (4) | |
288 | // dy = yc - y0 - (dy/dx + h*dz/dx)*x + h*(zc-z0) | |
289 | // the last term introduces the correction on y direction due to tilting pads. There are 2 ways to account for this: | |
290 | // 1. use tilting correction for calculating the y | |
291 | // 2. neglect tilting correction here and account for it in the error parametrization of the tracklet. | |
0906e73e | 292 | |
75bd7f81 | 293 | const Float_t kRatio = 0.8; |
e4f2f73d | 294 | const Int_t kClmin = 5; |
75bd7f81 | 295 | const Float_t kmaxtan = 2; |
296 | ||
e4f2f73d | 297 | if (TMath::Abs(fYref[1]) > kmaxtan){ |
298 | //printf("Exit: Abs(fYref[1]) = %3.3f, kmaxtan = %3.3f\n", TMath::Abs(fYref[1]), kmaxtan); | |
299 | return; // Track inclined too much | |
300 | } | |
ad670fba | 301 | |
302 | Float_t sigmaexp = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction | |
303 | Float_t ycrosscor = fPadLength * fTilt * 0.5; // Y correction for crossing | |
304 | fNChange = 0; | |
305 | ||
306 | Double_t sumw; | |
307 | Double_t sumwx; | |
308 | Double_t sumwx2; | |
309 | Double_t sumwy; | |
310 | Double_t sumwxy; | |
311 | Double_t sumwz; | |
312 | Double_t sumwxz; | |
313 | ||
e4f2f73d | 314 | // Buffering: Leave it constant fot Performance issues |
315 | Int_t zints[knTimebins]; // Histograming of the z coordinate | |
ad670fba | 316 | // Get 1 and second max probable coodinates in z |
e4f2f73d | 317 | Int_t zouts[2*knTimebins]; |
318 | Float_t allowedz[knTimebins]; // Allowed z for given time bin | |
87f70a90 | 319 | memset(allowedz, 0, knTimebins*sizeof(Float_t)); |
e4f2f73d | 320 | Float_t yres[knTimebins]; // Residuals from reference |
0849f128 | 321 | //Float_t anglecor = fTilt * fZref[1]; // Correction to the angle |
75bd7f81 | 322 | |
323 | ||
ad670fba | 324 | fN = 0; |
325 | fN2 = 0; | |
2985ffcb | 326 | for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) { |
ad670fba | 327 | yres[i] = 10000.0; |
75bd7f81 | 328 | if (!fClusters[i]) continue; |
bcb6fb78 | 329 | if(!fClusters[i]->IsInChamber()) continue; |
0849f128 | 330 | // Residual y |
331 | //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] + fTilt*(fZ[i] - fZref[0]); | |
332 | yres[i] = fY[i] - fTilt*(fZ[i] - (fZref[0] - fX[i]*fZref[1])); | |
75bd7f81 | 333 | zints[fN] = Int_t(fZ[i]); |
334 | fN++; | |
335 | } | |
336 | ||
e4f2f73d | 337 | if (fN < kClmin){ |
338 | //printf("Exit fN < kClmin: fN = %d\n", fN); | |
339 | return; | |
340 | } | |
bcb6fb78 | 341 | Int_t nz = AliTRDtracker::Freq(fN, zints, zouts, kFALSE); |
75bd7f81 | 342 | fZProb = zouts[0]; |
ad670fba | 343 | if (nz <= 1) zouts[3] = 0; |
e4f2f73d | 344 | if (zouts[1] + zouts[3] < kClmin) { |
345 | //printf("Exit zouts[1] = %d, zouts[3] = %d\n",zouts[1],zouts[3]); | |
346 | return; | |
347 | } | |
75bd7f81 | 348 | |
ad670fba | 349 | // Z distance bigger than pad - length |
bcb6fb78 | 350 | if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) zouts[3] = 0; |
75bd7f81 | 351 | |
352 | Int_t breaktime = -1; | |
353 | Bool_t mbefore = kFALSE; | |
e4f2f73d | 354 | Int_t cumul[knTimebins][2]; |
ad670fba | 355 | Int_t counts[2] = { 0, 0 }; |
75bd7f81 | 356 | |
ad670fba | 357 | if (zouts[3] >= 3) { |
35c24814 | 358 | |
359 | // | |
360 | // Find the break time allowing one chage on pad-rows | |
361 | // with maximal number of accepted clusters | |
362 | // | |
ad670fba | 363 | fNChange = 1; |
2985ffcb | 364 | for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) { |
75bd7f81 | 365 | cumul[i][0] = counts[0]; |
366 | cumul[i][1] = counts[1]; | |
ad670fba | 367 | if (TMath::Abs(fZ[i]-zouts[0]) < 2) counts[0]++; |
368 | if (TMath::Abs(fZ[i]-zouts[2]) < 2) counts[1]++; | |
75bd7f81 | 369 | } |
ad670fba | 370 | Int_t maxcount = 0; |
2985ffcb | 371 | for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) { |
35c24814 | 372 | Int_t after = cumul[AliTRDtrackerV1::GetNTimeBins()][0] - cumul[i][0]; |
75bd7f81 | 373 | Int_t before = cumul[i][1]; |
ad670fba | 374 | if (after + before > maxcount) { |
35c24814 | 375 | maxcount = after + before; |
376 | breaktime = i; | |
377 | mbefore = kFALSE; | |
75bd7f81 | 378 | } |
35c24814 | 379 | after = cumul[AliTRDtrackerV1::GetNTimeBins()-1][1] - cumul[i][1]; |
75bd7f81 | 380 | before = cumul[i][0]; |
ad670fba | 381 | if (after + before > maxcount) { |
35c24814 | 382 | maxcount = after + before; |
383 | breaktime = i; | |
384 | mbefore = kTRUE; | |
75bd7f81 | 385 | } |
386 | } | |
35c24814 | 387 | |
ad670fba | 388 | breaktime -= 1; |
35c24814 | 389 | |
75bd7f81 | 390 | } |
ad670fba | 391 | |
2985ffcb | 392 | for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) { |
ad670fba | 393 | if (i > breaktime) allowedz[i] = mbefore ? zouts[2] : zouts[0]; |
394 | if (i <= breaktime) allowedz[i] = (!mbefore) ? zouts[2] : zouts[0]; | |
75bd7f81 | 395 | } |
ad670fba | 396 | |
2985ffcb | 397 | if (((allowedz[0] > allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] < 0)) || |
398 | ((allowedz[0] < allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] > 0))) { | |
75bd7f81 | 399 | // |
ad670fba | 400 | // Tracklet z-direction not in correspondance with track z direction |
75bd7f81 | 401 | // |
ad670fba | 402 | fNChange = 0; |
35c24814 | 403 | for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) { |
404 | allowedz[i] = zouts[0]; // Only longest taken | |
405 | } | |
75bd7f81 | 406 | } |
407 | ||
35c24814 | 408 | if (fNChange > 0) { |
75bd7f81 | 409 | // |
ad670fba | 410 | // Cross pad -row tracklet - take the step change into account |
75bd7f81 | 411 | // |
2985ffcb | 412 | for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) { |
75bd7f81 | 413 | if (!fClusters[i]) continue; |
bcb6fb78 | 414 | if(!fClusters[i]->IsInChamber()) continue; |
ad670fba | 415 | if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue; |
0849f128 | 416 | // Residual y |
417 | //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] /*+ fTilt*(fZ[i] - fZref[0])*/; | |
418 | yres[i] = fY[i] - fTilt*(fZ[i] - (fZref[0] - fX[i]*fZref[1])); | |
419 | /* if (TMath::Abs(fZ[i] - fZProb) > 2) { | |
420 | if (fZ[i] > fZProb) yres[i] += fTilt * fPadLength; | |
421 | if (fZ[i] < fZProb) yres[i] -= fTilt * fPadLength; | |
422 | }*/ | |
75bd7f81 | 423 | } |
35c24814 | 424 | } |
75bd7f81 | 425 | |
e4f2f73d | 426 | Double_t yres2[knTimebins]; |
ad670fba | 427 | Double_t mean; |
428 | Double_t sigma; | |
2985ffcb | 429 | for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) { |
75bd7f81 | 430 | if (!fClusters[i]) continue; |
bcb6fb78 | 431 | if(!fClusters[i]->IsInChamber()) continue; |
ad670fba | 432 | if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue; |
433 | yres2[fN2] = yres[i]; | |
75bd7f81 | 434 | fN2++; |
435 | } | |
ad670fba | 436 | if (fN2 < kClmin) { |
e4f2f73d | 437 | //printf("Exit fN2 < kClmin: fN2 = %d\n", fN2); |
75bd7f81 | 438 | fN2 = 0; |
439 | return; | |
440 | } | |
e4f2f73d | 441 | AliMathBase::EvaluateUni(fN2,yres2,mean,sigma, Int_t(fN2*kRatio-2.)); |
ad670fba | 442 | if (sigma < sigmaexp * 0.8) { |
443 | sigma = sigmaexp; | |
444 | } | |
75bd7f81 | 445 | fSigmaY = sigma; |
ad670fba | 446 | |
447 | // Reset sums | |
448 | sumw = 0; | |
449 | sumwx = 0; | |
450 | sumwx2 = 0; | |
451 | sumwy = 0; | |
452 | sumwxy = 0; | |
453 | sumwz = 0; | |
454 | sumwxz = 0; | |
455 | ||
456 | fN2 = 0; | |
457 | fMeanz = 0; | |
458 | fMPads = 0; | |
459 | ||
2985ffcb | 460 | for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) { |
ad670fba | 461 | |
462 | fUsable[i] = kFALSE; | |
75bd7f81 | 463 | if (!fClusters[i]) continue; |
bcb6fb78 | 464 | if (!fClusters[i]->IsInChamber()) continue; |
4d6aee34 | 465 | if (TMath::Abs(fZ[i] - allowedz[i]) > 2){fClusters[i] = NULL; continue;} |
466 | if (TMath::Abs(yres[i] - mean) > 4.0 * sigma){fClusters[i] = NULL; continue;} | |
75bd7f81 | 467 | fUsable[i] = kTRUE; |
468 | fN2++; | |
ad670fba | 469 | fMPads += fClusters[i]->GetNPads(); |
470 | Float_t weight = 1.0; | |
471 | if (fClusters[i]->GetNPads() > 4) weight = 0.5; | |
472 | if (fClusters[i]->GetNPads() > 5) weight = 0.2; | |
473 | ||
0906e73e | 474 | |
75bd7f81 | 475 | Double_t x = fX[i]; |
0906e73e | 476 | //printf("x = %7.3f dy = %7.3f fit %7.3f\n", x, yres[i], fY[i]-yres[i]); |
477 | ||
ad670fba | 478 | sumw += weight; |
479 | sumwx += x * weight; | |
480 | sumwx2 += x*x * weight; | |
481 | sumwy += weight * yres[i]; | |
482 | sumwxy += weight * (yres[i]) * x; | |
483 | sumwz += weight * fZ[i]; | |
484 | sumwxz += weight * fZ[i] * x; | |
485 | ||
75bd7f81 | 486 | } |
ad670fba | 487 | |
488 | if (fN2 < kClmin){ | |
e4f2f73d | 489 | //printf("Exit fN2 < kClmin(2): fN2 = %d\n",fN2); |
75bd7f81 | 490 | fN2 = 0; |
491 | return; | |
492 | } | |
ad670fba | 493 | fMeanz = sumwz / sumw; |
494 | Float_t correction = 0; | |
495 | if (fNChange > 0) { | |
496 | // Tracklet on boundary | |
497 | if (fMeanz < fZProb) correction = ycrosscor; | |
498 | if (fMeanz > fZProb) correction = -ycrosscor; | |
75bd7f81 | 499 | } |
ad670fba | 500 | |
501 | Double_t det = sumw * sumwx2 - sumwx * sumwx; | |
502 | fYfitR[0] = (sumwx2 * sumwy - sumwx * sumwxy) / det; | |
503 | fYfitR[1] = (sumw * sumwxy - sumwx * sumwy) / det; | |
75bd7f81 | 504 | |
ad670fba | 505 | fSigmaY2 = 0; |
2985ffcb | 506 | for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) { |
75bd7f81 | 507 | if (!fUsable[i]) continue; |
ad670fba | 508 | Float_t delta = yres[i] - fYfitR[0] - fYfitR[1] * fX[i]; |
509 | fSigmaY2 += delta*delta; | |
75bd7f81 | 510 | } |
ad670fba | 511 | fSigmaY2 = TMath::Sqrt(fSigmaY2 / Float_t(fN2-2)); |
0906e73e | 512 | // TEMPORARY UNTIL covariance properly calculated |
513 | fSigmaY2 = TMath::Max(fSigmaY2, Float_t(.1)); | |
ad670fba | 514 | |
515 | fZfitR[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det; | |
516 | fZfitR[1] = (sumw * sumwxz - sumwx * sumwz) / det; | |
517 | fZfit[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det; | |
518 | fZfit[1] = (sumw * sumwxz - sumwx * sumwz) / det; | |
0849f128 | 519 | // fYfitR[0] += fYref[0] + correction; |
520 | // fYfitR[1] += fYref[1]; | |
ad670fba | 521 | fYfit[0] = fYfitR[0]; |
0849f128 | 522 | fYfit[1] = -fYfitR[1]; |
0906e73e | 523 | |
524 | //printf("y0 = %7.3f tgy = %7.3f z0 = %7.3f tgz = %7.3f \n", fYfitR[0], fYfitR[1], fZfitR[0], fZfitR[1]); | |
525 | ||
75bd7f81 | 526 | UpdateUsed(); |
527 | ||
528 | } | |
529 | ||
530 | //_____________________________________________________________________________ | |
531 | void AliTRDseed::UpdateUsed() | |
532 | { | |
533 | // | |
534 | // Update used seed | |
535 | // | |
536 | ||
537 | fNUsed = 0; | |
2985ffcb | 538 | for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) { |
0906e73e | 539 | if (!fClusters[i]) continue; |
540 | if(!fUsable[i]) continue; | |
541 | if ((fClusters[i]->IsUsed())) fNUsed++; | |
75bd7f81 | 542 | } |
543 | ||
544 | } | |
545 | ||
75bd7f81 | 546 | //_____________________________________________________________________________ |
547 | Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror) | |
548 | { | |
549 | // | |
550 | // Fit the Rieman tilt | |
551 | // | |
552 | ||
ad670fba | 553 | // Fitting with tilting pads - kz not fixed |
554 | TLinearFitter fitterT2(4,"hyp4"); | |
75bd7f81 | 555 | fitterT2.StoreData(kTRUE); |
e4f2f73d | 556 | |
ad670fba | 557 | Float_t xref2 = (cseed[2].fX0 + cseed[3].fX0) * 0.5; // Reference x0 for z |
558 | ||
559 | Int_t npointsT = 0; | |
75bd7f81 | 560 | fitterT2.ClearPoints(); |
ad670fba | 561 | |
562 | for (Int_t iLayer = 0; iLayer < 6; iLayer++) { | |
563 | ||
75bd7f81 | 564 | if (!cseed[iLayer].IsOK()) continue; |
565 | Double_t tilt = cseed[iLayer].fTilt; | |
566 | ||
2985ffcb | 567 | for (Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins()+1; itime++) { |
ad670fba | 568 | |
75bd7f81 | 569 | if (!cseed[iLayer].fUsable[itime]) continue; |
570 | // x relative to the midle chamber | |
ad670fba | 571 | Double_t x = cseed[iLayer].fX[itime] + cseed[iLayer].fX0 - xref2; |
572 | Double_t y = cseed[iLayer].fY[itime]; | |
573 | Double_t z = cseed[iLayer].fZ[itime]; | |
75bd7f81 | 574 | |
575 | // | |
ad670fba | 576 | // Tilted rieman |
75bd7f81 | 577 | // |
578 | Double_t uvt[6]; | |
ad670fba | 579 | Double_t x2 = cseed[iLayer].fX[itime] + cseed[iLayer].fX0; // Global x |
580 | Double_t t = 1.0 / (x2*x2 + y*y); | |
581 | uvt[1] = t; | |
582 | uvt[0] = 2.0 * x2 * uvt[1]; | |
583 | uvt[2] = 2.0 * tilt * uvt[1]; | |
584 | uvt[3] = 2.0 * tilt *uvt[1] * x; | |
585 | uvt[4] = 2.0 * (y + tilt * z) * uvt[1]; | |
586 | ||
587 | Double_t error = 2.0 * uvt[1]; | |
588 | if (terror) { | |
589 | error *= cseed[iLayer].fSigmaY; | |
590 | } | |
591 | else { | |
592 | error *= 0.2; //Default error | |
593 | } | |
75bd7f81 | 594 | fitterT2.AddPoint(uvt,uvt[4],error); |
595 | npointsT++; | |
ad670fba | 596 | |
75bd7f81 | 597 | } |
ad670fba | 598 | |
75bd7f81 | 599 | } |
ad670fba | 600 | |
75bd7f81 | 601 | fitterT2.Eval(); |
602 | Double_t rpolz0 = fitterT2.GetParameter(3); | |
603 | Double_t rpolz1 = fitterT2.GetParameter(4); | |
ad670fba | 604 | |
75bd7f81 | 605 | // |
ad670fba | 606 | // Linear fitter - not possible to make boundaries |
75bd7f81 | 607 | // non accept non possible z and dzdx combination |
608 | // | |
ad670fba | 609 | Bool_t acceptablez = kTRUE; |
610 | for (Int_t iLayer = 0; iLayer < 6; iLayer++) { | |
e4f2f73d | 611 | if (!cseed[iLayer].IsOK()) continue; |
612 | Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].fX0 - xref2); | |
613 | if (TMath::Abs(cseed[iLayer].fZProb - zT2) > cseed[iLayer].fPadLength * 0.5 + 1.0) acceptablez = kFALSE; | |
75bd7f81 | 614 | } |
ad670fba | 615 | if (!acceptablez) { |
616 | Double_t zmf = cseed[2].fZref[0] + cseed[2].fZref[1] * (xref2 - cseed[2].fX0); | |
617 | Double_t dzmf = (cseed[2].fZref[1] + cseed[3].fZref[1]) * 0.5; | |
75bd7f81 | 618 | fitterT2.FixParameter(3,zmf); |
619 | fitterT2.FixParameter(4,dzmf); | |
620 | fitterT2.Eval(); | |
621 | fitterT2.ReleaseParameter(3); | |
622 | fitterT2.ReleaseParameter(4); | |
623 | rpolz0 = fitterT2.GetParameter(3); | |
624 | rpolz1 = fitterT2.GetParameter(4); | |
625 | } | |
626 | ||
ad670fba | 627 | Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT); |
75bd7f81 | 628 | Double_t params[3]; |
ad670fba | 629 | params[0] = fitterT2.GetParameter(0); |
630 | params[1] = fitterT2.GetParameter(1); | |
631 | params[2] = fitterT2.GetParameter(2); | |
632 | Double_t curvature = 1.0 + params[1] * params[1] - params[2] * params[0]; | |
633 | ||
e4f2f73d | 634 | |
ad670fba | 635 | for (Int_t iLayer = 0; iLayer < 6; iLayer++) { |
ad670fba | 636 | Double_t x = cseed[iLayer].fX0; |
637 | Double_t y = 0; | |
638 | Double_t dy = 0; | |
639 | Double_t z = 0; | |
640 | Double_t dz = 0; | |
641 | ||
75bd7f81 | 642 | // y |
ad670fba | 643 | Double_t res2 = (x * params[0] + params[1]); |
644 | res2 *= res2; | |
645 | res2 = 1.0 - params[2]*params[0] + params[1]*params[1] - res2; | |
646 | if (res2 >= 0) { | |
75bd7f81 | 647 | res2 = TMath::Sqrt(res2); |
ad670fba | 648 | y = (1.0 - res2) / params[0]; |
75bd7f81 | 649 | } |
ad670fba | 650 | |
75bd7f81 | 651 | //dy |
ad670fba | 652 | Double_t x0 = -params[1] / params[0]; |
653 | if (-params[2]*params[0] + params[1]*params[1] + 1 > 0) { | |
654 | Double_t rm1 = params[0] / TMath::Sqrt(-params[2]*params[0] + params[1]*params[1] + 1); | |
655 | if (1.0/(rm1*rm1) - (x-x0) * (x-x0) > 0.0) { | |
60e55aee | 656 | Double_t res = (x - x0) / TMath::Sqrt((1./rm1-(x-x0))*(1./rm1+(x-x0))); |
e4f2f73d | 657 | if (params[0] < 0) res *= -1.0; |
658 | dy = res; | |
75bd7f81 | 659 | } |
660 | } | |
ad670fba | 661 | z = rpolz0 + rpolz1 * (x - xref2); |
75bd7f81 | 662 | dz = rpolz1; |
663 | cseed[iLayer].fYref[0] = y; | |
664 | cseed[iLayer].fYref[1] = dy; | |
665 | cseed[iLayer].fZref[0] = z; | |
666 | cseed[iLayer].fZref[1] = dz; | |
ad670fba | 667 | cseed[iLayer].fC = curvature; |
75bd7f81 | 668 | |
669 | } | |
670 | ||
671 | return chi2TR; | |
672 | ||
673 | } |