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