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