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