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