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
  25 #include "AliTRDseed.h"
  26 #include "AliTRDcluster.h"
  27 #include "AliTRDtracker.h"
  28
  29 #include "TMath.h"
  30 #include "TLinearFitter.h"
  31
  32 ClassImp(AliTRDseed)
  33
  34 //_____________________________________________________________________________
  35 AliTRDseed::AliTRDseed() :
  36   TObject(),
  37   fTilt(0),
  38   fPadLength(0),
  39   fX0(0),
  40   fSigmaY(0),
  41   fSigmaY2(0),
  42   fMeanz(0),
  43   fZProb(0),
  44   fN(0),
  45   fN2(0),
  46   fNUsed(0),
  47   fFreq(0),
  48   fNChange(0),
  49   fMPads(0),
  50   fC(0),
  51   fCC(0),
  52   fChi2(1e10),
  53   fChi2Z(1e10)
  54 {
  55   //
  56   // Default constructor
  57   //
  58
  59   for (Int_t i = 0; i < 25; i++) {
  60     fX[i]        = 0; // x position
  61     fY[i]        = 0; // y position
  62     fZ[i]        = 0; // z position
  63     fIndexes[i]  = 0; // indexes
  64     fClusters[i] = 0x0; // clusters
  65     fUsable[i]   = 0; // indication  - usable cluster
  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 void AliTRDseed::Reset()
  81 {
  82   //
  83   // Reset seed
  84   //
  85
  86   for (Int_t i = 0; i < 25; i++) {
  87     fX[i]        = 0;  // !x position
  88     fY[i]        = 0;  // !y position
  89     fZ[i]        = 0;  // !z position
  90     fIndexes[i]  = 0;  // !indexes
  91     fClusters[i] = 0;  // !clusters
  92     fUsable[i]   = kFALSE;
  93   }
  94   for (Int_t i = 0; i < 2; i++) {
  95     fYref[i]   = 0;    // reference y
  96     fZref[i]   = 0;    // reference z
  97     fYfit[i]   = 0;    // y fit position +derivation
  98     fYfitR[i]  = 0;    // y fit position +derivation
  99     fZfit[i]   = 0;    // z fit position
 100     fZfitR[i]  = 0;    // z fit position
 101     fLabels[i] =-1;    // labels
 102   }
 103   fSigmaY  = 0;        //"robust" sigma in y
 104   fSigmaY2 = 0;        //"robust" sigma in y
 105   fMeanz   = 0;        // mean vaue of z
 106   fZProb   = 0;        // max probbable z
 107   fMPads   = 0;
 108   fN       = 0;        // number of associated clusters
 109   fN2      = 0;        // number of not crossed
 110   fNUsed   = 0;        // number of used clusters
 111   fNChange = 0;        // change z counter
 112
 113 }
 114
 115 //_____________________________________________________________________________
 116 void AliTRDseed::CookLabels()
 117 {
 118   //
 119   // Cook 2 labels for seed
 120   //
 121
 122   Int_t labels[200];
 123   Int_t out[200];
 124   Int_t nlab =0;
 125   for (Int_t i=0;i<25;i++){
 126     if (!fClusters[i]) continue;
 127     for (Int_t ilab=0;ilab<3;ilab++){
 128       if (fClusters[i]->GetLabel(ilab)>=0){
 129         labels[nlab] = fClusters[i]->GetLabel(ilab);
 130         nlab++;
 131       }
 132     }
 133   }
 134   Int_t nlab2 = AliTRDtracker::Freq(nlab,labels,out,kTRUE);
 135   fLabels[0] = out[0];
 136   if (nlab2>1 && out[3]>1) fLabels[1] =out[2];
 137
 138 }
 139
 140 //_____________________________________________________________________________
 141 void AliTRDseed::UseClusters()
 142 {
 143   //
 144   // Use clusters
 145   //
 146
 147   for (Int_t i=0;i<25;i++){
 148     if (!fClusters[i]) continue;
 149     if (!(fClusters[i]->IsUsed())) fClusters[i]->Use();
 150   }
 151
 152 }
 153
 154 //_____________________________________________________________________________
 155 void AliTRDseed::Update()
 156 {
 157   //
 158   // Update the seed.
 159   //
 160
 161   const Float_t kRatio  = 0.8;
 162   const Int_t   kClmin  = 6;
 163   const Float_t kmaxtan = 2;
 164
 165   if (TMath::Abs(fYref[1])>kmaxtan) return;             // too much inclined track
 166
 167   Float_t  sigmaexp = 0.05+TMath::Abs(fYref[1]*0.25);   // expected r.m.s in y direction
 168   Float_t  ycrosscor = fPadLength*fTilt*0.5;             // y correction for crossing
 169   fNChange =0;
 170
 171   Double_t sumw, sumwx,sumwx2;
 172   Double_t sumwy, sumwxy, sumwz,sumwxz;
 173   Int_t    zints[25];        // histograming of the z coordinate
 174                              // get 1 and second max probable coodinates in z
 175   Int_t    zouts[50];        //
 176   Float_t  allowedz[25];     // allowed z for given time bin
 177   Float_t  yres[25];         // residuals from reference
 178   Float_t  anglecor = fTilt*fZref[1];  //correction to the angle
 179
 180
 181   fN=0; fN2 =0;
 182   for (Int_t i=0;i<25;i++){
 183     yres[i] =10000;
 184     if (!fClusters[i]) continue;
 185     yres[i] = fY[i]-fYref[0]-(fYref[1]+anglecor)*fX[i];   // residual y
 186     zints[fN] = Int_t(fZ[i]);
 187     fN++;
 188   }
 189
 190   if (fN<kClmin) return;
 191   Int_t nz = AliTRDtracker::Freq(fN,zints,zouts,kFALSE);
 192   fZProb   = zouts[0];
 193   if (nz<=1) zouts[3]=0;
 194   if (zouts[1]+zouts[3]<kClmin) return;
 195
 196   if (TMath::Abs(zouts[0]-zouts[2])>12.) zouts[3]=0;   // z distance bigger than pad - length
 197
 198   Int_t  breaktime = -1;
 199   Bool_t mbefore   = kFALSE;
 200   Int_t  cumul[25][2];
 201   Int_t  counts[2]={0,0};
 202
 203   if (zouts[3]>=3){
 204     //
 205     // find the break time allowing one chage on pad-rows with maximal numebr of accepted clusters
 206     //
 207     fNChange=1;
 208     for (Int_t i=0;i<25;i++){
 209       cumul[i][0] = counts[0];
 210       cumul[i][1] = counts[1];
 211       if (TMath::Abs(fZ[i]-zouts[0])<2) counts[0]++;
 212       if (TMath::Abs(fZ[i]-zouts[2])<2) counts[1]++;
 213     }
 214     Int_t  maxcount  = 0;
 215     for (Int_t i=0;i<24;i++) {
 216       Int_t after  = cumul[24][0]-cumul[i][0];
 217       Int_t before = cumul[i][1];
 218       if (after+before>maxcount) {
 219         maxcount=after+before;
 220         breaktime=i;
 221         mbefore=kFALSE;
 222       }
 223       after  = cumul[24][1]-cumul[i][1];
 224       before = cumul[i][0];
 225       if (after+before>maxcount) {
 226         maxcount=after+before;
 227         breaktime=i;
 228         mbefore=kTRUE;
 229       }
 230     }
 231     breaktime-=1;
 232   }
 233   for (Int_t i=0;i<25;i++){
 234     if (i>breaktime)  allowedz[i] =   mbefore  ? zouts[2]:zouts[0];
 235     if (i<=breaktime) allowedz[i] = (!mbefore) ? zouts[2]:zouts[0];
 236   }
 237   if ( (allowedz[0]>allowedz[24] && fZref[1]<0) || (allowedz[0]<allowedz[24] &&  fZref[1]>0)){
 238     //
 239     // tracklet z-direction not in correspondance with track z direction
 240     //
 241     fNChange =0;
 242     for (Int_t i=0;i<25;i++){
 243       allowedz[i] =  zouts[0];  //only longest taken
 244     }
 245   }
 246
 247   if (fNChange>0){
 248     //
 249     // cross pad -row tracklet  - take the step change into account
 250     //
 251     for (Int_t i=0;i<25;i++){
 252       if (!fClusters[i]) continue;
 253       if (TMath::Abs(fZ[i]-allowedz[i])>2) continue;
 254       yres[i] = fY[i]-fYref[0]-(fYref[1]+anglecor)*fX[i];   // residual y
 255       if (TMath::Abs(fZ[i]-fZProb)>2){
 256         if (fZ[i]>fZProb) yres[i]+=fTilt*fPadLength;
 257         if (fZ[i]<fZProb) yres[i]-=fTilt*fPadLength;
 258       }
 259     }
 260   }
 261
 262   Double_t yres2[25];
 263   Double_t mean,sigma;
 264   for (Int_t i=0;i<25;i++){
 265     if (!fClusters[i]) continue;
 266     if (TMath::Abs(fZ[i]-allowedz[i])>2) continue;
 267     yres2[fN2] =  yres[i];
 268     fN2++;
 269   }
 270   if (fN2<kClmin){
 271     fN2 = 0;
 272     return;
 273   }
 274   EvaluateUni(fN2,yres2,mean,sigma,Int_t(fN2*kRatio-2));
 275   if (sigma<sigmaexp*0.8) sigma=sigmaexp;
 276   fSigmaY = sigma;
 277   //
 278   //
 279   // reset sums
 280   sumw=0; sumwx=0; sumwx2=0;
 281   sumwy=0; sumwxy=0; sumwz=0;sumwxz=0;
 282   fN2 =0;
 283   fMeanz =0;
 284   fMPads =0;
 285   //
 286   for (Int_t i=0;i<25;i++){
 287     fUsable[i]=kFALSE;
 288     if (!fClusters[i]) continue;
 289     if (TMath::Abs(fZ[i]-allowedz[i])>2)  continue;
 290     if (TMath::Abs(yres[i]-mean)>4.*sigma) continue;
 291     fUsable[i] = kTRUE;
 292     fN2++;
 293     fMPads+=fClusters[i]->GetNPads();
 294     Float_t weight =1;
 295     if (fClusters[i]->GetNPads()>4) weight=0.5;
 296     if (fClusters[i]->GetNPads()>5) weight=0.2;
 297     //
 298     Double_t x = fX[i];
 299     sumw+=weight; sumwx+=x*weight; sumwx2+=x*x*weight;
 300     sumwy+=weight*yres[i];  sumwxy+=weight*(yres[i])*x;
 301     sumwz+=weight*fZ[i];    sumwxz+=weight*fZ[i]*x;
 302   }
 303   if (fN2<kClmin){
 304     fN2 = 0;
 305     return;
 306   }
 307   fMeanz       = sumwz/sumw;
 308   Float_t correction =0;
 309   if (fNChange>0){
 310     // tracklet on boundary
 311     if (fMeanz<fZProb) correction =   ycrosscor;
 312     if (fMeanz>fZProb) correction =  -ycrosscor;
 313   }
 314   Double_t det = sumw*sumwx2-sumwx*sumwx;
 315   fYfitR[0]    = (sumwx2*sumwy-sumwx*sumwxy)/det;
 316   fYfitR[1]    = (sumw*sumwxy-sumwx*sumwy)/det;
 317
 318   fSigmaY2     =0;
 319   for (Int_t i=0;i<25;i++){
 320     if (!fUsable[i]) continue;
 321     Float_t delta = yres[i]-fYfitR[0]-fYfitR[1]*fX[i];
 322     fSigmaY2+=delta*delta;
 323   }
 324   fSigmaY2 = TMath::Sqrt(fSigmaY2/Float_t(fN2-2));
 325   //
 326   fZfitR[0]    = (sumwx2*sumwz-sumwx*sumwxz)/det;
 327   fZfitR[1]    = (sumw*sumwxz-sumwx*sumwz)/det;
 328   fZfit[0]     = (sumwx2*sumwz-sumwx*sumwxz)/det;
 329   fZfit[1]     = (sumw*sumwxz-sumwx*sumwz)/det;
 330   fYfitR[0]   += fYref[0]+correction;
 331   fYfitR[1]   += fYref[1];
 332   fYfit[0]     = fYfitR[0];
 333   fYfit[1]     = fYfitR[1];
 334
 335   UpdateUsed();
 336
 337 }
 338
 339 //_____________________________________________________________________________
 340 void AliTRDseed::UpdateUsed()
 341 {
 342   //
 343   // Update used seed
 344   //
 345
 346   fNUsed = 0;
 347   for (Int_t i=0;i<25;i++){
 348      if (!fClusters[i]) continue;
 349      if ((fClusters[i]->IsUsed())) fNUsed++;
 350   }
 351
 352 }
 353
 354 //_____________________________________________________________________________
 355 void AliTRDseed::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean
 356                            , Double_t &sigma, Int_t hh)
 357 {
 358   //
 359   // Robust estimator in 1D case MI version
 360   //
 361   // For the univariate case
 362   // estimates of location and scatter are returned in mean and sigma parameters
 363   // the algorithm works on the same principle as in multivariate case -
 364   // it finds a subset of size hh with smallest sigma, and then returns mean and
 365   // sigma of this subset
 366   //
 367
 368   if (hh==0)
 369     hh=(nvectors+2)/2;
 370   Double_t faclts[]={2.6477,2.5092,2.3826,2.2662,2.1587,2.0589,1.9660,1.879,1.7973,1.7203,1.6473};
 371   Int_t *index=new Int_t[nvectors];
 372   TMath::Sort(nvectors, data, index, kFALSE);
 373
 374   Int_t    nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
 375   Double_t factor = faclts[nquant-1];
 376
 377   Double_t sumx  =0;
 378   Double_t sumx2 =0;
 379   Int_t    bestindex = -1;
 380   Double_t bestmean  = 0;
 381   Double_t bestsigma = data[index[nvectors-1]]-data[index[0]];   // maximal possible sigma
 382   for (Int_t i=0; i<hh; i++){
 383     sumx  += data[index[i]];
 384     sumx2 += data[index[i]]*data[index[i]];
 385   }
 386
 387   Double_t norm = 1./Double_t(hh);
 388   Double_t norm2 = 1./Double_t(hh-1);
 389   for (Int_t i=hh; i<nvectors; i++){
 390     Double_t cmean  = sumx*norm;
 391     Double_t csigma = (sumx2 - hh*cmean*cmean)*norm2;
 392     if (csigma<bestsigma){
 393       bestmean  = cmean;
 394       bestsigma = csigma;
 395       bestindex = i-hh;
 396     }
 397
 398     sumx  += data[index[i]]-data[index[i-hh]];
 399     sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]];
 400   }
 401
 402   Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
 403   mean  = bestmean;
 404   sigma = bstd;
 405   delete [] index;
 406
 407 }
 408
 409 //_____________________________________________________________________________
 410 Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror)
 411 {
 412   //
 413   // Fit the Rieman tilt
 414   //
 415
 416   TLinearFitter fitterT2(4,"hyp4");  // fitting with tilting pads - kz not fixed
 417   fitterT2.StoreData(kTRUE);
 418   Float_t xref2 = (cseed[2].fX0+cseed[3].fX0)*0.5; // reference x0 for z
 419   //
 420   Int_t npointsT =0;
 421   fitterT2.ClearPoints();
 422   for (Int_t iLayer=0; iLayer<6;iLayer++){
 423     if (!cseed[iLayer].IsOK()) continue;
 424     Double_t tilt = cseed[iLayer].fTilt;
 425
 426     for (Int_t itime=0;itime<25;itime++){
 427       if (!cseed[iLayer].fUsable[itime]) continue;
 428       // x relative to the midle chamber
 429       Double_t x   = cseed[iLayer].fX[itime]+cseed[iLayer].fX0-xref2;
 430       Double_t y   = cseed[iLayer].fY[itime];
 431       Double_t z   = cseed[iLayer].fZ[itime];
 432
 433       //
 434       // tilted rieman
 435       //
 436       Double_t uvt[6];
 437       Double_t x2 = cseed[iLayer].fX[itime]+cseed[iLayer].fX0;      // global x
 438       Double_t t = 1./(x2*x2+y*y);
 439       uvt[1]  = t;    // t
 440       uvt[0]  = 2.*x2*uvt[1];      // u
 441       uvt[2]  = 2.0*tilt*uvt[1];
 442       uvt[3]  = 2.0*tilt*x*uvt[1];
 443       uvt[4]  = 2.0*(y+tilt*z)*uvt[1];
 444       //
 445       Double_t error = 2*uvt[1];
 446       if (terror) error*=cseed[iLayer].fSigmaY;
 447       else {error *=0.2;} //default error
 448       fitterT2.AddPoint(uvt,uvt[4],error);
 449       npointsT++;
 450     }
 451   }
 452   fitterT2.Eval();
 453   Double_t rpolz0 = fitterT2.GetParameter(3);
 454   Double_t rpolz1 = fitterT2.GetParameter(4);
 455   //
 456   // linear fitter  - not possible to make boundaries
 457   // non accept non possible z and dzdx combination
 458   //
 459   Bool_t   acceptablez =kTRUE;
 460   for (Int_t iLayer=0; iLayer<6;iLayer++){
 461     if (cseed[iLayer].IsOK()){
 462       Double_t zT2 =  rpolz0+rpolz1*(cseed[iLayer].fX0 - xref2);
 463       if (TMath::Abs(cseed[iLayer].fZProb-zT2)>cseed[iLayer].fPadLength*0.5+1)
 464         acceptablez = kFALSE;
 465     }
 466   }
 467   if (!acceptablez){
 468     Double_t zmf  = cseed[2].fZref[0]+cseed[2].fZref[1]*(xref2-cseed[2].fX0);
 469     Double_t dzmf = (cseed[2].fZref[1]+ cseed[3].fZref[1])*0.5;
 470     fitterT2.FixParameter(3,zmf);
 471     fitterT2.FixParameter(4,dzmf);
 472     fitterT2.Eval();
 473     fitterT2.ReleaseParameter(3);
 474     fitterT2.ReleaseParameter(4);
 475     rpolz0 = fitterT2.GetParameter(3);
 476     rpolz1 = fitterT2.GetParameter(4);
 477   }
 478
 479   Double_t chi2TR = fitterT2.GetChisquare()/Float_t(npointsT);
 480   Double_t params[3];
 481   params[0]     =  fitterT2.GetParameter(0);
 482   params[1]     =  fitterT2.GetParameter(1);
 483   params[2]     =  fitterT2.GetParameter(2);
 484   Double_t curvature     =  1+params[1]*params[1]-params[2]*params[0];
 485   for (Int_t iLayer = 0; iLayer<6;iLayer++){
 486     Double_t  x = cseed[iLayer].fX0;
 487     Double_t  y=0,dy=0, z=0, dz=0;
 488     // y
 489     Double_t res2 = (x*params[0]+params[1]);
 490     res2*=res2;
 491     res2 = 1.-params[2]*params[0]+params[1]*params[1]-res2;
 492     if (res2>=0){
 493       res2 = TMath::Sqrt(res2);
 494       y    = (1-res2)/params[0];
 495     }
 496     //dy
 497     Double_t x0 = -params[1]/params[0];
 498     if (-params[2]*params[0]+params[1]*params[1]+1>0){
 499       Double_t rm1  = params[0]/TMath::Sqrt(-params[2]*params[0]+params[1]*params[1]+1);
 500       if ( 1./(rm1*rm1)-(x-x0)*(x-x0)>0){
 501         Double_t res = (x-x0)/TMath::Sqrt(1./(rm1*rm1)-(x-x0)*(x-x0));
 502         if (params[0]<0) res*=-1.;
 503         dy = res;
 504       }
 505     }
 506     z  = rpolz0+rpolz1*(x-xref2);
 507     dz = rpolz1;
 508     cseed[iLayer].fYref[0] = y;
 509     cseed[iLayer].fYref[1] = dy;
 510     cseed[iLayer].fZref[0] = z;
 511     cseed[iLayer].fZref[1] = dz;
 512     cseed[iLayer].fC  = curvature;
 513
 514   }
 515
 516   return chi2TR;
 517
 518 }