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