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