Bug fix in the HoughTransformerRow PeakFinder.
[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 <strings.h>
7 #include "AliL3StandardIncludes.h"
8
9 #ifndef no_root
10 #include <TNtuple.h>
11 #include <TFile.h>
12 #endif
13
14 #include "AliL3Logging.h"
15 #include "AliL3HoughMaxFinder.h"
16 #include "AliL3Histogram.h"
17 #include "AliL3TrackArray.h"
18 #include "AliL3HoughTrack.h"
19 #include "AliL3Transform.h"
20 #include "AliL3HoughTransformerRow.h"
21
22 #if __GNUC__ == 3
23 using namespace std;
24 #endif
25
26 /** \class AliL3HoughMaxFinder
27 <pre>
28 //_____________________________________________________________
29 // AliL3HoughMaxFinder
30 //
31 // Maximum finder
32 //
33 </pre>
34 */
35
36 ClassImp(AliL3HoughMaxFinder)
37
38   
39 AliL3HoughMaxFinder::AliL3HoughMaxFinder()
40 {
41   //Default constructor
42   fThreshold = 0;
43   fCurrentEtaSlice = -1;
44   fHistoType=0;
45   fXPeaks=0;
46   fYPeaks=0;
47   fWeight=0;
48   fNPeaks=0;
49   fN1PeaksPrevEtaSlice=0;
50   fN2PeaksPrevEtaSlice=0;
51   fSTARTXPeaks=0;
52   fSTARTYPeaks=0;
53   fENDXPeaks=0;
54   fENDYPeaks=0;
55   fSTARTETAPeaks=0;
56   fENDETAPeaks=0;
57   fNMax=0;
58   fGradX=1;
59   fGradY=1;
60 #ifndef no_root
61   fNtuppel = 0;
62 #endif
63 }
64
65 AliL3HoughMaxFinder::AliL3HoughMaxFinder(Char_t *histotype,Int_t nmax,AliL3Histogram *hist)
66 {
67   //Constructor
68
69   //fTracks = new AliL3TrackArray("AliL3HoughTrack");
70   if(strcmp(histotype,"KappaPhi")==0) fHistoType='c';
71   if(strcmp(histotype,"DPsi")==0) fHistoType='l';
72
73   fCurrentEtaSlice = -1;
74   
75   if(hist)
76     fCurrentHisto = hist;
77   
78   fGradX=1;
79   fGradY=1;
80   fNMax=nmax;
81   fXPeaks = new Float_t[fNMax];
82   fYPeaks = new Float_t[fNMax];
83   fWeight = new Int_t[fNMax];
84   fSTARTXPeaks = new Int_t[fNMax];
85   fSTARTYPeaks = new Int_t[fNMax];
86   fENDXPeaks = new Int_t[fNMax];
87   fENDYPeaks = new Int_t[fNMax];
88   fSTARTETAPeaks = new Int_t[fNMax];
89   fENDETAPeaks = new Int_t[fNMax];
90 #ifndef no_root
91   fNtuppel = 0;
92 #endif
93   fThreshold=0;
94 }
95
96 AliL3HoughMaxFinder::~AliL3HoughMaxFinder()
97 {
98   //Destructor
99   if(fXPeaks)
100     delete [] fXPeaks;
101   if(fYPeaks)
102     delete [] fYPeaks;
103   if(fWeight)
104     delete [] fWeight;
105   if(fSTARTXPeaks)
106     delete [] fSTARTXPeaks;
107   if(fSTARTYPeaks)
108     delete [] fSTARTYPeaks;
109   if(fENDXPeaks)
110     delete [] fENDXPeaks;
111   if(fENDYPeaks)
112     delete [] fENDYPeaks;
113   if(fSTARTETAPeaks)
114     delete [] fSTARTETAPeaks;
115   if(fENDETAPeaks)
116     delete [] fENDETAPeaks;
117 #ifndef no_root
118   if(fNtuppel)
119     delete fNtuppel;
120 #endif
121 }
122
123 void AliL3HoughMaxFinder::Reset()
124 {
125   // Method to reinit the Peak Finder
126   for(Int_t i=0; i<fNMax; i++)
127     {
128       fXPeaks[i]=0;
129       fYPeaks[i]=0;
130       fWeight[i]=0;
131       fSTARTXPeaks[i]=0;
132       fSTARTYPeaks[i]=0;
133       fENDXPeaks[i]=0;
134       fENDYPeaks[i]=0;
135       fSTARTETAPeaks[i]=0;
136       fENDETAPeaks[i]=0;
137     }
138   fNPeaks=0;
139   fN1PeaksPrevEtaSlice=0;
140   fN2PeaksPrevEtaSlice=0;
141 }
142
143 void AliL3HoughMaxFinder::CreateNtuppel()
144 {
145   // Fill a NTuple with the peak parameters
146 #ifndef no_root
147   //content#; neighbouring bins of the peak.
148   fNtuppel = new TNtuple("ntuppel","Peak charateristics","kappa:phi0:weigth:content3:content5:content1:content7");
149   fNtuppel->SetDirectory(0);
150 #endif  
151 }
152
153 void AliL3HoughMaxFinder::WriteNtuppel(Char_t *filename)
154 {
155   // Write the NTuple with the peak parameters
156 #ifndef no_root
157   TFile *file = TFile::Open(filename,"RECREATE");
158   if(!file)
159     {
160       cerr<<"AliL3HoughMaxFinder::WriteNtuppel : Error opening file "<<filename<<endl;
161       return;
162     }
163   fNtuppel->Write();
164   file->Close();
165 #endif
166 }
167
168 void AliL3HoughMaxFinder::FindAbsMaxima()
169 {
170   // Simple Peak Finder in the Hough space
171   if(!fCurrentHisto)
172     {
173       cerr<<"AliL3HoughMaxFinder::FindAbsMaxima : No histogram"<<endl;
174       return;
175     }
176   AliL3Histogram *hist = fCurrentHisto;
177   
178   if(hist->GetNEntries() == 0)
179     return;
180   
181   Int_t xmin = hist->GetFirstXbin();
182   Int_t xmax = hist->GetLastXbin();
183   Int_t ymin = hist->GetFirstYbin();
184   Int_t ymax = hist->GetLastYbin();  
185   Int_t bin;
186   Double_t value,maxvalue=0;
187   
188   Int_t maxxbin=0,maxybin=0;
189   for(Int_t xbin=xmin; xbin<=xmax; xbin++)
190     {
191       for(Int_t ybin=ymin; ybin<=ymax; ybin++)
192         {
193           bin = hist->GetBin(xbin,ybin);
194           value = hist->GetBinContent(bin);
195           if(value>maxvalue)
196             {
197               maxvalue = value;
198               maxxbin = xbin;
199               maxybin = ybin;
200             }
201         }
202     }
203   
204   if(maxvalue == 0)
205     return;
206   
207   if(fNPeaks > fNMax)
208     {
209       cerr<<"AliL3HoughMaxFinder::FindAbsMaxima : Array out of range : "<<fNPeaks<<endl;
210       return;
211     }
212   
213   Double_t maxx = hist->GetBinCenterX(maxxbin);
214   Double_t maxy = hist->GetBinCenterY(maxybin);
215   fXPeaks[fNPeaks] = maxx;
216   fYPeaks[fNPeaks] = maxy;
217   fWeight[fNPeaks] = (Int_t)maxvalue;
218
219   fNPeaks++;
220 #ifndef no_root
221   if(fNtuppel)
222     {
223       Int_t bin3 = hist->GetBin(maxxbin-1,maxybin);
224       Int_t bin5 = hist->GetBin(maxxbin+1,maxybin);
225       Int_t bin1 = hist->GetBin(maxxbin,maxybin-1);
226       Int_t bin7 = hist->GetBin(maxxbin,maxybin+1);
227       
228       fNtuppel->Fill(maxx,maxy,maxvalue,hist->GetBinContent(bin3),hist->GetBinContent(bin5),hist->GetBinContent(bin1),hist->GetBinContent(bin7));
229     }
230 #endif  
231 }
232
233 void AliL3HoughMaxFinder::FindBigMaxima()
234 {
235   // Another Peak finder  
236   AliL3Histogram *hist = fCurrentHisto;
237   
238   if(hist->GetNEntries() == 0)
239     return;
240   
241   Int_t xmin = hist->GetFirstXbin();
242   Int_t xmax = hist->GetLastXbin();
243   Int_t ymin = hist->GetFirstYbin();
244   Int_t ymax = hist->GetLastYbin();
245   Int_t bin[25],binindex;
246   Double_t value[25];
247   
248   for(Int_t xbin=xmin+2; xbin<xmax-3; xbin++)
249     {
250       for(Int_t ybin=ymin+2; ybin<ymax-3; ybin++)
251         {
252           binindex=0;
253           for(Int_t xb=xbin-2; xb<=xbin+2; xb++)
254             {
255               for(Int_t yb=ybin-2; yb<=ybin+2; yb++)
256                 {
257                   bin[binindex]=hist->GetBin(xb,yb);
258                   value[binindex]=hist->GetBinContent(bin[binindex]);
259                   binindex++;
260                 }
261             }
262           if(value[12]==0) continue;
263           Int_t b=0;
264           while(1)
265             {
266               if(value[b] > value[12] || b==binindex) break;
267               b++;
268               //printf("b %d\n",b);
269             }
270           if(b == binindex)
271             {
272               //Found maxima
273               if(fNPeaks > fNMax)
274                 {
275                   cerr<<"AliL3HoughMaxFinder::FindBigMaxima : Array out of range "<<fNPeaks<<endl;
276                   return;
277                 }
278               
279               Double_t maxx = hist->GetBinCenterX(xbin);
280               Double_t maxy = hist->GetBinCenterY(ybin);
281               fXPeaks[fNPeaks] = maxx;
282               fYPeaks[fNPeaks] = maxy;
283               fNPeaks++;
284             }
285         }
286     }
287 }
288
289 void AliL3HoughMaxFinder::FindMaxima(Int_t threshold)
290 {
291   //Locate all the maxima in input histogram.
292   //Maxima is defined as bins with more entries than the
293   //immediately neighbouring bins. 
294   
295   if(fCurrentHisto->GetNEntries() == 0)
296     return;
297   
298   Int_t xmin = fCurrentHisto->GetFirstXbin();
299   Int_t xmax = fCurrentHisto->GetLastXbin();
300   Int_t ymin = fCurrentHisto->GetFirstYbin();
301   Int_t ymax = fCurrentHisto->GetLastYbin();
302   Int_t bin[9];
303   Double_t value[9];
304   
305   //Float_t max_kappa = 0.001;
306   //Float_t max_phi0 = 0.08;
307
308   for(Int_t xbin=xmin+1; xbin<=xmax-1; xbin++)
309     {
310       for(Int_t ybin=ymin+1; ybin<=ymax-1; ybin++)
311         {
312           bin[0] = fCurrentHisto->GetBin(xbin-1,ybin-1);
313           bin[1] = fCurrentHisto->GetBin(xbin,ybin-1);
314           bin[2] = fCurrentHisto->GetBin(xbin+1,ybin-1);
315           bin[3] = fCurrentHisto->GetBin(xbin-1,ybin);
316           bin[4] = fCurrentHisto->GetBin(xbin,ybin);
317           bin[5] = fCurrentHisto->GetBin(xbin+1,ybin);
318           bin[6] = fCurrentHisto->GetBin(xbin-1,ybin+1);
319           bin[7] = fCurrentHisto->GetBin(xbin,ybin+1);
320           bin[8] = fCurrentHisto->GetBin(xbin+1,ybin+1);
321           value[0] = fCurrentHisto->GetBinContent(bin[0]);
322           value[1] = fCurrentHisto->GetBinContent(bin[1]);
323           value[2] = fCurrentHisto->GetBinContent(bin[2]);
324           value[3] = fCurrentHisto->GetBinContent(bin[3]);
325           value[4] = fCurrentHisto->GetBinContent(bin[4]);
326           value[5] = fCurrentHisto->GetBinContent(bin[5]);
327           value[6] = fCurrentHisto->GetBinContent(bin[6]);
328           value[7] = fCurrentHisto->GetBinContent(bin[7]);
329           value[8] = fCurrentHisto->GetBinContent(bin[8]);
330           
331           
332           if(value[4]>value[0] && value[4]>value[1] && value[4]>value[2]
333              && value[4]>value[3] && value[4]>value[5] && value[4]>value[6]
334              && value[4]>value[7] && value[4]>value[8])
335             {
336               //Found a local maxima
337               Float_t maxx = fCurrentHisto->GetBinCenterX(xbin);
338               Float_t maxy = fCurrentHisto->GetBinCenterY(ybin);
339               
340               if((Int_t)value[4] <= threshold) continue;//central bin below threshold
341               if(fNPeaks >= fNMax)
342                 {
343                   cout<<"AliL3HoughMaxFinder::FindMaxima : Array out of range "<<fNPeaks<<endl;
344                   return;
345                 }
346               
347               //Check the gradient:
348               if(value[3]/value[4] > fGradX && value[5]/value[4] > fGradX)
349                 continue;
350
351               if(value[1]/value[4] > fGradY && value[7]/value[4] > fGradY)
352                 continue;
353
354               fXPeaks[fNPeaks] = maxx;
355               fYPeaks[fNPeaks] = maxy;
356               fWeight[fNPeaks] = (Int_t)value[4];
357               fNPeaks++;
358
359               /*
360               //Check if the peak is overlapping with a previous:
361               Bool_t bigger = kFALSE;
362               for(Int_t p=0; p<fNPeaks; p++)
363                 {
364                   if(fabs(maxx - fXPeaks[p]) < max_kappa && fabs(maxy - fYPeaks[p]) < max_phi0)
365                     {
366                       bigger = kTRUE;
367                       if(value[4] > fWeight[p]) //this peak is bigger.
368                         {
369                           fXPeaks[p] = maxx;
370                           fYPeaks[p] = maxy;
371                           fWeight[p] = (Int_t)value[4];
372                         }
373                       else
374                         continue; //previous peak is bigger.
375                     }
376                 }
377               if(!bigger) //there were no overlapping peaks.
378                 {
379                   fXPeaks[fNPeaks] = maxx;
380                   fYPeaks[fNPeaks] = maxy;
381                   fWeight[fNPeaks] = (Int_t)value[4];
382                   fNPeaks++;
383                 }
384               */
385             }
386         }
387     }
388   
389 }
390
391 struct AliL3Window 
392 {
393   Int_t fStart; // Start
394   Int_t fSum; // Sum
395 };
396
397 void AliL3HoughMaxFinder::FindAdaptedPeaks(Int_t kappawindow,Float_t cutratio)
398 {
399   //Peak finder which looks for peaks with a certain shape.
400   //The first step involves a pre-peak finder, which looks for peaks
401   //in windows (size controlled by kappawindow) summing over each psi-bin.
402   //These pre-preaks are then matched between neighbouring kappa-bins to
403   //look for real 2D peaks exhbiting the typical cross-shape in the Hough circle transform.
404   //The maximum bin within this region is marked as the peak itself, and
405   //a few checks is performed to avoid the clear fake peaks (asymmetry check etc.)
406   
407   
408   AliL3Histogram *hist = fCurrentHisto;
409   
410   if(!hist)
411     {
412       cerr<<"AliL3HoughMaxFinder : No histogram!"<<endl;
413       return;
414     }
415   
416   if(hist->GetNEntries() == 0)
417     return;
418
419   Int_t xmin = hist->GetFirstXbin();
420   Int_t xmax = hist->GetLastXbin();
421   Int_t ymin = hist->GetFirstYbin();
422   Int_t ymax = hist->GetLastYbin();
423
424
425   //Start by looking for pre-peaks:
426   
427   AliL3Window **localmaxima = new AliL3Window*[hist->GetNbinsY()];
428   
429   Short_t *nmaxs = new Short_t[hist->GetNbinsY()];
430   Int_t n,lastsum,sum;
431   Bool_t sumwasrising;
432   for(Int_t ybin=ymin; ybin<=ymax; ybin++)
433     {
434       localmaxima[ybin-ymin] = new AliL3Window[hist->GetNbinsX()];
435       nmaxs[ybin-ymin] = 0;
436       sumwasrising=0;
437       lastsum=0;
438       n=0;
439       for(Int_t xbin=xmin; xbin<=xmax-kappawindow; xbin++)
440         {
441           sum=0;
442           for(Int_t lbin=xbin; lbin<xbin+kappawindow; lbin++)
443             sum += hist->GetBinContent(hist->GetBin(lbin,ybin));
444           
445           if(sum < lastsum)
446             {
447               if(sum > fThreshold)
448                 if(sumwasrising)//Previous sum was a local maxima
449                   {
450                     localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fStart = xbin-1;
451                     localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fSum = lastsum;
452                     nmaxs[ybin-ymin]++;
453                   }
454               
455               sumwasrising=0;
456             }
457           else if(sum > 0) 
458             sumwasrising=1;
459           lastsum=sum;
460         }
461     }
462   
463   Int_t match=0;
464   Int_t *starts = new Int_t[hist->GetNbinsY()+1];
465   Int_t *maxs = new Int_t[hist->GetNbinsY()+1];
466   
467   for(Int_t ybin=ymax; ybin >= ymin+1; ybin--)
468     {
469       for(Int_t i=0; i<nmaxs[ybin-ymin]; i++)
470         {
471           Int_t lw = localmaxima[ybin-ymin][i].fSum;
472
473           if(lw<0)
474             continue; //already used
475
476           Int_t maxvalue=0,maxybin=0,maxxbin=0,maxwindow=0;
477           for(Int_t k=localmaxima[ybin-ymin][i].fStart; k<localmaxima[ybin-ymin][i].fStart + kappawindow; k++)
478             if(hist->GetBinContent(hist->GetBin(k,ybin)) > maxvalue)
479               {
480                 maxvalue = hist->GetBinContent(hist->GetBin(k,ybin));
481                 maxybin = ybin;
482                 maxxbin = k;
483               }
484           
485           //start expanding in the psi-direction:
486
487           Int_t lb = localmaxima[ybin-ymin][i].fStart;
488           //Int_t ystart=ybin;
489           starts[ybin] = localmaxima[ybin-ymin][i].fStart;
490           maxs[ybin] = maxxbin;
491           Int_t yl=ybin-1,nybins=1;
492           
493           //cout<<"Starting search at ybin "<<ybin<<" start "<<lb<<" with sum "<<localmaxima[ybin-ymin][i].sum<<endl;
494           while(yl >= ymin)
495             {
496               Bool_t found=0;
497               for(Int_t j=0; j<nmaxs[yl-ymin]; j++)
498                 {
499                   if( localmaxima[yl-ymin][j].fStart - lb < 0) continue;
500                   if( localmaxima[yl-ymin][j].fStart < lb + kappawindow + match &&
501                       localmaxima[yl-ymin][j].fStart >= lb && localmaxima[yl-ymin][j].fSum > 0)
502                     {
503                       
504                       //cout<<"match at ybin "<<yl<<" yvalue "<<hist->GetBinCenterY(yl)<<" start "<<localmaxima[yl-ymin][j].start<<" sum "<<localmaxima[yl-ymin][j].sum<<endl;
505                       
506                       Int_t lmaxvalue=0,lmaxxbin=0;
507                       for(Int_t k=localmaxima[yl-ymin][j].fStart; k<localmaxima[yl-ymin][j].fStart + kappawindow; k++)
508                         {
509                           if(hist->GetBinContent(hist->GetBin(k,yl)) > maxvalue)
510                             {
511                               maxvalue = hist->GetBinContent(hist->GetBin(k,yl));
512                               maxxbin = k;
513                               maxybin = yl;
514                               maxwindow = j;
515                             }
516                           if(hist->GetBinContent(hist->GetBin(k,yl)) > lmaxvalue)//local maxima value
517                             {
518                               lmaxvalue=hist->GetBinContent(hist->GetBin(k,yl));
519                               lmaxxbin=k;
520                             }
521                         }
522                       nybins++;
523                       starts[yl] = localmaxima[yl-ymin][j].fStart;
524                       maxs[yl] = lmaxxbin;
525                       localmaxima[yl-ymin][j].fSum=-1; //Mark as used
526                       found=1;
527                       lb = localmaxima[yl-ymin][j].fStart;
528                       break;//Since we found a match in this bin, we dont have to search it anymore, goto next bin.
529                     }
530                 }
531               if(!found || yl == ymin)//no more local maximas to be matched, so write the final peak and break the expansion:
532                 {
533                   if(nybins > 4)
534                     {
535                       //cout<<"ystart "<<ystart<<" and nybins "<<nybins<<endl;
536
537                       Bool_t truepeak=kTRUE;
538                       
539                       //cout<<"Maxima found at xbin "<<maxxbin<<" ybin "<<maxybin<<" value "<<maxvalue<<endl;
540                       //cout<<"Starting to sum at xbin "<<starts[maxybin-ymin]<<endl;
541                       
542                       
543                       //Look in a window on both sides to probe the asymmetry
544                       Float_t right=0,left=0;
545                       for(Int_t w=maxxbin+1; w<=maxxbin+3; w++)
546                         {
547                           for(Int_t r=maxybin+1; r<=maxybin+3; r++)
548                             {
549                               right += (Float_t)hist->GetBinContent(hist->GetBin(w,r));
550                             }
551                         }
552                       
553                       for(Int_t w=maxxbin-1; w>=maxxbin-3; w--)
554                         {
555                           for(Int_t r=maxybin+1; r<=maxybin+3; r++)
556                             {
557                               left += (Float_t)hist->GetBinContent(hist->GetBin(w,r));
558                             }
559                         }
560                       
561                       //cout<<"ratio "<<right/left<<endl;
562                       
563                       Float_t upperratio=1,lowerratio=1;
564                       if(left)
565                         upperratio = right/left;
566                       
567                       right=left=0;
568                       for(Int_t w=maxxbin+1; w<=maxxbin+3; w++)
569                         {
570                           for(Int_t r=maxybin-1; r>=maxybin-3; r--)
571                             {
572                               right += (Float_t)hist->GetBinContent(hist->GetBin(w,r));
573                             }
574                         }
575                       
576                       for(Int_t w=maxxbin-1; w>=maxxbin-3; w--)
577                         {
578                           for(Int_t r=maxybin-1; r>=maxybin-3; r--)
579                             {
580                               left += (Float_t)hist->GetBinContent(hist->GetBin(w,r));
581                             }
582                         }
583                       
584                       //cout<<"ratio "<<left/right<<endl;
585                       
586                       if(right)
587                         lowerratio = left/right;
588                       
589                       if(upperratio > cutratio || lowerratio > cutratio)
590                         truepeak=kFALSE;
591                       
592                       if(truepeak)
593                         {
594                           
595                           fXPeaks[fNPeaks] = hist->GetBinCenterX(maxxbin);
596                           fYPeaks[fNPeaks] = hist->GetBinCenterY(maxybin);
597                           fWeight[fNPeaks] = maxvalue;
598                           fNPeaks++;
599                           
600                           /*
601                           //Calculate the peak using weigthed means:
602                           Float_t sum=0;
603                           fYPeaks[fNPeaks]=0;
604                           for(Int_t k=maxybin-1; k<=maxybin+1; k++)
605                             {
606                               Float_t lsum = 0;
607                               for(Int_t l=starts[k]; l<starts[k]+kappawindow; l++)
608                                 {
609                                   lsum += (Float_t)hist->GetBinContent(hist->GetBin(l,k));
610                                   sum += (Float_t)hist->GetBinContent(hist->GetBin(l,k));
611                                 }
612                               fYPeaks[fNPeaks] += lsum*hist->GetBinCenterY(k);
613                             }
614                           fYPeaks[fNPeaks] /= sum;
615                           Int_t ybin1,ybin2;
616                           if(fYPeaks[fNPeaks] < hist->GetBinCenterY(hist->FindYbin(fYPeaks[fNPeaks])))
617                             {
618                               ybin1 = hist->FindYbin(fYPeaks[fNPeaks])-1;
619                               ybin2 = ybin1+1;
620                             }
621                           else
622                             {
623                               ybin1 = hist->FindYbin(fYPeaks[fNPeaks]);
624                               ybin2 = ybin1+1;
625                             }
626
627                           Float_t kappa1=0,kappa2=0;
628                           sum=0;
629                           for(Int_t k=starts[ybin1]; k<starts[ybin1] + kappawindow; k++)
630                             {
631                               kappa1 += hist->GetBinCenterX(k)*hist->GetBinContent(hist->GetBin(k,ybin1));
632                               sum += (Float_t)hist->GetBinContent(hist->GetBin(k,ybin1));
633                             }
634                           kappa1 /= sum;
635                           sum=0;
636                           for(Int_t k=starts[ybin2]; k<starts[ybin2] + kappawindow; k++)
637                             {
638                               kappa2 += hist->GetBinCenterX(k)*hist->GetBinContent(hist->GetBin(k,ybin2));
639                               sum += (Float_t)hist->GetBinContent(hist->GetBin(k,ybin2));
640                             }
641                           kappa2 /= sum;
642                           
643                           fXPeaks[fNPeaks] = ( kappa1*( hist->GetBinCenterY(ybin2) - fYPeaks[fNPeaks] ) + 
644                                                kappa2*( fYPeaks[fNPeaks] - hist->GetBinCenterY(ybin1) ) )  / 
645                             (hist->GetBinCenterY(ybin2) - hist->GetBinCenterY(ybin1));
646
647                           fNPeaks++;
648                           */
649                         }
650                     }
651                   break;
652                 }
653               else
654                 yl--;//Search continues...
655             }
656         }
657     }
658
659   for(Int_t i=0; i<hist->GetNbinsY(); i++)
660     delete localmaxima[i];
661
662   delete [] localmaxima;
663   delete [] nmaxs;
664   delete [] starts;
665   delete [] maxs;
666 }
667
668 struct AliL3PreYPeak 
669 {
670   Int_t fStartPosition; // Start position in X
671   Int_t fEndPosition; // End position in X
672   Int_t fMinValue; // Minimum value inside the prepeak
673   Int_t fMaxValue; // Maximum value inside the prepeak
674   Int_t fPrevValue; // Neighbour values
675   Int_t fLeftValue; // Neighbour values
676   Int_t fRightValue; // Neighbour values
677 };
678
679 struct AliL3Pre2DPeak
680 {
681   Float_t fX; // X coordinate of the preak
682   Float_t fY; // Y coordinate of the preak
683   Float_t fSizeX; // Size of the peak
684   Float_t fSizeY; // Size of the peak
685   Int_t fStartX; // Start position of the peak
686   Int_t fStartY; // Start position of the peak
687   Int_t fEndX; // End position of the peak
688   Int_t fEndY; // End position of the peak
689   Float_t fWeight; // Weight assigned to the peak
690 };
691
692 void AliL3HoughMaxFinder::FindAdaptedRowPeaks(Int_t kappawindow,Int_t xsize,Int_t ysize)
693 {
694   // Peak finder which is working over the Hough Space provided by the AliL3HoughTransformerRow class
695   AliL3Histogram *hist = fCurrentHisto;
696   
697   if(!hist)
698     {
699       cerr<<"AliL3HoughMaxFinder : No histogram!"<<endl;
700       return;
701     }
702   
703   if(hist->GetNEntries() == 0)
704     return;
705   
706   Int_t xmin = hist->GetFirstXbin();
707   Int_t xmax = hist->GetLastXbin();
708   Int_t ymin = hist->GetFirstYbin();
709   Int_t ymax = hist->GetLastYbin();
710   
711   //Start by looking for pre-peaks:
712   
713   AliL3PreYPeak **localmaxima = new AliL3PreYPeak*[hist->GetNbinsY()];
714   
715   Short_t *nmaxs = new Short_t[hist->GetNbinsY()];
716   Int_t lastvalue=0,value=0;
717   for(Int_t ybin=ymin; ybin<=ymax; ybin++)
718     {
719       localmaxima[ybin-ymin] = new AliL3PreYPeak[hist->GetNbinsX()];
720       nmaxs[ybin-ymin] = 0;
721       lastvalue = 0;
722       Bool_t found = 0;
723       for(Int_t xbin=xmin; xbin<=xmax; xbin++)
724         {
725           value = hist->GetBinContent(hist->GetBin(xbin,ybin));
726           if(value > 0)
727             {
728               if((value - lastvalue) > 1)
729                 {
730                   localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fStartPosition = xbin;
731                   localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fEndPosition = xbin;
732                   localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMinValue = value;
733                   localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMaxValue = value;
734                   localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fPrevValue = 0;
735                   localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fLeftValue = lastvalue;
736                   found = 1;
737                 }
738               if(abs(value - lastvalue) <= 1)
739                 {
740                   if(found) {
741                     localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fEndPosition = xbin;
742                     if(value>localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMaxValue)
743                       localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMaxValue = value;
744                     if(value<localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMinValue)
745                       localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMinValue = value;
746                   }
747                 }
748               if((value - lastvalue) < -1)
749                 {
750                   if(found) {
751                     localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fRightValue = value;
752                     nmaxs[ybin-ymin]++;
753                     found = 0;
754                   }
755                 }
756             }
757           else
758             {
759               if(found) {
760                 localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fRightValue = value;
761                 nmaxs[ybin-ymin]++;
762                 found = 0;
763               }
764             }
765           lastvalue = value;
766               
767         }
768       if(found) {
769         localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fRightValue = 0;
770         nmaxs[ybin-ymin]++;
771       }
772     }
773   
774   AliL3Pre2DPeak maxima[500];
775   Int_t nmaxima = 0;
776
777   for(Int_t ybin=ymax; ybin >= ymin; ybin--)
778     {
779       for(Int_t i=0; i<nmaxs[ybin-ymin]; i++)
780         {
781           Int_t localminvalue = localmaxima[ybin-ymin][i].fMinValue;
782           Int_t localmaxvalue = localmaxima[ybin-ymin][i].fMaxValue;
783           Int_t localprevvalue = localmaxima[ybin-ymin][i].fPrevValue;
784           Int_t localnextvalue = 0;
785           Int_t localleftvalue = localmaxima[ybin-ymin][i].fLeftValue;
786           Int_t localrightvalue = localmaxima[ybin-ymin][i].fRightValue;
787
788           if(localminvalue<0)
789             continue; //already used
790
791           //start expanding in the psi-direction:
792
793           Int_t localxstart = localmaxima[ybin-ymin][i].fStartPosition;
794           Int_t localxend = localmaxima[ybin-ymin][i].fEndPosition;
795           Int_t tempxstart = localmaxima[ybin-ymin][i].fStartPosition;
796           Int_t tempxend = localmaxima[ybin-ymin][i].fEndPosition;
797
798           Int_t localy=ybin-1,nybins=1;
799
800           while(localy >= ymin)
801             {
802               Bool_t found=0;
803               for(Int_t j=0; j<nmaxs[localy-ymin]; j++)
804                 {
805                   if( (localmaxima[localy-ymin][j].fStartPosition <= (tempxend + kappawindow)) && (localmaxima[localy-ymin][j].fEndPosition >= (tempxstart - kappawindow))) 
806                     {
807                       if((localmaxima[localy-ymin][j].fMinValue <= localmaxvalue) && (localmaxima[localy-ymin][j].fMaxValue >= localminvalue))
808                         {
809                           if(localmaxima[localy-ymin][j].fEndPosition > localxend)
810                             localxend = localmaxima[localy-ymin][j].fEndPosition;
811                           if(localmaxima[localy-ymin][j].fStartPosition < localxstart)
812                             localxstart = localmaxima[localy-ymin][j].fStartPosition;
813                           tempxstart = localmaxima[localy-ymin][j].fStartPosition;
814                           tempxend = localmaxima[localy-ymin][j].fEndPosition;
815                           if(localmaxima[localy-ymin][j].fMinValue < localminvalue)
816                             localminvalue = localmaxima[localy-ymin][j].fMinValue;
817                           if(localmaxima[localy-ymin][j].fMaxValue > localmaxvalue)
818                             localmaxvalue = localmaxima[localy-ymin][j].fMaxValue;
819                           if(localmaxima[localy-ymin][j].fRightValue > localrightvalue)
820                             localrightvalue = localmaxima[localy-ymin][j].fRightValue;
821                           if(localmaxima[localy-ymin][j].fLeftValue > localleftvalue)
822                             localleftvalue = localmaxima[localy-ymin][j].fLeftValue;
823                           localmaxima[localy-ymin][j].fMinValue = -1;
824                           found = 1;
825                           nybins++;
826                           break;
827                         }
828                       else
829                         {
830                           if(localmaxvalue > localmaxima[localy-ymin][j].fPrevValue)
831                             localmaxima[localy-ymin][j].fPrevValue = localmaxvalue;
832                           if(localmaxima[localy-ymin][j].fMaxValue > localnextvalue)
833                             localnextvalue = localmaxima[localy-ymin][j].fMaxValue;
834                         }
835                     }
836                 }
837               if(!found || localy == ymin)//no more local maximas to be matched, so write the final peak and break the expansion:
838                 {
839                   if((nybins > ysize) && ((localxend-localxstart+1) > xsize) && (localmaxvalue > localprevvalue) && (localmaxvalue > localnextvalue) && (localmaxvalue > localleftvalue) && (localmaxvalue > localrightvalue))
840                   //              if((nybins > ysize) && ((localxend-localxstart+1) > xsize))
841                     {
842                       maxima[nmaxima].fX = ((Float_t)localxstart+(Float_t)localxend)/2.0;
843                       maxima[nmaxima].fY = ((Float_t)ybin+(Float_t)(localy+1))/2.0;
844                       maxima[nmaxima].fSizeX = localxend-localxstart+1;
845                       maxima[nmaxima].fSizeY = nybins;
846                       maxima[nmaxima].fWeight = (localminvalue+localmaxvalue)/2;
847                       maxima[nmaxima].fStartX = localxstart;
848                       maxima[nmaxima].fEndX = localxend;
849                       maxima[nmaxima].fStartY = localy +1;
850                       maxima[nmaxima].fEndY = ybin;
851 #ifdef do_mc
852                       //                      cout<<"Peak found at: "<<((Float_t)localxstart+(Float_t)localxend)/2.0<<" "<<((Float_t)ybin+(Float_t)(localy+1))/2.0<<" "<<localminvalue<<" "<<localmaxvalue<<" "<<" with weight "<<(localminvalue+localmaxvalue)/2<<" and size "<<localxend-localxstart+1<<" by "<<nybins<<endl;
853 #endif
854                       nmaxima++;
855                     }
856                   break;
857                 }
858               else
859                 localy--;//Search continues...
860             }
861         }
862     }
863
864   //remove fake tracks
865   Int_t nxbins = hist->GetNbinsX()+2;
866   
867   for(Int_t i = 0; i < (nmaxima - 1); i++)
868     {
869       if(maxima[i].fWeight < 0) continue;
870       for(Int_t j = i + 1; j < nmaxima; j++)
871         {
872           if(maxima[j].fWeight < 0) continue;
873           Int_t xtrack1=0,xtrack2=0,ytrack1=0,ytrack2=0;
874           Int_t deltax = 9999;
875           for(Int_t ix1 = maxima[i].fStartX; ix1 <= maxima[i].fEndX; ix1++) {
876             for(Int_t ix2 = maxima[j].fStartX; ix2 <= maxima[j].fEndX; ix2++) {
877               if(abs(ix1 - ix2) < deltax) {
878                 deltax = abs(ix1 - ix2);
879                 xtrack1 = ix1;
880                 xtrack2 = ix2;
881               }
882             }
883           }
884           Int_t deltay = 9999;
885           for(Int_t iy1 = maxima[i].fStartY; iy1 <= maxima[i].fEndY; iy1++) {
886             for(Int_t iy2 = maxima[j].fStartY; iy2 <= maxima[j].fEndY; iy2++) {
887               if(abs(iy1 - iy2) < deltay) {
888                 deltay = abs(iy1 - iy2);
889                 ytrack1 = iy1;
890                 ytrack2 = iy2;
891               }
892             }
893           }
894           Int_t firstrow1 = AliL3HoughTransformerRow::GetTrackFirstRow()[xtrack1 + nxbins*ytrack1];
895           Int_t lastrow1 = AliL3HoughTransformerRow::GetTrackLastRow()[xtrack1 + nxbins*ytrack1];
896           Int_t firstrow2 = AliL3HoughTransformerRow::GetTrackFirstRow()[xtrack1 + nxbins*ytrack1];
897           Int_t lastrow2 = AliL3HoughTransformerRow::GetTrackLastRow()[xtrack1 + nxbins*ytrack1];
898           Int_t firstrow,lastrow;
899           if(firstrow1 < firstrow2)
900             firstrow = firstrow2;
901           else
902             firstrow = firstrow1;
903
904           if(lastrow1 > lastrow2)
905             lastrow = lastrow2;
906           else
907             lastrow = lastrow1;
908          
909           AliL3HoughTrack track1;
910           Float_t x1 = hist->GetPreciseBinCenterX(xtrack1);
911           Float_t y1 = hist->GetPreciseBinCenterY(ytrack1);
912           Float_t psi1 = atan((x1-y1)/(AliL3HoughTransformerRow::GetBeta1()-AliL3HoughTransformerRow::GetBeta2()));
913           Float_t kappa1 = 2.0*(x1*cos(psi1)-AliL3HoughTransformerRow::GetBeta1()*sin(psi1));
914           track1.SetTrackParameters(kappa1,psi1,1);
915           Float_t firsthit1[3];
916           if(!track1.GetCrossingPoint(firstrow,firsthit1)) continue;
917           Float_t lasthit1[3];
918           if(!track1.GetCrossingPoint(lastrow,lasthit1)) continue;
919
920           AliL3HoughTrack track2;
921           Float_t x2 = hist->GetPreciseBinCenterX(xtrack2);
922           Float_t y2 = hist->GetPreciseBinCenterY(ytrack2);
923           Float_t psi2 = atan((x2-y2)/(AliL3HoughTransformerRow::GetBeta1()-AliL3HoughTransformerRow::GetBeta2()));
924           Float_t kappa2 = 2.0*(x2*cos(psi2)-AliL3HoughTransformerRow::GetBeta1()*sin(psi2));
925           track2.SetTrackParameters(kappa2,psi2,1);
926           Float_t firsthit2[3];
927           if(!track2.GetCrossingPoint(firstrow,firsthit2)) continue;
928           Float_t lasthit2[3];
929           if(!track2.GetCrossingPoint(lastrow,lasthit2)) continue;
930           
931           Float_t padpitchlow = AliL3Transform::GetPadPitchWidth(AliL3Transform::GetPatch(firstrow));
932           Float_t padpitchup = AliL3Transform::GetPadPitchWidth(AliL3Transform::GetPatch(lastrow));
933           // check the distance between tracks at the edges
934 #ifdef do_mc
935           //      cout<<"DEBUG Merge peaks "<<i<<" "<<j<<" "<<firsthit1[1]<<" "<<firsthit2[1]<<"     "<<lasthit1[1]<<" "<<lasthit2[1]<<endl;
936 #endif
937           //cvetan please check!!! I added a cast to Int_t
938           if((fabs(firsthit1[1]-firsthit2[1]) < 3.0*padpitchlow) && (fabs(lasthit1[1]-lasthit2[1]) < 3.0*padpitchup)) {
939             if(maxima[i].fSizeX*maxima[i].fSizeY > maxima[j].fSizeX*maxima[j].fSizeY)
940               maxima[j].fWeight = -maxima[j].fWeight;
941             if(maxima[i].fSizeX*maxima[i].fSizeY < maxima[j].fSizeX*maxima[j].fSizeY)
942               maxima[i].fWeight = -maxima[i].fWeight;
943 #ifdef do_mc
944             //      cout<<"Merge peaks "<<i<<" "<<j<<" "<<maxima[i].fWeight<<" "<<maxima[j].fWeight<<endl;
945 #endif
946           }
947         }
948     }
949
950   //merge tracks in neighbour eta slices
951   /*
952   for(Int_t i = 0; i < nmaxima; i++) {
953     if(maxima[i].fWeight > 0) {
954       fXPeaks[fNPeaks] = hist->GetPreciseBinCenterX(maxima[i].fX);
955       fYPeaks[fNPeaks] = hist->GetPreciseBinCenterY(maxima[i].fY);
956       fWeight[fNPeaks] = (Int_t)maxima[i].fWeight;
957 #ifdef do_mc
958       cout<<"Final Peak found at: "<<maxima[i].fX<<" "<<maxima[i].fY<<" "<<" "<<fXPeaks[fNPeaks]<<" "<<fYPeaks[fNPeaks]<<" with weight "<<fWeight[fNPeaks]<<" and size "<<maxima[i].fSizeX<<" by "<<maxima[i].fSizeY<<endl;
959 #endif
960       fNPeaks++;
961     }
962   }
963   */
964
965   Int_t currentnpeaks = fNPeaks;
966   for(Int_t i = 0; i < nmaxima; i++) {
967     if(maxima[i].fWeight < 0) continue;
968     Bool_t merged = kFALSE;
969     for(Int_t j = fN1PeaksPrevEtaSlice; j < fN2PeaksPrevEtaSlice; j++) {
970       if(fWeight[j] < 0) continue;
971       if((fENDETAPeaks[j]-fSTARTETAPeaks[j]) >= 1) continue;
972       if((maxima[i].fStartX <= fENDXPeaks[j]+1) && (maxima[i].fEndX >= fSTARTXPeaks[j]-1)) {
973         if((maxima[i].fStartY <= fENDYPeaks[j]+1) && (maxima[i].fEndY >= fSTARTYPeaks[j]-1)) {
974           //merge
975           merged = kTRUE;
976           fXPeaks[fNPeaks] = (hist->GetPreciseBinCenterX(maxima[i].fX)+(fENDETAPeaks[j]-fSTARTETAPeaks[j]+1)*fXPeaks[j])/(fENDETAPeaks[j]-fSTARTETAPeaks[j]+2);
977           fYPeaks[fNPeaks] = (hist->GetPreciseBinCenterY(maxima[i].fY)+(fENDETAPeaks[j]-fSTARTETAPeaks[j]+1)*fYPeaks[j])/(fENDETAPeaks[j]-fSTARTETAPeaks[j]+2);
978           fWeight[fNPeaks] = (Int_t)maxima[i].fWeight + fWeight[j];
979           fSTARTXPeaks[fNPeaks] = maxima[i].fStartX;
980           fSTARTYPeaks[fNPeaks] = maxima[i].fStartY;
981           fENDXPeaks[fNPeaks] = maxima[i].fEndX;
982           fENDYPeaks[fNPeaks] = maxima[i].fEndY;
983           fSTARTETAPeaks[fNPeaks] = fSTARTETAPeaks[j];
984           fENDETAPeaks[fNPeaks] = fCurrentEtaSlice;
985           fNPeaks++;
986           fWeight[j] = -fWeight[j];
987         }
988       }
989     }
990     fXPeaks[fNPeaks] = hist->GetPreciseBinCenterX(maxima[i].fX);
991     fYPeaks[fNPeaks] = hist->GetPreciseBinCenterY(maxima[i].fY);
992     if(!merged)
993       fWeight[fNPeaks] = (Int_t)maxima[i].fWeight;
994     else
995       fWeight[fNPeaks] = -(Int_t)maxima[i].fWeight;
996     fSTARTXPeaks[fNPeaks] = maxima[i].fStartX;
997     fSTARTYPeaks[fNPeaks] = maxima[i].fStartY;
998     fENDXPeaks[fNPeaks] = maxima[i].fEndX;
999     fENDYPeaks[fNPeaks] = maxima[i].fEndY;
1000     fSTARTETAPeaks[fNPeaks] = fCurrentEtaSlice;
1001     fENDETAPeaks[fNPeaks] = fCurrentEtaSlice;
1002     fNPeaks++;
1003   }  
1004   fN1PeaksPrevEtaSlice = currentnpeaks;    
1005   fN2PeaksPrevEtaSlice = fNPeaks;
1006
1007   for(Int_t i=0; i<hist->GetNbinsY(); i++)
1008     delete localmaxima[i];
1009
1010   delete [] localmaxima;
1011   delete [] nmaxs;
1012 }
1013
1014 void AliL3HoughMaxFinder::FindPeak1(Int_t ywindow,Int_t xbinsides)
1015 {
1016   //Testing mutliple peakfinding.
1017   //The algorithm searches the histogram for prepreaks by looking in windows
1018   //for each bin on the xaxis. The size of these windows is controlled by ywindow.
1019   //Then the prepreaks are sorted according to their weight (sum inside window),
1020   //and the peak positions are calculated by taking the weighted mean in both
1021   //x and y direction. The size of the peak in x-direction is controlled by xbinsides.
1022
1023   if(!fCurrentHisto)
1024     {
1025       printf("AliL3HoughMaxFinder::FindPeak1 : No input histogram\n");
1026       return;
1027     }  
1028   if(fCurrentHisto->GetNEntries()==0)
1029     return;
1030   
1031   //Int_t ywindow=2;
1032   //Int_t xbinsides=1;
1033   
1034   //Float_t max_kappa = 0.001;
1035   //Float_t max_phi0 = 0.08;
1036   
1037   Int_t maxsum=0;
1038   
1039   Int_t xmin = fCurrentHisto->GetFirstXbin();
1040   Int_t xmax = fCurrentHisto->GetLastXbin();
1041   Int_t ymin = fCurrentHisto->GetFirstYbin();
1042   Int_t ymax = fCurrentHisto->GetLastYbin();
1043   Int_t nbinsx = fCurrentHisto->GetNbinsX()+1;
1044   
1045   AliL3AxisWindow **windowPt = new AliL3AxisWindow*[nbinsx];
1046   AliL3AxisWindow **anotherPt = new AliL3AxisWindow*[nbinsx];
1047   
1048   for(Int_t i=0; i<nbinsx; i++)
1049     {
1050       windowPt[i] = new AliL3AxisWindow;
1051 #if defined(__DECCXX)
1052       bzero((char *)windowPt[i],sizeof(AliL3AxisWindow));
1053 #else
1054       bzero((void*)windowPt[i],sizeof(AliL3AxisWindow));
1055 #endif
1056       anotherPt[i] = windowPt[i];
1057     }
1058   
1059   for(Int_t xbin=xmin; xbin<=xmax; xbin++)
1060     {
1061       maxsum = 0;
1062       for(Int_t ybin=ymin; ybin<=ymax-ywindow; ybin++)
1063         {
1064           Int_t suminwindow=0;
1065           for(Int_t b=ybin; b<ybin+ywindow; b++)
1066             {
1067               //inside window
1068               Int_t bin = fCurrentHisto->GetBin(xbin,b);
1069               suminwindow += (Int_t)fCurrentHisto->GetBinContent(bin);
1070             }
1071           
1072           if(suminwindow > maxsum)
1073             {
1074               maxsum = suminwindow;
1075               windowPt[xbin]->fYmin = ybin;
1076               windowPt[xbin]->fYmax = ybin + ywindow;
1077               windowPt[xbin]->fWeight = suminwindow;
1078               windowPt[xbin]->fXbin = xbin;
1079             }
1080         }
1081     }
1082
1083   //Sort the windows according to the weight
1084   SortPeaks(windowPt,0,nbinsx);
1085   
1086   Float_t top,butt;
1087   for(Int_t i=0; i<nbinsx; i++)
1088     {
1089       top=butt=0;
1090       Int_t xbin = windowPt[i]->fXbin;
1091       
1092       if(xbin<xmin || xbin > xmax-1) continue;
1093       
1094       //Check if this is really a local maxima
1095       if(anotherPt[xbin-1]->fWeight > anotherPt[xbin]->fWeight ||
1096          anotherPt[xbin+1]->fWeight > anotherPt[xbin]->fWeight)
1097         continue;
1098
1099       for(Int_t j=windowPt[i]->fYmin; j<windowPt[i]->fYmax; j++)
1100         {
1101           //Calculate the mean in y direction:
1102           Int_t bin = fCurrentHisto->GetBin(windowPt[i]->fXbin,j);
1103           top += (fCurrentHisto->GetBinCenterY(j))*(fCurrentHisto->GetBinContent(bin));
1104           butt += fCurrentHisto->GetBinContent(bin);
1105         }
1106       
1107       if(butt < fThreshold)
1108         continue;
1109       
1110       fXPeaks[fNPeaks] = fCurrentHisto->GetBinCenterX(windowPt[i]->fXbin);
1111       fYPeaks[fNPeaks] = top/butt;
1112       fWeight[fNPeaks] = (Int_t)butt;
1113       //cout<<"mean in y "<<ypeaks[n]<<" on x "<<windowPt[i]->xbin<<" content "<<butt<<endl;
1114       fNPeaks++;
1115       if(fNPeaks==fNMax) 
1116         {
1117           cerr<<"AliL3HoughMaxFinder::FindPeak1 : Peak array out of range!!!"<<endl;
1118           break;
1119         }
1120     }
1121
1122   
1123   //Improve the peaks by including the region around in x.
1124   Float_t ytop,ybutt;
1125   Int_t prev;
1126   Int_t w;
1127   for(Int_t i=0; i<fNPeaks; i++)
1128     {
1129       Int_t xbin = fCurrentHisto->FindXbin(fXPeaks[i]);
1130       if(xbin - xbinsides < xmin || xbin + xbinsides > xmax) continue;
1131       top=butt=0;
1132       ytop=0,ybutt=0;     
1133       w=0;
1134       prev = xbin - xbinsides+1;
1135       for(Int_t j=xbin-xbinsides; j<=xbin+xbinsides; j++)
1136         {
1137           /*
1138           //Check if the windows are overlapping:
1139           if(anotherPt[j]->ymin > anotherPt[prev]->ymax) {prev=j; continue;}
1140           if(anotherPt[j]->ymax < anotherPt[prev]->ymin) {prev=j; continue;}
1141           prev = j;
1142           */
1143           
1144           top += fCurrentHisto->GetBinCenterX(j)*anotherPt[j]->fWeight;
1145           butt += anotherPt[j]->fWeight;
1146           
1147           for(Int_t k=anotherPt[j]->fYmin; k<anotherPt[j]->fYmax; k++)
1148             {
1149               Int_t bin = fCurrentHisto->GetBin(j,k);
1150               ytop += (fCurrentHisto->GetBinCenterY(k))*(fCurrentHisto->GetBinContent(bin));
1151               ybutt += fCurrentHisto->GetBinContent(bin);
1152               w+=(Int_t)fCurrentHisto->GetBinContent(bin);
1153             }
1154         }
1155       
1156       fXPeaks[i] = top/butt;
1157       fYPeaks[i] = ytop/ybutt;
1158       fWeight[i] = w;
1159       //cout<<"Setting weight "<<w<<" kappa "<<fXPeaks[i]<<" phi0 "<<fYPeaks[i]<<endl;
1160       
1161       /*
1162       //Check if this peak is overlapping with a previous:
1163       for(Int_t p=0; p<i; p++)
1164         {
1165           //cout<<fabs(fXPeaks[p] - fXPeaks[i])<<" "<<fabs(fYPeaks[p] - fYPeaks[i])<<endl;
1166           if(fabs(fXPeaks[p] - fXPeaks[i]) < max_kappa &&
1167              fabs(fYPeaks[p] - fYPeaks[i]) < max_phi0)
1168             {
1169               fWeight[i]=0;
1170               //break;
1171             }
1172         }
1173       */
1174     }
1175   
1176   for(Int_t i=0; i<nbinsx; i++)
1177     delete windowPt[i];
1178   delete [] windowPt;
1179   delete [] anotherPt;
1180 }
1181
1182 void AliL3HoughMaxFinder::SortPeaks(struct AliL3AxisWindow **a,Int_t first,Int_t last)
1183 {
1184   //General sorting routine
1185   //Sort according to PeakCompare()
1186
1187   static struct AliL3AxisWindow *tmp;
1188   static int i;           // "static" to save stack space
1189   int j;
1190   
1191   while (last - first > 1) {
1192     i = first;
1193     j = last;
1194     for (;;) {
1195       while (++i < last && PeakCompare(a[i], a[first]) < 0)
1196         ;
1197       while (--j > first && PeakCompare(a[j], a[first]) > 0)
1198         ;
1199       if (i >= j)
1200         break;
1201       
1202       tmp  = a[i];
1203       a[i] = a[j];
1204       a[j] = tmp;
1205     }
1206     if (j == first) {
1207       ++first;
1208       continue;
1209     }
1210     tmp = a[first];
1211     a[first] = a[j];
1212     a[j] = tmp;
1213     if (j - first < last - (j + 1)) {
1214       SortPeaks(a, first, j);
1215       first = j + 1;   // QSort(j + 1, last);
1216     } else {
1217       SortPeaks(a, j + 1, last);
1218       last = j;        // QSort(first, j);
1219     }
1220   }
1221   
1222 }
1223
1224 Int_t AliL3HoughMaxFinder::PeakCompare(struct AliL3AxisWindow *a,struct AliL3AxisWindow *b) const
1225 {
1226   // Peak comparison based on peaks weight
1227   if(a->fWeight < b->fWeight) return 1;
1228   if(a->fWeight > b->fWeight) return -1;
1229   return 0;
1230 }
1231
1232 void AliL3HoughMaxFinder::FindPeak(Int_t t1,Double_t t2,Int_t t3)
1233 {
1234   //Attempt of a more sophisticated peak finder.
1235   //Finds the best peak in the histogram, and returns the corresponding
1236   //track object.
1237
1238   if(!fCurrentHisto)
1239     {
1240       printf("AliL3HoughMaxFinder::FindPeak : No histogram!!\n");
1241       return;
1242     }
1243   AliL3Histogram *hist = fCurrentHisto;
1244   if(hist->GetNEntries()==0)
1245     return;
1246
1247   Int_t xmin = hist->GetFirstXbin();
1248   Int_t xmax = hist->GetLastXbin();
1249   Int_t ymin = hist->GetFirstYbin();
1250   Int_t ymax = hist->GetLastYbin();
1251   Int_t nbinsx = hist->GetNbinsX()+1;
1252   
1253   Int_t *m = new Int_t[nbinsx];
1254   Int_t *mlow = new Int_t[nbinsx];
1255   Int_t *mup = new Int_t[nbinsx];
1256   
1257   
1258  recompute:  //this is a goto.
1259   
1260   for(Int_t i=0; i<nbinsx; i++)
1261     {
1262       m[i]=0;
1263       mlow[i]=0;
1264       mup[i]=0;
1265     }
1266
1267   Int_t maxx=0,sum=0,maxxbin=0,bin=0;
1268
1269   for(Int_t xbin=xmin; xbin<=xmax; xbin++)
1270     {
1271       for(Int_t ybin=ymin; ybin <= ymax - t1; ybin++)
1272         {
1273           sum = 0;
1274           for(Int_t y=ybin; y <= ybin+t1; y++)
1275             {
1276               if(y>ymax) break;
1277               //Inside window
1278               bin = hist->GetBin(xbin,y);
1279               sum += (Int_t)hist->GetBinContent(bin);
1280               
1281             }
1282           if(sum > m[xbin]) //Max value locally in this xbin
1283             {
1284               m[xbin]=sum;
1285               mlow[xbin]=ybin;
1286               mup[xbin]=ybin + t1;
1287             }
1288           
1289         }
1290       
1291       if(m[xbin] > maxx) //Max value globally in x-direction
1292         {
1293           maxxbin = xbin;
1294           maxx = m[xbin];//sum;
1295         }
1296     }
1297   //printf("maxxbin %d maxx %d mlow %d mup %d\n",maxxbin,maxx,mlow[maxxbin],mup[maxxbin]);
1298   //printf("ylow %f yup %f\n",hist->GetBinCenterY(mlow[maxxbin]),hist->GetBinCenterY(mup[maxxbin]));
1299
1300   //Determine a width in the x-direction
1301   Int_t xlow=0,xup=0;
1302   
1303   for(Int_t xbin=maxxbin-1; xbin >= xmin; xbin--)
1304     {
1305       if(m[xbin] < maxx*t2)
1306         {
1307           xlow = xbin+1;
1308           break;
1309         }
1310     }
1311   for(Int_t xbin = maxxbin+1; xbin <=xmax; xbin++)
1312     {
1313       if(m[xbin] < maxx*t2)
1314         {
1315           xup = xbin-1;
1316           break;
1317         }
1318     }
1319   
1320   Double_t top=0,butt=0,value,xpeak;
1321   if(xup - xlow + 1 > t3)
1322     {
1323       t1 -= 1;
1324       printf("\nxrange out if limit xup %d xlow %d t1 %d\n\n",xlow,xup,t1);
1325       if(t1 > 1)
1326         goto recompute;
1327       else
1328         {
1329           xpeak = hist->GetBinCenterX(maxxbin);
1330           goto moveon;
1331         }
1332     }
1333   
1334   //printf("xlow %f xup %f\n",hist->GetBinCenterX(xlow),hist->GetBinCenterX(xup));
1335   //printf("Spread in x %d\n",xup-xlow +1);
1336
1337   //Now, calculate the center of mass in x-direction
1338   for(Int_t xbin=xlow; xbin <= xup; xbin++)
1339     {
1340       value = hist->GetBinCenterX(xbin);
1341       top += value*m[xbin];
1342       butt += m[xbin];
1343     }
1344   xpeak = top/butt;
1345   
1346  moveon:
1347   
1348   //Find the peak in y direction:
1349   Int_t xl = hist->FindXbin(xpeak);
1350   if(hist->GetBinCenterX(xl) > xpeak)
1351     xl--;
1352
1353   Int_t xu = xl + 1;
1354   
1355   if(hist->GetBinCenterX(xl) > xpeak || hist->GetBinCenterX(xu) <= xpeak)
1356     printf("\nAliL3HoughMaxFinder::FindPeak : Wrong xrange %f %f %f\n\n",hist->GetBinCenterX(xl),xpeak,hist->GetBinCenterX(xu));
1357     
1358     //printf("\nxlow %f xup %f\n",hist->GetBinCenterX(xl),hist->GetBinCenterX(xu));
1359
1360   value=top=butt=0;
1361   
1362   //printf("ylow %f yup %f\n",hist->GetBinCenterY(mlow[xl]),hist->GetBinCenterY(mup[xl]));
1363   //printf("ylow %f yup %f\n",hist->GetBinCenterY(mlow[xu]),hist->GetBinCenterY(mup[xu]));
1364   
1365   for(Int_t ybin=mlow[xl]; ybin <= mup[xl]; ybin++)
1366     {
1367       value = hist->GetBinCenterY(ybin);
1368       bin = hist->GetBin(xl,ybin);
1369       top += value*hist->GetBinContent(bin);
1370       butt += hist->GetBinContent(bin);
1371     }
1372   Double_t ypeaklow = top/butt;
1373   
1374   //printf("ypeaklow %f\n",ypeaklow);
1375
1376   value=top=butt=0;
1377   for(Int_t ybin=mlow[xu]; ybin <= mup[xu]; ybin++)
1378     {
1379       value = hist->GetBinCenterY(ybin);
1380       bin = hist->GetBin(xu,ybin);
1381       top += value*hist->GetBinContent(bin);
1382       butt += hist->GetBinContent(bin);
1383     }
1384   Double_t ypeakup = top/butt;
1385   
1386   //printf("ypeakup %f\n",ypeakup);
1387
1388   Double_t xvalueup = hist->GetBinCenterX(xu);
1389   Double_t xvaluelow = hist->GetBinCenterX(xl);
1390
1391   Double_t ypeak = (ypeaklow*(xvalueup - xpeak) + ypeakup*(xpeak - xvaluelow))/(xvalueup - xvaluelow);
1392
1393
1394   //Find the weight:
1395   //bin = hist->FindBin(xpeak,ypeak);
1396   //Int_t weight = (Int_t)hist->GetBinContent(bin);
1397
1398   //AliL3HoughTrack *track = new AliL3HoughTrack();
1399   //track->SetTrackParameters(xpeak,ypeak,weight);
1400   fXPeaks[fNPeaks]=xpeak;
1401   fYPeaks[fNPeaks]=ypeak;
1402   fWeight[fNPeaks]=(Int_t)hist->GetBinContent(bin);
1403   fNPeaks++;
1404   
1405   delete [] m;
1406   delete [] mlow;
1407   delete [] mup;
1408   
1409   //return track;
1410 }
1411
1412 Float_t AliL3HoughMaxFinder::GetXPeakSize(Int_t i) const
1413 {
1414   // Get X size of a peak
1415   if(i<0 || i>fNMax)
1416     {
1417       STDCERR<<"AliL3HoughMaxFinder::GetXPeakSize : Invalid index "<<i<<STDENDL;
1418       return 0;
1419     }
1420   Float_t binwidth = fCurrentHisto->GetBinWidthX();
1421   return binwidth*(fENDXPeaks[i]-fSTARTXPeaks[i]+1);
1422 }
1423
1424 Float_t AliL3HoughMaxFinder::GetYPeakSize(Int_t i) const
1425 {
1426   // Get Y size of a peak
1427   if(i<0 || i>fNMax)
1428     {
1429       STDCERR<<"AliL3HoughMaxFinder::GetYPeak : Invalid index "<<i<<STDENDL;
1430       return 0;
1431     }
1432   Float_t binwidth = fCurrentHisto->GetBinWidthY();
1433   return binwidth*(fENDYPeaks[i]-fSTARTYPeaks[i]+1);
1434 }