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