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