]> 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
move warning to debug message
[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
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//_____________________________________________________________________________
531void 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//_____________________________________________________________________________
547Float_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}