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