TPC/AliTPCClusterFinder.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 /*
  17 $Log$
  18 Revision 1.7  2002/02/05 09:12:26  hristov
  19 Small mods for gcc 3.02
  20
  21 Revision 1.6  2001/10/21 19:07:24  hristov
  22 Several pointers were set to zero in the default constructors to avoid memory management problems
  23
  24 Revision 1.5  2001/01/26 19:57:22  hristov
  25 Major upgrade of AliRoot code
  26
  27 Revision 1.4  2000/10/05 16:08:15  kowal2
  28 Changes due to a new class AliComplexCluster. Forward declarations.
  29
  30 Revision 1.3  2000/07/10 20:57:39  hristov
  31 Update of TPC code and macros by M.Kowalski
  32
  33 Revision 1.2  2000/06/30 12:07:49  kowal2
  34 Updated from the TPC-PreRelease branch
  35
  36 Revision 1.1.2.1  2000/06/25 08:52:51  kowal2
  37 replacing AliClusterFinder
  38
  39 */
  40
  41 //-----------------------------------------------------------------------------
  42 //
  43 //  Implementation of class ALITPCCLUSTERFINDER
  44 //
  45 //Class for cluster finding in two dimension.
  46 //In the present there are implemented two algorithm
  47 //primitive recursion algorithm. (FindPeaks)
  48 //Algorithm is not working in case of overlaping clusters
  49 //Maximum - minimum in direction algoritm (Find clusters)
  50 //In this algoritm we suppose that each cluster has local
  51 //maximum. From this local maximum I mus see each point
  52 //of cluster.
  53 //From maximum i can accept every point in radial
  54 //direction which is before border in direction
  55 //Border in direction occur if we have next in
  56 //direction nder threshold or response begin
  57 //to increase in given radial direction
  58 //-----------------------------------------------------------------------------
  59
  60 #include "TMinuit.h"
  61 #include "AliArrayI.h"
  62 #include "TClonesArray.h"
  63 #include "AliTPC.h"
  64 #include "TRandom.h"
  65 #include "AliH2F.h"
  66 #include "TMarker.h"
  67 #include "AliComplexCluster.h"
  68 #include "AliTPCClusterFinder.h"
  69 #include <Riostream.h>
  70 #include <Riostream.h>
  71
  72 //direction constants possible direction in 8 different sectors
  73 //
  74
  75
  76 const Int_t kClStackSize =1000;
  77
  78
  79
  80
  81 static AliTPCClusterFinder * gClusterFinder; //for fitting routine
  82
  83 void gauss(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
  84 {
  85   AliArrayI * points = gClusterFinder->GetStack();
  86   const Int_t nbins = gClusterFinder->GetStackIndex();
  87   Int_t i;
  88   //calculate chisquare
  89    Double_t chisq = 0;
  90    Double_t delta;
  91    for (i=0;i<nbins; i++) {
  92      Float_t x = points->At(i*3);
  93      Float_t y = points->At(i*3+1);
  94      Float_t z = points->At(i*3+2);
  95      Float_t deltax2 = (x-par[1]);
  96      deltax2*=deltax2;
  97      deltax2*=par[3];
  98      Float_t deltay2 = (y-par[2]);
  99      deltay2*=deltay2;
 100      deltay2*=par[4];
 101
 102      delta  = z-par[0]*TMath::Exp(-deltax2-deltay2);
 103      chisq += delta*delta;
 104    }
 105    f = chisq;
 106 }
 107
 108
 109 ClassImp(AliTPCClusterFinder)
 110   //ClassImp(AliCell)
 111
 112 AliTPCClusterFinder::AliTPCClusterFinder()
 113 {
 114   fDigits =0;
 115   fCells = 0;
 116   fDimX = 0;
 117   fDimY = 0;
 118   fNoiseTh = 3;
 119   fMulSigma2 = 16; //4 sigma
 120   fDirSigmaFac = 1.4;
 121   fDirAmpFac =1.3;
 122   fNType=8;
 123   fThreshold = 2;
 124   fStack = new AliArrayI;
 125   fStack->Set(kClStackSize);
 126   fClustersArray =0;
 127   SetSigmaX(1,0,0);
 128   SetSigmaY(1,0,0);
 129
 130
 131   fDetectorParam = 0;
 132   fDetectorIndex = 0;
 133   ResetStatus();
 134   fBFit = kFALSE;
 135   fMinuit= new TMinuit(5);
 136   fMinuit->SetFCN(gauss);
 137   gClusterFinder = this;
 138
 139 }
 140
 141
 142 AliTPCClusterFinder::~AliTPCClusterFinder()
 143 {
 144  if (fDigits  != 0) delete fDigits;
 145 }
 146
 147 void AliTPCClusterFinder::SetSigmaX(Float_t s0, Float_t s1x, Float_t s1y)
 148 {
 149   fSigmaX[0]=s0;
 150   fSigmaX[1]=s1x;
 151   fSigmaX[2]=s1y;
 152
 153 }
 154 void AliTPCClusterFinder::SetSigmaY(Float_t s0, Float_t s1x, Float_t s1y)
 155 {
 156   fSigmaY[0]=s0;
 157   fSigmaY[1]=s1x;
 158   fSigmaY[2]=s1y;
 159 }
 160
 161
 162
 163 Bool_t AliTPCClusterFinder::SetSigma2(Int_t i, Int_t j, Float_t & sigmax2, Float_t &sigmay2)
 164 {
 165   //
 166   //set sigmax2 and sigma y2  accordig i and j position of cell
 167   //
 168
 169   //  Float_t x[3] = {ItoX(i),JtoY(j),0};
 170   Float_t x= ItoX(i);
 171   Float_t y= JtoY(j);
 172
 173   sigmax2= fSigmaX[0]+fSigmaX[1]*x+fSigmaX[2]*y;
 174   sigmay2= fSigmaY[0]+fSigmaY[1]*x+fSigmaY[2]*y;
 175   return kTRUE;
 176 }
 177
 178 /*
 179 Bool_t AliTPCClusterFinder::SetSigma2(Int_t i, Int_t j, Float_t & sigmax2, Float_t &sigmay2)
 180 {
 181   //
 182   //set sigmax2 and sigma y2  accordig i and j position of cell
 183   //
 184   if (fDetectorParam==0) {
 185     sigmax2=1;
 186     sigmay2=1;
 187     return kFALSE;
 188   }
 189   Float_t x[3] = {ItoX(i),JtoY(j),0};
 190   Float_t sigma[2];
 191   fDetectorParam->GetClusterSize(x,fDetectorIndex,0,0,sigma);
 192   sigmax2=sigma[0]*(fX2-fX1)*(fX2-fX1)/(fDimX*fDimX);
 193   sigmay2=sigma[1]*(fY2-fY1)*(fY2-fY1)/(fDimY*fDimY);
 194   return kTRUE;
 195 }
 196 */
 197
 198
 199 void AliTPCClusterFinder::GetHisto(TH2F * his2)
 200 {
 201
 202   UInt_t idim =his2->GetNbinsX();
 203   UInt_t jdim =his2->GetNbinsY();
 204   fX1 = his2->GetXaxis()->GetXmin();
 205   fX2 = his2->GetXaxis()->GetXmax();
 206   fY1 = his2->GetYaxis()->GetXmin();
 207   fY2 = his2->GetYaxis()->GetXmax();
 208
 209   if ( (idim>0) && (jdim>0))
 210     {
 211       rOK = kTRUE;
 212       fDimX = idim;
 213       fDimY = jdim;
 214       Int_t size =idim*jdim;
 215       if (fDigits !=0) delete fDigits;
 216       fDigits  = (Int_t*) new Int_t[size];
 217       fCells  = (AliCell*) new AliCell[size];
 218
 219     }  else
 220       rOK=kFALSE;
 221   for (Int_t i = 0; i<(Int_t)idim;i++)
 222     for (Int_t j = 0; j<(Int_t)jdim;j++)
 223       {
 224         Int_t index = his2->GetBin(i+1,j+1);
 225         //AliCell * cell = GetCell(i,j);
 226         //if (cell!=0) cell->SetSignal(his2->GetBinContent(index));
 227     SetSignal(static_cast<int>(his2->GetBinContent(index)),i,j);
 228       }
 229
 230 }
 231
 232
 233
 234
 235 void AliTPCClusterFinder::FindMaxima()
 236 {
 237   for (Int_t i=0; i<fDimX; i++)
 238     for (Int_t j=0;j<fDimY; j++)
 239       if (IsMaximum(i,j)) cout<<i<<"   "<<j<<"\n";
 240 }
 241
 242
 243 void  AliTPCClusterFinder::Transform(AliDigitCluster * c)
 244 {
 245   //transform coordinata from bin coordinata to "normal coordinata"
 246   //for example if we initialize finder with histogram
 247   //it transform values from bin coordinata to the histogram coordinata
 248   c->fX=ItoX(c->fX);
 249   c->fY=JtoY(c->fY);
 250   c->fMaxX=ItoX(c->fMaxX);
 251   c->fMaxY=JtoY(c->fMaxY);
 252
 253   c->fSigmaX2=c->fSigmaX2*(fX2-fX1)*(fX2-fX1)/(fDimX*fDimX);
 254   c->fSigmaY2=c->fSigmaY2*(fY2-fY1)*(fY2-fY1)/(fDimY*fDimY);
 255   c->fArea   =c->fArea*(fX2-fX1)*(fY2-fY1)/(fDimX*fDimY);
 256 }
 257 void  AliTPCClusterFinder::AddToStack(Int_t i, Int_t j, Int_t signal)
 258 {
 259   //
 260   //add digit to stack
 261   //
 262   if ( ((fStackIndex+2)>=kClStackSize) || (fStackIndex<0) ) return;
 263   fStack->AddAt(i,fStackIndex);
 264   fStack->AddAt(j,fStackIndex+1);
 265   fStack->AddAt(signal,fStackIndex+2);
 266   fStackIndex+=3;
 267 }
 268
 269 void AliTPCClusterFinder::GetClusterStatistic(AliDigitCluster & cluster)
 270 {
 271   //
 272   //calculate statistic of cluster
 273   //
 274   Double_t sumxw,sumyw,sumx2w,sumy2w,sumxyw,sumw;
 275   Int_t minx,maxx,miny,maxy;
 276   sumxw=sumyw=sumx2w=sumy2w=sumxyw=sumw=0;
 277   minx=fDimX;
 278   maxx=-fDimX;
 279   miny=fDimY;
 280   maxy=-fDimY;
 281   Int_t x0=fStack->At(0);
 282   Int_t y0=fStack->At(1);
 283   Int_t maxQx =x0;
 284   Int_t maxQy =y0;
 285   Int_t maxQ=fStack->At(2);
 286
 287
 288   for (Int_t i = 0; i<fStackIndex;i+=3){
 289     Int_t x = fStack->At(i);
 290     Int_t y = fStack->At(i+1);
 291     Int_t dx=x-x0;
 292     Int_t dy=y-y0;
 293     Int_t w = fStack->At(i+2);
 294     if (w>maxQ){
 295       maxQ = w;
 296       maxQx = x;
 297       maxQy=y;
 298     }
 299     if (x<minx) minx=x;
 300     if (y<miny) miny=y;
 301     if (x>maxx) maxx=x;
 302     if (y>maxy) maxy=y;
 303     sumxw+=dx*w;
 304     sumyw+=dy*w;
 305     sumx2w+=dx*dx*w;
 306     sumy2w+=dy*dy*w;
 307     sumxyw+=dx*dy*w;
 308     sumw+=w;
 309   }
 310   cluster.fQ = sumw;
 311   if (sumw>0){
 312     cluster.fX = sumxw/sumw;
 313     cluster.fY = sumyw/sumw;
 314     cluster.fQ = sumw;
 315     cluster.fSigmaX2 = sumx2w/sumw-cluster.fX*cluster.fX;
 316     cluster.fSigmaY2 = sumy2w/sumw-cluster.fY*cluster.fY;
 317     cluster.fSigmaXY = sumxyw/sumw-cluster.fX*cluster.fY;
 318     cluster.fMaxX = maxQx;
 319     cluster.fMaxY = maxQy;
 320     cluster.fMax = maxQ;
 321     cluster.fArea = fStackIndex/3;
 322     cluster.fNx = maxx-minx+1;
 323     cluster.fNy = maxy-miny+1;
 324     cluster.fX +=x0;
 325     cluster.fY +=y0;
 326   }
 327 }
 328 void AliTPCClusterFinder::GetClusterFit(AliDigitCluster & cluster)
 329 {
 330   //
 331   //calculate statistic of cluster
 332   //
 333   Double_t arglist[10];
 334   Int_t ierflg = 0;
 335
 336   arglist[0] = 1;
 337   fMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
 338
 339   //fistly find starting parameters
 340   Int_t minx,maxx,miny,maxy,maxQ,maxQx,maxQy;
 341   Int_t over =0;
 342   Float_t sumxw,sumyw,sumw;
 343   sumxw=sumyw=sumw=0;
 344   minx=fDimX;
 345   maxx=-fDimX;
 346   miny=fDimY;
 347   maxy=-fDimY;
 348   maxQx=fStack->At(0);
 349   maxQy=fStack->At(1);
 350   maxQ=fStack->At(2);
 351
 352   for (Int_t i = 0; i<fStackIndex;i+=3){
 353     Int_t x = fStack->At(i);
 354     Int_t y = fStack->At(i+1);
 355     Int_t w = fStack->At(i+2);
 356     if (w>fThreshold) {
 357       over++;
 358       sumw+=w;
 359       sumxw+=x*w;
 360       sumyw+=y*w;
 361       if (x<minx) minx=x;
 362       if (y<miny) miny=y;
 363       if (x>maxx) maxx=x;
 364       if (y>maxy) maxy=y;
 365       if (w>maxQ) {
 366         maxQ=w;
 367         maxQx=x;
 368         maxQy=y;
 369       }
 370     }
 371   }
 372   Int_t nx = maxx-minx+1;
 373   Int_t ny = maxy-miny+1;
 374
 375   SetSigma2(maxQx,maxQy,fCurrentSigmaX2,fCurrentSigmaY2);
 376   Double_t vstart[5]={maxQ,sumxw/sumw,sumyw/sumw,1/(2*fCurrentSigmaX2),1/(2*fCurrentSigmaY2)};
 377   Double_t step[5]={1.,0.01,0.01,0.01,0.01};
 378   fMinuit->mnparm(0, "amp", vstart[0], step[0], 0,0,ierflg);
 379   fMinuit->mnparm(1, "x0", vstart[1], step[1], 0,0,ierflg);
 380   fMinuit->mnparm(2, "y0", vstart[2], step[2], 0,0,ierflg);
 381   fMinuit->mnparm(3, "sx2", vstart[3], step[3], 0,0,ierflg);
 382   fMinuit->mnparm(4, "sy2", vstart[4], step[4], 0,0,ierflg);
 383   arglist[0] = 500;
 384   arglist[1] = 1.;
 385
 386   fMinuit->mnfree(0);  //set unfixed all parameters
 387   //if we have area less then
 388   if (over<=21) {  //if we dont't have more then  7  points
 389     fMinuit->FixParameter(3);
 390     fMinuit->FixParameter(4);
 391   }
 392   else {
 393     if (nx<3)  fMinuit->FixParameter(3); //fix sigma x if no data in x direction
 394     if (ny<3)  fMinuit->FixParameter(4);  //fix sigma y if no data in y direction
 395   }
 396   fMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);
 397
 398   if (sumw>0){
 399     Double_t  x[5];
 400     Double_t  error[5];
 401     fMinuit->GetParameter(0,x[0],error[0]);
 402     fMinuit->GetParameter(1,x[1],error[1]);
 403     fMinuit->GetParameter(2,x[2],error[2]);
 404     fMinuit->GetParameter(3,x[3],error[3]);
 405     fMinuit->GetParameter(4,x[4],error[4]);
 406
 407     cluster.fX = x[1];
 408     cluster.fY = x[2];
 409     cluster.fMaxX = maxQx;
 410     cluster.fMaxY = maxQy;
 411
 412     cluster.fQ = sumw;
 413     cluster.fSigmaX2 = 1/TMath::Sqrt(2*x[3]);
 414     cluster.fSigmaY2 = 1/TMath::Sqrt(2*x[4]);
 415     cluster.fSigmaXY = 0;
 416     cluster.fMax = x[0];
 417     cluster.fArea = over;
 418     cluster.fNx = nx;
 419     cluster.fNy = ny;
 420   }
 421 }
 422
 423 Bool_t   AliTPCClusterFinder::CheckIfDirBorder(Float_t x, Float_t y,
 424                                           Int_t i,Int_t j)
 425 {
 426   //
 427   //function which control if given cell with index i, j is the
 428   //minimum in direction
 429   // x and y are estimate of local maximum
 430   //direction is given by the
 431   Float_t virtualcell;
 432   AliCell * cellor= GetCell(i,j);
 433   Int_t     sigor = GetSignal(i,j);
 434
 435   //control derivation in direction
 436   //if function grows up in direction then there is border
 437   Float_t dx = i-x;
 438   Float_t dy = j-y;
 439   Float_t dd = TMath::Sqrt(dx*dx+dy*dy);
 440   Float_t ddx = TMath::Abs(dx);
 441   ddx =   (ddx>0.5) ? ddx-0.5: 0;
 442   ddx*=ddx;
 443   Float_t ddy = TMath::Abs(dy);
 444   ddy = (ddy>0.5) ? ddy-0.5: 0;
 445   ddy*=ddy;
 446   Float_t d2 = ddx/(2*fDirSigmaFac*fCurrentSigmaX2)+ddy/(2*fDirSigmaFac*fCurrentSigmaY2); //safety factor
 447   //I accept sigmax and sigma y bigge by factor sqrt(fDirsigmaFac)
 448   Float_t amp = TMath::Exp(-d2)*fCurrentMaxAmp*fDirAmpFac; //safety factor fDirFac>1
 449
 450   if (sigor>amp) return kTRUE;
 451   if (dd==0) return kFALSE;
 452
 453   dx/=dd;
 454   dy/=dd;
 455   virtualcell = GetVirtualSignal(i+dx,j+dy);
 456   if (virtualcell <=fThreshold) return kFALSE;
 457   if (virtualcell>sigor)
 458     if (virtualcell>(sigor+fNoiseTh))
 459       {cellor->SetDirBorder(fIndex); return kTRUE;}
 460     else
 461       {
 462         virtualcell = GetVirtualSignal(i+2*dx,j+2*dy);
 463         if (virtualcell>sigor)
 464           { cellor->SetDirBorder(fIndex); return kTRUE;}
 465       };
 466   return kFALSE;
 467 }
 468
 469
 470
 471
 472
 473 Bool_t  AliTPCClusterFinder::IsMaximum(Int_t i, Int_t  j)
 474 {
 475   //there is maximum if given digits is 1 sigma over all adjacent
 476   //in 8 neighborow
 477   //or ther exist virual maximum
 478   //is maximum on 24 points neighboring
 479   //  Bool_t res = kFALSE;
 480   Int_t over =0;
 481   Int_t overth=0;
 482   Int_t oversigma =0;
 483   AliCell * cell = GetCell(i,j);
 484   Int_t signal = GetSignal(i,j);
 485   if (cell == 0) return kFALSE;
 486   for ( Int_t di=-1;di<=1;di++)
 487     for ( Int_t dj=-1;dj<=1;dj++){
 488       if ( (di!=0) || (dj!=0))
 489         {
 490           AliCell * cell2=GetCell(i+di,j+dj);
 491           Int_t signal2 = GetSignal(i+di,j+dj);
 492           if (cell2 == 0) {
 493             over+=1;
 494             oversigma+=1;
 495           }
 496           else
 497             {
 498               if (signal2>signal) return kFALSE;
 499               if (signal2>fThreshold) overth++;
 500               if (signal2==signal) {
 501                 if (di<0) return kFALSE;
 502                 if ( (di+dj)<0) return kFALSE;
 503               }
 504               //              if (signal>=signal2){
 505               over+=1;
 506               if (signal>fNoiseTh+signal2)
 507                 oversigma+=1;
 508               //}
 509             }
 510         }
 511     }
 512   //if I have only one neighborough over threshold
 513   if (overth<2) return kFALSE;
 514   if (over<8)   return kFALSE;
 515   if (oversigma==8) {
 516     fCurrentMaxX = i;
 517     fCurrentMaxY = j;
 518     fCurrentMaxAmp =signal;
 519     SetMaximum(fIndex,i,j);
 520     return kTRUE;
 521   }
 522   //check if there exist virtual maximum
 523   for (Float_t ddi=0.;(ddi<1.);ddi+=0.5)
 524     for (Float_t ddj=0.;(ddj<1.);ddj+=0.5)
 525       if (IsVirtualMaximum(Float_t(i)+ddi,Float_t(j)+ddj)){
 526         fCurrentMaxX = i+ddi;
 527         fCurrentMaxY = j+ddj;
 528         fCurrentMaxAmp =signal;
 529         SetMaximum(fIndex,i,j);
 530         return kTRUE;
 531       }
 532   return kFALSE;
 533 }
 534
 535 Bool_t  AliTPCClusterFinder::IsVirtualMaximum(Float_t x, Float_t  y)
 536 {
 537   //there is maximum if given digits is 1 sigma over all adjacent
 538   //in 8 neighborow or
 539   //is maximum on 24 points neighboring
 540   Bool_t res = kFALSE;
 541   Int_t over =0;
 542   Int_t overth=0;
 543   Int_t oversigma =0;
 544   Float_t virtualcell = GetVirtualSignal(x,y);
 545   if (virtualcell < 0) return kFALSE;
 546   for ( Int_t di=-1;di<=1;di++)
 547     for ( Int_t dj=-1;dj<=1;dj++)
 548       if ( (di!=0) || (dj!=0))
 549         {
 550           Float_t virtualcell2=GetVirtualSignal(x+di,y+dj);
 551           if (virtualcell2 < 0) {
 552             over+=1;
 553             oversigma+=1;
 554           }
 555           else
 556             {
 557               if (virtualcell2>fThreshold) overth++;
 558               if (virtualcell>=virtualcell2){
 559                 over+=1;
 560                 if (virtualcell>fNoiseTh+virtualcell2)
 561                   oversigma+=1;
 562               }
 563             }
 564         }
 565   if (overth<2) return kFALSE;
 566   //if there exist only one or less neighboring above threshold
 567   if (oversigma==8)  res = kTRUE;
 568   else if ((over==8)&&(GetNType()==8)) res=kTRUE;
 569   else if (over ==8 )
 570     for ( Int_t di=-2;di<=2;di++)
 571       for ( Int_t dj=-2;dj<=2;dj++)
 572         if ( (di==2)||(di==-2) || (dj==2)|| (dj==-2) )
 573           {
 574             Float_t virtualcell2=GetVirtualSignal(x+di,y+dj);
 575             if (virtualcell2 < 0) {
 576               over+=1;
 577               oversigma+=1;
 578             }
 579             else
 580               {
 581                 if (virtualcell>=virtualcell2) over+=1;
 582               }
 583           }
 584   if (over == 24) res=kTRUE;
 585   return res;
 586
 587 }
 588
 589
 590 void AliTPCClusterFinder::ResetSignal()
 591 {
 592    //reset dignals to 0
 593   Int_t size = fDimX*fDimY;
 594   AliCell *dig=fCells;
 595   if (rOK==kTRUE) for (Int_t i=0 ; i<size;i++) dig[i] = 0;
 596 }
 597
 598
 599
 600 void AliTPCClusterFinder::ResetStatus()
 601 {
 602    //reset status of signals to not used
 603   Int_t size = fDimX*fDimY;
 604   AliCell *dig=fCells;
 605   if (rOK==kTRUE) for (Int_t i=0 ; i<size;i++)
 606       dig[i].SetStatus(0);
 607 }
 608
 609
 610 AliCell  *  AliTPCClusterFinder::GetCell(Int_t i, Int_t j)
 611 {
 612   //return reference to the cell with index i,j
 613   if (rOK == kTRUE)
 614     if ( (i>=0) && (i<fDimX) && (j>=0) && (j<fDimY) )
 615       return &fCells[i+j*fDimX];
 616   return 0;
 617 }
 618
 619 Float_t   AliTPCClusterFinder::GetVirtualSignal(Float_t ri, Float_t rj)
 620 {
 621   //it generate virtual cell as mean value from different cels
 622   //after using it must be destructed !!!
 623   Int_t i=(Int_t)ri;
 624   Int_t j=(Int_t)rj;
 625   Int_t ddi = (ri>i)? 1:0;
 626   Int_t ddj = (rj>j)? 1:0;
 627   Float_t sum = 0;
 628   Float_t sumw= 0;
 629   for (Int_t di=0;di<=ddi;di++)
 630     for (Int_t dj=0;dj<=ddj;dj++)
 631       {
 632         Float_t w = (ri-i-di)*(ri-i-di)+(rj-j-dj)*(rj-j-dj);
 633         if (w>0) w=1/TMath::Sqrt(w);
 634         else w=9999999;
 635         AliCell * cel2 =GetCell(i+di,j+dj);
 636         Int_t signal2 = GetSignal(i+di,j+dj);
 637         if (cel2!=0) {
 638           sumw+=w;
 639           sum+= signal2*w;
 640         }
 641       }
 642   if (sumw>0)  return (sum/sumw);
 643   else
 644     return -1;
 645 }
 646
 647
 648
 649 void AliTPCClusterFinder::SetBlockIndex(Int_t * index)
 650 {
 651   //
 652   //calculate which indexes we must check for border
 653   //
 654   if (TMath::Abs(index[0])<2) index[2] = 0;
 655   else {
 656     index[2] = TMath::Abs(index[0])-1;
 657     if (index[0]<0) index[2]*=-1;   //first x block
 658   }
 659   if (TMath::Abs(index[1])<2) index[3] = 0;
 660   else {
 661     index[3] = TMath::Abs(index[1])-1;
 662     if (index[1]<0) index[3]*=-1;   //first y block
 663   }
 664   if (TMath::Abs(index[0])<TMath::Abs(index[1])){
 665     index[4]=index[0];
 666     index[5]=index[3];
 667   }
 668   else
 669     if (index[0]==index[1]) {
 670       index[4]=0;
 671       index[5]=0;
 672     }
 673     else{
 674       index[4]=index[2];
 675       index[5]=index[1];
 676     }
 677   return;
 678 }
 679
 680 //***********************************************************************
 681 //***********************************************************************
 682
 683 TClonesArray * AliTPCClusterFinder::FindPeaks1(TClonesArray *arr)
 684 {
 685   //find peaks and write it in form of AliTPCcluster to array
 686   if (arr==0){
 687     fClustersArray=new TClonesArray("AliDigitCluster",300);
 688     fIndex=1;
 689   }
 690   else {
 691     fClustersArray = arr;
 692     fIndex = fClustersArray->GetEntriesFast();
 693   }
 694
 695   AliDigitCluster c;
 696   ResetStatus();
 697    for (Int_t i=0; i<fDimX; i++)
 698      for (Int_t j=0;j<fDimY; j++)
 699        {
 700          fStackIndex=0;
 701          fBDistType = kFALSE;
 702          AliCell * cell = GetCell(i,j);
 703          if (!(cell->IsChecked()))  Adjacent(i,j);
 704          //if there exists more then 2 digits cluster
 705          if (fStackIndex >2 ){
 706            if (fBFit==kFALSE) GetClusterStatistic(c);
 707              else GetClusterFit(c);
 708            //write some important chracteristic area of cluster
 709            //
 710            Transform(&c);
 711            //write cluster information to array
 712            TClonesArray &lclusters = *fClustersArray;
 713            new (lclusters[fIndex++])  AliDigitCluster(c);
 714            //             cout<<"fx="<<c.fX<<"   fy"<<c.fY<<"\n";
 715          }
 716        }
 717    return fClustersArray;
 718 }
 719
 720
 721 TClonesArray * AliTPCClusterFinder::FindPeaks2(TClonesArray *arr)
 722 {
 723   //find peaks and write it in form of AliTPCcluster to array
 724   if (arr==0){
 725     fClustersArray=new TClonesArray("AliDigitCluster",300);
 726     fIndex=1;
 727   }
 728   else {
 729     fClustersArray = arr;
 730     fIndex = fClustersArray->GetEntriesFast();
 731   }
 732
 733   AliDigitCluster c;
 734   ResetStatus();
 735
 736    for (Int_t i=0; i<fDimX; i++)
 737      for (Int_t j=0;j<fDimY; j++)
 738        {
 739          fStackIndex=0;
 740          if (IsMaximum(i,j) == kTRUE){
 741            SetSigma2(i,j,fCurrentSigmaX2,fCurrentSigmaY2);
 742            fBDistType = kTRUE;
 743            Adjacent(i,j);
 744            //if there exists more then 2 digits cluster
 745            if (fStackIndex >2 ){
 746              if (fBFit==kFALSE) GetClusterStatistic(c);
 747              else GetClusterFit(c);
 748              //write some important chracteristic area of cluster
 749              //
 750              Transform(&c);
 751              //write cluster information to array
 752              TClonesArray &lclusters = *fClustersArray;
 753              new(lclusters[fIndex++]) AliDigitCluster(c);
 754              //             cout<<"fx="<<c.fX<<"   fy"<<c.fY<<"\n";
 755            }
 756          }
 757        }
 758    return fClustersArray;
 759 }
 760
 761
 762 TClonesArray * AliTPCClusterFinder::FindPeaks3(TClonesArray *arr)
 763 {
 764   //find peaks and write it in form of AliTPCcluster to array
 765   if (arr==0){
 766     fClustersArray=new TClonesArray("AliDigitCluster",300);
 767     fIndex=1;
 768   }
 769   else {
 770     fClustersArray = arr;
 771     fIndex = fClustersArray->GetEntriesFast();
 772   }
 773
 774   AliDigitCluster c;
 775   ResetStatus();
 776
 777   Int_t dmax=5;
 778   Int_t naccepted =1;
 779    for (Int_t i=0; i<fDimX; i++)
 780      for (Int_t j=0;j<fDimY; j++)
 781        {
 782          fStackIndex=0;
 783          if (IsMaximum(i,j) == kTRUE){
 784            SetSigma2(i,j,fCurrentSigmaX2,fCurrentSigmaY2);
 785            AddToStack(i,j,GetSignal(i,j));
 786
 787            //loop over different distance
 788            naccepted =1;
 789            for ( Int_t dd =1;((dd<=dmax) && (naccepted>0));dd++){
 790              naccepted=0;
 791              for (Int_t di = -dd;di<=dd;di++){
 792                Int_t ddj = dd-TMath::Abs(di);
 793                Int_t sigstart = (ddj>0) ?  -1 : 0;
 794                for (Int_t sig = sigstart;sig<=1;sig+=2){
 795                  Int_t dj= sig*ddj;
 796                  AliCell *cell= GetCell(i+di,j+dj);
 797                  Int_t signal = GetSignal(i+di,j+dj);
 798                  if (cell==0) continue;
 799                  Int_t index[6];
 800                  index[0]=di;
 801                  index[1]=dj;
 802                  if (dd>2) {
 803                    SetBlockIndex(index);  //adjust index to control
 804                    if ( IsBorder(fIndex,i+index[2],j+index[3]) ||
 805                         IsBorder(fIndex,i+index[4],j+index[5])) {
 806                      cell->SetBorder(fIndex);
 807                      continue;
 808                    }
 809                  }
 810                  if ( signal<=fThreshold ){
 811                    //if under threshold
 812                    cell->SetThBorder(fIndex);
 813                    if (fBFit==kTRUE)  AddToStack(i+di,j+dj,signal);
 814                    continue;
 815                  }
 816                  naccepted++;
 817                  if (CheckIfDirBorder(fCurrentMaxX,fCurrentMaxY,i+di,j+dj) == kTRUE) {
 818                    if (fBFit==kFALSE) AddToStack(i+di,j+dj,signal/2);
 819                    continue;
 820                  }
 821                  AddToStack(i+di,j+dj,signal);
 822
 823                } //loop over sig dj
 824              } //loop over di
 825
 826            }//loop over dd
 827          } //if there is maximum
 828          //if there exists more then 2 digits cluster
 829          if (fStackIndex >2 ){
 830            if (fBFit==kFALSE) GetClusterStatistic(c);
 831            else GetClusterFit(c);
 832            //write some important chracteristic area of cluster
 833            //
 834            Transform(&c);
 835            //write cluster information to array
 836            TClonesArray &lclusters = *fClustersArray;
 837            new(lclusters[fIndex++]) AliDigitCluster(c);
 838            //             cout<<"fx="<<c.fX<<"   fy"<<c.fY<<"\n";
 839          }
 840        } //lopp over all digits
 841
 842    return fClustersArray;
 843 }
 844
 845
 846
 847
 848
 849
 850 void AliTPCClusterFinder::Adjacent(Int_t i,Int_t j)
 851 {
 852   //
 853   //recursive agorithm program
 854   //
 855   if (fBDistType==kTRUE) {
 856     Float_t delta = (i-fCurrentMaxX)*(i-fCurrentMaxX)/fCurrentSigmaX2;
 857     delta+=(j-fCurrentMaxY)*(j-fCurrentMaxY)/fCurrentSigmaY2;
 858     if (delta > fMulSigma2) {
 859        SetDirBorder(fIndex,i,j);
 860       return;
 861     }
 862   }
 863   AliCell *cell = GetCell(i,j);
 864   Int_t signal = GetSignal(i,j);
 865   Int_t q=signal;
 866   cell->SetChecked(fIndex);
 867   if ( (q>fThreshold) || (fBFit==kTRUE))   AddToStack(i,j,q);
 868   if ( q >fThreshold )
 869     {
 870
 871       AliCell * newcel;
 872       newcel = GetCell(i-1,j);
 873       if (newcel !=0) if (!newcel->IsChecked(fIndex) ) Adjacent(i-1,j);
 874       newcel = GetCell(i,j-1);
 875       if (newcel !=0) if (!newcel->IsChecked(fIndex) ) Adjacent(i,j-1);
 876       newcel = GetCell(i+1,j);
 877       if (newcel !=0) if (!newcel->IsChecked(fIndex) ) Adjacent(i+1,j);
 878       newcel = GetCell(i,j+1);
 879       if (newcel !=0) if (!newcel->IsChecked(fIndex) ) Adjacent(i,j+1);
 880     }
 881   else cell->SetThBorder(fIndex);
 882 }
 883
 884
 885
 886 AliH2F *  AliTPCClusterFinder::DrawHisto( const char *option=0,
 887                               Float_t x1, Float_t x2, Float_t y1, Float_t y2)
 888 {
 889   //
 890   //draw digits in given array
 891   //
 892   //make digits histo
 893   char ch[30];
 894   sprintf(ch,"Cluster finder digits ");
 895   if ( (fDimX<1)|| (fDimY<1)) {
 896     return 0;
 897   }
 898   AliH2F * his  = new AliH2F(ch,ch,fDimX,fX1,fX2,fDimY,fY1,fY2);
 899   //set histogram  values
 900   for (Int_t i = 0; i<fDimX;i++)
 901     for (Int_t j = 0; j<fDimY;j++){
 902       Float_t x = ItoX(i);
 903       Float_t y= JtoY(j);
 904       his->Fill(x,y,GetSignal(i,j));
 905     }
 906   if (x1>=0) {
 907       AliH2F *h2fsub = his->GetSubrange2d(x1,x2,y1,y2);
 908       delete his;
 909       his=h2fsub;
 910   }
 911   if (his==0) return 0;
 912   if (option!=0) his->Draw(option);
 913   else his->Draw();
 914   return his;
 915 }
 916
 917
 918 void AliTPCClusterFinder::DrawCluster(
 919                                   Int_t color, Int_t size, Int_t style)
 920 {
 921
 922   if (fClustersArray==0) return;
 923   //draw marker for each of cluster
 924   Int_t ncl=fClustersArray->GetEntriesFast();
 925   for (Int_t i=0;i<ncl;i++){
 926     AliComplexCluster *cl = (AliComplexCluster*)fClustersArray->UncheckedAt(i);
 927     TMarker * marker = new TMarker;
 928     marker->SetX(cl->fX);
 929     marker->SetY(cl->fY);
 930     marker->SetMarkerSize(size);
 931     marker->SetMarkerStyle(style);
 932     marker->SetMarkerColor(color);
 933     marker->Draw();
 934   }
 935 }
 936
 937
 938
 939 AliH2F *  AliTPCClusterFinder::DrawBorders( const char *option,  AliH2F *h, Int_t type ,
 940                               Float_t x1, Float_t x2, Float_t y1, Float_t y2)
 941 {
 942   //
 943   //draw digits in given array
 944   //
 945   //make digits histo
 946   char ch[30];
 947   sprintf(ch,"Cluster finder digits borders");
 948   if ( (fDimX<1)|| (fDimY<1)) {
 949     return 0;
 950   }
 951   AliH2F * his;
 952   if (h!=0) his =h;
 953   else his  = new AliH2F(ch,ch,fDimX,fX1,fX2,fDimY,fY1,fY2);
 954   //set histogram  values
 955   for (Int_t i = 0; i<fDimX;i++)
 956     for (Int_t j = 0; j<fDimY;j++){
 957       Float_t x = ItoX(i);
 958       Float_t y= JtoY(j);
 959       if (((type==1)||(type==0))&&IsMaximum(0,i,j)) his->Fill(x,y,16);
 960       if (((type==3)||(type==0))&&(IsDirBorder(0,i,j))) his->Fill(x,y,8);
 961       if (((type==4)||(type==0))&&(IsThBorder(0,i,j))) his->Fill(x,y,4);
 962       if (((type==2)||(type==0))&&IsBorder(0,i,j)) his->Fill(x,y,1);
 963
 964     }
 965
 966   if (x1>=0) {
 967       AliH2F *h2fsub = his->GetSubrange2d(x1,x2,y1,y2);
 968       delete his;
 969       his=h2fsub;
 970   }
 971   if (his==0) return 0;
 972   if (option!=0) his->Draw(option);
 973   else his->Draw();
 974   return his;
 975 }