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