Added support for NEWIO, merged cern-hlt tree, updated to latest track candidate...
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughMaxFinder.cxx
1 // @(#) $Id$
2
3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright &copy ALICE HLT Group
5
6 #include "AliL3StandardIncludes.h"
7
8 #ifndef no_root
9 #include <TNtuple.h>
10 #include <TFile.h>
11 #endif
12
13 #include "AliL3Logging.h"
14 #include "AliL3HoughMaxFinder.h"
15 #include "AliL3Histogram.h"
16 #include "AliL3TrackArray.h"
17 #include "AliL3HoughTrack.h"
18
19 #if GCCVERSION == 3
20 using namespace std;
21 #endif
22
23 //_____________________________________________________________
24 // AliL3HoughMaxFinder
25 //
26 // Maximum finder
27
28 ClassImp(AliL3HoughMaxFinder)
29
30   
31 AliL3HoughMaxFinder::AliL3HoughMaxFinder()
32 {
33   //Default constructor
34   fThreshold = 0;
35   fHistoType=0;
36   fXPeaks=0;
37   fYPeaks=0;
38   fNPeaks=0;
39   fNMax=0;
40   fGradX=1;
41   fGradY=1;
42 #ifndef no_root
43   fNtuppel = 0;
44 #endif
45 }
46
47
48 AliL3HoughMaxFinder::AliL3HoughMaxFinder(Char_t *histotype,Int_t nmax,AliL3Histogram *hist)
49 {
50   //Constructor
51
52   //fTracks = new AliL3TrackArray("AliL3HoughTrack");
53   if(strcmp(histotype,"KappaPhi")==0) fHistoType='c';
54   if(strcmp(histotype,"DPsi")==0) fHistoType='l';
55   
56   if(hist)
57     fCurrentHisto = hist;
58   
59   fGradX=1;
60   fGradY=1;
61   fNMax=nmax;
62   fXPeaks = new Float_t[fNMax];
63   fYPeaks = new Float_t[fNMax];
64   fWeight = new Int_t[fNMax];
65 #ifndef no_root
66   fNtuppel = 0;
67 #endif
68   fThreshold=0;
69 }
70
71
72 AliL3HoughMaxFinder::~AliL3HoughMaxFinder()
73 {
74   //Destructor
75   if(fXPeaks)
76     delete [] fXPeaks;
77   if(fYPeaks)
78     delete [] fYPeaks;
79   if(fWeight)
80     delete [] fWeight;
81 #ifndef no_root
82   if(fNtuppel)
83     delete fNtuppel;
84 #endif
85 }
86
87 void AliL3HoughMaxFinder::Reset()
88 {
89   for(Int_t i=0; i<fNMax; i++)
90     {
91       fXPeaks[i]=0;
92       fYPeaks[i]=0;
93       fWeight[i]=0;
94     }
95   fNPeaks=0;
96 }
97
98 void AliL3HoughMaxFinder::CreateNtuppel()
99 {
100 #ifndef no_root
101   //content#; neighbouring bins of the peak.
102   fNtuppel = new TNtuple("ntuppel","Peak charateristics","kappa:phi0:weigth:content3:content5:content1:content7");
103   fNtuppel->SetDirectory(0);
104 #endif  
105 }
106
107 void AliL3HoughMaxFinder::WriteNtuppel(Char_t *filename)
108 {
109 #ifndef no_root
110   TFile *file = TFile::Open(filename,"RECREATE");
111   if(!file)
112     {
113       cerr<<"AliL3HoughMaxFinder::WriteNtuppel : Error opening file "<<filename<<endl;
114       return;
115     }
116   fNtuppel->Write();
117   file->Close();
118 #endif
119 }
120
121 void AliL3HoughMaxFinder::FindAbsMaxima()
122 {
123   
124   if(!fCurrentHisto)
125     {
126       cerr<<"AliL3HoughMaxFinder::FindAbsMaxima : No histogram"<<endl;
127       return;
128     }
129   AliL3Histogram *hist = fCurrentHisto;
130   
131   if(hist->GetNEntries() == 0)
132     return;
133   
134   Int_t xmin = hist->GetFirstXbin();
135   Int_t xmax = hist->GetLastXbin();
136   Int_t ymin = hist->GetFirstYbin();
137   Int_t ymax = hist->GetLastYbin();  
138   Int_t bin;
139   Double_t value,max_value=0;
140   
141   Int_t max_xbin=0,max_ybin=0;
142   for(Int_t xbin=xmin; xbin<=xmax; xbin++)
143     {
144       for(Int_t ybin=ymin; ybin<=ymax; ybin++)
145         {
146           bin = hist->GetBin(xbin,ybin);
147           value = hist->GetBinContent(bin);
148           if(value>max_value)
149             {
150               max_value = value;
151               max_xbin = xbin;
152               max_ybin = ybin;
153             }
154         }
155     }
156   
157   if(max_value == 0)
158     return;
159   
160   if(fNPeaks > fNMax)
161     {
162       cerr<<"AliL3HoughMaxFinder::FindAbsMaxima : Array out of range : "<<fNPeaks<<endl;
163       return;
164     }
165   
166   Double_t max_x = hist->GetBinCenterX(max_xbin);
167   Double_t max_y = hist->GetBinCenterY(max_ybin);
168   fXPeaks[fNPeaks] = max_x;
169   fYPeaks[fNPeaks] = max_y;
170   fWeight[fNPeaks] = (Int_t)max_value;
171
172   fNPeaks++;
173 #ifndef no_root
174   if(fNtuppel)
175     {
176       Int_t bin3 = hist->GetBin(max_xbin-1,max_ybin);
177       Int_t bin5 = hist->GetBin(max_xbin+1,max_ybin);
178       Int_t bin1 = hist->GetBin(max_xbin,max_ybin-1);
179       Int_t bin7 = hist->GetBin(max_xbin,max_ybin+1);
180       
181       fNtuppel->Fill(max_x,max_y,max_value,hist->GetBinContent(bin3),hist->GetBinContent(bin5),hist->GetBinContent(bin1),hist->GetBinContent(bin7));
182     }
183 #endif  
184 }
185
186 void AliL3HoughMaxFinder::FindBigMaxima()
187 {
188   
189   AliL3Histogram *hist = fCurrentHisto;
190   
191   if(hist->GetNEntries() == 0)
192     return;
193   
194   Int_t xmin = hist->GetFirstXbin();
195   Int_t xmax = hist->GetLastXbin();
196   Int_t ymin = hist->GetFirstYbin();
197   Int_t ymax = hist->GetLastYbin();
198   Int_t bin[25],bin_index;
199   Double_t value[25];
200   
201   for(Int_t xbin=xmin+2; xbin<xmax-3; xbin++)
202     {
203       for(Int_t ybin=ymin+2; ybin<ymax-3; ybin++)
204         {
205           bin_index=0;
206           for(Int_t xb=xbin-2; xb<=xbin+2; xb++)
207             {
208               for(Int_t yb=ybin-2; yb<=ybin+2; yb++)
209                 {
210                   bin[bin_index]=hist->GetBin(xb,yb);
211                   value[bin_index]=hist->GetBinContent(bin[bin_index]);
212                   bin_index++;
213                 }
214               
215             }
216           if(value[12]==0) continue;
217           Int_t b=0;
218           while(1)
219             {
220               if(value[b] > value[12] || b==bin_index) break;
221               b++;
222               //printf("b %d\n",b);
223             }
224           if(b == bin_index)
225             {
226               //Found maxima
227               if(fNPeaks > fNMax)
228                 {
229                   cerr<<"AliL3HoughMaxFinder::FindBigMaxima : Array out of range "<<fNPeaks<<endl;
230                   return;
231                 }
232               
233               Double_t max_x = hist->GetBinCenterX(xbin);
234               Double_t max_y = hist->GetBinCenterY(ybin);
235               fXPeaks[fNPeaks] = max_x;
236               fYPeaks[fNPeaks] = max_y;
237               fNPeaks++;
238             }
239         }
240     }
241 }
242
243 void AliL3HoughMaxFinder::FindMaxima(Int_t threshold)
244 {
245   //Locate all the maxima in input histogram.
246   //Maxima is defined as bins with more entries than the
247   //immediately neighbouring bins. 
248   
249   if(fCurrentHisto->GetNEntries() == 0)
250     return;
251   
252   Int_t xmin = fCurrentHisto->GetFirstXbin();
253   Int_t xmax = fCurrentHisto->GetLastXbin();
254   Int_t ymin = fCurrentHisto->GetFirstYbin();
255   Int_t ymax = fCurrentHisto->GetLastYbin();
256   Int_t bin[9];
257   Double_t value[9];
258   
259   //Float_t max_kappa = 0.001;
260   //Float_t max_phi0 = 0.08;
261
262   for(Int_t xbin=xmin+1; xbin<=xmax-1; xbin++)
263     {
264       for(Int_t ybin=ymin+1; ybin<=ymax-1; ybin++)
265         {
266           bin[0] = fCurrentHisto->GetBin(xbin-1,ybin-1);
267           bin[1] = fCurrentHisto->GetBin(xbin,ybin-1);
268           bin[2] = fCurrentHisto->GetBin(xbin+1,ybin-1);
269           bin[3] = fCurrentHisto->GetBin(xbin-1,ybin);
270           bin[4] = fCurrentHisto->GetBin(xbin,ybin);
271           bin[5] = fCurrentHisto->GetBin(xbin+1,ybin);
272           bin[6] = fCurrentHisto->GetBin(xbin-1,ybin+1);
273           bin[7] = fCurrentHisto->GetBin(xbin,ybin+1);
274           bin[8] = fCurrentHisto->GetBin(xbin+1,ybin+1);
275           value[0] = fCurrentHisto->GetBinContent(bin[0]);
276           value[1] = fCurrentHisto->GetBinContent(bin[1]);
277           value[2] = fCurrentHisto->GetBinContent(bin[2]);
278           value[3] = fCurrentHisto->GetBinContent(bin[3]);
279           value[4] = fCurrentHisto->GetBinContent(bin[4]);
280           value[5] = fCurrentHisto->GetBinContent(bin[5]);
281           value[6] = fCurrentHisto->GetBinContent(bin[6]);
282           value[7] = fCurrentHisto->GetBinContent(bin[7]);
283           value[8] = fCurrentHisto->GetBinContent(bin[8]);
284           
285           
286           if(value[4]>value[0] && value[4]>value[1] && value[4]>value[2]
287              && value[4]>value[3] && value[4]>value[5] && value[4]>value[6]
288              && value[4]>value[7] && value[4]>value[8])
289             {
290               //Found a local maxima
291               Float_t max_x = fCurrentHisto->GetBinCenterX(xbin);
292               Float_t max_y = fCurrentHisto->GetBinCenterY(ybin);
293               
294               if((Int_t)value[4] <= threshold) continue;//central bin below threshold
295               if(fNPeaks >= fNMax)
296                 {
297                   cout<<"AliL3HoughMaxFinder::FindMaxima : Array out of range "<<fNPeaks<<endl;
298                   return;
299                 }
300               
301               //Check the gradient:
302               if(value[3]/value[4] > fGradX && value[5]/value[4] > fGradX)
303                 continue;
304
305               if(value[1]/value[4] > fGradY && value[7]/value[4] > fGradY)
306                 continue;
307
308               fXPeaks[fNPeaks] = max_x;
309               fYPeaks[fNPeaks] = max_y;
310               fWeight[fNPeaks] = (Int_t)value[4];
311               fNPeaks++;
312
313               /*
314               //Check if the peak is overlapping with a previous:
315               Bool_t bigger = kFALSE;
316               for(Int_t p=0; p<fNPeaks; p++)
317                 {
318                   if(fabs(max_x - fXPeaks[p]) < max_kappa && fabs(max_y - fYPeaks[p]) < max_phi0)
319                     {
320                       bigger = kTRUE;
321                       if(value[4] > fWeight[p]) //this peak is bigger.
322                         {
323                           fXPeaks[p] = max_x;
324                           fYPeaks[p] = max_y;
325                           fWeight[p] = (Int_t)value[4];
326                         }
327                       else
328                         continue; //previous peak is bigger.
329                     }
330                 }
331               if(!bigger) //there were no overlapping peaks.
332                 {
333                   fXPeaks[fNPeaks] = max_x;
334                   fYPeaks[fNPeaks] = max_y;
335                   fWeight[fNPeaks] = (Int_t)value[4];
336                   fNPeaks++;
337                 }
338               */
339             }
340         }
341     }
342   
343 }
344
345 struct Window 
346 {
347   Int_t start;
348   Int_t sum;
349 };
350
351 void AliL3HoughMaxFinder::FindAdaptedPeaks(Int_t kappawindow)
352 {
353   //Peak finder which looks for peaks with a certain shape.
354   //The first step involves a pre-peak finder, which looks for peaks
355   //in windows (size controlled by kappawindow) summing over each psi-bin.
356   //These pre-preaks are then matched between neighbouring kappa-bins to
357   //look for real 2D peaks exhbiting the typical cross-shape in the Hough circle transform.
358   //The maximum bin within this region is marked as the peak itself, and
359   //a few checks is performed to avoid the clear fake peaks (asymmetry check etc.)
360   
361   
362   AliL3Histogram *hist = fCurrentHisto;
363   
364   if(!hist)
365     {
366       cerr<<"AliL3HoughMaxFinder : No histogram!"<<endl;
367       return;
368     }
369   
370   if(hist->GetNEntries() == 0)
371     return;
372   
373   Int_t xmin = hist->GetFirstXbin();
374   Int_t xmax = hist->GetLastXbin();
375   Int_t ymin = hist->GetFirstYbin();
376   Int_t ymax = hist->GetLastYbin();
377   
378   //Start by looking for pre-peaks:
379   
380   Window **local_maxima = new Window*[hist->GetNbinsY()];
381   
382   Short_t *nmaxs = new Short_t[hist->GetNbinsY()];
383   Int_t n,last_sum,sum;
384   Bool_t sum_was_rising;
385   for(Int_t ybin=ymin; ybin<=ymax; ybin++)
386     {
387       local_maxima[ybin-ymin] = new Window[hist->GetNbinsX()];
388       nmaxs[ybin-ymin] = 0;
389       sum_was_rising=0;
390       last_sum=0;
391       n=0;
392       for(Int_t xbin=xmin; xbin<=xmax-kappawindow; xbin++)
393         {
394           sum=0;
395           for(Int_t lbin=xbin; lbin<xbin+kappawindow; lbin++)
396             sum += hist->GetBinContent(hist->GetBin(lbin,ybin));
397           
398           if(sum < last_sum)
399             {
400               if(sum > fThreshold)
401                 if(sum_was_rising)//Previous sum was a local maxima
402                   {
403                     local_maxima[ybin-ymin][nmaxs[ybin-ymin]].start = xbin-1;
404                     local_maxima[ybin-ymin][nmaxs[ybin-ymin]].sum = last_sum;
405                     nmaxs[ybin-ymin]++;
406                   }
407               
408               sum_was_rising=0;
409             }
410           else if(sum > 0) 
411             sum_was_rising=1;
412           last_sum=sum;
413         }
414     }
415   
416   Float_t cut_ratio=0.5;
417   Int_t match=0;
418   Int_t *starts = new Int_t[hist->GetNbinsY()+1];
419   Int_t *maxs = new Int_t[hist->GetNbinsY()+1];
420   
421   for(Int_t ybin=ymax; ybin >= ymin+1; ybin--)
422     {
423       for(Int_t i=0; i<nmaxs[ybin-ymin]; i++)
424         {
425           Int_t lw = local_maxima[ybin-ymin][i].sum;
426
427           if(lw<0)
428             continue; //already used
429
430           Int_t maxvalue=0,maxybin=0,maxxbin=0,maxwindow=0;
431           for(Int_t k=local_maxima[ybin-ymin][i].start; k<local_maxima[ybin-ymin][i].start + kappawindow; k++)
432             if(hist->GetBinContent(hist->GetBin(k,ybin)) > maxvalue)
433               {
434                 maxvalue = hist->GetBinContent(hist->GetBin(k,ybin));
435                 maxybin = ybin;
436                 maxxbin = k;
437               }
438           
439           //start expanding in the psi-direction:
440
441           Int_t lb = local_maxima[ybin-ymin][i].start;
442           //Int_t ystart=ybin;
443           starts[ybin] = local_maxima[ybin-ymin][i].start;
444           maxs[ybin] = maxxbin;
445           Int_t yl=ybin-1,nybins=1;
446           
447           //cout<<"Starting search at ybin "<<ybin<<" start "<<lb<<" with sum "<<local_maxima[ybin-ymin][i].sum<<endl;
448           while(yl >= ymin)
449             {
450               Bool_t found=0;
451               for(Int_t j=0; j<nmaxs[yl-ymin]; j++)
452                 {
453                   if( local_maxima[yl-ymin][j].start - lb < 0) continue;
454                   if( local_maxima[yl-ymin][j].start < lb + kappawindow + match &&
455                       local_maxima[yl-ymin][j].start >= lb && local_maxima[yl-ymin][j].sum > 0)
456                     {
457                       
458                       //cout<<"match at ybin "<<yl<<" yvalue "<<hist->GetBinCenterY(yl)<<" start "<<local_maxima[yl-ymin][j].start<<" sum "<<local_maxima[yl-ymin][j].sum<<endl;
459                       
460                       Int_t lmaxvalue=0,lmaxxbin=0;
461                       for(Int_t k=local_maxima[yl-ymin][j].start; k<local_maxima[yl-ymin][j].start + kappawindow; k++)
462                         {
463                           if(hist->GetBinContent(hist->GetBin(k,yl)) > maxvalue)
464                             {
465                               maxvalue = hist->GetBinContent(hist->GetBin(k,yl));
466                               maxxbin = k;
467                               maxybin = yl;
468                               maxwindow = j;
469                             }
470                           if(hist->GetBinContent(hist->GetBin(k,yl)) > lmaxvalue)//local maxima value
471                             {
472                               lmaxvalue=hist->GetBinContent(hist->GetBin(k,yl));
473                               lmaxxbin=k;
474                             }
475                         }
476                       nybins++;
477                       starts[yl] = local_maxima[yl-ymin][j].start;
478                       maxs[yl] = lmaxxbin;
479                       local_maxima[yl-ymin][j].sum=-1; //Mark as used
480                       found=1;
481                       lb = local_maxima[yl-ymin][j].start;
482                       break;//Since we found a match in this bin, we dont have to search it anymore, goto next bin.
483                     }
484                 }
485               if(!found || yl == ymin)//no more local maximas to be matched, so write the final peak and break the expansion:
486                 {
487                   if(nybins > 4)
488                     {
489                       //cout<<"ystart "<<ystart<<" and nybins "<<nybins<<endl;
490
491                       Bool_t truepeak=kTRUE;
492                       
493                       //cout<<"Maxima found at xbin "<<maxxbin<<" ybin "<<maxybin<<" value "<<maxvalue<<endl;
494                       //cout<<"Starting to sum at xbin "<<starts[maxybin-ymin]<<endl;
495                       
496                       
497                       //Look in a window on both sides to probe the asymmetry
498                       Float_t right=0,left=0;
499                       for(Int_t w=maxxbin+1; w<=maxxbin+3; w++)
500                         {
501                           for(Int_t r=maxybin+1; r<=maxybin+3; r++)
502                             {
503                               right += (Float_t)hist->GetBinContent(hist->GetBin(w,r));
504                             }
505                         }
506                       
507                       for(Int_t w=maxxbin-1; w>=maxxbin-3; w--)
508                         {
509                           for(Int_t r=maxybin+1; r<=maxybin+3; r++)
510                             {
511                               left += (Float_t)hist->GetBinContent(hist->GetBin(w,r));
512                             }
513                         }
514                       
515                       //cout<<"ratio "<<right/left<<endl;
516                       
517                       Float_t upper_ratio=1,lower_ratio=1;
518                       if(left)
519                         upper_ratio = right/left;
520                       
521                       right=left=0;
522                       for(Int_t w=maxxbin+1; w<=maxxbin+3; w++)
523                         {
524                           for(Int_t r=maxybin-1; r>=maxybin-3; r--)
525                             {
526                               right += (Float_t)hist->GetBinContent(hist->GetBin(w,r));
527                             }
528                         }
529                       
530                       for(Int_t w=maxxbin-1; w>=maxxbin-3; w--)
531                         {
532                           for(Int_t r=maxybin-1; r>=maxybin-3; r--)
533                             {
534                               left += (Float_t)hist->GetBinContent(hist->GetBin(w,r));
535                             }
536                         }
537                       
538                       //cout<<"ratio "<<left/right<<endl;
539                       
540                       if(right)
541                         lower_ratio = left/right;
542                       
543                       if(upper_ratio > cut_ratio || lower_ratio > cut_ratio)
544                         truepeak=kFALSE;
545                       
546                       if(truepeak)
547                         {
548                           
549                           fXPeaks[fNPeaks] = hist->GetBinCenterX(maxxbin);
550                           fYPeaks[fNPeaks] = hist->GetBinCenterY(maxybin);
551                           fWeight[fNPeaks] = maxvalue;
552                           fNPeaks++;
553                           
554                           /*
555                           //Calculate the peak using weigthed means:
556                           Float_t sum=0;
557                           fYPeaks[fNPeaks]=0;
558                           for(Int_t k=maxybin-1; k<=maxybin+1; k++)
559                             {
560                               Float_t lsum = 0;
561                               for(Int_t l=starts[k]; l<starts[k]+kappawindow; l++)
562                                 {
563                                   lsum += (Float_t)hist->GetBinContent(hist->GetBin(l,k));
564                                   sum += (Float_t)hist->GetBinContent(hist->GetBin(l,k));
565                                 }
566                               fYPeaks[fNPeaks] += lsum*hist->GetBinCenterY(k);
567                             }
568                           fYPeaks[fNPeaks] /= sum;
569                           Int_t ybin1,ybin2;
570                           if(fYPeaks[fNPeaks] < hist->GetBinCenterY(hist->FindYbin(fYPeaks[fNPeaks])))
571                             {
572                               ybin1 = hist->FindYbin(fYPeaks[fNPeaks])-1;
573                               ybin2 = ybin1+1;
574                             }
575                           else
576                             {
577                               ybin1 = hist->FindYbin(fYPeaks[fNPeaks]);
578                               ybin2 = ybin1+1;
579                             }
580
581                           Float_t kappa1=0,kappa2=0;
582                           sum=0;
583                           for(Int_t k=starts[ybin1]; k<starts[ybin1] + kappawindow; k++)
584                             {
585                               kappa1 += hist->GetBinCenterX(k)*hist->GetBinContent(hist->GetBin(k,ybin1));
586                               sum += (Float_t)hist->GetBinContent(hist->GetBin(k,ybin1));
587                             }
588                           kappa1 /= sum;
589                           sum=0;
590                           for(Int_t k=starts[ybin2]; k<starts[ybin2] + kappawindow; k++)
591                             {
592                               kappa2 += hist->GetBinCenterX(k)*hist->GetBinContent(hist->GetBin(k,ybin2));
593                               sum += (Float_t)hist->GetBinContent(hist->GetBin(k,ybin2));
594                             }
595                           kappa2 /= sum;
596                           
597                           fXPeaks[fNPeaks] = ( kappa1*( hist->GetBinCenterY(ybin2) - fYPeaks[fNPeaks] ) + 
598                                                kappa2*( fYPeaks[fNPeaks] - hist->GetBinCenterY(ybin1) ) )  / 
599                             (hist->GetBinCenterY(ybin2) - hist->GetBinCenterY(ybin1));
600
601                           fNPeaks++;
602                           */
603                         }
604                     }
605                   break;
606                 }
607               else
608                 yl--;//Search continues...
609             }
610         }
611     }
612
613   for(Int_t i=0; i<hist->GetNbinsY(); i++)
614     delete local_maxima[i];
615
616   delete [] local_maxima;
617   delete [] nmaxs;
618   delete [] starts;
619   delete [] maxs;
620 }
621
622
623 void AliL3HoughMaxFinder::FindPeak1(Int_t y_window,Int_t x_bin_sides)
624 {
625   //Testing mutliple peakfinding.
626   //The algorithm searches the histogram for prepreaks by looking in windows
627   //for each bin on the xaxis. The size of these windows is controlled by y_window.
628   //Then the prepreaks are sorted according to their weight (sum inside window),
629   //and the peak positions are calculated by taking the weighted mean in both
630   //x and y direction. The size of the peak in x-direction is controlled by x_bin_sides.
631
632   if(!fCurrentHisto)
633     {
634       printf("AliL3HoughMaxFinder::FindPeak1 : No input histogram\n");
635       return;
636     }  
637   if(fCurrentHisto->GetNEntries()==0)
638     return;
639   
640   //Int_t y_window=2;
641   //Int_t x_bin_sides=1;
642   
643   //Float_t max_kappa = 0.001;
644   //Float_t max_phi0 = 0.08;
645   
646   Int_t max_sum=0;
647   
648   Int_t xmin = fCurrentHisto->GetFirstXbin();
649   Int_t xmax = fCurrentHisto->GetLastXbin();
650   Int_t ymin = fCurrentHisto->GetFirstYbin();
651   Int_t ymax = fCurrentHisto->GetLastYbin();
652   Int_t nbinsx = fCurrentHisto->GetNbinsX()+1;
653   
654   AxisWindow **windowPt = new AxisWindow*[nbinsx];
655   AxisWindow **anotherPt = new AxisWindow*[nbinsx];
656   
657   for(Int_t i=0; i<nbinsx; i++)
658     {
659       windowPt[i] = new AxisWindow;
660       bzero((void*)windowPt[i],sizeof(AxisWindow));
661       anotherPt[i] = windowPt[i];
662     }
663   
664   for(Int_t xbin=xmin; xbin<=xmax; xbin++)
665     {
666       max_sum = 0;
667       for(Int_t ybin=ymin; ybin<=ymax-y_window; ybin++)
668         {
669           Int_t sum_in_window=0;
670           for(Int_t b=ybin; b<ybin+y_window; b++)
671             {
672               //inside window
673               Int_t bin = fCurrentHisto->GetBin(xbin,b);
674               sum_in_window += (Int_t)fCurrentHisto->GetBinContent(bin);
675             }
676           
677           if(sum_in_window > max_sum)
678             {
679               max_sum = sum_in_window;
680               windowPt[xbin]->ymin = ybin;
681               windowPt[xbin]->ymax = ybin + y_window;
682               windowPt[xbin]->weight = sum_in_window;
683               windowPt[xbin]->xbin = xbin;
684             }
685         }
686     }
687
688   //Sort the windows according to the weight
689   SortPeaks(windowPt,0,nbinsx);
690   
691   Float_t top,butt;
692   for(Int_t i=0; i<nbinsx; i++)
693     {
694       top=butt=0;
695       Int_t xbin = windowPt[i]->xbin;
696       
697       if(xbin<xmin || xbin > xmax-1) continue;
698       
699       //Check if this is really a local maxima
700       if(anotherPt[xbin-1]->weight > anotherPt[xbin]->weight ||
701          anotherPt[xbin+1]->weight > anotherPt[xbin]->weight)
702         continue;
703
704       for(Int_t j=windowPt[i]->ymin; j<windowPt[i]->ymax; j++)
705         {
706           //Calculate the mean in y direction:
707           Int_t bin = fCurrentHisto->GetBin(windowPt[i]->xbin,j);
708           top += (fCurrentHisto->GetBinCenterY(j))*(fCurrentHisto->GetBinContent(bin));
709           butt += fCurrentHisto->GetBinContent(bin);
710         }
711       
712       if(butt < fThreshold)
713         continue;
714       
715       fXPeaks[fNPeaks] = fCurrentHisto->GetBinCenterX(windowPt[i]->xbin);
716       fYPeaks[fNPeaks] = top/butt;
717       fWeight[fNPeaks] = (Int_t)butt;
718       //cout<<"mean in y "<<ypeaks[n]<<" on x "<<windowPt[i]->xbin<<" content "<<butt<<endl;
719       fNPeaks++;
720       if(fNPeaks==fNMax) 
721         {
722           cerr<<"AliL3HoughMaxFinder::FindPeak1 : Peak array out of range!!!"<<endl;
723           break;
724         }
725     }
726
727   
728   //Improve the peaks by including the region around in x.
729   Float_t ytop,ybutt;
730   Int_t prev;
731   Int_t w;
732   for(Int_t i=0; i<fNPeaks; i++)
733     {
734       Int_t xbin = fCurrentHisto->FindXbin(fXPeaks[i]);
735       if(xbin - x_bin_sides < xmin || xbin + x_bin_sides > xmax) continue;
736       top=butt=0;
737       ytop=0,ybutt=0;     
738       w=0;
739       prev = xbin - x_bin_sides+1;
740       for(Int_t j=xbin-x_bin_sides; j<=xbin+x_bin_sides; j++)
741         {
742           /*
743           //Check if the windows are overlapping:
744           if(anotherPt[j]->ymin > anotherPt[prev]->ymax) {prev=j; continue;}
745           if(anotherPt[j]->ymax < anotherPt[prev]->ymin) {prev=j; continue;}
746           prev = j;
747           */
748           
749           top += fCurrentHisto->GetBinCenterX(j)*anotherPt[j]->weight;
750           butt += anotherPt[j]->weight;
751           
752           for(Int_t k=anotherPt[j]->ymin; k<anotherPt[j]->ymax; k++)
753             {
754               Int_t bin = fCurrentHisto->GetBin(j,k);
755               ytop += (fCurrentHisto->GetBinCenterY(k))*(fCurrentHisto->GetBinContent(bin));
756               ybutt += fCurrentHisto->GetBinContent(bin);
757               w+=(Int_t)fCurrentHisto->GetBinContent(bin);
758             }
759         }
760       
761       fXPeaks[i] = top/butt;
762       fYPeaks[i] = ytop/ybutt;
763       fWeight[i] = w;
764       //cout<<"Setting weight "<<w<<" kappa "<<fXPeaks[i]<<" phi0 "<<fYPeaks[i]<<endl;
765       
766       /*
767       //Check if this peak is overlapping with a previous:
768       for(Int_t p=0; p<i; p++)
769         {
770           //cout<<fabs(fXPeaks[p] - fXPeaks[i])<<" "<<fabs(fYPeaks[p] - fYPeaks[i])<<endl;
771           if(fabs(fXPeaks[p] - fXPeaks[i]) < max_kappa &&
772              fabs(fYPeaks[p] - fYPeaks[i]) < max_phi0)
773             {
774               fWeight[i]=0;
775               //break;
776             }
777         }
778       */
779     }
780   
781   for(Int_t i=0; i<nbinsx; i++)
782     delete windowPt[i];
783   delete [] windowPt;
784   delete [] anotherPt;
785   
786 }
787
788 void AliL3HoughMaxFinder::SortPeaks(struct AxisWindow **a,Int_t first,Int_t last)
789 {
790   //General sorting routine
791   //Sort according to PeakCompare()
792
793   static struct AxisWindow *tmp;
794   static int i;           // "static" to save stack space
795   int j;
796   
797   while (last - first > 1) {
798     i = first;
799     j = last;
800     for (;;) {
801       while (++i < last && PeakCompare(a[i], a[first]) < 0)
802         ;
803       while (--j > first && PeakCompare(a[j], a[first]) > 0)
804         ;
805       if (i >= j)
806         break;
807       
808       tmp  = a[i];
809       a[i] = a[j];
810       a[j] = tmp;
811     }
812     if (j == first) {
813       ++first;
814       continue;
815     }
816     tmp = a[first];
817     a[first] = a[j];
818     a[j] = tmp;
819     if (j - first < last - (j + 1)) {
820       SortPeaks(a, first, j);
821       first = j + 1;   // QSort(j + 1, last);
822     } else {
823       SortPeaks(a, j + 1, last);
824       last = j;        // QSort(first, j);
825     }
826   }
827   
828 }
829
830 Int_t AliL3HoughMaxFinder::PeakCompare(struct AxisWindow *a,struct AxisWindow *b)
831 {
832   if(a->weight < b->weight) return 1;
833   if(a->weight > b->weight) return -1;
834   return 0;
835
836 }
837
838 void AliL3HoughMaxFinder::FindPeak(Int_t t1,Double_t t2,Int_t t3)
839 {
840   //Attempt of a more sophisticated peak finder.
841   //Finds the best peak in the histogram, and returns the corresponding
842   //track object.
843
844   if(!fCurrentHisto)
845     {
846       printf("AliL3HoughMaxFinder::FindPeak : No histogram!!\n");
847       return;
848     }
849   AliL3Histogram *hist = fCurrentHisto;
850   if(hist->GetNEntries()==0)
851     return;
852
853   Int_t xmin = hist->GetFirstXbin();
854   Int_t xmax = hist->GetLastXbin();
855   Int_t ymin = hist->GetFirstYbin();
856   Int_t ymax = hist->GetLastYbin();
857   Int_t nbinsx = hist->GetNbinsX()+1;
858   
859   Int_t *m = new Int_t[nbinsx];
860   Int_t *m_low = new Int_t[nbinsx];
861   Int_t *m_up = new Int_t[nbinsx];
862   
863   
864  recompute:  //this is a goto.
865   
866   for(Int_t i=0; i<nbinsx; i++)
867     {
868       m[i]=0;
869       m_low[i]=0;
870       m_up[i]=0;
871     }
872
873   Int_t max_x=0,sum=0,max_xbin=0,bin=0;
874
875   for(Int_t xbin=xmin; xbin<=xmax; xbin++)
876     {
877       for(Int_t ybin=ymin; ybin <= ymax - t1; ybin++)
878         {
879           sum = 0;
880           for(Int_t y=ybin; y <= ybin+t1; y++)
881             {
882               if(y>ymax) break;
883               //Inside window
884               bin = hist->GetBin(xbin,y);
885               sum += (Int_t)hist->GetBinContent(bin);
886               
887             }
888           if(sum > m[xbin]) //Max value locally in this xbin
889             {
890               m[xbin]=sum;
891               m_low[xbin]=ybin;
892               m_up[xbin]=ybin + t1;
893             }
894           
895         }
896       
897       if(m[xbin] > max_x) //Max value globally in x-direction
898         {
899           max_xbin = xbin;
900           max_x = m[xbin];//sum;
901         }
902     }
903   //printf("max_xbin %d max_x %d m_low %d m_up %d\n",max_xbin,max_x,m_low[max_xbin],m_up[max_xbin]);
904   //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[max_xbin]),hist->GetBinCenterY(m_up[max_xbin]));
905
906   //Determine a width in the x-direction
907   Int_t x_low=0,x_up=0;
908   
909   for(Int_t xbin=max_xbin-1; xbin >= xmin; xbin--)
910     {
911       if(m[xbin] < max_x*t2)
912         {
913           x_low = xbin+1;
914           break;
915         }
916     }
917   for(Int_t xbin = max_xbin+1; xbin <=xmax; xbin++)
918     {
919       if(m[xbin] < max_x*t2)
920         {
921           x_up = xbin-1;
922           break;
923         }
924     }
925   
926   Double_t top=0,butt=0,value,x_peak;
927   if(x_up - x_low + 1 > t3)
928     {
929       t1 -= 1;
930       printf("\nxrange out if limit x_up %d x_low %d t1 %d\n\n",x_low,x_up,t1);
931       if(t1 > 1)
932         goto recompute;
933       else
934         {
935           x_peak = hist->GetBinCenterX(max_xbin);
936           goto moveon;
937         }
938     }
939   
940   //printf("xlow %f xup %f\n",hist->GetBinCenterX(x_low),hist->GetBinCenterX(x_up));
941   //printf("Spread in x %d\n",x_up-x_low +1);
942
943   //Now, calculate the center of mass in x-direction
944   for(Int_t xbin=x_low; xbin <= x_up; xbin++)
945     {
946       value = hist->GetBinCenterX(xbin);
947       top += value*m[xbin];
948       butt += m[xbin];
949     }
950   x_peak = top/butt;
951   
952  moveon:
953   
954   //Find the peak in y direction:
955   Int_t x_l = hist->FindXbin(x_peak);
956   if(hist->GetBinCenterX(x_l) > x_peak)
957     x_l--;
958
959   Int_t x_u = x_l + 1;
960   
961   if(hist->GetBinCenterX(x_l) > x_peak || hist->GetBinCenterX(x_u) <= x_peak)
962     printf("\nAliL3HoughMaxFinder::FindPeak : Wrong xrange %f %f %f\n\n",hist->GetBinCenterX(x_l),x_peak,hist->GetBinCenterX(x_u));
963     
964     //printf("\nxlow %f xup %f\n",hist->GetBinCenterX(x_l),hist->GetBinCenterX(x_u));
965
966   value=top=butt=0;
967   
968   //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_l]),hist->GetBinCenterY(m_up[x_l]));
969   //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_u]),hist->GetBinCenterY(m_up[x_u]));
970   
971   for(Int_t ybin=m_low[x_l]; ybin <= m_up[x_l]; ybin++)
972     {
973       value = hist->GetBinCenterY(ybin);
974       bin = hist->GetBin(x_l,ybin);
975       top += value*hist->GetBinContent(bin);
976       butt += hist->GetBinContent(bin);
977     }
978   Double_t y_peak_low = top/butt;
979   
980   //printf("y_peak_low %f\n",y_peak_low);
981
982   value=top=butt=0;
983   for(Int_t ybin=m_low[x_u]; ybin <= m_up[x_u]; ybin++)
984     {
985       value = hist->GetBinCenterY(ybin);
986       bin = hist->GetBin(x_u,ybin);
987       top += value*hist->GetBinContent(bin);
988       butt += hist->GetBinContent(bin);
989     }
990   Double_t y_peak_up = top/butt;
991   
992   //printf("y_peak_up %f\n",y_peak_up);
993
994   Double_t x_value_up = hist->GetBinCenterX(x_u);
995   Double_t x_value_low = hist->GetBinCenterX(x_l);
996
997   Double_t y_peak = (y_peak_low*(x_value_up - x_peak) + y_peak_up*(x_peak - x_value_low))/(x_value_up - x_value_low);
998
999
1000   //Find the weight:
1001   //bin = hist->FindBin(x_peak,y_peak);
1002   //Int_t weight = (Int_t)hist->GetBinContent(bin);
1003
1004   //AliL3HoughTrack *track = new AliL3HoughTrack();
1005   //track->SetTrackParameters(x_peak,y_peak,weight);
1006   fXPeaks[fNPeaks]=x_peak;
1007   fYPeaks[fNPeaks]=y_peak;
1008   fWeight[fNPeaks]=(Int_t)hist->GetBinContent(bin);
1009   fNPeaks++;
1010   
1011   delete [] m;
1012   delete [] m_low;
1013   delete [] m_up;
1014   
1015   //return track;
1016
1017     
1018 }
1019