]> 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
A new method DrawPMDModule is added
[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
284 if (TMath::Abs(fYref[1]) > kmaxtan){
285 //printf("Exit: Abs(fYref[1]) = %3.3f, kmaxtan = %3.3f\n", TMath::Abs(fYref[1]), kmaxtan);
286 return; // Track inclined too much
287 }
ad670fba 288
289 Float_t sigmaexp = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction
290 Float_t ycrosscor = fPadLength * fTilt * 0.5; // Y correction for crossing
291 fNChange = 0;
292
293 Double_t sumw;
294 Double_t sumwx;
295 Double_t sumwx2;
296 Double_t sumwy;
297 Double_t sumwxy;
298 Double_t sumwz;
299 Double_t sumwxz;
300
e4f2f73d 301 // Buffering: Leave it constant fot Performance issues
302 Int_t zints[knTimebins]; // Histograming of the z coordinate
ad670fba 303 // Get 1 and second max probable coodinates in z
e4f2f73d 304 Int_t zouts[2*knTimebins];
305 Float_t allowedz[knTimebins]; // Allowed z for given time bin
306 Float_t yres[knTimebins]; // Residuals from reference
ad670fba 307 Float_t anglecor = fTilt * fZref[1]; // Correction to the angle
75bd7f81 308
309
ad670fba 310 fN = 0;
311 fN2 = 0;
2985ffcb 312 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
ad670fba 313 yres[i] = 10000.0;
75bd7f81 314 if (!fClusters[i]) continue;
bcb6fb78 315 if(!fClusters[i]->IsInChamber()) continue;
316 yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] + fTilt*(fZ[i] - fZref[0]); // Residual y
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) {
342
75bd7f81 343 //
ad670fba 344 // Find the break time allowing one chage on pad-rows
bcb6fb78 345 // with maximal number of accepted clusters
75bd7f81 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++) {
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) {
359 maxcount = after + before;
360 breaktime = i;
361 mbefore = kFALSE;
75bd7f81 362 }
2985ffcb 363 after = cumul[AliTRDtrackerV1::GetNTimeBins()-1][1] - cumul[i][1];
75bd7f81 364 before = cumul[i][0];
ad670fba 365 if (after + before > maxcount) {
366 maxcount = after + before;
367 breaktime = i;
368 mbefore = kTRUE;
75bd7f81 369 }
370 }
ad670fba 371
372 breaktime -= 1;
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;
2985ffcb 387 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
ad670fba 388 allowedz[i] = zouts[0]; // Only longest taken
75bd7f81 389 }
390 }
391
ad670fba 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;
bcb6fb78 400 yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] /*+ fTilt*(fZ[i] - fZref[0])*/; // Residual y
ad670fba 401 if (TMath::Abs(fZ[i] - fZProb) > 2) {
402 if (fZ[i] > fZProb) yres[i] += fTilt * fPadLength;
403 if (fZ[i] < fZProb) yres[i] -= fTilt * fPadLength;
75bd7f81 404 }
405 }
406 }
407
e4f2f73d 408 Double_t yres2[knTimebins];
ad670fba 409 Double_t mean;
410 Double_t sigma;
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;
415 yres2[fN2] = yres[i];
75bd7f81 416 fN2++;
417 }
ad670fba 418 if (fN2 < kClmin) {
e4f2f73d 419 //printf("Exit fN2 < kClmin: fN2 = %d\n", fN2);
75bd7f81 420 fN2 = 0;
421 return;
422 }
e4f2f73d 423 AliMathBase::EvaluateUni(fN2,yres2,mean,sigma, Int_t(fN2*kRatio-2.));
ad670fba 424 if (sigma < sigmaexp * 0.8) {
425 sigma = sigmaexp;
426 }
75bd7f81 427 fSigmaY = sigma;
ad670fba 428
429 // Reset sums
430 sumw = 0;
431 sumwx = 0;
432 sumwx2 = 0;
433 sumwy = 0;
434 sumwxy = 0;
435 sumwz = 0;
436 sumwxz = 0;
437
438 fN2 = 0;
439 fMeanz = 0;
440 fMPads = 0;
441
2985ffcb 442 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
ad670fba 443
444 fUsable[i] = kFALSE;
75bd7f81 445 if (!fClusters[i]) continue;
bcb6fb78 446 if (!fClusters[i]->IsInChamber()) continue;
0906e73e 447 if (TMath::Abs(fZ[i] - allowedz[i]) > 2){fClusters[i] = 0x0; continue;}
448 if (TMath::Abs(yres[i] - mean) > 4.0 * sigma){fClusters[i] = 0x0; continue;}
75bd7f81 449 fUsable[i] = kTRUE;
450 fN2++;
ad670fba 451 fMPads += fClusters[i]->GetNPads();
452 Float_t weight = 1.0;
453 if (fClusters[i]->GetNPads() > 4) weight = 0.5;
454 if (fClusters[i]->GetNPads() > 5) weight = 0.2;
455
0906e73e 456
75bd7f81 457 Double_t x = fX[i];
0906e73e 458 //printf("x = %7.3f dy = %7.3f fit %7.3f\n", x, yres[i], fY[i]-yres[i]);
459
ad670fba 460 sumw += weight;
461 sumwx += x * weight;
462 sumwx2 += x*x * weight;
463 sumwy += weight * yres[i];
464 sumwxy += weight * (yres[i]) * x;
465 sumwz += weight * fZ[i];
466 sumwxz += weight * fZ[i] * x;
467
75bd7f81 468 }
ad670fba 469
470 if (fN2 < kClmin){
e4f2f73d 471 //printf("Exit fN2 < kClmin(2): fN2 = %d\n",fN2);
75bd7f81 472 fN2 = 0;
473 return;
474 }
ad670fba 475 fMeanz = sumwz / sumw;
476 Float_t correction = 0;
477 if (fNChange > 0) {
478 // Tracklet on boundary
479 if (fMeanz < fZProb) correction = ycrosscor;
480 if (fMeanz > fZProb) correction = -ycrosscor;
75bd7f81 481 }
ad670fba 482
483 Double_t det = sumw * sumwx2 - sumwx * sumwx;
484 fYfitR[0] = (sumwx2 * sumwy - sumwx * sumwxy) / det;
485 fYfitR[1] = (sumw * sumwxy - sumwx * sumwy) / det;
75bd7f81 486
ad670fba 487 fSigmaY2 = 0;
2985ffcb 488 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
75bd7f81 489 if (!fUsable[i]) continue;
ad670fba 490 Float_t delta = yres[i] - fYfitR[0] - fYfitR[1] * fX[i];
491 fSigmaY2 += delta*delta;
75bd7f81 492 }
ad670fba 493 fSigmaY2 = TMath::Sqrt(fSigmaY2 / Float_t(fN2-2));
0906e73e 494 // TEMPORARY UNTIL covariance properly calculated
495 fSigmaY2 = TMath::Max(fSigmaY2, Float_t(.1));
ad670fba 496
497 fZfitR[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det;
498 fZfitR[1] = (sumw * sumwxz - sumwx * sumwz) / det;
499 fZfit[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det;
500 fZfit[1] = (sumw * sumwxz - sumwx * sumwz) / det;
501 fYfitR[0] += fYref[0] + correction;
502 fYfitR[1] += fYref[1];
503 fYfit[0] = fYfitR[0];
504 fYfit[1] = fYfitR[1];
0906e73e 505
506 //printf("y0 = %7.3f tgy = %7.3f z0 = %7.3f tgz = %7.3f \n", fYfitR[0], fYfitR[1], fZfitR[0], fZfitR[1]);
507
75bd7f81 508 UpdateUsed();
509
510}
511
512//_____________________________________________________________________________
513void AliTRDseed::UpdateUsed()
514{
515 //
516 // Update used seed
517 //
518
519 fNUsed = 0;
2985ffcb 520 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
0906e73e 521 if (!fClusters[i]) continue;
522 if(!fUsable[i]) continue;
523 if ((fClusters[i]->IsUsed())) fNUsed++;
75bd7f81 524 }
525
526}
527
75bd7f81 528//_____________________________________________________________________________
529Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror)
530{
531 //
532 // Fit the Rieman tilt
533 //
534
ad670fba 535 // Fitting with tilting pads - kz not fixed
536 TLinearFitter fitterT2(4,"hyp4");
75bd7f81 537 fitterT2.StoreData(kTRUE);
e4f2f73d 538
ad670fba 539 Float_t xref2 = (cseed[2].fX0 + cseed[3].fX0) * 0.5; // Reference x0 for z
540
541 Int_t npointsT = 0;
75bd7f81 542 fitterT2.ClearPoints();
ad670fba 543
544 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
545
75bd7f81 546 if (!cseed[iLayer].IsOK()) continue;
547 Double_t tilt = cseed[iLayer].fTilt;
548
2985ffcb 549 for (Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins()+1; itime++) {
ad670fba 550
75bd7f81 551 if (!cseed[iLayer].fUsable[itime]) continue;
552 // x relative to the midle chamber
ad670fba 553 Double_t x = cseed[iLayer].fX[itime] + cseed[iLayer].fX0 - xref2;
554 Double_t y = cseed[iLayer].fY[itime];
555 Double_t z = cseed[iLayer].fZ[itime];
75bd7f81 556
557 //
ad670fba 558 // Tilted rieman
75bd7f81 559 //
560 Double_t uvt[6];
ad670fba 561 Double_t x2 = cseed[iLayer].fX[itime] + cseed[iLayer].fX0; // Global x
562 Double_t t = 1.0 / (x2*x2 + y*y);
563 uvt[1] = t;
564 uvt[0] = 2.0 * x2 * uvt[1];
565 uvt[2] = 2.0 * tilt * uvt[1];
566 uvt[3] = 2.0 * tilt *uvt[1] * x;
567 uvt[4] = 2.0 * (y + tilt * z) * uvt[1];
568
569 Double_t error = 2.0 * uvt[1];
570 if (terror) {
571 error *= cseed[iLayer].fSigmaY;
572 }
573 else {
574 error *= 0.2; //Default error
575 }
75bd7f81 576 fitterT2.AddPoint(uvt,uvt[4],error);
577 npointsT++;
ad670fba 578
75bd7f81 579 }
ad670fba 580
75bd7f81 581 }
ad670fba 582
75bd7f81 583 fitterT2.Eval();
584 Double_t rpolz0 = fitterT2.GetParameter(3);
585 Double_t rpolz1 = fitterT2.GetParameter(4);
ad670fba 586
75bd7f81 587 //
ad670fba 588 // Linear fitter - not possible to make boundaries
75bd7f81 589 // non accept non possible z and dzdx combination
590 //
ad670fba 591 Bool_t acceptablez = kTRUE;
592 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
e4f2f73d 593 if (!cseed[iLayer].IsOK()) continue;
594 Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].fX0 - xref2);
595 if (TMath::Abs(cseed[iLayer].fZProb - zT2) > cseed[iLayer].fPadLength * 0.5 + 1.0) acceptablez = kFALSE;
75bd7f81 596 }
ad670fba 597 if (!acceptablez) {
598 Double_t zmf = cseed[2].fZref[0] + cseed[2].fZref[1] * (xref2 - cseed[2].fX0);
599 Double_t dzmf = (cseed[2].fZref[1] + cseed[3].fZref[1]) * 0.5;
75bd7f81 600 fitterT2.FixParameter(3,zmf);
601 fitterT2.FixParameter(4,dzmf);
602 fitterT2.Eval();
603 fitterT2.ReleaseParameter(3);
604 fitterT2.ReleaseParameter(4);
605 rpolz0 = fitterT2.GetParameter(3);
606 rpolz1 = fitterT2.GetParameter(4);
607 }
608
ad670fba 609 Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT);
75bd7f81 610 Double_t params[3];
ad670fba 611 params[0] = fitterT2.GetParameter(0);
612 params[1] = fitterT2.GetParameter(1);
613 params[2] = fitterT2.GetParameter(2);
614 Double_t curvature = 1.0 + params[1] * params[1] - params[2] * params[0];
615
e4f2f73d 616
ad670fba 617 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
ad670fba 618 Double_t x = cseed[iLayer].fX0;
619 Double_t y = 0;
620 Double_t dy = 0;
621 Double_t z = 0;
622 Double_t dz = 0;
623
75bd7f81 624 // y
ad670fba 625 Double_t res2 = (x * params[0] + params[1]);
626 res2 *= res2;
627 res2 = 1.0 - params[2]*params[0] + params[1]*params[1] - res2;
628 if (res2 >= 0) {
75bd7f81 629 res2 = TMath::Sqrt(res2);
ad670fba 630 y = (1.0 - res2) / params[0];
75bd7f81 631 }
ad670fba 632
75bd7f81 633 //dy
ad670fba 634 Double_t x0 = -params[1] / params[0];
635 if (-params[2]*params[0] + params[1]*params[1] + 1 > 0) {
636 Double_t rm1 = params[0] / TMath::Sqrt(-params[2]*params[0] + params[1]*params[1] + 1);
637 if (1.0/(rm1*rm1) - (x-x0) * (x-x0) > 0.0) {
e4f2f73d 638 Double_t res = (x - x0) / TMath::Sqrt(1.0 / (rm1*rm1) - (x-x0)*(x-x0));
639 if (params[0] < 0) res *= -1.0;
640 dy = res;
75bd7f81 641 }
642 }
ad670fba 643 z = rpolz0 + rpolz1 * (x - xref2);
75bd7f81 644 dz = rpolz1;
645 cseed[iLayer].fYref[0] = y;
646 cseed[iLayer].fYref[1] = dy;
647 cseed[iLayer].fZref[0] = z;
648 cseed[iLayer].fZref[1] = dz;
ad670fba 649 cseed[iLayer].fC = curvature;
75bd7f81 650
651 }
652
653 return chi2TR;
654
655}