3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright © ALICE HLT Group
7 #include "AliL3StandardIncludes.h"
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"
26 /** \class AliL3HoughMaxFinder
28 //_____________________________________________________________
29 // AliL3HoughMaxFinder
36 ClassImp(AliL3HoughMaxFinder)
39 AliL3HoughMaxFinder::AliL3HoughMaxFinder()
43 fCurrentEtaSlice = -1;
49 fN1PeaksPrevEtaSlice=0;
50 fN2PeaksPrevEtaSlice=0;
65 AliL3HoughMaxFinder::AliL3HoughMaxFinder(Char_t *histotype,Int_t nmax,AliL3Histogram *hist)
69 //fTracks = new AliL3TrackArray("AliL3HoughTrack");
70 if(strcmp(histotype,"KappaPhi")==0) fHistoType='c';
71 if(strcmp(histotype,"DPsi")==0) fHistoType='l';
73 fCurrentEtaSlice = -1;
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];
96 AliL3HoughMaxFinder::~AliL3HoughMaxFinder()
106 delete [] fSTARTXPeaks;
108 delete [] fSTARTYPeaks;
110 delete [] fENDXPeaks;
112 delete [] fENDYPeaks;
114 delete [] fSTARTETAPeaks;
116 delete [] fENDETAPeaks;
123 void AliL3HoughMaxFinder::Reset()
125 // Method to reinit the Peak Finder
126 for(Int_t i=0; i<fNMax; i++)
139 fN1PeaksPrevEtaSlice=0;
140 fN2PeaksPrevEtaSlice=0;
143 void AliL3HoughMaxFinder::CreateNtuppel()
145 // Fill a NTuple with the peak parameters
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);
153 void AliL3HoughMaxFinder::WriteNtuppel(Char_t *filename)
155 // Write the NTuple with the peak parameters
157 TFile *file = TFile::Open(filename,"RECREATE");
160 cerr<<"AliL3HoughMaxFinder::WriteNtuppel : Error opening file "<<filename<<endl;
168 void AliL3HoughMaxFinder::FindAbsMaxima()
170 // Simple Peak Finder in the Hough space
173 cerr<<"AliL3HoughMaxFinder::FindAbsMaxima : No histogram"<<endl;
176 AliL3Histogram *hist = fCurrentHisto;
178 if(hist->GetNEntries() == 0)
181 Int_t xmin = hist->GetFirstXbin();
182 Int_t xmax = hist->GetLastXbin();
183 Int_t ymin = hist->GetFirstYbin();
184 Int_t ymax = hist->GetLastYbin();
186 Double_t value,maxvalue=0;
188 Int_t maxxbin=0,maxybin=0;
189 for(Int_t xbin=xmin; xbin<=xmax; xbin++)
191 for(Int_t ybin=ymin; ybin<=ymax; ybin++)
193 bin = hist->GetBin(xbin,ybin);
194 value = hist->GetBinContent(bin);
209 cerr<<"AliL3HoughMaxFinder::FindAbsMaxima : Array out of range : "<<fNPeaks<<endl;
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;
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);
228 fNtuppel->Fill(maxx,maxy,maxvalue,hist->GetBinContent(bin3),hist->GetBinContent(bin5),hist->GetBinContent(bin1),hist->GetBinContent(bin7));
233 void AliL3HoughMaxFinder::FindBigMaxima()
235 // Another Peak finder
236 AliL3Histogram *hist = fCurrentHisto;
238 if(hist->GetNEntries() == 0)
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;
248 for(Int_t xbin=xmin+2; xbin<xmax-3; xbin++)
250 for(Int_t ybin=ymin+2; ybin<ymax-3; ybin++)
253 for(Int_t xb=xbin-2; xb<=xbin+2; xb++)
255 for(Int_t yb=ybin-2; yb<=ybin+2; yb++)
257 bin[binindex]=hist->GetBin(xb,yb);
258 value[binindex]=hist->GetBinContent(bin[binindex]);
262 if(value[12]==0) continue;
266 if(value[b] > value[12] || b==binindex) break;
268 //printf("b %d\n",b);
275 cerr<<"AliL3HoughMaxFinder::FindBigMaxima : Array out of range "<<fNPeaks<<endl;
279 Double_t maxx = hist->GetBinCenterX(xbin);
280 Double_t maxy = hist->GetBinCenterY(ybin);
281 fXPeaks[fNPeaks] = maxx;
282 fYPeaks[fNPeaks] = maxy;
289 void AliL3HoughMaxFinder::FindMaxima(Int_t threshold)
291 //Locate all the maxima in input histogram.
292 //Maxima is defined as bins with more entries than the
293 //immediately neighbouring bins.
295 if(fCurrentHisto->GetNEntries() == 0)
298 Int_t xmin = fCurrentHisto->GetFirstXbin();
299 Int_t xmax = fCurrentHisto->GetLastXbin();
300 Int_t ymin = fCurrentHisto->GetFirstYbin();
301 Int_t ymax = fCurrentHisto->GetLastYbin();
305 //Float_t max_kappa = 0.001;
306 //Float_t max_phi0 = 0.08;
308 for(Int_t xbin=xmin+1; xbin<=xmax-1; xbin++)
310 for(Int_t ybin=ymin+1; ybin<=ymax-1; ybin++)
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]);
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])
336 //Found a local maxima
337 Float_t maxx = fCurrentHisto->GetBinCenterX(xbin);
338 Float_t maxy = fCurrentHisto->GetBinCenterY(ybin);
340 if((Int_t)value[4] <= threshold) continue;//central bin below threshold
343 cout<<"AliL3HoughMaxFinder::FindMaxima : Array out of range "<<fNPeaks<<endl;
347 //Check the gradient:
348 if(value[3]/value[4] > fGradX && value[5]/value[4] > fGradX)
351 if(value[1]/value[4] > fGradY && value[7]/value[4] > fGradY)
354 fXPeaks[fNPeaks] = maxx;
355 fYPeaks[fNPeaks] = maxy;
356 fWeight[fNPeaks] = (Int_t)value[4];
360 //Check if the peak is overlapping with a previous:
361 Bool_t bigger = kFALSE;
362 for(Int_t p=0; p<fNPeaks; p++)
364 if(fabs(maxx - fXPeaks[p]) < max_kappa && fabs(maxy - fYPeaks[p]) < max_phi0)
367 if(value[4] > fWeight[p]) //this peak is bigger.
371 fWeight[p] = (Int_t)value[4];
374 continue; //previous peak is bigger.
377 if(!bigger) //there were no overlapping peaks.
379 fXPeaks[fNPeaks] = maxx;
380 fYPeaks[fNPeaks] = maxy;
381 fWeight[fNPeaks] = (Int_t)value[4];
393 Int_t fStart; // Start
397 void AliL3HoughMaxFinder::FindAdaptedPeaks(Int_t kappawindow,Float_t cutratio)
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.)
408 AliL3Histogram *hist = fCurrentHisto;
412 cerr<<"AliL3HoughMaxFinder : No histogram!"<<endl;
416 if(hist->GetNEntries() == 0)
419 Int_t xmin = hist->GetFirstXbin();
420 Int_t xmax = hist->GetLastXbin();
421 Int_t ymin = hist->GetFirstYbin();
422 Int_t ymax = hist->GetLastYbin();
425 //Start by looking for pre-peaks:
427 AliL3Window **localmaxima = new AliL3Window*[hist->GetNbinsY()];
429 Short_t *nmaxs = new Short_t[hist->GetNbinsY()];
432 for(Int_t ybin=ymin; ybin<=ymax; ybin++)
434 localmaxima[ybin-ymin] = new AliL3Window[hist->GetNbinsX()];
435 nmaxs[ybin-ymin] = 0;
439 for(Int_t xbin=xmin; xbin<=xmax-kappawindow; xbin++)
442 for(Int_t lbin=xbin; lbin<xbin+kappawindow; lbin++)
443 sum += hist->GetBinContent(hist->GetBin(lbin,ybin));
448 if(sumwasrising)//Previous sum was a local maxima
450 localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fStart = xbin-1;
451 localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fSum = lastsum;
464 Int_t *starts = new Int_t[hist->GetNbinsY()+1];
465 Int_t *maxs = new Int_t[hist->GetNbinsY()+1];
467 for(Int_t ybin=ymax; ybin >= ymin+1; ybin--)
469 for(Int_t i=0; i<nmaxs[ybin-ymin]; i++)
471 Int_t lw = localmaxima[ybin-ymin][i].fSum;
474 continue; //already used
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)
480 maxvalue = hist->GetBinContent(hist->GetBin(k,ybin));
485 //start expanding in the psi-direction:
487 Int_t lb = localmaxima[ybin-ymin][i].fStart;
489 starts[ybin] = localmaxima[ybin-ymin][i].fStart;
490 maxs[ybin] = maxxbin;
491 Int_t yl=ybin-1,nybins=1;
493 //cout<<"Starting search at ybin "<<ybin<<" start "<<lb<<" with sum "<<localmaxima[ybin-ymin][i].sum<<endl;
497 for(Int_t j=0; j<nmaxs[yl-ymin]; j++)
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)
504 //cout<<"match at ybin "<<yl<<" yvalue "<<hist->GetBinCenterY(yl)<<" start "<<localmaxima[yl-ymin][j].start<<" sum "<<localmaxima[yl-ymin][j].sum<<endl;
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++)
509 if(hist->GetBinContent(hist->GetBin(k,yl)) > maxvalue)
511 maxvalue = hist->GetBinContent(hist->GetBin(k,yl));
516 if(hist->GetBinContent(hist->GetBin(k,yl)) > lmaxvalue)//local maxima value
518 lmaxvalue=hist->GetBinContent(hist->GetBin(k,yl));
523 starts[yl] = localmaxima[yl-ymin][j].fStart;
525 localmaxima[yl-ymin][j].fSum=-1; //Mark as used
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.
531 if(!found || yl == ymin)//no more local maximas to be matched, so write the final peak and break the expansion:
535 //cout<<"ystart "<<ystart<<" and nybins "<<nybins<<endl;
537 Bool_t truepeak=kTRUE;
539 //cout<<"Maxima found at xbin "<<maxxbin<<" ybin "<<maxybin<<" value "<<maxvalue<<endl;
540 //cout<<"Starting to sum at xbin "<<starts[maxybin-ymin]<<endl;
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++)
547 for(Int_t r=maxybin+1; r<=maxybin+3; r++)
549 right += (Float_t)hist->GetBinContent(hist->GetBin(w,r));
553 for(Int_t w=maxxbin-1; w>=maxxbin-3; w--)
555 for(Int_t r=maxybin+1; r<=maxybin+3; r++)
557 left += (Float_t)hist->GetBinContent(hist->GetBin(w,r));
561 //cout<<"ratio "<<right/left<<endl;
563 Float_t upperratio=1,lowerratio=1;
565 upperratio = right/left;
568 for(Int_t w=maxxbin+1; w<=maxxbin+3; w++)
570 for(Int_t r=maxybin-1; r>=maxybin-3; r--)
572 right += (Float_t)hist->GetBinContent(hist->GetBin(w,r));
576 for(Int_t w=maxxbin-1; w>=maxxbin-3; w--)
578 for(Int_t r=maxybin-1; r>=maxybin-3; r--)
580 left += (Float_t)hist->GetBinContent(hist->GetBin(w,r));
584 //cout<<"ratio "<<left/right<<endl;
587 lowerratio = left/right;
589 if(upperratio > cutratio || lowerratio > cutratio)
595 fXPeaks[fNPeaks] = hist->GetBinCenterX(maxxbin);
596 fYPeaks[fNPeaks] = hist->GetBinCenterY(maxybin);
597 fWeight[fNPeaks] = maxvalue;
601 //Calculate the peak using weigthed means:
604 for(Int_t k=maxybin-1; k<=maxybin+1; k++)
607 for(Int_t l=starts[k]; l<starts[k]+kappawindow; l++)
609 lsum += (Float_t)hist->GetBinContent(hist->GetBin(l,k));
610 sum += (Float_t)hist->GetBinContent(hist->GetBin(l,k));
612 fYPeaks[fNPeaks] += lsum*hist->GetBinCenterY(k);
614 fYPeaks[fNPeaks] /= sum;
616 if(fYPeaks[fNPeaks] < hist->GetBinCenterY(hist->FindYbin(fYPeaks[fNPeaks])))
618 ybin1 = hist->FindYbin(fYPeaks[fNPeaks])-1;
623 ybin1 = hist->FindYbin(fYPeaks[fNPeaks]);
627 Float_t kappa1=0,kappa2=0;
629 for(Int_t k=starts[ybin1]; k<starts[ybin1] + kappawindow; k++)
631 kappa1 += hist->GetBinCenterX(k)*hist->GetBinContent(hist->GetBin(k,ybin1));
632 sum += (Float_t)hist->GetBinContent(hist->GetBin(k,ybin1));
636 for(Int_t k=starts[ybin2]; k<starts[ybin2] + kappawindow; k++)
638 kappa2 += hist->GetBinCenterX(k)*hist->GetBinContent(hist->GetBin(k,ybin2));
639 sum += (Float_t)hist->GetBinContent(hist->GetBin(k,ybin2));
643 fXPeaks[fNPeaks] = ( kappa1*( hist->GetBinCenterY(ybin2) - fYPeaks[fNPeaks] ) +
644 kappa2*( fYPeaks[fNPeaks] - hist->GetBinCenterY(ybin1) ) ) /
645 (hist->GetBinCenterY(ybin2) - hist->GetBinCenterY(ybin1));
654 yl--;//Search continues...
659 for(Int_t i=0; i<hist->GetNbinsY(); i++)
660 delete localmaxima[i];
662 delete [] localmaxima;
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
679 struct AliL3Pre2DPeak
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
692 void AliL3HoughMaxFinder::FindAdaptedRowPeaks(Int_t kappawindow,Int_t xsize,Int_t ysize)
694 // Peak finder which is working over the Hough Space provided by the AliL3HoughTransformerRow class
695 AliL3Histogram *hist = fCurrentHisto;
699 cerr<<"AliL3HoughMaxFinder : No histogram!"<<endl;
703 if(hist->GetNEntries() == 0)
706 Int_t xmin = hist->GetFirstXbin();
707 Int_t xmax = hist->GetLastXbin();
708 Int_t ymin = hist->GetFirstYbin();
709 Int_t ymax = hist->GetLastYbin();
711 //Start by looking for pre-peaks:
713 AliL3PreYPeak **localmaxima = new AliL3PreYPeak*[hist->GetNbinsY()];
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++)
719 localmaxima[ybin-ymin] = new AliL3PreYPeak[hist->GetNbinsX()];
720 nmaxs[ybin-ymin] = 0;
723 for(Int_t xbin=xmin; xbin<=xmax; xbin++)
725 value = hist->GetBinContent(hist->GetBin(xbin,ybin));
728 if(abs(value - lastvalue) > 1)
731 localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fRightValue = value;
734 localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fStartPosition = xbin;
735 localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fEndPosition = xbin;
736 localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMinValue = value;
737 localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMaxValue = value;
738 localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fPrevValue = 0;
739 localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fLeftValue = lastvalue;
742 if(abs(value - lastvalue) <= 1)
744 localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fEndPosition = xbin;
745 if(value>localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMaxValue)
746 localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMaxValue = value;
747 if(value<localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMinValue)
748 localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMinValue = value;
755 localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fRightValue = value;
760 AliL3Pre2DPeak maxima[500];
763 for(Int_t ybin=ymax; ybin >= ymin; ybin--)
765 for(Int_t i=0; i<nmaxs[ybin-ymin]; i++)
767 Int_t localminvalue = localmaxima[ybin-ymin][i].fMinValue;
768 Int_t localmaxvalue = localmaxima[ybin-ymin][i].fMaxValue;
769 Int_t localprevvalue = localmaxima[ybin-ymin][i].fPrevValue;
770 Int_t localnextvalue = 0;
771 Int_t localleftvalue = localmaxima[ybin-ymin][i].fLeftValue;
772 Int_t localrightvalue = localmaxima[ybin-ymin][i].fRightValue;
775 continue; //already used
777 //start expanding in the psi-direction:
779 Int_t localxstart = localmaxima[ybin-ymin][i].fStartPosition;
780 Int_t localxend = localmaxima[ybin-ymin][i].fEndPosition;
781 Int_t tempxstart = localmaxima[ybin-ymin][i].fStartPosition;
782 Int_t tempxend = localmaxima[ybin-ymin][i].fEndPosition;
784 Int_t localy=ybin-1,nybins=1;
786 while(localy >= ymin)
789 for(Int_t j=0; j<nmaxs[localy-ymin]; j++)
791 if( (localmaxima[localy-ymin][j].fStartPosition <= (tempxend + kappawindow)) && (localmaxima[localy-ymin][j].fEndPosition >= (tempxstart - kappawindow)))
793 if(((localmaxima[localy-ymin][j].fMinValue <= localmaxvalue) && (localmaxima[localy-ymin][j].fMinValue >= localminvalue)) ||
794 ((localmaxima[localy-ymin][j].fMaxValue >= localminvalue) && (localmaxima[localy-ymin][j].fMaxValue <= localmaxvalue)))
796 if(localmaxima[localy-ymin][j].fEndPosition > localxend)
797 localxend = localmaxima[localy-ymin][j].fEndPosition;
798 if(localmaxima[localy-ymin][j].fStartPosition < localxstart)
799 localxstart = localmaxima[localy-ymin][j].fStartPosition;
800 tempxstart = localmaxima[localy-ymin][j].fStartPosition;
801 tempxend = localmaxima[localy-ymin][j].fEndPosition;
802 if(localmaxima[localy-ymin][j].fMinValue < localminvalue)
803 localminvalue = localmaxima[localy-ymin][j].fMinValue;
804 if(localmaxima[localy-ymin][j].fMaxValue > localmaxvalue)
805 localmaxvalue = localmaxima[localy-ymin][j].fMaxValue;
806 if(localmaxima[localy-ymin][j].fRightValue > localrightvalue)
807 localrightvalue = localmaxima[localy-ymin][j].fRightValue;
808 if(localmaxima[localy-ymin][j].fLeftValue > localleftvalue)
809 localleftvalue = localmaxima[localy-ymin][j].fLeftValue;
810 localmaxima[localy-ymin][j].fMinValue = -1;
817 if(localmaxvalue > localmaxima[localy-ymin][j].fPrevValue)
818 localmaxima[localy-ymin][j].fPrevValue = localmaxvalue;
819 if(localmaxima[localy-ymin][j].fMaxValue > localnextvalue)
820 localnextvalue = localmaxima[localy-ymin][j].fMaxValue;
824 if(!found || localy == ymin)//no more local maximas to be matched, so write the final peak and break the expansion:
826 if((nybins > ysize) && ((localxend-localxstart+1) > xsize) && (localmaxvalue > localprevvalue) && (localmaxvalue > localnextvalue) && (localmaxvalue > localleftvalue) && (localmaxvalue > localrightvalue))
827 // if((nybins > ysize) && ((localxend-localxstart+1) > xsize))
829 maxima[nmaxima].fX = ((Float_t)localxstart+(Float_t)localxend)/2.0;
830 maxima[nmaxima].fY = ((Float_t)ybin+(Float_t)(localy+1))/2.0;
831 maxima[nmaxima].fSizeX = localxend-localxstart+1;
832 maxima[nmaxima].fSizeY = nybins;
833 maxima[nmaxima].fWeight = (localminvalue+localmaxvalue)/2;
834 maxima[nmaxima].fStartX = localxstart;
835 maxima[nmaxima].fEndX = localxend;
836 maxima[nmaxima].fStartY = localy +1;
837 maxima[nmaxima].fEndY = ybin;
839 // 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;
846 localy--;//Search continues...
852 Int_t nxbins = hist->GetNbinsX()+2;
854 for(Int_t i = 0; i < (nmaxima - 1); i++)
856 if(maxima[i].fWeight < 0) continue;
857 for(Int_t j = i + 1; j < nmaxima; j++)
859 if(maxima[j].fWeight < 0) continue;
860 Int_t xtrack1=0,xtrack2=0,ytrack1=0,ytrack2=0;
862 for(Int_t ix1 = maxima[i].fStartX; ix1 <= maxima[i].fEndX; ix1++) {
863 for(Int_t ix2 = maxima[j].fStartX; ix2 <= maxima[j].fEndX; ix2++) {
864 if(abs(ix1 - ix2) < deltax) {
865 deltax = abs(ix1 - ix2);
872 for(Int_t iy1 = maxima[i].fStartY; iy1 <= maxima[i].fEndY; iy1++) {
873 for(Int_t iy2 = maxima[j].fStartY; iy2 <= maxima[j].fEndY; iy2++) {
874 if(abs(iy1 - iy2) < deltay) {
875 deltay = abs(iy1 - iy2);
881 Int_t firstrow1 = AliL3HoughTransformerRow::GetTrackFirstRow()[xtrack1 + nxbins*ytrack1];
882 Int_t lastrow1 = AliL3HoughTransformerRow::GetTrackLastRow()[xtrack1 + nxbins*ytrack1];
883 Int_t firstrow2 = AliL3HoughTransformerRow::GetTrackFirstRow()[xtrack1 + nxbins*ytrack1];
884 Int_t lastrow2 = AliL3HoughTransformerRow::GetTrackLastRow()[xtrack1 + nxbins*ytrack1];
885 Int_t firstrow,lastrow;
886 if(firstrow1 < firstrow2)
887 firstrow = firstrow2;
889 firstrow = firstrow1;
891 if(lastrow1 > lastrow2)
896 AliL3HoughTrack track1;
897 Float_t x1 = hist->GetPreciseBinCenterX(xtrack1);
898 Float_t y1 = hist->GetPreciseBinCenterY(ytrack1);
899 Float_t psi1 = atan((x1-y1)/(AliL3HoughTransformerRow::GetBeta1()-AliL3HoughTransformerRow::GetBeta2()));
900 Float_t kappa1 = 2.0*(x1*cos(psi1)-AliL3HoughTransformerRow::GetBeta1()*sin(psi1));
901 track1.SetTrackParameters(kappa1,psi1,1);
902 Float_t firsthit1[3];
903 if(!track1.GetCrossingPoint(firstrow,firsthit1)) continue;
905 if(!track1.GetCrossingPoint(lastrow,lasthit1)) continue;
907 AliL3HoughTrack track2;
908 Float_t x2 = hist->GetPreciseBinCenterX(xtrack2);
909 Float_t y2 = hist->GetPreciseBinCenterY(ytrack2);
910 Float_t psi2 = atan((x2-y2)/(AliL3HoughTransformerRow::GetBeta1()-AliL3HoughTransformerRow::GetBeta2()));
911 Float_t kappa2 = 2.0*(x2*cos(psi2)-AliL3HoughTransformerRow::GetBeta1()*sin(psi2));
912 track2.SetTrackParameters(kappa2,psi2,1);
913 Float_t firsthit2[3];
914 if(!track2.GetCrossingPoint(firstrow,firsthit2)) continue;
916 if(!track2.GetCrossingPoint(lastrow,lasthit2)) continue;
918 Float_t padpitchlow = AliL3Transform::GetPadPitchWidth(AliL3Transform::GetPatch(firstrow));
919 Float_t padpitchup = AliL3Transform::GetPadPitchWidth(AliL3Transform::GetPatch(lastrow));
920 // check the distance between tracks at the edges
922 // cout<<"DEBUG Merge peaks "<<i<<" "<<j<<" "<<firsthit1[1]<<" "<<firsthit2[1]<<" "<<lasthit1[1]<<" "<<lasthit2[1]<<endl;
924 //cvetan please check!!! I added a cast to Int_t
925 if((fabs(firsthit1[1]-firsthit2[1]) < 3.0*padpitchlow) && (fabs(lasthit1[1]-lasthit2[1]) < 3.0*padpitchup)) {
926 if(maxima[i].fSizeX*maxima[i].fSizeY > maxima[j].fSizeX*maxima[j].fSizeY)
927 maxima[j].fWeight = -maxima[j].fWeight;
928 if(maxima[i].fSizeX*maxima[i].fSizeY < maxima[j].fSizeX*maxima[j].fSizeY)
929 maxima[i].fWeight = -maxima[i].fWeight;
931 // cout<<"Merge peaks "<<i<<" "<<j<<" "<<maxima[i].fWeight<<" "<<maxima[j].fWeight<<endl;
937 //merge tracks in neighbour eta slices
939 for(Int_t i = 0; i < nmaxima; i++) {
940 if(maxima[i].fWeight > 0) {
941 fXPeaks[fNPeaks] = hist->GetPreciseBinCenterX(maxima[i].fX);
942 fYPeaks[fNPeaks] = hist->GetPreciseBinCenterY(maxima[i].fY);
943 fWeight[fNPeaks] = (Int_t)maxima[i].fWeight;
945 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;
952 Int_t currentnpeaks = fNPeaks;
953 for(Int_t i = 0; i < nmaxima; i++) {
954 if(maxima[i].fWeight < 0) continue;
955 Bool_t merged = kFALSE;
956 for(Int_t j = fN1PeaksPrevEtaSlice; j < fN2PeaksPrevEtaSlice; j++) {
957 if(fWeight[j] < 0) continue;
958 if((fENDETAPeaks[j]-fSTARTETAPeaks[j]) >= 1) continue;
959 if((maxima[i].fStartX <= fENDXPeaks[j]+1) && (maxima[i].fEndX >= fSTARTXPeaks[j]-1)) {
960 if((maxima[i].fStartY <= fENDYPeaks[j]+1) && (maxima[i].fEndY >= fSTARTYPeaks[j]-1)) {
963 fXPeaks[fNPeaks] = (hist->GetPreciseBinCenterX(maxima[i].fX)+(fENDETAPeaks[j]-fSTARTETAPeaks[j]+1)*fXPeaks[j])/(fENDETAPeaks[j]-fSTARTETAPeaks[j]+2);
964 fYPeaks[fNPeaks] = (hist->GetPreciseBinCenterY(maxima[i].fY)+(fENDETAPeaks[j]-fSTARTETAPeaks[j]+1)*fYPeaks[j])/(fENDETAPeaks[j]-fSTARTETAPeaks[j]+2);
965 fWeight[fNPeaks] = (Int_t)maxima[i].fWeight + fWeight[j];
966 fSTARTXPeaks[fNPeaks] = maxima[i].fStartX;
967 fSTARTYPeaks[fNPeaks] = maxima[i].fStartY;
968 fENDXPeaks[fNPeaks] = maxima[i].fEndX;
969 fENDYPeaks[fNPeaks] = maxima[i].fEndY;
970 fSTARTETAPeaks[fNPeaks] = fSTARTETAPeaks[j];
971 fENDETAPeaks[fNPeaks] = fCurrentEtaSlice;
973 fWeight[j] = -fWeight[j];
977 fXPeaks[fNPeaks] = hist->GetPreciseBinCenterX(maxima[i].fX);
978 fYPeaks[fNPeaks] = hist->GetPreciseBinCenterY(maxima[i].fY);
980 fWeight[fNPeaks] = (Int_t)maxima[i].fWeight;
982 fWeight[fNPeaks] = -(Int_t)maxima[i].fWeight;
983 fSTARTXPeaks[fNPeaks] = maxima[i].fStartX;
984 fSTARTYPeaks[fNPeaks] = maxima[i].fStartY;
985 fENDXPeaks[fNPeaks] = maxima[i].fEndX;
986 fENDYPeaks[fNPeaks] = maxima[i].fEndY;
987 fSTARTETAPeaks[fNPeaks] = fCurrentEtaSlice;
988 fENDETAPeaks[fNPeaks] = fCurrentEtaSlice;
991 fN1PeaksPrevEtaSlice = currentnpeaks;
992 fN2PeaksPrevEtaSlice = fNPeaks;
994 for(Int_t i=0; i<hist->GetNbinsY(); i++)
995 delete localmaxima[i];
997 delete [] localmaxima;
1001 void AliL3HoughMaxFinder::FindPeak1(Int_t ywindow,Int_t xbinsides)
1003 //Testing mutliple peakfinding.
1004 //The algorithm searches the histogram for prepreaks by looking in windows
1005 //for each bin on the xaxis. The size of these windows is controlled by ywindow.
1006 //Then the prepreaks are sorted according to their weight (sum inside window),
1007 //and the peak positions are calculated by taking the weighted mean in both
1008 //x and y direction. The size of the peak in x-direction is controlled by xbinsides.
1012 printf("AliL3HoughMaxFinder::FindPeak1 : No input histogram\n");
1015 if(fCurrentHisto->GetNEntries()==0)
1019 //Int_t xbinsides=1;
1021 //Float_t max_kappa = 0.001;
1022 //Float_t max_phi0 = 0.08;
1026 Int_t xmin = fCurrentHisto->GetFirstXbin();
1027 Int_t xmax = fCurrentHisto->GetLastXbin();
1028 Int_t ymin = fCurrentHisto->GetFirstYbin();
1029 Int_t ymax = fCurrentHisto->GetLastYbin();
1030 Int_t nbinsx = fCurrentHisto->GetNbinsX()+1;
1032 AliL3AxisWindow **windowPt = new AliL3AxisWindow*[nbinsx];
1033 AliL3AxisWindow **anotherPt = new AliL3AxisWindow*[nbinsx];
1035 for(Int_t i=0; i<nbinsx; i++)
1037 windowPt[i] = new AliL3AxisWindow;
1038 #if defined(__DECCXX)
1039 bzero((char *)windowPt[i],sizeof(AliL3AxisWindow));
1041 bzero((void*)windowPt[i],sizeof(AliL3AxisWindow));
1043 anotherPt[i] = windowPt[i];
1046 for(Int_t xbin=xmin; xbin<=xmax; xbin++)
1049 for(Int_t ybin=ymin; ybin<=ymax-ywindow; ybin++)
1051 Int_t suminwindow=0;
1052 for(Int_t b=ybin; b<ybin+ywindow; b++)
1055 Int_t bin = fCurrentHisto->GetBin(xbin,b);
1056 suminwindow += (Int_t)fCurrentHisto->GetBinContent(bin);
1059 if(suminwindow > maxsum)
1061 maxsum = suminwindow;
1062 windowPt[xbin]->fYmin = ybin;
1063 windowPt[xbin]->fYmax = ybin + ywindow;
1064 windowPt[xbin]->fWeight = suminwindow;
1065 windowPt[xbin]->fXbin = xbin;
1070 //Sort the windows according to the weight
1071 SortPeaks(windowPt,0,nbinsx);
1074 for(Int_t i=0; i<nbinsx; i++)
1077 Int_t xbin = windowPt[i]->fXbin;
1079 if(xbin<xmin || xbin > xmax-1) continue;
1081 //Check if this is really a local maxima
1082 if(anotherPt[xbin-1]->fWeight > anotherPt[xbin]->fWeight ||
1083 anotherPt[xbin+1]->fWeight > anotherPt[xbin]->fWeight)
1086 for(Int_t j=windowPt[i]->fYmin; j<windowPt[i]->fYmax; j++)
1088 //Calculate the mean in y direction:
1089 Int_t bin = fCurrentHisto->GetBin(windowPt[i]->fXbin,j);
1090 top += (fCurrentHisto->GetBinCenterY(j))*(fCurrentHisto->GetBinContent(bin));
1091 butt += fCurrentHisto->GetBinContent(bin);
1094 if(butt < fThreshold)
1097 fXPeaks[fNPeaks] = fCurrentHisto->GetBinCenterX(windowPt[i]->fXbin);
1098 fYPeaks[fNPeaks] = top/butt;
1099 fWeight[fNPeaks] = (Int_t)butt;
1100 //cout<<"mean in y "<<ypeaks[n]<<" on x "<<windowPt[i]->xbin<<" content "<<butt<<endl;
1104 cerr<<"AliL3HoughMaxFinder::FindPeak1 : Peak array out of range!!!"<<endl;
1110 //Improve the peaks by including the region around in x.
1114 for(Int_t i=0; i<fNPeaks; i++)
1116 Int_t xbin = fCurrentHisto->FindXbin(fXPeaks[i]);
1117 if(xbin - xbinsides < xmin || xbin + xbinsides > xmax) continue;
1121 prev = xbin - xbinsides+1;
1122 for(Int_t j=xbin-xbinsides; j<=xbin+xbinsides; j++)
1125 //Check if the windows are overlapping:
1126 if(anotherPt[j]->ymin > anotherPt[prev]->ymax) {prev=j; continue;}
1127 if(anotherPt[j]->ymax < anotherPt[prev]->ymin) {prev=j; continue;}
1131 top += fCurrentHisto->GetBinCenterX(j)*anotherPt[j]->fWeight;
1132 butt += anotherPt[j]->fWeight;
1134 for(Int_t k=anotherPt[j]->fYmin; k<anotherPt[j]->fYmax; k++)
1136 Int_t bin = fCurrentHisto->GetBin(j,k);
1137 ytop += (fCurrentHisto->GetBinCenterY(k))*(fCurrentHisto->GetBinContent(bin));
1138 ybutt += fCurrentHisto->GetBinContent(bin);
1139 w+=(Int_t)fCurrentHisto->GetBinContent(bin);
1143 fXPeaks[i] = top/butt;
1144 fYPeaks[i] = ytop/ybutt;
1146 //cout<<"Setting weight "<<w<<" kappa "<<fXPeaks[i]<<" phi0 "<<fYPeaks[i]<<endl;
1149 //Check if this peak is overlapping with a previous:
1150 for(Int_t p=0; p<i; p++)
1152 //cout<<fabs(fXPeaks[p] - fXPeaks[i])<<" "<<fabs(fYPeaks[p] - fYPeaks[i])<<endl;
1153 if(fabs(fXPeaks[p] - fXPeaks[i]) < max_kappa &&
1154 fabs(fYPeaks[p] - fYPeaks[i]) < max_phi0)
1163 for(Int_t i=0; i<nbinsx; i++)
1166 delete [] anotherPt;
1169 void AliL3HoughMaxFinder::SortPeaks(struct AliL3AxisWindow **a,Int_t first,Int_t last)
1171 //General sorting routine
1172 //Sort according to PeakCompare()
1174 static struct AliL3AxisWindow *tmp;
1175 static int i; // "static" to save stack space
1178 while (last - first > 1) {
1182 while (++i < last && PeakCompare(a[i], a[first]) < 0)
1184 while (--j > first && PeakCompare(a[j], a[first]) > 0)
1200 if (j - first < last - (j + 1)) {
1201 SortPeaks(a, first, j);
1202 first = j + 1; // QSort(j + 1, last);
1204 SortPeaks(a, j + 1, last);
1205 last = j; // QSort(first, j);
1211 Int_t AliL3HoughMaxFinder::PeakCompare(struct AliL3AxisWindow *a,struct AliL3AxisWindow *b) const
1213 // Peak comparison based on peaks weight
1214 if(a->fWeight < b->fWeight) return 1;
1215 if(a->fWeight > b->fWeight) return -1;
1219 void AliL3HoughMaxFinder::FindPeak(Int_t t1,Double_t t2,Int_t t3)
1221 //Attempt of a more sophisticated peak finder.
1222 //Finds the best peak in the histogram, and returns the corresponding
1227 printf("AliL3HoughMaxFinder::FindPeak : No histogram!!\n");
1230 AliL3Histogram *hist = fCurrentHisto;
1231 if(hist->GetNEntries()==0)
1234 Int_t xmin = hist->GetFirstXbin();
1235 Int_t xmax = hist->GetLastXbin();
1236 Int_t ymin = hist->GetFirstYbin();
1237 Int_t ymax = hist->GetLastYbin();
1238 Int_t nbinsx = hist->GetNbinsX()+1;
1240 Int_t *m = new Int_t[nbinsx];
1241 Int_t *mlow = new Int_t[nbinsx];
1242 Int_t *mup = new Int_t[nbinsx];
1245 recompute: //this is a goto.
1247 for(Int_t i=0; i<nbinsx; i++)
1254 Int_t maxx=0,sum=0,maxxbin=0,bin=0;
1256 for(Int_t xbin=xmin; xbin<=xmax; xbin++)
1258 for(Int_t ybin=ymin; ybin <= ymax - t1; ybin++)
1261 for(Int_t y=ybin; y <= ybin+t1; y++)
1265 bin = hist->GetBin(xbin,y);
1266 sum += (Int_t)hist->GetBinContent(bin);
1269 if(sum > m[xbin]) //Max value locally in this xbin
1273 mup[xbin]=ybin + t1;
1278 if(m[xbin] > maxx) //Max value globally in x-direction
1281 maxx = m[xbin];//sum;
1284 //printf("maxxbin %d maxx %d mlow %d mup %d\n",maxxbin,maxx,mlow[maxxbin],mup[maxxbin]);
1285 //printf("ylow %f yup %f\n",hist->GetBinCenterY(mlow[maxxbin]),hist->GetBinCenterY(mup[maxxbin]));
1287 //Determine a width in the x-direction
1290 for(Int_t xbin=maxxbin-1; xbin >= xmin; xbin--)
1292 if(m[xbin] < maxx*t2)
1298 for(Int_t xbin = maxxbin+1; xbin <=xmax; xbin++)
1300 if(m[xbin] < maxx*t2)
1307 Double_t top=0,butt=0,value,xpeak;
1308 if(xup - xlow + 1 > t3)
1311 printf("\nxrange out if limit xup %d xlow %d t1 %d\n\n",xlow,xup,t1);
1316 xpeak = hist->GetBinCenterX(maxxbin);
1321 //printf("xlow %f xup %f\n",hist->GetBinCenterX(xlow),hist->GetBinCenterX(xup));
1322 //printf("Spread in x %d\n",xup-xlow +1);
1324 //Now, calculate the center of mass in x-direction
1325 for(Int_t xbin=xlow; xbin <= xup; xbin++)
1327 value = hist->GetBinCenterX(xbin);
1328 top += value*m[xbin];
1335 //Find the peak in y direction:
1336 Int_t xl = hist->FindXbin(xpeak);
1337 if(hist->GetBinCenterX(xl) > xpeak)
1342 if(hist->GetBinCenterX(xl) > xpeak || hist->GetBinCenterX(xu) <= xpeak)
1343 printf("\nAliL3HoughMaxFinder::FindPeak : Wrong xrange %f %f %f\n\n",hist->GetBinCenterX(xl),xpeak,hist->GetBinCenterX(xu));
1345 //printf("\nxlow %f xup %f\n",hist->GetBinCenterX(xl),hist->GetBinCenterX(xu));
1349 //printf("ylow %f yup %f\n",hist->GetBinCenterY(mlow[xl]),hist->GetBinCenterY(mup[xl]));
1350 //printf("ylow %f yup %f\n",hist->GetBinCenterY(mlow[xu]),hist->GetBinCenterY(mup[xu]));
1352 for(Int_t ybin=mlow[xl]; ybin <= mup[xl]; ybin++)
1354 value = hist->GetBinCenterY(ybin);
1355 bin = hist->GetBin(xl,ybin);
1356 top += value*hist->GetBinContent(bin);
1357 butt += hist->GetBinContent(bin);
1359 Double_t ypeaklow = top/butt;
1361 //printf("ypeaklow %f\n",ypeaklow);
1364 for(Int_t ybin=mlow[xu]; ybin <= mup[xu]; ybin++)
1366 value = hist->GetBinCenterY(ybin);
1367 bin = hist->GetBin(xu,ybin);
1368 top += value*hist->GetBinContent(bin);
1369 butt += hist->GetBinContent(bin);
1371 Double_t ypeakup = top/butt;
1373 //printf("ypeakup %f\n",ypeakup);
1375 Double_t xvalueup = hist->GetBinCenterX(xu);
1376 Double_t xvaluelow = hist->GetBinCenterX(xl);
1378 Double_t ypeak = (ypeaklow*(xvalueup - xpeak) + ypeakup*(xpeak - xvaluelow))/(xvalueup - xvaluelow);
1382 //bin = hist->FindBin(xpeak,ypeak);
1383 //Int_t weight = (Int_t)hist->GetBinContent(bin);
1385 //AliL3HoughTrack *track = new AliL3HoughTrack();
1386 //track->SetTrackParameters(xpeak,ypeak,weight);
1387 fXPeaks[fNPeaks]=xpeak;
1388 fYPeaks[fNPeaks]=ypeak;
1389 fWeight[fNPeaks]=(Int_t)hist->GetBinContent(bin);
1399 Float_t AliL3HoughMaxFinder::GetXPeakSize(Int_t i) const
1401 // Get X size of a peak
1404 STDCERR<<"AliL3HoughMaxFinder::GetXPeakSize : Invalid index "<<i<<STDENDL;
1407 Float_t binwidth = fCurrentHisto->GetBinWidthX();
1408 return binwidth*(fENDXPeaks[i]-fSTARTXPeaks[i]+1);
1411 Float_t AliL3HoughMaxFinder::GetYPeakSize(Int_t i) const
1413 // Get Y size of a peak
1416 STDCERR<<"AliL3HoughMaxFinder::GetYPeak : Invalid index "<<i<<STDENDL;
1419 Float_t binwidth = fCurrentHisto->GetBinWidthY();
1420 return binwidth*(fENDYPeaks[i]-fSTARTYPeaks[i]+1);