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