]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - TRD/AliTRDseed.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Coding rules
[u/mrichter/AliRoot.git] / TRD / AliTRDseed.cxx
... / ...
CommitLineData
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// //
22///////////////////////////////////////////////////////////////////////////////
23
24#include "TMath.h"
25#include "TLinearFitter.h"
26
27#include "AliMathBase.h"
28
29#include "AliTRDseed.h"
30#include "AliTRDcalibDB.h"
31#include "AliTRDcluster.h"
32#include "AliTRDtracker.h"
33#include "AliTRDtrackerV1.h"
34
35ClassImp(AliTRDseed)
36
37//_____________________________________________________________________________
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)
57{
58 //
59 // Default constructor
60 //
61
62 for (Int_t i = 0; i < knTimebins; i++) {
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)
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)
103{
104 //
105 // Copy constructor
106 //
107
108 for (Int_t i = 0; i < knTimebins; i++) {
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
115 }
116
117 for (Int_t i = 0; i < 2; i++) {
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
125 }
126
127}
128
129//_____________________________________________________________________________
130void AliTRDseed::Copy(TObject &o) const
131{
132 //printf("AliTRDseed::Copy()\n");
133
134 AliTRDseed &seed = (AliTRDseed &)o;
135
136 seed.fTilt = fTilt;
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;
153 for (Int_t i = 0; i < knTimebins; i++) {
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
176//_____________________________________________________________________________
177void AliTRDseed::Reset()
178{
179 //
180 // Reset seed
181 //
182
183 for (Int_t i = 0; i < knTimebins; i++) {
184 fX[i] = 0; // X position
185 fY[i] = 0; // Y position
186 fZ[i] = 0; // Z position
187 fIndexes[i] = 0; // Indexes
188 fClusters[i] = 0x0; // Clusters
189 fUsable[i] = kFALSE;
190 }
191
192 for (Int_t i = 0; i < 2; i++) {
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
200 }
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
205 fMPads = 0;
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
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];
222 Int_t nlab = 0;
223
224 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
225 if (!fClusters[i]) continue;
226 for (Int_t ilab = 0; ilab < 3; ilab++) {
227 if (fClusters[i]->GetLabel(ilab) >= 0) {
228 labels[nlab] = fClusters[i]->GetLabel(ilab);
229 nlab++;
230 }
231 }
232 }
233
234 Int_t nlab2 = AliTRDtracker::Freq(nlab,labels,out,kTRUE);
235 fLabels[0] = out[0];
236 if ((nlab2 > 1) &&
237 (out[3] > 1)) {
238 fLabels[1] = out[2];
239 }
240
241}
242
243//_____________________________________________________________________________
244void AliTRDseed::UseClusters()
245{
246 //
247 // Use clusters
248 //
249
250 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
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
264
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.
278
279 const Float_t kRatio = 0.8;
280 const Int_t kClmin = 5;
281 const Float_t kmaxtan = 2;
282
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 }
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
300 // Buffering: Leave it constant fot Performance issues
301 Int_t zints[knTimebins]; // Histograming of the z coordinate
302 // Get 1 and second max probable coodinates in z
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
306 //Float_t anglecor = fTilt * fZref[1]; // Correction to the angle
307
308
309 fN = 0;
310 fN2 = 0;
311 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
312 yres[i] = 10000.0;
313 if (!fClusters[i]) continue;
314 if(!fClusters[i]->IsInChamber()) continue;
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]));
318 zints[fN] = Int_t(fZ[i]);
319 fN++;
320 }
321
322 if (fN < kClmin){
323 //printf("Exit fN < kClmin: fN = %d\n", fN);
324 return;
325 }
326 Int_t nz = AliTRDtracker::Freq(fN, zints, zouts, kFALSE);
327 fZProb = zouts[0];
328 if (nz <= 1) zouts[3] = 0;
329 if (zouts[1] + zouts[3] < kClmin) {
330 //printf("Exit zouts[1] = %d, zouts[3] = %d\n",zouts[1],zouts[3]);
331 return;
332 }
333
334 // Z distance bigger than pad - length
335 if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) zouts[3] = 0;
336
337 Int_t breaktime = -1;
338 Bool_t mbefore = kFALSE;
339 Int_t cumul[knTimebins][2];
340 Int_t counts[2] = { 0, 0 };
341
342 if (zouts[3] >= 3) {
343
344 //
345 // Find the break time allowing one chage on pad-rows
346 // with maximal number of accepted clusters
347 //
348 fNChange = 1;
349 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
350 cumul[i][0] = counts[0];
351 cumul[i][1] = counts[1];
352 if (TMath::Abs(fZ[i]-zouts[0]) < 2) counts[0]++;
353 if (TMath::Abs(fZ[i]-zouts[2]) < 2) counts[1]++;
354 }
355 Int_t maxcount = 0;
356 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
357 Int_t after = cumul[AliTRDtrackerV1::GetNTimeBins()][0] - cumul[i][0];
358 Int_t before = cumul[i][1];
359 if (after + before > maxcount) {
360 maxcount = after + before;
361 breaktime = i;
362 mbefore = kFALSE;
363 }
364 after = cumul[AliTRDtrackerV1::GetNTimeBins()-1][1] - cumul[i][1];
365 before = cumul[i][0];
366 if (after + before > maxcount) {
367 maxcount = after + before;
368 breaktime = i;
369 mbefore = kTRUE;
370 }
371 }
372
373 breaktime -= 1;
374
375 }
376
377 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
378 if (i > breaktime) allowedz[i] = mbefore ? zouts[2] : zouts[0];
379 if (i <= breaktime) allowedz[i] = (!mbefore) ? zouts[2] : zouts[0];
380 }
381
382 if (((allowedz[0] > allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] < 0)) ||
383 ((allowedz[0] < allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] > 0))) {
384 //
385 // Tracklet z-direction not in correspondance with track z direction
386 //
387 fNChange = 0;
388 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
389 allowedz[i] = zouts[0]; // Only longest taken
390 }
391 }
392
393 if (fNChange > 0) {
394 //
395 // Cross pad -row tracklet - take the step change into account
396 //
397 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
398 if (!fClusters[i]) continue;
399 if(!fClusters[i]->IsInChamber()) continue;
400 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
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 }*/
408 }
409 }
410
411 Double_t yres2[knTimebins];
412 Double_t mean;
413 Double_t sigma;
414 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
415 if (!fClusters[i]) continue;
416 if(!fClusters[i]->IsInChamber()) continue;
417 if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
418 yres2[fN2] = yres[i];
419 fN2++;
420 }
421 if (fN2 < kClmin) {
422 //printf("Exit fN2 < kClmin: fN2 = %d\n", fN2);
423 fN2 = 0;
424 return;
425 }
426 AliMathBase::EvaluateUni(fN2,yres2,mean,sigma, Int_t(fN2*kRatio-2.));
427 if (sigma < sigmaexp * 0.8) {
428 sigma = sigmaexp;
429 }
430 fSigmaY = sigma;
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
445 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
446
447 fUsable[i] = kFALSE;
448 if (!fClusters[i]) continue;
449 if (!fClusters[i]->IsInChamber()) continue;
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;}
452 fUsable[i] = kTRUE;
453 fN2++;
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
459
460 Double_t x = fX[i];
461 //printf("x = %7.3f dy = %7.3f fit %7.3f\n", x, yres[i], fY[i]-yres[i]);
462
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
471 }
472
473 if (fN2 < kClmin){
474 //printf("Exit fN2 < kClmin(2): fN2 = %d\n",fN2);
475 fN2 = 0;
476 return;
477 }
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;
484 }
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;
489
490 fSigmaY2 = 0;
491 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
492 if (!fUsable[i]) continue;
493 Float_t delta = yres[i] - fYfitR[0] - fYfitR[1] * fX[i];
494 fSigmaY2 += delta*delta;
495 }
496 fSigmaY2 = TMath::Sqrt(fSigmaY2 / Float_t(fN2-2));
497 // TEMPORARY UNTIL covariance properly calculated
498 fSigmaY2 = TMath::Max(fSigmaY2, Float_t(.1));
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;
504// fYfitR[0] += fYref[0] + correction;
505// fYfitR[1] += fYref[1];
506 fYfit[0] = fYfitR[0];
507 fYfit[1] = -fYfitR[1];
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
511 UpdateUsed();
512
513}
514
515//_____________________________________________________________________________
516void AliTRDseed::UpdateUsed()
517{
518 //
519 // Update used seed
520 //
521
522 fNUsed = 0;
523 for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
524 if (!fClusters[i]) continue;
525 if(!fUsable[i]) continue;
526 if ((fClusters[i]->IsUsed())) fNUsed++;
527 }
528
529}
530
531//_____________________________________________________________________________
532Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror)
533{
534 //
535 // Fit the Rieman tilt
536 //
537
538 // Fitting with tilting pads - kz not fixed
539 TLinearFitter fitterT2(4,"hyp4");
540 fitterT2.StoreData(kTRUE);
541
542 Float_t xref2 = (cseed[2].fX0 + cseed[3].fX0) * 0.5; // Reference x0 for z
543
544 Int_t npointsT = 0;
545 fitterT2.ClearPoints();
546
547 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
548
549 if (!cseed[iLayer].IsOK()) continue;
550 Double_t tilt = cseed[iLayer].fTilt;
551
552 for (Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins()+1; itime++) {
553
554 if (!cseed[iLayer].fUsable[itime]) continue;
555 // x relative to the midle chamber
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];
559
560 //
561 // Tilted rieman
562 //
563 Double_t uvt[6];
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 }
579 fitterT2.AddPoint(uvt,uvt[4],error);
580 npointsT++;
581
582 }
583
584 }
585
586 fitterT2.Eval();
587 Double_t rpolz0 = fitterT2.GetParameter(3);
588 Double_t rpolz1 = fitterT2.GetParameter(4);
589
590 //
591 // Linear fitter - not possible to make boundaries
592 // non accept non possible z and dzdx combination
593 //
594 Bool_t acceptablez = kTRUE;
595 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
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;
599 }
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;
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
612 Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT);
613 Double_t params[3];
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
619
620 for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
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
627 // y
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) {
632 res2 = TMath::Sqrt(res2);
633 y = (1.0 - res2) / params[0];
634 }
635
636 //dy
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) {
641 Double_t res = (x - x0) / TMath::Sqrt((1./rm1-(x-x0))*(1./rm1+(x-x0)));
642 if (params[0] < 0) res *= -1.0;
643 dy = res;
644 }
645 }
646 z = rpolz0 + rpolz1 * (x - xref2);
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;
652 cseed[iLayer].fC = curvature;
653
654 }
655
656 return chi2TR;
657
658}