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