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