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