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