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