TRD/AliTRDseed.cxx

   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
  34 ClassImp(AliTRDseed)
  35
  36 Int_t AliTRDseed::fgTimeBins = 0;
  37
  38 //_____________________________________________________________________________
  39 AliTRDseed::AliTRDseed()
  40   :TObject()
  41   ,fTimeBinsRange(0)
  42   ,fTimeBin0(0)
  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)
  60 {
  61   //
  62   // Default constructor
  63   //
  64
  65   for (Int_t i = 0; i < knTimebins; i++) {
  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 //_____________________________________________________________________________
  87 AliTRDseed::AliTRDseed(const AliTRDseed &s)
  88   :TObject(s)
  89   ,fTimeBinsRange(s.fTimeBinsRange)
  90   ,fTimeBin0(s.fTimeBin0)
  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)
 108 {
 109   //
 110   // Copy constructor
 111   //
 112
 113   for (Int_t i = 0; i < knTimebins; i++) {
 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
 120   }
 121
 122   for (Int_t i = 0; i < 2; i++) {
 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
 130   }
 131
 132 }
 133
 134 //_____________________________________________________________________________
 135 void AliTRDseed::Copy(TObject &o) const
 136 {
 137         //printf("AliTRDseed::Copy()\n");
 138
 139         AliTRDseed &seed = (AliTRDseed &)o;
 140
 141         seed.fTimeBinsRange = fTimeBinsRange;
 142         seed.fTimeBin0 = fTimeBin0;
 143         seed.fTilt = fTilt;
 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;
 160         for (Int_t i = 0; i < knTimebins; i++) {
 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
 183 //_____________________________________________________________________________
 184 void AliTRDseed::Reset()
 185 {
 186   //
 187   // Reset seed
 188   //
 189
 190   for (Int_t i = 0; i < knTimebins; i++) {
 191     fX[i]        = 0;  // X position
 192     fY[i]        = 0;  // Y position
 193     fZ[i]        = 0;  // Z position
 194     fIndexes[i]  = 0;  // Indexes
 195     fClusters[i] = 0x0;  // Clusters
 196     fUsable[i]   = kFALSE;
 197   }
 198
 199   for (Int_t i = 0; i < 2; i++) {
 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
 207   }
 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
 212   fMPads   = 0;
 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
 217
 218 }
 219
 220 //_____________________________________________________________________________
 221 void AliTRDseed::CookLabels()
 222 {
 223   //
 224   // Cook 2 labels for seed
 225   //
 226
 227   Int_t labels[200];
 228   Int_t out[200];
 229   Int_t nlab = 0;
 230
 231   for (Int_t i = 0; i < fgTimeBins+1; i++) {
 232     if (!fClusters[i]) continue;
 233     for (Int_t ilab = 0; ilab < 3; ilab++) {
 234       if (fClusters[i]->GetLabel(ilab) >= 0) {
 235         labels[nlab] = fClusters[i]->GetLabel(ilab);
 236         nlab++;
 237       }
 238     }
 239   }
 240
 241   Int_t nlab2 = AliTRDtracker::Freq(nlab,labels,out,kTRUE);
 242   fLabels[0] = out[0];
 243   if ((nlab2  > 1) &&
 244       (out[3] > 1)) {
 245     fLabels[1] = out[2];
 246   }
 247
 248 }
 249
 250 //_____________________________________________________________________________
 251 void AliTRDseed::UseClusters()
 252 {
 253   //
 254   // Use clusters
 255   //
 256
 257   for (Int_t i = 0; i < fgTimeBins+1; i++) {
 258     if (!fClusters[i]) continue;
 259     if (!(fClusters[i]->IsUsed())) fClusters[i]->Use();
 260   }
 261
 262 }
 263
 264 //_____________________________________________________________________________
 265 void AliTRDseed::Update()
 266 {
 267   //
 268   // Update the seed.
 269   //
 270
 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
 276   const Float_t kRatio  = 0.8;
 277   const Int_t   kClmin  = 5;
 278   const Float_t kmaxtan = 2;
 279
 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         }
 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
 298         // Buffering: Leave it constant fot Performance issues
 299   Int_t    zints[knTimebins];            // Histograming of the z coordinate
 300                                          // Get 1 and second max probable coodinates in z
 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
 304   Float_t  anglecor = fTilt * fZref[1];  // Correction to the angle
 305
 306
 307   fN  = 0;
 308   fN2 = 0;
 309   for (Int_t i = 0; i < fTimeBinsRange; i++) {
 310     yres[i] = 10000.0;
 311     if (!fClusters[i]) continue;
 312     yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i];   // Residual y
 313     zints[fN] = Int_t(fZ[i]);
 314     fN++;
 315   }
 316
 317   if (fN < kClmin){
 318                 //printf("Exit fN < kClmin: fN = %d\n", fN);
 319                 return;
 320         }
 321   Int_t nz = AliTRDtracker::Freq(fN,zints,zouts,kFALSE);
 322   fZProb   = zouts[0];
 323   if (nz <= 1) zouts[3] = 0;
 324   if (zouts[1] + zouts[3] < kClmin) {
 325                 //printf("Exit zouts[1] = %d, zouts[3] = %d\n",zouts[1],zouts[3]);
 326                 return;
 327         }
 328
 329   // Z distance bigger than pad - length
 330   if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) {
 331     zouts[3]=0;
 332   }
 333
 334   Int_t  breaktime = -1;
 335   Bool_t mbefore   = kFALSE;
 336   Int_t  cumul[knTimebins][2];
 337   Int_t  counts[2] = { 0, 0 };
 338
 339   if (zouts[3] >= 3) {
 340
 341     //
 342     // Find the break time allowing one chage on pad-rows
 343     // with maximal numebr of accepted clusters
 344     //
 345     fNChange = 1;
 346     for (Int_t i = 0; i < fTimeBinsRange; i++) {
 347       cumul[i][0] = counts[0];
 348       cumul[i][1] = counts[1];
 349       if (TMath::Abs(fZ[i]-zouts[0]) < 2) counts[0]++;
 350       if (TMath::Abs(fZ[i]-zouts[2]) < 2) counts[1]++;
 351     }
 352     Int_t  maxcount = 0;
 353     for (Int_t i = 0; i < fTimeBinsRange; i++) {
 354       Int_t after  = cumul[fTimeBinsRange][0] - cumul[i][0];
 355       Int_t before = cumul[i][1];
 356       if (after + before > maxcount) {
 357         maxcount  = after + before;
 358         breaktime = i;
 359         mbefore   = kFALSE;
 360       }
 361       after  = cumul[fTimeBinsRange-1][1] - cumul[i][1];
 362       before = cumul[i][0];
 363       if (after + before > maxcount) {
 364         maxcount  = after + before;
 365         breaktime = i;
 366         mbefore   = kTRUE;
 367       }
 368     }
 369
 370     breaktime -= 1;
 371
 372   }
 373
 374   for (Int_t i = 0; i < fTimeBinsRange+1; i++) {
 375     if (i >  breaktime) allowedz[i] =   mbefore  ? zouts[2] : zouts[0];
 376     if (i <= breaktime) allowedz[i] = (!mbefore) ? zouts[2] : zouts[0];
 377   }
 378
 379   if (((allowedz[0] > allowedz[fTimeBinsRange]) && (fZref[1] < 0)) ||
 380       ((allowedz[0] < allowedz[fTimeBinsRange]) && (fZref[1] > 0))) {
 381     //
 382     // Tracklet z-direction not in correspondance with track z direction
 383     //
 384     fNChange = 0;
 385     for (Int_t i = 0; i < fTimeBinsRange+1; i++) {
 386       allowedz[i] = zouts[0];  // Only longest taken
 387     }
 388   }
 389
 390   if (fNChange > 0) {
 391     //
 392     // Cross pad -row tracklet  - take the step change into account
 393     //
 394     for (Int_t i = 0; i < fTimeBinsRange+1; i++) {
 395       if (!fClusters[i]) continue;
 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;
 401       }
 402     }
 403   }
 404
 405   Double_t yres2[knTimebins];
 406   Double_t mean;
 407   Double_t sigma;
 408   for (Int_t i = 0; i < fTimeBinsRange+1; i++) {
 409     if (!fClusters[i]) continue;
 410     if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
 411     yres2[fN2] = yres[i];
 412     fN2++;
 413   }
 414   if (fN2 < kClmin) {
 415                 //printf("Exit fN2 < kClmin: fN2 = %d\n", fN2);
 416     fN2 = 0;
 417     return;
 418   }
 419   AliMathBase::EvaluateUni(fN2,yres2,mean,sigma, Int_t(fN2*kRatio-2.));
 420   if (sigma < sigmaexp * 0.8) {
 421     sigma = sigmaexp;
 422   }
 423   fSigmaY = sigma;
 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
 438   for (Int_t i = 0; i < fTimeBinsRange+1; i++) {
 439
 440     fUsable[i] = kFALSE;
 441     if (!fClusters[i]) continue;
 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;}
 444     fUsable[i] = kTRUE;
 445     fN2++;
 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
 451
 452     Double_t x = fX[i];
 453     //printf("x = %7.3f dy = %7.3f fit %7.3f\n", x, yres[i], fY[i]-yres[i]);
 454
 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
 463   }
 464
 465   if (fN2 < kClmin){
 466                 //printf("Exit fN2 < kClmin(2): fN2 = %d\n",fN2);
 467     fN2 = 0;
 468     return;
 469   }
 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;
 476   }
 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;
 481
 482   fSigmaY2 = 0;
 483   for (Int_t i = 0; i < fTimeBinsRange+1; i++) {
 484     if (!fUsable[i]) continue;
 485     Float_t delta = yres[i] - fYfitR[0] - fYfitR[1] * fX[i];
 486     fSigmaY2 += delta*delta;
 487   }
 488   fSigmaY2 = TMath::Sqrt(fSigmaY2 / Float_t(fN2-2));
 489         // TEMPORARY UNTIL covariance properly calculated
 490         fSigmaY2 = TMath::Max(fSigmaY2, Float_t(.1));
 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];
 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
 503   UpdateUsed();
 504
 505 }
 506
 507 //_____________________________________________________________________________
 508 void AliTRDseed::UpdateUsed()
 509 {
 510   //
 511   // Update used seed
 512   //
 513
 514   fNUsed = 0;
 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++;
 519   }
 520
 521 }
 522
 523 //_____________________________________________________________________________
 524 Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror)
 525 {
 526   //
 527   // Fit the Rieman tilt
 528   //
 529
 530   // Fitting with tilting pads - kz not fixed
 531   TLinearFitter fitterT2(4,"hyp4");
 532   fitterT2.StoreData(kTRUE);
 533
 534   Float_t xref2 = (cseed[2].fX0 + cseed[3].fX0) * 0.5; // Reference x0 for z
 535
 536   Int_t npointsT = 0;
 537   fitterT2.ClearPoints();
 538
 539   for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
 540
 541     if (!cseed[iLayer].IsOK()) continue;
 542     Double_t tilt = cseed[iLayer].fTilt;
 543
 544     for (Int_t itime = 0; itime < fgTimeBins+1; itime++) {
 545
 546       if (!cseed[iLayer].fUsable[itime]) continue;
 547       // x relative to the midle chamber
 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];
 551
 552       //
 553       // Tilted rieman
 554       //
 555       Double_t uvt[6];
 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       }
 571       fitterT2.AddPoint(uvt,uvt[4],error);
 572       npointsT++;
 573
 574     }
 575
 576   }
 577
 578   fitterT2.Eval();
 579   Double_t rpolz0 = fitterT2.GetParameter(3);
 580   Double_t rpolz1 = fitterT2.GetParameter(4);
 581
 582   //
 583   // Linear fitter  - not possible to make boundaries
 584   // non accept non possible z and dzdx combination
 585   //
 586   Bool_t acceptablez = kTRUE;
 587   for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
 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;
 591   }
 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;
 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
 604   Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT);
 605   Double_t params[3];
 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
 611
 612   for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
 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
 619     // y
 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) {
 624       res2 = TMath::Sqrt(res2);
 625       y    = (1.0 - res2) / params[0];
 626     }
 627
 628     //dy
 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) {
 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;
 636       }
 637     }
 638     z  = rpolz0 + rpolz1 * (x - xref2);
 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;
 644     cseed[iLayer].fC       = curvature;
 645
 646   }
 647
 648   return chi2TR;
 649
 650 }