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