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