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