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 "AliTRDseed.h"
  28 #include "AliTRDcluster.h"
  29 #include "AliTRDtracker.h"
  30
  31 ClassImp(AliTRDseed)
  32
  33 //_____________________________________________________________________________
  34 AliTRDseed::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)
  53 {
  54   //
  55   // Default constructor
  56   //
  57
  58   for (Int_t i = 0; i < 25; i++) {
  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 //_____________________________________________________________________________
  80 AliTRDseed::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
 111   }
 112
 113   for (Int_t i = 0; i < 2; i++) {
 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
 121   }
 122
 123 }
 124
 125 //_____________________________________________________________________________
 126 void AliTRDseed::Reset()
 127 {
 128   //
 129   // Reset seed
 130   //
 131
 132   for (Int_t i = 0; i < 25; i++) {
 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
 138     fUsable[i]   = kFALSE;
 139   }
 140
 141   for (Int_t i = 0; i < 2; i++) {
 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
 149   }
 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
 154   fMPads   = 0;
 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
 159
 160 }
 161
 162 //_____________________________________________________________________________
 163 void AliTRDseed::CookLabels()
 164 {
 165   //
 166   // Cook 2 labels for seed
 167   //
 168
 169   Int_t labels[200];
 170   Int_t out[200];
 171   Int_t nlab = 0;
 172
 173   for (Int_t i = 0; i < 25; i++) {
 174     if (!fClusters[i]) continue;
 175     for (Int_t ilab = 0; ilab < 3; ilab++) {
 176       if (fClusters[i]->GetLabel(ilab) >= 0) {
 177         labels[nlab] = fClusters[i]->GetLabel(ilab);
 178         nlab++;
 179       }
 180     }
 181   }
 182
 183   Int_t nlab2 = AliTRDtracker::Freq(nlab,labels,out,kTRUE);
 184   fLabels[0] = out[0];
 185   if ((nlab2  > 1) &&
 186       (out[3] > 1)) {
 187     fLabels[1] = out[2];
 188   }
 189
 190 }
 191
 192 //_____________________________________________________________________________
 193 void AliTRDseed::UseClusters()
 194 {
 195   //
 196   // Use clusters
 197   //
 198
 199   for (Int_t i = 0; i < 25; i++) {
 200     if (!fClusters[i]) continue;
 201     if (!(fClusters[i]->IsUsed())) fClusters[i]->Use();
 202   }
 203
 204 }
 205
 206 //_____________________________________________________________________________
 207 void 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
 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
 237
 238
 239   fN  = 0;
 240   fN2 = 0;
 241   for (Int_t i = 0; i < 25; i++) {
 242     yres[i] = 10000.0;
 243     if (!fClusters[i]) continue;
 244     yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i];   // Residual y
 245     zints[fN] = Int_t(fZ[i]);
 246     fN++;
 247   }
 248
 249   if (fN < kClmin) return;
 250   Int_t nz = AliTRDtracker::Freq(fN,zints,zouts,kFALSE);
 251   fZProb   = zouts[0];
 252   if (nz <= 1) zouts[3] = 0;
 253   if (zouts[1] + zouts[3] < kClmin) return;
 254
 255   // Z distance bigger than pad - length
 256   if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) {
 257     zouts[3]=0;
 258   }
 259
 260   Int_t  breaktime = -1;
 261   Bool_t mbefore   = kFALSE;
 262   Int_t  cumul[25][2];
 263   Int_t  counts[2] = { 0, 0 };
 264
 265   if (zouts[3] >= 3) {
 266
 267     //
 268     // Find the break time allowing one chage on pad-rows
 269     // with maximal numebr of accepted clusters
 270     //
 271     fNChange = 1;
 272     for (Int_t i = 0; i < 25; i++) {
 273       cumul[i][0] = counts[0];
 274       cumul[i][1] = counts[1];
 275       if (TMath::Abs(fZ[i]-zouts[0]) < 2) counts[0]++;
 276       if (TMath::Abs(fZ[i]-zouts[2]) < 2) counts[1]++;
 277     }
 278     Int_t  maxcount = 0;
 279     for (Int_t i = 0; i < 24; i++) {
 280       Int_t after  = cumul[24][0] - cumul[i][0];
 281       Int_t before = cumul[i][1];
 282       if (after + before > maxcount) {
 283         maxcount  = after + before;
 284         breaktime = i;
 285         mbefore   = kFALSE;
 286       }
 287       after  = cumul[24][1] - cumul[i][1];
 288       before = cumul[i][0];
 289       if (after + before > maxcount) {
 290         maxcount  = after + before;
 291         breaktime = i;
 292         mbefore   = kTRUE;
 293       }
 294     }
 295
 296     breaktime -= 1;
 297
 298   }
 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];
 303   }
 304
 305   if (((allowedz[0] > allowedz[24]) && (fZref[1] < 0)) ||
 306       ((allowedz[0] < allowedz[24]) && (fZref[1] > 0))) {
 307     //
 308     // Tracklet z-direction not in correspondance with track z direction
 309     //
 310     fNChange = 0;
 311     for (Int_t i = 0; i < 25; i++) {
 312       allowedz[i] = zouts[0];  // Only longest taken
 313     }
 314   }
 315
 316   if (fNChange > 0) {
 317     //
 318     // Cross pad -row tracklet  - take the step change into account
 319     //
 320     for (Int_t i = 0; i < 25; i++) {
 321       if (!fClusters[i]) continue;
 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;
 327       }
 328     }
 329   }
 330
 331   Double_t yres2[25];
 332   Double_t mean;
 333   Double_t sigma;
 334   for (Int_t i = 0; i < 25; i++) {
 335     if (!fClusters[i]) continue;
 336     if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
 337     yres2[fN2] = yres[i];
 338     fN2++;
 339   }
 340   if (fN2 < kClmin) {
 341     fN2 = 0;
 342     return;
 343   }
 344   EvaluateUni(fN2,yres2,mean,sigma,Int_t(fN2*kRatio-2));
 345   if (sigma < sigmaexp * 0.8) {
 346     sigma = sigmaexp;
 347   }
 348   fSigmaY = sigma;
 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;
 366     if (!fClusters[i]) continue;
 367     if (TMath::Abs(fZ[i] - allowedz[i]) > 2)  continue;
 368     if (TMath::Abs(yres[i] - mean) > 4.0 * sigma) continue;
 369     fUsable[i] = kTRUE;
 370     fN2++;
 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
 376     Double_t x = fX[i];
 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
 385   }
 386
 387   if (fN2 < kClmin){
 388     fN2 = 0;
 389     return;
 390   }
 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;
 397   }
 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;
 402
 403   fSigmaY2 = 0;
 404   for (Int_t i = 0; i < 25; i++) {
 405     if (!fUsable[i]) continue;
 406     Float_t delta = yres[i] - fYfitR[0] - fYfitR[1] * fX[i];
 407     fSigmaY2 += delta*delta;
 408   }
 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];
 419
 420   UpdateUsed();
 421
 422 }
 423
 424 //_____________________________________________________________________________
 425 void AliTRDseed::UpdateUsed()
 426 {
 427   //
 428   // Update used seed
 429   //
 430
 431   fNUsed = 0;
 432   for (Int_t i = 0; i < 25; i++) {
 433     if (!fClusters[i]) continue;
 434     if ((fClusters[i]->IsUsed())) fNUsed++;
 435   }
 436
 437 }
 438
 439 //_____________________________________________________________________________
 440 void 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
 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];
 461   TMath::Sort(nvectors, data, index, kFALSE);
 462
 463   Int_t    nquant = TMath::Min(Int_t(Double_t(((hh * 1.0 / nvectors) - 0.5) * 40)) + 1,11);
 464   Double_t factor = faclts[nquant-1];
 465
 466   Double_t sumx      = 0.0;
 467   Double_t sumx2     = 0.0;
 468   Int_t    bestindex = -1;
 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++) {
 472     sumx  += data[index[i]];
 473     sumx2 += data[index[i]]*data[index[i]];
 474   }
 475
 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
 480     Double_t cmean  = sumx*norm;
 481     Double_t csigma = (sumx2 - hh*cmean*cmean) * norm2;
 482     if (csigma < bestsigma) {
 483       bestmean  = cmean;
 484       bestsigma = csigma;
 485       bestindex = i - hh;
 486     }
 487
 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
 491   }
 492
 493   Double_t bstd = factor * TMath::Sqrt(TMath::Abs(bestsigma));
 494   mean  = bestmean;
 495   sigma = bstd;
 496
 497   delete [] index;
 498
 499 }
 500
 501 //_____________________________________________________________________________
 502 Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror)
 503 {
 504   //
 505   // Fit the Rieman tilt
 506   //
 507
 508   // Fitting with tilting pads - kz not fixed
 509   TLinearFitter fitterT2(4,"hyp4");
 510   fitterT2.StoreData(kTRUE);
 511   Float_t xref2 = (cseed[2].fX0 + cseed[3].fX0) * 0.5; // Reference x0 for z
 512
 513   Int_t npointsT = 0;
 514   fitterT2.ClearPoints();
 515
 516   for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
 517
 518     if (!cseed[iLayer].IsOK()) continue;
 519     Double_t tilt = cseed[iLayer].fTilt;
 520
 521     for (Int_t itime = 0; itime < 25; itime++) {
 522
 523       if (!cseed[iLayer].fUsable[itime]) continue;
 524       // x relative to the midle chamber
 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];
 528
 529       //
 530       // Tilted rieman
 531       //
 532       Double_t uvt[6];
 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       }
 548       fitterT2.AddPoint(uvt,uvt[4],error);
 549       npointsT++;
 550
 551     }
 552
 553   }
 554
 555   fitterT2.Eval();
 556   Double_t rpolz0 = fitterT2.GetParameter(3);
 557   Double_t rpolz1 = fitterT2.GetParameter(4);
 558
 559   //
 560   // Linear fitter  - not possible to make boundaries
 561   // non accept non possible z and dzdx combination
 562   //
 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) {
 568         acceptablez = kFALSE;
 569       }
 570     }
 571   }
 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;
 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
 584   Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT);
 585   Double_t params[3];
 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
 599     // y
 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) {
 604       res2 = TMath::Sqrt(res2);
 605       y    = (1.0 - res2) / params[0];
 606     }
 607
 608     //dy
 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;
 615         dy = res;
 616       }
 617     }
 618     z  = rpolz0 + rpolz1 * (x - xref2);
 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;
 624     cseed[iLayer].fC       = curvature;
 625
 626   }
 627
 628   return chi2TR;
 629
 630 }