]> 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
New stand alone tracking by Alexadru and Martin
[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)
35
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
66 fClusters[i] = 0x0; // Clusters
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
e4f2f73d 128//_____________________________________________________________________________
129void AliTRDseed::Copy(TObject &o) const {
130 AliTRDseed &seed = (AliTRDseed &)o;
131 seed.fTilt = fTilt;
132 seed.fPadLength = fPadLength;
133 seed.fX0 = fX0;
134 seed.fSigmaY = fSigmaY;
135 seed.fSigmaY2 = fSigmaY2;
136 seed.fMeanz = fMeanz;
137 seed.fZProb = fZProb;
138 seed.fN = fN;
139 seed.fN2 = fN2;
140 seed.fNUsed = fNUsed;
141 seed.fFreq = fFreq;
142 seed.fNChange = fNChange;
143 seed.fMPads = fMPads;
144 seed.fC = fC;
145 seed.fCC = fCC;
146 seed.fChi2 = fChi2;
147 seed.fChi2Z = fChi2Z;
148 for (Int_t i = 0; i < knTimebins; i++) {
149 seed.fX[i] = fX[i];
150 seed.fY[i] = fY[i];
151 seed.fZ[i] = fZ[i];
152 seed.fIndexes[i] = fIndexes[i];
153 seed.fClusters[i] = fClusters[i];
154 seed.fUsable[i] = fUsable[i];
155 }
156
157 for (Int_t i = 0; i < 2; i++) {
158 seed.fYref[i] = fYref[i];
159 seed.fZref[i] = fZref[i];
160 seed.fYfit[i] = fYfit[i];
161 seed.fYfitR[i] = fYfitR[i];
162 seed.fZfit[i] = fZfit[i];
163 seed.fZfitR[i] = fZfitR[i];
164 seed.fLabels[i] = fLabels[i];
165 }
166
167 TObject::Copy(seed);
168
169}
170
75bd7f81 171//_____________________________________________________________________________
172void AliTRDseed::Reset()
173{
174 //
175 // Reset seed
176 //
177
e4f2f73d 178 for (Int_t i = 0; i < knTimebins; i++) {
ad670fba 179 fX[i] = 0; // X position
180 fY[i] = 0; // Y position
181 fZ[i] = 0; // Z position
182 fIndexes[i] = 0; // Indexes
e4f2f73d 183 fClusters[i] = 0x0; // Clusters
75bd7f81 184 fUsable[i] = kFALSE;
185 }
ad670fba 186
75bd7f81 187 for (Int_t i = 0; i < 2; i++) {
ad670fba 188 fYref[i] = 0; // Reference y
189 fZref[i] = 0; // Reference z
190 fYfit[i] = 0; // Y fit position +derivation
191 fYfitR[i] = 0; // Y fit position +derivation
192 fZfit[i] = 0; // Z fit position
193 fZfitR[i] = 0; // Z fit position
194 fLabels[i] = -1; // Labels
75bd7f81 195 }
ad670fba 196 fSigmaY = 0; // "Robust" sigma in y
197 fSigmaY2 = 0; // "Robust" sigma in y
198 fMeanz = 0; // Mean vaue of z
199 fZProb = 0; // Max probbable z
75bd7f81 200 fMPads = 0;
ad670fba 201 fN = 0; // Number of associated clusters
202 fN2 = 0; // Number of not crossed
203 fNUsed = 0; // Number of used clusters
204 fNChange = 0; // Change z counter
75bd7f81 205
206}
207
208//_____________________________________________________________________________
209void AliTRDseed::CookLabels()
210{
211 //
212 // Cook 2 labels for seed
213 //
214
e4f2f73d 215 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
216 Int_t nTimeBins = cal->GetNumberOfTimeBins();
217
75bd7f81 218 Int_t labels[200];
219 Int_t out[200];
ad670fba 220 Int_t nlab = 0;
221
e4f2f73d 222 for (Int_t i = 0; i < nTimeBins+1; i++) {
75bd7f81 223 if (!fClusters[i]) continue;
ad670fba 224 for (Int_t ilab = 0; ilab < 3; ilab++) {
225 if (fClusters[i]->GetLabel(ilab) >= 0) {
75bd7f81 226 labels[nlab] = fClusters[i]->GetLabel(ilab);
227 nlab++;
228 }
229 }
230 }
ad670fba 231
75bd7f81 232 Int_t nlab2 = AliTRDtracker::Freq(nlab,labels,out,kTRUE);
233 fLabels[0] = out[0];
ad670fba 234 if ((nlab2 > 1) &&
235 (out[3] > 1)) {
236 fLabels[1] = out[2];
237 }
75bd7f81 238
239}
240
241//_____________________________________________________________________________
242void AliTRDseed::UseClusters()
243{
244 //
245 // Use clusters
246 //
247
e4f2f73d 248 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
249 Int_t nTimeBins = cal->GetNumberOfTimeBins();
250
251 for (Int_t i = 0; i < nTimeBins+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
265 const Float_t kRatio = 0.8;
e4f2f73d 266 const Int_t kClmin = 5;
75bd7f81 267 const Float_t kmaxtan = 2;
268
e4f2f73d 269 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
270 Int_t nTimeBins = cal->GetNumberOfTimeBins();
271
272 if (TMath::Abs(fYref[1]) > kmaxtan){
273 //printf("Exit: Abs(fYref[1]) = %3.3f, kmaxtan = %3.3f\n", TMath::Abs(fYref[1]), kmaxtan);
274 return; // Track inclined too much
275 }
ad670fba 276
277 Float_t sigmaexp = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction
278 Float_t ycrosscor = fPadLength * fTilt * 0.5; // Y correction for crossing
279 fNChange = 0;
280
281 Double_t sumw;
282 Double_t sumwx;
283 Double_t sumwx2;
284 Double_t sumwy;
285 Double_t sumwxy;
286 Double_t sumwz;
287 Double_t sumwxz;
288
e4f2f73d 289 // Buffering: Leave it constant fot Performance issues
290 Int_t zints[knTimebins]; // Histograming of the z coordinate
ad670fba 291 // Get 1 and second max probable coodinates in z
e4f2f73d 292 Int_t zouts[2*knTimebins];
293 Float_t allowedz[knTimebins]; // Allowed z for given time bin
294 Float_t yres[knTimebins]; // Residuals from reference
ad670fba 295 Float_t anglecor = fTilt * fZref[1]; // Correction to the angle
75bd7f81 296
297
ad670fba 298 fN = 0;
299 fN2 = 0;
e4f2f73d 300 for (Int_t i = 0; i < nTimeBins; i++) {
ad670fba 301 yres[i] = 10000.0;
75bd7f81 302 if (!fClusters[i]) continue;
ad670fba 303 yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i]; // Residual y
75bd7f81 304 zints[fN] = Int_t(fZ[i]);
305 fN++;
306 }
307
e4f2f73d 308 if (fN < kClmin){
309 //printf("Exit fN < kClmin: fN = %d\n", fN);
310 return;
311 }
75bd7f81 312 Int_t nz = AliTRDtracker::Freq(fN,zints,zouts,kFALSE);
313 fZProb = zouts[0];
ad670fba 314 if (nz <= 1) zouts[3] = 0;
e4f2f73d 315 if (zouts[1] + zouts[3] < kClmin) {
316 //printf("Exit zouts[1] = %d, zouts[3] = %d\n",zouts[1],zouts[3]);
317 return;
318 }
75bd7f81 319
ad670fba 320 // Z distance bigger than pad - length
321 if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) {
322 zouts[3]=0;
323 }
75bd7f81 324
325 Int_t breaktime = -1;
326 Bool_t mbefore = kFALSE;
e4f2f73d 327 Int_t cumul[knTimebins][2];
ad670fba 328 Int_t counts[2] = { 0, 0 };
75bd7f81 329
ad670fba 330 if (zouts[3] >= 3) {
331
75bd7f81 332 //
ad670fba 333 // Find the break time allowing one chage on pad-rows
334 // with maximal numebr of accepted clusters
75bd7f81 335 //
ad670fba 336 fNChange = 1;
e4f2f73d 337 for (Int_t i = 0; i < nTimeBins; i++) {
75bd7f81 338 cumul[i][0] = counts[0];
339 cumul[i][1] = counts[1];
ad670fba 340 if (TMath::Abs(fZ[i]-zouts[0]) < 2) counts[0]++;
341 if (TMath::Abs(fZ[i]-zouts[2]) < 2) counts[1]++;
75bd7f81 342 }
ad670fba 343 Int_t maxcount = 0;
e4f2f73d 344 for (Int_t i = 0; i < nTimeBins; i++) {
345 Int_t after = cumul[nTimeBins][0] - cumul[i][0];
75bd7f81 346 Int_t before = cumul[i][1];
ad670fba 347 if (after + before > maxcount) {
348 maxcount = after + before;
349 breaktime = i;
350 mbefore = kFALSE;
75bd7f81 351 }
e4f2f73d 352 after = cumul[nTimeBins-1][1] - cumul[i][1];
75bd7f81 353 before = cumul[i][0];
ad670fba 354 if (after + before > maxcount) {
355 maxcount = after + before;
356 breaktime = i;
357 mbefore = kTRUE;
75bd7f81 358 }
359 }
ad670fba 360
361 breaktime -= 1;
362
75bd7f81 363 }
ad670fba 364
e4f2f73d 365 for (Int_t i = 0; i < nTimeBins+1; i++) {
ad670fba 366 if (i > breaktime) allowedz[i] = mbefore ? zouts[2] : zouts[0];
367 if (i <= breaktime) allowedz[i] = (!mbefore) ? zouts[2] : zouts[0];
75bd7f81 368 }
ad670fba 369
e4f2f73d 370 if (((allowedz[0] > allowedz[nTimeBins]) && (fZref[1] < 0)) ||
371 ((allowedz[0] < allowedz[nTimeBins]) && (fZref[1] > 0))) {
75bd7f81 372 //
ad670fba 373 // Tracklet z-direction not in correspondance with track z direction
75bd7f81 374 //
ad670fba 375 fNChange = 0;
e4f2f73d 376 for (Int_t i = 0; i < nTimeBins+1; i++) {
ad670fba 377 allowedz[i] = zouts[0]; // Only longest taken
75bd7f81 378 }
379 }
380
ad670fba 381 if (fNChange > 0) {
75bd7f81 382 //
ad670fba 383 // Cross pad -row tracklet - take the step change into account
75bd7f81 384 //
e4f2f73d 385 for (Int_t i = 0; i < nTimeBins+1; i++) {
75bd7f81 386 if (!fClusters[i]) continue;
ad670fba 387 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
388 yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i]; // Residual y
389 if (TMath::Abs(fZ[i] - fZProb) > 2) {
390 if (fZ[i] > fZProb) yres[i] += fTilt * fPadLength;
391 if (fZ[i] < fZProb) yres[i] -= fTilt * fPadLength;
75bd7f81 392 }
393 }
394 }
395
e4f2f73d 396 Double_t yres2[knTimebins];
ad670fba 397 Double_t mean;
398 Double_t sigma;
e4f2f73d 399 for (Int_t i = 0; i < nTimeBins+1; i++) {
75bd7f81 400 if (!fClusters[i]) continue;
ad670fba 401 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
402 yres2[fN2] = yres[i];
75bd7f81 403 fN2++;
404 }
ad670fba 405 if (fN2 < kClmin) {
e4f2f73d 406 //printf("Exit fN2 < kClmin: fN2 = %d\n", fN2);
75bd7f81 407 fN2 = 0;
408 return;
409 }
e4f2f73d 410 AliMathBase::EvaluateUni(fN2,yres2,mean,sigma, Int_t(fN2*kRatio-2.));
ad670fba 411 if (sigma < sigmaexp * 0.8) {
412 sigma = sigmaexp;
413 }
75bd7f81 414 fSigmaY = sigma;
ad670fba 415
416 // Reset sums
417 sumw = 0;
418 sumwx = 0;
419 sumwx2 = 0;
420 sumwy = 0;
421 sumwxy = 0;
422 sumwz = 0;
423 sumwxz = 0;
424
425 fN2 = 0;
426 fMeanz = 0;
427 fMPads = 0;
428
e4f2f73d 429 for (Int_t i = 0; i < nTimeBins+1; i++) {
ad670fba 430
431 fUsable[i] = kFALSE;
75bd7f81 432 if (!fClusters[i]) continue;
ad670fba 433 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
434 if (TMath::Abs(yres[i] - mean) > 4.0 * sigma) continue;
75bd7f81 435 fUsable[i] = kTRUE;
436 fN2++;
ad670fba 437 fMPads += fClusters[i]->GetNPads();
438 Float_t weight = 1.0;
439 if (fClusters[i]->GetNPads() > 4) weight = 0.5;
440 if (fClusters[i]->GetNPads() > 5) weight = 0.2;
441
75bd7f81 442 Double_t x = fX[i];
ad670fba 443 sumw += weight;
444 sumwx += x * weight;
445 sumwx2 += x*x * weight;
446 sumwy += weight * yres[i];
447 sumwxy += weight * (yres[i]) * x;
448 sumwz += weight * fZ[i];
449 sumwxz += weight * fZ[i] * x;
450
75bd7f81 451 }
ad670fba 452
453 if (fN2 < kClmin){
e4f2f73d 454 //printf("Exit fN2 < kClmin(2): fN2 = %d\n",fN2);
75bd7f81 455 fN2 = 0;
456 return;
457 }
ad670fba 458 fMeanz = sumwz / sumw;
459 Float_t correction = 0;
460 if (fNChange > 0) {
461 // Tracklet on boundary
462 if (fMeanz < fZProb) correction = ycrosscor;
463 if (fMeanz > fZProb) correction = -ycrosscor;
75bd7f81 464 }
ad670fba 465
466 Double_t det = sumw * sumwx2 - sumwx * sumwx;
467 fYfitR[0] = (sumwx2 * sumwy - sumwx * sumwxy) / det;
468 fYfitR[1] = (sumw * sumwxy - sumwx * sumwy) / det;
75bd7f81 469
ad670fba 470 fSigmaY2 = 0;
e4f2f73d 471 for (Int_t i = 0; i < nTimeBins+1; i++) {
75bd7f81 472 if (!fUsable[i]) continue;
ad670fba 473 Float_t delta = yres[i] - fYfitR[0] - fYfitR[1] * fX[i];
474 fSigmaY2 += delta*delta;
75bd7f81 475 }
ad670fba 476 fSigmaY2 = TMath::Sqrt(fSigmaY2 / Float_t(fN2-2));
477
478 fZfitR[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det;
479 fZfitR[1] = (sumw * sumwxz - sumwx * sumwz) / det;
480 fZfit[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det;
481 fZfit[1] = (sumw * sumwxz - sumwx * sumwz) / det;
482 fYfitR[0] += fYref[0] + correction;
483 fYfitR[1] += fYref[1];
484 fYfit[0] = fYfitR[0];
485 fYfit[1] = fYfitR[1];
75bd7f81 486
487 UpdateUsed();
488
489}
490
491//_____________________________________________________________________________
492void AliTRDseed::UpdateUsed()
493{
494 //
495 // Update used seed
496 //
497
e4f2f73d 498 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
499 Int_t nTimeBins = cal->GetNumberOfTimeBins();
500
75bd7f81 501 fNUsed = 0;
e4f2f73d 502 for (Int_t i = 0; i < nTimeBins; i++) {
95a5eb17 503 if (!fClusters[i]) {
504 continue;
505 }
506 if ((fClusters[i]->IsUsed())) {
507 fNUsed++;
508 }
75bd7f81 509 }
510
511}
512
75bd7f81 513//_____________________________________________________________________________
514Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror)
515{
516 //
517 // Fit the Rieman tilt
518 //
519
ad670fba 520 // Fitting with tilting pads - kz not fixed
521 TLinearFitter fitterT2(4,"hyp4");
75bd7f81 522 fitterT2.StoreData(kTRUE);
e4f2f73d 523
524 AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
525 Int_t nTimeBins = cal->GetNumberOfTimeBins();
526
ad670fba 527 Float_t xref2 = (cseed[2].fX0 + cseed[3].fX0) * 0.5; // Reference x0 for z
528
529 Int_t npointsT = 0;
75bd7f81 530 fitterT2.ClearPoints();
ad670fba 531
532 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
533
75bd7f81 534 if (!cseed[iLayer].IsOK()) continue;
535 Double_t tilt = cseed[iLayer].fTilt;
536
e4f2f73d 537 for (Int_t itime = 0; itime < nTimeBins+1; itime++) {
ad670fba 538
75bd7f81 539 if (!cseed[iLayer].fUsable[itime]) continue;
540 // x relative to the midle chamber
ad670fba 541 Double_t x = cseed[iLayer].fX[itime] + cseed[iLayer].fX0 - xref2;
542 Double_t y = cseed[iLayer].fY[itime];
543 Double_t z = cseed[iLayer].fZ[itime];
75bd7f81 544
545 //
ad670fba 546 // Tilted rieman
75bd7f81 547 //
548 Double_t uvt[6];
ad670fba 549 Double_t x2 = cseed[iLayer].fX[itime] + cseed[iLayer].fX0; // Global x
550 Double_t t = 1.0 / (x2*x2 + y*y);
551 uvt[1] = t;
552 uvt[0] = 2.0 * x2 * uvt[1];
553 uvt[2] = 2.0 * tilt * uvt[1];
554 uvt[3] = 2.0 * tilt *uvt[1] * x;
555 uvt[4] = 2.0 * (y + tilt * z) * uvt[1];
556
557 Double_t error = 2.0 * uvt[1];
558 if (terror) {
559 error *= cseed[iLayer].fSigmaY;
560 }
561 else {
562 error *= 0.2; //Default error
563 }
75bd7f81 564 fitterT2.AddPoint(uvt,uvt[4],error);
565 npointsT++;
ad670fba 566
75bd7f81 567 }
ad670fba 568
75bd7f81 569 }
ad670fba 570
75bd7f81 571 fitterT2.Eval();
572 Double_t rpolz0 = fitterT2.GetParameter(3);
573 Double_t rpolz1 = fitterT2.GetParameter(4);
ad670fba 574
75bd7f81 575 //
ad670fba 576 // Linear fitter - not possible to make boundaries
75bd7f81 577 // non accept non possible z and dzdx combination
578 //
ad670fba 579 Bool_t acceptablez = kTRUE;
580 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
e4f2f73d 581 if (!cseed[iLayer].IsOK()) continue;
582 Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].fX0 - xref2);
583 if (TMath::Abs(cseed[iLayer].fZProb - zT2) > cseed[iLayer].fPadLength * 0.5 + 1.0) acceptablez = kFALSE;
75bd7f81 584 }
ad670fba 585 if (!acceptablez) {
586 Double_t zmf = cseed[2].fZref[0] + cseed[2].fZref[1] * (xref2 - cseed[2].fX0);
587 Double_t dzmf = (cseed[2].fZref[1] + cseed[3].fZref[1]) * 0.5;
75bd7f81 588 fitterT2.FixParameter(3,zmf);
589 fitterT2.FixParameter(4,dzmf);
590 fitterT2.Eval();
591 fitterT2.ReleaseParameter(3);
592 fitterT2.ReleaseParameter(4);
593 rpolz0 = fitterT2.GetParameter(3);
594 rpolz1 = fitterT2.GetParameter(4);
595 }
596
ad670fba 597 Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT);
75bd7f81 598 Double_t params[3];
ad670fba 599 params[0] = fitterT2.GetParameter(0);
600 params[1] = fitterT2.GetParameter(1);
601 params[2] = fitterT2.GetParameter(2);
602 Double_t curvature = 1.0 + params[1] * params[1] - params[2] * params[0];
603
e4f2f73d 604
ad670fba 605 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
ad670fba 606 Double_t x = cseed[iLayer].fX0;
607 Double_t y = 0;
608 Double_t dy = 0;
609 Double_t z = 0;
610 Double_t dz = 0;
611
75bd7f81 612 // y
ad670fba 613 Double_t res2 = (x * params[0] + params[1]);
614 res2 *= res2;
615 res2 = 1.0 - params[2]*params[0] + params[1]*params[1] - res2;
616 if (res2 >= 0) {
75bd7f81 617 res2 = TMath::Sqrt(res2);
ad670fba 618 y = (1.0 - res2) / params[0];
75bd7f81 619 }
ad670fba 620
75bd7f81 621 //dy
ad670fba 622 Double_t x0 = -params[1] / params[0];
623 if (-params[2]*params[0] + params[1]*params[1] + 1 > 0) {
624 Double_t rm1 = params[0] / TMath::Sqrt(-params[2]*params[0] + params[1]*params[1] + 1);
625 if (1.0/(rm1*rm1) - (x-x0) * (x-x0) > 0.0) {
e4f2f73d 626 Double_t res = (x - x0) / TMath::Sqrt(1.0 / (rm1*rm1) - (x-x0)*(x-x0));
627 if (params[0] < 0) res *= -1.0;
628 dy = res;
75bd7f81 629 }
630 }
ad670fba 631 z = rpolz0 + rpolz1 * (x - xref2);
75bd7f81 632 dz = rpolz1;
633 cseed[iLayer].fYref[0] = y;
634 cseed[iLayer].fYref[1] = dy;
635 cseed[iLayer].fZref[0] = z;
636 cseed[iLayer].fZref[1] = dz;
ad670fba 637 cseed[iLayer].fC = curvature;
75bd7f81 638
639 }
640
641 return chi2TR;
642
643}