-//Author: Anders Strand Vestbo
-//Last Modified: 28.6.01
+//$Id$
+
+// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
+//*-- Copyright © ASV
-#include <string.h>
#include <math.h>
#include <stdlib.h>
+#include <string.h>
+
+#include <TNtuple.h>
+#include <TFile.h>
#include "AliL3Histogram.h"
#include "AliL3TrackArray.h"
#include "AliL3HoughTrack.h"
#include "AliL3HoughMaxFinder.h"
+//_____________________________________________________________
+// AliL3HoughMaxFinder
+//
+// Maximum finder
+
ClassImp(AliL3HoughMaxFinder)
{
//Default constructor
fThreshold = 0;
- //fTracks = 0;
fHistoType=0;
+ fXPeaks=0;
+ fYPeaks=0;
+ fNPeaks=0;
+ fNMax=0;
+ fNtuppel = 0;
}
-AliL3HoughMaxFinder::AliL3HoughMaxFinder(Char_t *histotype,AliL3Histogram *hist)
+AliL3HoughMaxFinder::AliL3HoughMaxFinder(Char_t *histotype,Int_t nmax,AliL3Histogram *hist)
{
//Constructor
if(hist)
fCurrentHisto = hist;
+
+ fNMax=nmax;
+ fXPeaks = new Float_t[fNMax];
+ fYPeaks = new Float_t[fNMax];
+ fWeight = new Int_t[fNMax];
+ fNtuppel = 0;
+ fThreshold=0;
}
AliL3HoughMaxFinder::~AliL3HoughMaxFinder()
{
//Destructor
+ if(fXPeaks)
+ delete [] fXPeaks;
+ if(fYPeaks)
+ delete [] fYPeaks;
+ if(fWeight)
+ delete [] fWeight;
+ if(fNtuppel)
+ delete fNtuppel;
+}
+
+void AliL3HoughMaxFinder::Reset()
+{
+ for(Int_t i=0; i<fNMax; i++)
+ {
+ fXPeaks[i]=0;
+ fYPeaks[i]=0;
+ fWeight[i]=0;
+ }
+ fNPeaks=0;
+}
+void AliL3HoughMaxFinder::CreateNtuppel()
+{
+ fNtuppel = new TNtuple("ntuppel","Peak charateristics","kappa:phi0:weigth:content3:content5:content1:content7");
+ //content#; neighbouring bins of the peak.
+
}
-void AliL3HoughMaxFinder::FindAbsMaxima(Int_t &max_xbin,Int_t &max_ybin)
+void AliL3HoughMaxFinder::WriteNtuppel(Char_t *filename)
{
+ TFile *file = TFile::Open(filename,"RECREATE");
+ if(!file)
+ {
+ cerr<<"AliL3HoughMaxFinder::WriteNtuppel : Error opening file "<<filename<<endl;
+ return;
+ }
+ fNtuppel->Write();
+ file->Close();
+}
+
+void AliL3HoughMaxFinder::FindAbsMaxima()
+{
+
if(!fCurrentHisto)
{
- printf("AliL3HoughMaxFinder::FindAbsMaxima : No histogram!\n");
+ cerr<<"AliL3HoughMaxFinder::FindAbsMaxima : No histogram"<<endl;
return;
}
AliL3Histogram *hist = fCurrentHisto;
Int_t ymin = hist->GetFirstYbin();
Int_t ymax = hist->GetLastYbin();
Int_t bin;
- Stat_t value,max_value=0;
-
+ Double_t value,max_value=0;
+
+ Int_t max_xbin=0,max_ybin=0;
for(Int_t xbin=xmin; xbin<=xmax; xbin++)
{
for(Int_t ybin=ymin; ybin<=ymax; ybin++)
}
}
+ if(fNPeaks > fNMax)
+ {
+ cerr<<"AliL3HoughMaxFinder::FindAbsMaxima : Array out of range : "<<fNPeaks<<endl;
+ return;
+ }
+
+
+ Double_t max_x = hist->GetBinCenterX(max_xbin);
+ Double_t max_y = hist->GetBinCenterY(max_ybin);
+ fXPeaks[fNPeaks] = max_x;
+ fYPeaks[fNPeaks] = max_y;
+ fWeight[fNPeaks] = (Int_t)max_value;
+ fNPeaks++;
+
+ if(fNtuppel)
+ {
+ Int_t bin3 = hist->GetBin(max_xbin-1,max_ybin);
+ Int_t bin5 = hist->GetBin(max_xbin+1,max_ybin);
+ Int_t bin1 = hist->GetBin(max_xbin,max_ybin-1);
+ Int_t bin7 = hist->GetBin(max_xbin,max_ybin+1);
+
+ fNtuppel->Fill(max_x,max_y,max_value,hist->GetBinContent(bin3),hist->GetBinContent(bin5),hist->GetBinContent(bin1),hist->GetBinContent(bin7));
+ }
+
}
-AliL3TrackArray *AliL3HoughMaxFinder::FindBigMaxima(AliL3Histogram *hist)
+void AliL3HoughMaxFinder::FindBigMaxima()
{
+ AliL3Histogram *hist = fCurrentHisto;
Int_t xmin = hist->GetFirstXbin();
Int_t xmax = hist->GetLastXbin();
Int_t ymin = hist->GetFirstYbin();
Int_t ymax = hist->GetLastYbin();
Int_t bin[25],bin_index;
- Stat_t value[25];
-
- AliL3TrackArray *tracks = new AliL3TrackArray("AliL3HoughTrack");
- AliL3HoughTrack *track;
+ Double_t value[25];
for(Int_t xbin=xmin+2; xbin<xmax-3; xbin++)
{
{
if(value[b] > value[12] || b==bin_index) break;
b++;
- printf("b %d\n",b);
+ //printf("b %d\n",b);
}
if(b == bin_index)
{
//Found maxima
+ if(fNPeaks > fNMax)
+ {
+ cerr<<"AliL3HoughMaxFinder::FindBigMaxima : Array out of range "<<fNPeaks<<endl;
+ return;
+ }
+
Double_t max_x = hist->GetBinCenterX(xbin);
Double_t max_y = hist->GetBinCenterY(ybin);
- track = (AliL3HoughTrack*)tracks->NextTrack();
- track->SetTrackParameters(max_x,max_y,(Int_t)value[12]);
+ fXPeaks[fNPeaks] = max_x;
+ fYPeaks[fNPeaks] = max_y;
+ fNPeaks++;
}
}
}
-
- tracks->QSort();
- return tracks;
}
-AliL3TrackArray *AliL3HoughMaxFinder::FindMaxima(AliL3Histogram *hist,Int_t *rowrange,Int_t ref_row)
+void AliL3HoughMaxFinder::FindMaxima(Double_t grad_x,Double_t grad_y)
{
//Locate all the maxima in input histogram.
//Maxima is defined as bins with more entries than the
//immediately neighbouring bins.
- Int_t xmin = hist->GetFirstXbin();
- Int_t xmax = hist->GetLastXbin();
- Int_t ymin = hist->GetFirstYbin();
- Int_t ymax = hist->GetLastYbin();
- Int_t bin[9],track_counter=0;
- Stat_t value[9];
+ Int_t xmin = fCurrentHisto->GetFirstXbin();
+ Int_t xmax = fCurrentHisto->GetLastXbin();
+ Int_t ymin = fCurrentHisto->GetFirstYbin();
+ Int_t ymax = fCurrentHisto->GetLastYbin();
+ Int_t bin[9];
+ Double_t value[9];
- AliL3TrackArray *tracks = new AliL3TrackArray("AliL3HoughTrack");
- AliL3HoughTrack *track;
-
for(Int_t xbin=xmin+1; xbin<xmax-1; xbin++)
{
for(Int_t ybin=ymin+1; ybin<ymax-1; ybin++)
{
- bin[0] = hist->GetBin(xbin-1,ybin-1);
- bin[1] = hist->GetBin(xbin,ybin-1);
- bin[2] = hist->GetBin(xbin+1,ybin-1);
- bin[3] = hist->GetBin(xbin-1,ybin);
- bin[4] = hist->GetBin(xbin,ybin);
- bin[5] = hist->GetBin(xbin+1,ybin);
- bin[6] = hist->GetBin(xbin-1,ybin+1);
- bin[7] = hist->GetBin(xbin,ybin+1);
- bin[8] = hist->GetBin(xbin+1,ybin+1);
- value[0] = hist->GetBinContent(bin[0]);
- value[1] = hist->GetBinContent(bin[1]);
- value[2] = hist->GetBinContent(bin[2]);
- value[3] = hist->GetBinContent(bin[3]);
- value[4] = hist->GetBinContent(bin[4]);
- value[5] = hist->GetBinContent(bin[5]);
- value[6] = hist->GetBinContent(bin[6]);
- value[7] = hist->GetBinContent(bin[7]);
- value[8] = hist->GetBinContent(bin[8]);
+ bin[0] = fCurrentHisto->GetBin(xbin-1,ybin-1);
+ bin[1] = fCurrentHisto->GetBin(xbin,ybin-1);
+ bin[2] = fCurrentHisto->GetBin(xbin+1,ybin-1);
+ bin[3] = fCurrentHisto->GetBin(xbin-1,ybin);
+ bin[4] = fCurrentHisto->GetBin(xbin,ybin);
+ bin[5] = fCurrentHisto->GetBin(xbin+1,ybin);
+ bin[6] = fCurrentHisto->GetBin(xbin-1,ybin+1);
+ bin[7] = fCurrentHisto->GetBin(xbin,ybin+1);
+ bin[8] = fCurrentHisto->GetBin(xbin+1,ybin+1);
+ value[0] = fCurrentHisto->GetBinContent(bin[0]);
+ value[1] = fCurrentHisto->GetBinContent(bin[1]);
+ value[2] = fCurrentHisto->GetBinContent(bin[2]);
+ value[3] = fCurrentHisto->GetBinContent(bin[3]);
+ value[4] = fCurrentHisto->GetBinContent(bin[4]);
+ value[5] = fCurrentHisto->GetBinContent(bin[5]);
+ value[6] = fCurrentHisto->GetBinContent(bin[6]);
+ value[7] = fCurrentHisto->GetBinContent(bin[7]);
+ value[8] = fCurrentHisto->GetBinContent(bin[8]);
+
- if(value[4] <= fThreshold) continue;//central bin below threshold
if(value[4]>value[0] && value[4]>value[1] && value[4]>value[2]
&& value[4]>value[3] && value[4]>value[5] && value[4]>value[6]
&& value[4]>value[7] && value[4]>value[8])
{
//Found a local maxima
- Float_t max_x = hist->GetBinCenterX(xbin);
- Float_t max_y = hist->GetBinCenterY(ybin);
-
- track = (AliL3HoughTrack*)tracks->NextTrack();
+ Float_t max_x = fCurrentHisto->GetBinCenterX(xbin);
+ Float_t max_y = fCurrentHisto->GetBinCenterY(ybin);
+ cout<<"Checking for threshols "<<value[4]<<" "<<fThreshold<<endl;
+ if((Int_t)value[4] <= fThreshold) continue;//central bin below threshold
- switch(fHistoType)
+ if(fNPeaks > fNMax)
{
- case 'c':
- track->SetTrackParameters(max_x,max_y,(Int_t)value[4]);
- break;
- case 'l':
- track->SetLineParameters(max_x,max_y,(Int_t)value[4],rowrange,ref_row);
- break;
- default:
- printf("AliL3HoughMaxFinder: Error in tracktype\n");
+ cerr<<"AliL3HoughMaxFinder::FindMaxima : Array out of range "<<fNPeaks<<endl;
+ return;
}
- track_counter++;
+ //Check the gradient:
+ if(value[4]/value[3] < grad_x || value[4]/value[5] < grad_x ||
+ value[4]/value[1] < grad_y || value[4]/value[7] < grad_y)
+ continue;
-
+ fXPeaks[fNPeaks] = max_x;
+ fYPeaks[fNPeaks] = max_y;
+ fWeight[fNPeaks] = (Int_t)value[4];
+ fNPeaks++;
+
+ /*
+ //Check if the peak is overlapping with a previous:
+ Bool_t bigger = kFALSE;
+ for(Int_t p=0; p<entries; p++)
+ {
+ if(fabs(max_x - xpeaks[p]) < kappa_overlap && fabs(max_y - ypeaks[p]) < phi_overlap)
+ {
+ bigger = kTRUE;
+ if(value[4] > weight[p]) //this peak is bigger.
+ {
+ xpeaks[p] = max_x;
+ ypeaks[p] = max_y;
+ weight[p] = (Int_t)value[4];
+ }
+ else
+ continue; //previous peak is bigger.
+ }
+ }
+ if(!bigger) //there were no overlapping peaks.
+ {
+ xpeaks[entries] = max_x;
+ ypeaks[entries] = max_y;
+ weight[entries] = (Int_t)value[4];
+ entries++;
+ }
+ */
}
else
continue; //not a maxima
}
}
- tracks->QSort();
- return tracks;
+
}
-void AliL3HoughMaxFinder::FindPeak1(Float_t *xpeaks,Float_t *ypeaks,Int_t &n)
+void AliL3HoughMaxFinder::FindPeak1(Float_t *xpeaks,Float_t *ypeaks,Int_t *weight,Int_t &n,
+ Int_t y_window,Int_t x_bin_sides)
{
- //testing mutliple peakfinding
-
+ //Testing mutliple peakfinding.
+ //The algorithm searches the histogram for prepreaks by looking in windows
+ //for each bin on the xaxis. The size of these windows is controlled by y_window.
+ //Then the prepreaks are sorted according to their weight (sum inside window),
+ //and the peak positions are calculated by taking the weighted mean in both
+ //x and y direction. The size of the peak in x-direction is controlled by x_bin_sides.
+
Int_t max_entries = n;
n=0;
if(!fCurrentHisto)
printf("AliL3HoughMaxFinder::FindPeak1 : No input histogram\n");
return;
}
- Int_t y_window=2;
- Int_t x_bin_sides=2;
+ //Int_t y_window=2;
+ //Int_t x_bin_sides=1;
+
+ Float_t max_kappa = 0.001;
+ Float_t max_phi0 = 0.05;
+
Int_t max_sum=0;
Int_t xmin = fCurrentHisto->GetFirstXbin();
//Sort the windows according to the weight
SortPeaks(windowPt,0,nbinsx);
- //for(Int_t i=0; i<nbinsx; i++)
- //printf("xbin %f weight %d\n",fCurrentHisto->GetBinCenterX(windowPt[i]->xbin),windowPt[i]->weight);
-
Float_t top,butt;
for(Int_t i=0; i<nbinsx; i++)
{
if(xbin<xmin || xbin > xmax-1) continue;
- //Check if this is really a peak
- if(anotherPt[xbin-1]->weight > windowPt[i]->weight ||
- anotherPt[xbin+1]->weight > windowPt[i]->weight)
+ //Check if this is really a local maxima
+ if(anotherPt[xbin-1]->weight > anotherPt[xbin]->weight ||
+ anotherPt[xbin+1]->weight > anotherPt[xbin]->weight)
continue;
for(Int_t j=windowPt[i]->ymin; j<windowPt[i]->ymax; j++)
}
xpeaks[n] = fCurrentHisto->GetBinCenterX(windowPt[i]->xbin);
ypeaks[n] = top/butt;
+ weight[n] = (Int_t)butt;
+ //cout<<"mean in y "<<ypeaks[n]<<" on x "<<windowPt[i]->xbin<<" content "<<butt<<endl;
n++;
if(n==max_entries) break;
}
//Improve the peaks by including the region around in x.
Float_t ytop,ybutt;
+ Int_t prev;
+ Int_t w;
for(Int_t i=0; i<n; i++)
{
Int_t xbin = fCurrentHisto->FindXbin(xpeaks[i]);
if(xbin - x_bin_sides < xmin || xbin + x_bin_sides > xmax) continue;
top=butt=0;
ytop=0,ybutt=0;
+ w=0;
+ prev = xbin - x_bin_sides+1;
for(Int_t j=xbin-x_bin_sides; j<=xbin+x_bin_sides; j++)
{
+ /*
+ //Check if the windows are overlapping:
+ if(anotherPt[j]->ymin > anotherPt[prev]->ymax) {prev=j; continue;}
+ if(anotherPt[j]->ymax < anotherPt[prev]->ymin) {prev=j; continue;}
+ prev = j;
+ */
+
top += fCurrentHisto->GetBinCenterX(j)*anotherPt[j]->weight;
butt += anotherPt[j]->weight;
+
for(Int_t k=anotherPt[j]->ymin; k<anotherPt[j]->ymax; k++)
{
Int_t bin = fCurrentHisto->GetBin(j,k);
ytop += (fCurrentHisto->GetBinCenterY(k))*(fCurrentHisto->GetBinContent(bin));
ybutt += fCurrentHisto->GetBinContent(bin);
+ w+=(Int_t)fCurrentHisto->GetBinContent(bin);
}
}
+
xpeaks[i] = top/butt;
ypeaks[i] = ytop/ybutt;
+ weight[i] = w;
+ //cout<<"Setting weight "<<w<<" kappa "<<xpeaks[i]<<" phi0 "<<ypeaks[i]<<endl;
+
+ //Check if this peak is overlapping with a previous:
+ for(Int_t p=0; p<i-1; p++)
+ {
+ if(fabs(xpeaks[p] - xpeaks[i]) < max_kappa ||
+ fabs(ypeaks[p] - ypeaks[i]) < max_phi0)
+ {
+ weight[i]=0;
+ break;
+ }
+ }
+
}
for(Int_t i=0; i<nbinsx; i++)
}
-AliL3HoughTrack *AliL3HoughMaxFinder::FindPeak(Int_t t1,Double_t t2,Int_t t3)
+void AliL3HoughMaxFinder::FindPeak(Int_t t1,Double_t t2,Int_t t3,Float_t &kappa,Float_t &phi0)
{
//Attempt of a more sophisticated peak finder.
//Finds the best peak in the histogram, and returns the corresponding
if(!fCurrentHisto)
{
printf("AliL3HoughMaxFinder::FindPeak : No histogram!!\n");
- return 0;
+ return;
}
AliL3Histogram *hist = fCurrentHisto;
//Find the weight:
- bin = hist->FindBin(x_peak,y_peak);
- Int_t weight = (Int_t)hist->GetBinContent(bin);
+ //bin = hist->FindBin(x_peak,y_peak);
+ //Int_t weight = (Int_t)hist->GetBinContent(bin);
- AliL3HoughTrack *track = new AliL3HoughTrack();
- track->SetTrackParameters(x_peak,y_peak,weight);
-
+ //AliL3HoughTrack *track = new AliL3HoughTrack();
+ //track->SetTrackParameters(x_peak,y_peak,weight);
+ kappa = x_peak;
+ phi0 = y_peak;
delete [] m;
delete [] m_low;
delete [] m_up;
- return track;
+ //return track;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff