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