]> 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 protection for friend files.
[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
7e5909fe 128//_____________________________________________________________________________
129AliTRDseed &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 144void 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//_____________________________________________________________________________
191void 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//_____________________________________________________________________________
228void 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//_____________________________________________________________________________
258void 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//_____________________________________________________________________________
272void 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
319 Float_t yres[knTimebins]; // Residuals from reference
0849f128 320 //Float_t anglecor = fTilt * fZref[1]; // Correction to the angle
75bd7f81 321
322
ad670fba 323 fN = 0;
324 fN2 = 0;
2985ffcb 325 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
ad670fba 326 yres[i] = 10000.0;
75bd7f81 327 if (!fClusters[i]) continue;
bcb6fb78 328 if(!fClusters[i]->IsInChamber()) continue;
0849f128 329 // Residual y
330 //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] + fTilt*(fZ[i] - fZref[0]);
331 yres[i] = fY[i] - fTilt*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
75bd7f81 332 zints[fN] = Int_t(fZ[i]);
333 fN++;
334 }
335
e4f2f73d 336 if (fN < kClmin){
337 //printf("Exit fN < kClmin: fN = %d\n", fN);
338 return;
339 }
bcb6fb78 340 Int_t nz = AliTRDtracker::Freq(fN, zints, zouts, kFALSE);
75bd7f81 341 fZProb = zouts[0];
ad670fba 342 if (nz <= 1) zouts[3] = 0;
e4f2f73d 343 if (zouts[1] + zouts[3] < kClmin) {
344 //printf("Exit zouts[1] = %d, zouts[3] = %d\n",zouts[1],zouts[3]);
345 return;
346 }
75bd7f81 347
ad670fba 348 // Z distance bigger than pad - length
bcb6fb78 349 if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) zouts[3] = 0;
75bd7f81 350
351 Int_t breaktime = -1;
352 Bool_t mbefore = kFALSE;
e4f2f73d 353 Int_t cumul[knTimebins][2];
ad670fba 354 Int_t counts[2] = { 0, 0 };
75bd7f81 355
ad670fba 356 if (zouts[3] >= 3) {
35c24814 357
358 //
359 // Find the break time allowing one chage on pad-rows
360 // with maximal number of accepted clusters
361 //
ad670fba 362 fNChange = 1;
2985ffcb 363 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
75bd7f81 364 cumul[i][0] = counts[0];
365 cumul[i][1] = counts[1];
ad670fba 366 if (TMath::Abs(fZ[i]-zouts[0]) < 2) counts[0]++;
367 if (TMath::Abs(fZ[i]-zouts[2]) < 2) counts[1]++;
75bd7f81 368 }
ad670fba 369 Int_t maxcount = 0;
2985ffcb 370 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
35c24814 371 Int_t after = cumul[AliTRDtrackerV1::GetNTimeBins()][0] - cumul[i][0];
75bd7f81 372 Int_t before = cumul[i][1];
ad670fba 373 if (after + before > maxcount) {
35c24814 374 maxcount = after + before;
375 breaktime = i;
376 mbefore = kFALSE;
75bd7f81 377 }
35c24814 378 after = cumul[AliTRDtrackerV1::GetNTimeBins()-1][1] - cumul[i][1];
75bd7f81 379 before = cumul[i][0];
ad670fba 380 if (after + before > maxcount) {
35c24814 381 maxcount = after + before;
382 breaktime = i;
383 mbefore = kTRUE;
75bd7f81 384 }
385 }
35c24814 386
ad670fba 387 breaktime -= 1;
35c24814 388
75bd7f81 389 }
ad670fba 390
2985ffcb 391 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
ad670fba 392 if (i > breaktime) allowedz[i] = mbefore ? zouts[2] : zouts[0];
393 if (i <= breaktime) allowedz[i] = (!mbefore) ? zouts[2] : zouts[0];
75bd7f81 394 }
ad670fba 395
2985ffcb 396 if (((allowedz[0] > allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] < 0)) ||
397 ((allowedz[0] < allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] > 0))) {
75bd7f81 398 //
ad670fba 399 // Tracklet z-direction not in correspondance with track z direction
75bd7f81 400 //
ad670fba 401 fNChange = 0;
35c24814 402 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
403 allowedz[i] = zouts[0]; // Only longest taken
404 }
75bd7f81 405 }
406
35c24814 407 if (fNChange > 0) {
75bd7f81 408 //
ad670fba 409 // Cross pad -row tracklet - take the step change into account
75bd7f81 410 //
2985ffcb 411 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
75bd7f81 412 if (!fClusters[i]) continue;
bcb6fb78 413 if(!fClusters[i]->IsInChamber()) continue;
ad670fba 414 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
0849f128 415 // Residual y
416 //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] /*+ fTilt*(fZ[i] - fZref[0])*/;
417 yres[i] = fY[i] - fTilt*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
418/* if (TMath::Abs(fZ[i] - fZProb) > 2) {
419 if (fZ[i] > fZProb) yres[i] += fTilt * fPadLength;
420 if (fZ[i] < fZProb) yres[i] -= fTilt * fPadLength;
421 }*/
75bd7f81 422 }
35c24814 423 }
75bd7f81 424
e4f2f73d 425 Double_t yres2[knTimebins];
ad670fba 426 Double_t mean;
427 Double_t sigma;
2985ffcb 428 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
75bd7f81 429 if (!fClusters[i]) continue;
bcb6fb78 430 if(!fClusters[i]->IsInChamber()) continue;
ad670fba 431 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
432 yres2[fN2] = yres[i];
75bd7f81 433 fN2++;
434 }
ad670fba 435 if (fN2 < kClmin) {
e4f2f73d 436 //printf("Exit fN2 < kClmin: fN2 = %d\n", fN2);
75bd7f81 437 fN2 = 0;
438 return;
439 }
e4f2f73d 440 AliMathBase::EvaluateUni(fN2,yres2,mean,sigma, Int_t(fN2*kRatio-2.));
ad670fba 441 if (sigma < sigmaexp * 0.8) {
442 sigma = sigmaexp;
443 }
75bd7f81 444 fSigmaY = sigma;
ad670fba 445
446 // Reset sums
447 sumw = 0;
448 sumwx = 0;
449 sumwx2 = 0;
450 sumwy = 0;
451 sumwxy = 0;
452 sumwz = 0;
453 sumwxz = 0;
454
455 fN2 = 0;
456 fMeanz = 0;
457 fMPads = 0;
458
2985ffcb 459 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
ad670fba 460
461 fUsable[i] = kFALSE;
75bd7f81 462 if (!fClusters[i]) continue;
bcb6fb78 463 if (!fClusters[i]->IsInChamber()) continue;
4d6aee34 464 if (TMath::Abs(fZ[i] - allowedz[i]) > 2){fClusters[i] = NULL; continue;}
465 if (TMath::Abs(yres[i] - mean) > 4.0 * sigma){fClusters[i] = NULL; continue;}
75bd7f81 466 fUsable[i] = kTRUE;
467 fN2++;
ad670fba 468 fMPads += fClusters[i]->GetNPads();
469 Float_t weight = 1.0;
470 if (fClusters[i]->GetNPads() > 4) weight = 0.5;
471 if (fClusters[i]->GetNPads() > 5) weight = 0.2;
472
0906e73e 473
75bd7f81 474 Double_t x = fX[i];
0906e73e 475 //printf("x = %7.3f dy = %7.3f fit %7.3f\n", x, yres[i], fY[i]-yres[i]);
476
ad670fba 477 sumw += weight;
478 sumwx += x * weight;
479 sumwx2 += x*x * weight;
480 sumwy += weight * yres[i];
481 sumwxy += weight * (yres[i]) * x;
482 sumwz += weight * fZ[i];
483 sumwxz += weight * fZ[i] * x;
484
75bd7f81 485 }
ad670fba 486
487 if (fN2 < kClmin){
e4f2f73d 488 //printf("Exit fN2 < kClmin(2): fN2 = %d\n",fN2);
75bd7f81 489 fN2 = 0;
490 return;
491 }
ad670fba 492 fMeanz = sumwz / sumw;
493 Float_t correction = 0;
494 if (fNChange > 0) {
495 // Tracklet on boundary
496 if (fMeanz < fZProb) correction = ycrosscor;
497 if (fMeanz > fZProb) correction = -ycrosscor;
75bd7f81 498 }
ad670fba 499
500 Double_t det = sumw * sumwx2 - sumwx * sumwx;
501 fYfitR[0] = (sumwx2 * sumwy - sumwx * sumwxy) / det;
502 fYfitR[1] = (sumw * sumwxy - sumwx * sumwy) / det;
75bd7f81 503
ad670fba 504 fSigmaY2 = 0;
2985ffcb 505 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
75bd7f81 506 if (!fUsable[i]) continue;
ad670fba 507 Float_t delta = yres[i] - fYfitR[0] - fYfitR[1] * fX[i];
508 fSigmaY2 += delta*delta;
75bd7f81 509 }
ad670fba 510 fSigmaY2 = TMath::Sqrt(fSigmaY2 / Float_t(fN2-2));
0906e73e 511 // TEMPORARY UNTIL covariance properly calculated
512 fSigmaY2 = TMath::Max(fSigmaY2, Float_t(.1));
ad670fba 513
514 fZfitR[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det;
515 fZfitR[1] = (sumw * sumwxz - sumwx * sumwz) / det;
516 fZfit[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det;
517 fZfit[1] = (sumw * sumwxz - sumwx * sumwz) / det;
0849f128 518// fYfitR[0] += fYref[0] + correction;
519// fYfitR[1] += fYref[1];
ad670fba 520 fYfit[0] = fYfitR[0];
0849f128 521 fYfit[1] = -fYfitR[1];
0906e73e 522
523 //printf("y0 = %7.3f tgy = %7.3f z0 = %7.3f tgz = %7.3f \n", fYfitR[0], fYfitR[1], fZfitR[0], fZfitR[1]);
524
75bd7f81 525 UpdateUsed();
526
527}
528
529//_____________________________________________________________________________
530void AliTRDseed::UpdateUsed()
531{
532 //
533 // Update used seed
534 //
535
536 fNUsed = 0;
2985ffcb 537 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
0906e73e 538 if (!fClusters[i]) continue;
539 if(!fUsable[i]) continue;
540 if ((fClusters[i]->IsUsed())) fNUsed++;
75bd7f81 541 }
542
543}
544
75bd7f81 545//_____________________________________________________________________________
546Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror)
547{
548 //
549 // Fit the Rieman tilt
550 //
551
ad670fba 552 // Fitting with tilting pads - kz not fixed
553 TLinearFitter fitterT2(4,"hyp4");
75bd7f81 554 fitterT2.StoreData(kTRUE);
e4f2f73d 555
ad670fba 556 Float_t xref2 = (cseed[2].fX0 + cseed[3].fX0) * 0.5; // Reference x0 for z
557
558 Int_t npointsT = 0;
75bd7f81 559 fitterT2.ClearPoints();
ad670fba 560
561 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
562
75bd7f81 563 if (!cseed[iLayer].IsOK()) continue;
564 Double_t tilt = cseed[iLayer].fTilt;
565
2985ffcb 566 for (Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins()+1; itime++) {
ad670fba 567
75bd7f81 568 if (!cseed[iLayer].fUsable[itime]) continue;
569 // x relative to the midle chamber
ad670fba 570 Double_t x = cseed[iLayer].fX[itime] + cseed[iLayer].fX0 - xref2;
571 Double_t y = cseed[iLayer].fY[itime];
572 Double_t z = cseed[iLayer].fZ[itime];
75bd7f81 573
574 //
ad670fba 575 // Tilted rieman
75bd7f81 576 //
577 Double_t uvt[6];
ad670fba 578 Double_t x2 = cseed[iLayer].fX[itime] + cseed[iLayer].fX0; // Global x
579 Double_t t = 1.0 / (x2*x2 + y*y);
580 uvt[1] = t;
581 uvt[0] = 2.0 * x2 * uvt[1];
582 uvt[2] = 2.0 * tilt * uvt[1];
583 uvt[3] = 2.0 * tilt *uvt[1] * x;
584 uvt[4] = 2.0 * (y + tilt * z) * uvt[1];
585
586 Double_t error = 2.0 * uvt[1];
587 if (terror) {
588 error *= cseed[iLayer].fSigmaY;
589 }
590 else {
591 error *= 0.2; //Default error
592 }
75bd7f81 593 fitterT2.AddPoint(uvt,uvt[4],error);
594 npointsT++;
ad670fba 595
75bd7f81 596 }
ad670fba 597
75bd7f81 598 }
ad670fba 599
75bd7f81 600 fitterT2.Eval();
601 Double_t rpolz0 = fitterT2.GetParameter(3);
602 Double_t rpolz1 = fitterT2.GetParameter(4);
ad670fba 603
75bd7f81 604 //
ad670fba 605 // Linear fitter - not possible to make boundaries
75bd7f81 606 // non accept non possible z and dzdx combination
607 //
ad670fba 608 Bool_t acceptablez = kTRUE;
609 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
e4f2f73d 610 if (!cseed[iLayer].IsOK()) continue;
611 Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].fX0 - xref2);
612 if (TMath::Abs(cseed[iLayer].fZProb - zT2) > cseed[iLayer].fPadLength * 0.5 + 1.0) acceptablez = kFALSE;
75bd7f81 613 }
ad670fba 614 if (!acceptablez) {
615 Double_t zmf = cseed[2].fZref[0] + cseed[2].fZref[1] * (xref2 - cseed[2].fX0);
616 Double_t dzmf = (cseed[2].fZref[1] + cseed[3].fZref[1]) * 0.5;
75bd7f81 617 fitterT2.FixParameter(3,zmf);
618 fitterT2.FixParameter(4,dzmf);
619 fitterT2.Eval();
620 fitterT2.ReleaseParameter(3);
621 fitterT2.ReleaseParameter(4);
622 rpolz0 = fitterT2.GetParameter(3);
623 rpolz1 = fitterT2.GetParameter(4);
624 }
625
ad670fba 626 Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT);
75bd7f81 627 Double_t params[3];
ad670fba 628 params[0] = fitterT2.GetParameter(0);
629 params[1] = fitterT2.GetParameter(1);
630 params[2] = fitterT2.GetParameter(2);
631 Double_t curvature = 1.0 + params[1] * params[1] - params[2] * params[0];
632
e4f2f73d 633
ad670fba 634 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
ad670fba 635 Double_t x = cseed[iLayer].fX0;
636 Double_t y = 0;
637 Double_t dy = 0;
638 Double_t z = 0;
639 Double_t dz = 0;
640
75bd7f81 641 // y
ad670fba 642 Double_t res2 = (x * params[0] + params[1]);
643 res2 *= res2;
644 res2 = 1.0 - params[2]*params[0] + params[1]*params[1] - res2;
645 if (res2 >= 0) {
75bd7f81 646 res2 = TMath::Sqrt(res2);
ad670fba 647 y = (1.0 - res2) / params[0];
75bd7f81 648 }
ad670fba 649
75bd7f81 650 //dy
ad670fba 651 Double_t x0 = -params[1] / params[0];
652 if (-params[2]*params[0] + params[1]*params[1] + 1 > 0) {
653 Double_t rm1 = params[0] / TMath::Sqrt(-params[2]*params[0] + params[1]*params[1] + 1);
654 if (1.0/(rm1*rm1) - (x-x0) * (x-x0) > 0.0) {
60e55aee 655 Double_t res = (x - x0) / TMath::Sqrt((1./rm1-(x-x0))*(1./rm1+(x-x0)));
e4f2f73d 656 if (params[0] < 0) res *= -1.0;
657 dy = res;
75bd7f81 658 }
659 }
ad670fba 660 z = rpolz0 + rpolz1 * (x - xref2);
75bd7f81 661 dz = rpolz1;
662 cseed[iLayer].fYref[0] = y;
663 cseed[iLayer].fYref[1] = dy;
664 cseed[iLayer].fZref[0] = z;
665 cseed[iLayer].fZref[1] = dz;
ad670fba 666 cseed[iLayer].fC = curvature;
75bd7f81 667
668 }
669
670 return chi2TR;
671
672}