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