fYmin = 0;
fXmax = 0;
fYmax = 0;
+ fFirstXbin = 0;
+ fLastXbin = 0;
+ fFirstYbin = 0;
+ fLastYbin = 0;
fEntries = 0;
fContent = 0;
}
-AliL3Histogram::AliL3Histogram(Char_t *name,Char_t *id,Int_t nxbin,Double_t xmin,Double_t xmax,Int_t nybin,Double_t ymin,Double_t ymax) : TH2F(name,id,nxbin,xmin,xmax,nybin,ymin,ymax)
+AliL3Histogram::AliL3Histogram(Char_t *name,Char_t *id,Int_t nxbin,Double_t xmin,Double_t xmax,Int_t nybin,Double_t ymin,Double_t ymax)
{
strcpy(fName,name);
fXmax = xmax;
fYmax = ymax;
fEntries = 0;
-
+ fFirstXbin = 1;
+ fFirstYbin = 1;
+ fLastXbin = nxbin;
+ fLastYbin = nybin;
+ fRootHisto = 0;
+
fContent = new Double_t[fNcells];
Reset();
}
//Destructor
if(fContent)
delete [] fContent;
+ if(fRootHisto)
+ delete fRootHisto;
}
}
Int_t AliL3Histogram::FindBin(Double_t x,Double_t y)
+{
+
+ Int_t xbin = FindXbin(x);
+ Int_t ybin = FindYbin(y);
+
+ return GetBin(xbin,ybin);
+}
+
+Int_t AliL3Histogram::FindXbin(Double_t x)
{
if(x < fXmin || x > fXmax)
{
- LOG(AliL3Log::kError,"AliL3Histogram::Find
Bin","array")<<AliL3Log::kDec<<
+ LOG(AliL3Log::kError,"AliL3Histogram::Find
Xbin","array")<<AliL3Log::kDec<<
"x-value out of range "<<x<<ENDLOG;
return 0;
}
+
+ return 1 + (Int_t)(fNxbins*(x-fXmin)/(fXmax-fXmin));
+
+}
+
+Int_t AliL3Histogram::FindYbin(Double_t y)
+{
if(y < fYmin || y > fYmax)
{
- LOG(AliL3Log::kError,"AliL3Histogram::Find
Bin","array")<<AliL3Log::kDec<<
+ LOG(AliL3Log::kError,"AliL3Histogram::Find
Ybin","array")<<AliL3Log::kDec<<
"y-value out of range "<<y<<ENDLOG;
return 0;
}
- Int_t xbin = 1 + (Int_t)(fNxbins*(x-fXmin)/(fXmax-fXmin));
- Int_t ybin = 1 + (Int_t)(fNybins*(y-fYmin)/(fYmax-fYmin));
+ return 1 + (Int_t)(fNybins*(y-fYmin)/(fYmax-fYmin));
+}
+
+Int_t AliL3Histogram::GetBin(Int_t xbin,Int_t ybin)
+{
+ if(xbin < 0 || xbin > GetLastXbin())
+ {
+ LOG(AliL3Log::kError,"AliL3Histogram::GetBin","array")<<AliL3Log::kDec<<
+ "xbin out of range "<<xbin<<ENDLOG;
+ return 0;
+ }
+ if(ybin < 0 || ybin > GetLastYbin())
+ {
+ LOG(AliL3Log::kError,"AliL3Histogram::FindYbin","array")<<AliL3Log::kDec<<
+ "ybin out of range "<<xbin<<ENDLOG;
+ return 0;
+ }
+
return xbin + ybin*(fNxbins+2);
+}
+
+Double_t AliL3Histogram::GetBinContent(Int_t bin)
+{
+ if(bin >= fNcells)
+ {
+ LOG(AliL3Log::kError,"AliL3Histogram::GetBinContent","array")<<AliL3Log::kDec<<
+ "bin out of range "<<bin<<ENDLOG;
+ return 0;
+ }
+
+ return fContent[bin];
+}
+
+void AliL3Histogram::SetBinContent(Int_t xbin,Int_t ybin)
+{
+ Int_t bin = GetBin(xbin,ybin);
+ SetBinContent(bin);
+}
+
+void AliL3Histogram::SetBinContent(Int_t bin)
+{
+
+ if(bin >= fNcells)
+ {
+ LOG(AliL3Log::kError,"AliL3Histogram::SetBinContent","array")<<AliL3Log::kDec<<
+ "bin out of range "<<bin<<ENDLOG;
+ return;
+ }
+ fContent[bin]=0;
}
void AliL3Histogram::AddBinContent(Int_t xbin,Int_t ybin,Int_t weight)
{
- Int_t bin = xbin + ybin*(fNxbins+2);
+ Int_t bin = GetBin(xbin,ybin);
AddBinContent(bin,weight);
}
fContent[bin] += weight;
}
-Double_t AliL3Histogram::GetXBinCenter(Int_t xbin)
+Double_t AliL3Histogram::GetBinCenterX(Int_t xbin)
{
Double_t binwidth = (fXmax - fXmin) / fNxbins;
}
-Double_t AliL3Histogram::GetYBinCenter(Int_t ybin)
+Double_t AliL3Histogram::GetBinCenterY(Int_t ybin)
{
Double_t binwidth = (fYmax - fYmin) / fNybins;
}
-void AliL3Histogram::Draw()
+void AliL3Histogram::Draw(Char_t *option)
{
- TH2F *hist = new TH2F(fName,"",fNxbins,fXmin,fXmax,fNybins,fYmin,fYmax);
+ fRootHisto = new TH2F(fName,"",fNxbins,fXmin,fXmax,fNybins,fYmin,fYmax);
for(Int_t bin=0; bin<fNcells; bin++)
{
- hist->AddBinContent(bin,fContent[bin]);
- if(fContent[bin]!=0) printf("bin %d\n",bin);
+ fRootHisto->AddBinContent(bin,fContent[bin]);
}
//printf("ncells %d %d\n",(hist->GetNbinsX()+2)*(hist->GetNbinsY()+2),fNcells);
- printf("maxbin %d\n",hist->GetMaximumBin());
- hist->Draw();
+ fRootHisto->Draw(option);
+
}
#include "AliL3RootTypes.h"
#include <TH2.h>
-class AliL3Histogram : public TH2F {
+
+class AliL3Histogram : public TObject {
private:
Int_t fNybins;
Int_t fNcells;
Int_t fEntries;
-
+ Int_t fFirstXbin;
+ Int_t fFirstYbin;
+ Int_t fLastXbin;
+ Int_t fLastYbin;
+
Double_t fXmin;
Double_t fYmin;
Double_t fXmax;
Double_t fYmax;
+ TH2F *fRootHisto;
public:
AliL3Histogram();
void Reset();
void Fill(Double_t x,Double_t y,Int_t weight);
Int_t FindBin(Double_t x,Double_t y);
+ Int_t FindXbin(Double_t x);
+ Int_t FindYbin(Double_t y);
+ Int_t GetBin(Int_t xbin,Int_t ybin);
+ Double_t GetBinContent(Int_t bin);
+ void SetBinContent(Int_t xbin,Int_t ybin);
+ void SetBinContent(Int_t bin);
void AddBinContent(Int_t xbin,Int_t ybin,Int_t weight);
void AddBinContent(Int_t bin,Int_t weight);
- void Draw();
+ void Draw(Char_t *option="hist");
+ TH2F *GetRootHisto() {return fRootHisto;}
Double_t GetXmin() {return fXmin;}
Double_t GetXmax() {return fXmax;}
Double_t GetYmin() {return fYmin;}
Double_t GetYmax() {return fXmax;}
- Double_t GetXBinCenter(Int_t xbin);
- Double_t GetYBinCenter(Int_t ybin);
- Int_t GetFirstXbin() {return 1 + (Int_t)(fNxbins*(fXmin-fXmin)/(fXmax-fXmin));}
- Int_t GetLastXbin() {return 1 + (Int_t)(fNxbins*(fXmax-fXmin)/(fXmax-fXmin));}
- Int_t GetFirstYbin() {return 1 + (Int_t)(fNxbins*(fXmin-fXmin)/(fXmax-fXmin));}
- Int_t GetLastYbin() {return 1 + (Int_t)(fNybins*(fYmax-fYmin)/(fYmax-fYmin));}
+ Double_t GetBinCenterX(Int_t xbin);
+ Double_t GetBinCenterY(Int_t ybin);
+ Int_t GetFirstXbin() {return fFirstXbin;}
+ Int_t GetLastXbin() {return fLastXbin;}
+ Int_t GetFirstYbin() {return fFirstYbin;}
+ Int_t GetLastYbin() {return fLastYbin;}
+ Int_t GetNbinsX() {return fNxbins;}
+ Int_t GetNbinsY() {return fNybins;}
Int_t GetNEntries() {return fEntries;}
ClassDef(AliL3Histogram,1)
-#include <string.h>
+//Author: Anders Strand Vestbo
+//Last Modified: 28.6.01
-#include <TH2.h>
+#include <string.h>
+#include "AliL3Histogram.h"
#include "AliL3TrackArray.h"
#include "AliL3HoughTrack.h"
#include "AliL3HoughMaxFinder.h"
}
-AliL3HoughMaxFinder::AliL3HoughMaxFinder(Char_t *histotype)
+AliL3HoughMaxFinder::AliL3HoughMaxFinder(Char_t *histotype,AliL3Histogram *hist)
{
//Constructor
if(strcmp(histotype,"KappaPhi")==0) fHistoType='c';
if(strcmp(histotype,"DPsi")==0) fHistoType='l';
+ if(hist)
+ fCurrentHisto = hist;
}
}
-AliL3TrackArray *AliL3HoughMaxFinder::FindBigMaxima(TH2F *hist)
+Int_t *AliL3HoughMaxFinder::FindAbsMaxima()
{
- Int_t xmin = hist->GetXaxis()->GetFirst();
- Int_t xmax = hist->GetXaxis()->GetLast();
- Int_t ymin = hist->GetYaxis()->GetFirst();
- Int_t ymax = hist->GetYaxis()->GetLast();
+ if(!fCurrentHisto)
+ {
+ printf("AliL3HoughMaxFinder::FindAbsMaxima : No histogram!\n");
+ return 0;
+ }
+ 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;
+ Int_t *max_bin = new Int_t[2];
+ Stat_t value,max_value=0;
+
+ for(Int_t xbin=xmin; xbin<=xmax; xbin++)
+ {
+ for(Int_t ybin=ymin; ybin<=ymax; ybin++)
+ {
+ bin = hist->GetBin(xbin,ybin);
+ value = hist->GetBinContent(bin);
+ if(value>max_value)
+ {
+ max_value = value;
+ max_bin[0] = xbin;
+ max_bin[1] = ybin;
+ }
+ }
+ }
+ return max_bin;
+}
+
+AliL3TrackArray *AliL3HoughMaxFinder::FindBigMaxima(AliL3Histogram *hist)
+{
+
+ 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];
for(Int_t ybin=ymin+2; ybin<ymax-3; ybin++)
{
bin_index=0;
- for(Int_t xb=xbin-2; xb<xbin+3; xb++)
+ for(Int_t xb=xbin-2; xb<=xbin+2; xb++)
{
- for(Int_t yb=ybin-2; yb<ybin+3; yb++)
+ for(Int_t yb=ybin-2; yb<=ybin+2; yb++)
{
bin[bin_index]=hist->GetBin(xb,yb);
value[bin_index]=hist->GetBinContent(bin[bin_index]);
if(b == bin_index)
{
//Found maxima
- Double_t max_x = hist->GetXaxis()->GetBinCenter(xbin);
- Double_t max_y = hist->GetYaxis()->GetBinCenter(ybin);
+ Double_t max_x = hist->GetBinCenterX(xbin);
+ Double_t max_y = hist->GetBinCenterY(ybin);
track = (AliL3HoughTrack*)tracks->NextTrack();
track->SetTrackParameters(max_x,max_y,value[12]);
}
}
-AliL3TrackArray *AliL3HoughMaxFinder::FindMaxima(TH2F *hist,Int_t *rowrange,Int_t ref_row)
+AliL3TrackArray *AliL3HoughMaxFinder::FindMaxima(AliL3Histogram *hist,Int_t *rowrange,Int_t ref_row)
{
//Locate all the maxima in input histogram.
//Maxima is defined as bins with more entries than the
//immediately neighbouring bins.
- Int_t xmin = hist->GetXaxis()->GetFirst();
- Int_t xmax = hist->GetXaxis()->GetLast();
- Int_t ymin = hist->GetYaxis()->GetFirst();
- Int_t ymax = hist->GetYaxis()->GetLast();
+ 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];
&& value[4]>value[7] && value[4]>value[8])
{
//Found a local maxima
- Float_t max_x = hist->GetXaxis()->GetBinCenter(xbin);
- Float_t max_y = hist->GetYaxis()->GetBinCenter(ybin);
+ Float_t max_x = hist->GetBinCenterX(xbin);
+ Float_t max_y = hist->GetBinCenterY(ybin);
track = (AliL3HoughTrack*)tracks->NextTrack();
}
-AliL3TrackArray *AliL3HoughMaxFinder::LookForPeaks(TH2F *hist,Int_t nbins)
+AliL3TrackArray *AliL3HoughMaxFinder::LookForPeaks(AliL3Histogram *hist,Int_t nbins)
{
AliL3TrackArray *tracks = new AliL3TrackArray("AliL3HoughTrack");
AliL3HoughTrack *track;
- Int_t xmin = hist->GetXaxis()->GetFirst();
- Int_t xmax = hist->GetXaxis()->GetLast();
- Int_t ymin = hist->GetYaxis()->GetFirst();
- Int_t ymax = hist->GetYaxis()->GetLast();
+ Int_t xmin = hist->GetFirstXbin();
+ Int_t xmax = hist->GetLastXbin();
+ Int_t ymin = hist->GetFirstYbin();
+ Int_t ymax = hist->GetLastYbin();
Int_t weight_loc;
Int_t candidate=0;
if(weight_loc > 0)
{
track = (AliL3HoughTrack*)tracks->NextTrack();
- track->SetTrackParameters(hist->GetXaxis()->GetBinCenter(xbin),hist->GetYaxis()->GetBinCenter(ybin),weight_loc);
+ track->SetTrackParameters(hist->GetBinCenterX(xbin),hist->GetBinCenterY(ybin),weight_loc);
}
}
{
track1 = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
if(!track1) continue;
- Int_t xbin1 = hist->GetXaxis()->FindBin(track1->GetKappa());
- Int_t ybin1 = hist->GetYaxis()->FindBin(track1->GetPhi0());
+ Int_t xbin1 = hist->FindXbin(track1->GetKappa());
+ Int_t ybin1 = hist->FindXbin(track1->GetPhi0());
for(Int_t j=0; j < i; j++)
{
track2 = (AliL3HoughTrack*)tracks->GetCheckedTrack(j);
if(!track2) continue;
- Int_t xbin2 = hist->GetXaxis()->FindBin(track2->GetKappa());
- Int_t ybin2 = hist->GetYaxis()->FindBin(track2->GetPhi0());
+ Int_t xbin2 = hist->FindXbin(track2->GetKappa());
+ Int_t ybin2 = hist->FindYbin(track2->GetPhi0());
if(abs(xbin1-xbin2) < 10 && abs(ybin1-ybin2) < 10)
{
tracks->Remove(i);
return tracks;
}
-AliL3TrackArray *AliL3HoughMaxFinder::LookInWindows(TH2F *hist,Int_t nbins,Int_t t1,Double_t t2,Int_t t3)
+AliL3TrackArray *AliL3HoughMaxFinder::LookInWindows(AliL3Histogram *hist,Int_t nbins,Int_t t1,Double_t t2,Int_t t3)
{
AliL3TrackArray *tracks = new AliL3TrackArray("AliL3HoughTrack");
AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->NextTrack();
- Int_t xmin = hist->GetXaxis()->GetFirst();
- Int_t xmax = hist->GetXaxis()->GetLast();
- Int_t ymin = hist->GetYaxis()->GetFirst();
- Int_t ymax = hist->GetYaxis()->GetLast();
+ Int_t xmin = hist->GetFirstXbin();
+ Int_t xmax = hist->GetLastXbin();
+ Int_t ymin = hist->GetFirstYbin();
+ Int_t ymax = hist->GetLastYbin();
for(Int_t xbin=xmin+nbins; xbin < xmax - nbins; xbin += nbins)
{
}
-Bool_t AliL3HoughMaxFinder::LocatePeak(TH2F *hist,AliL3HoughTrack *track,Int_t *xrange,Int_t *yrange,Int_t t1,Double_t t2,Int_t t3)
+AliL3HoughTrack *AliL3HoughMaxFinder::CalculatePeakInWindow(Int_t *maxbin,Int_t t0,Int_t t1,Double_t t2,Int_t t3)
{
+ //Try to expand the area around the maxbin +- t0
- /* Int_t xmin = hist->GetXaxis()->FindBin(track->GetKappa()) - nbins;
- Int_t xmax = hist->GetXaxis()->FindBin(track->GetKappa()) + nbins;
- Int_t ymin = hist->GetYaxis()->FindBin(track->GetPhi0()) - nbins;
- Int_t ymax = hist->GetYaxis()->FindBin(track->GetPhi0()) + nbins;
- Int_t nbinsx = nbins*2 + 1;
-
- if(xmin < 0 || xmax > hist->GetXaxis()->GetNbins() || ymin < 0 || ymax > hist->GetYaxis()->GetNbins())
- return false;
- */
- Int_t xmin = xrange[0];
- Int_t xmax = xrange[1];
- Int_t ymin = yrange[0];
- Int_t ymax = yrange[1];
-
- if(xmin < 0)
- xmin=0;
- if(xmax > hist->GetXaxis()->GetNbins())
- xmax = hist->GetXaxis()->GetNbins();
- if(ymin < 0)
- ymin=0;
- if(ymax > hist->GetYaxis()->GetNbins())
- ymax = hist->GetYaxis()->GetNbins();
-
- printf("Defined xrange %d %d yrange %d %d\n",xmin,xmax,ymin,ymax);
-
-
- Int_t totbins = hist->GetXaxis()->GetNbins();
- Int_t *m = new Int_t[totbins];
- Int_t *m_low = new Int_t[totbins];
- Int_t *m_up = new Int_t[totbins];
-
-
- for(Int_t i=0; i<totbins; i++)
+ if(!fCurrentHisto)
{
- m[i]=0;
- m_low[i]=0;
- m_up[i]=0;
+ printf("AliL3HoughMaxFinder::LocatePeak : No histogram\n");
+ return 0;
}
+ 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 xlow = maxbin[0]-t0;
+ if(xlow < xmin)
+ xlow = xmin;
+ Int_t xup = maxbin[0]+t0;
+ if(xup > xmax)
+ xup = xmax;
+ Int_t ylow = maxbin[1]-t0;
+ if(ylow < ymin)
+ ylow = ymin;
+ Int_t yup = maxbin[1]+t0;
+ if(yup > ymax)
+ yup = ymax;
+
+ Int_t nbinsx = hist->GetNbinsX()+1;
+
+ Int_t *m = new Int_t[nbinsx];
+ Int_t *m_low = new Int_t[nbinsx];
+ Int_t *m_up = new Int_t[nbinsx];
+
+ recompute:
Int_t max_x=0,sum=0,max_xbin=0,bin;
-
- for(Int_t xbin=xmin; xbin<=xmax; xbin++)
+
+ for(Int_t xbin=xlow; xbin<=xup; xbin++)
{
- for(Int_t ybin=ymin; ybin < ymax - t1; ybin++)
+ for(Int_t ybin=ylow; ybin <= yup; ybin++)
{
sum = 0;
- for(Int_t y=ybin; y < ybin+t1; y++)
+ for(Int_t y=ybin; y <= ybin+t1; y++)
{
+ if(y>yup) break; //reached the upper limit in y.
//Inside window
bin = hist->GetBin(xbin,y);
sum += (Int_t)hist->GetBinContent(bin);
{
m[xbin]=sum;
m_low[xbin]=ybin;
- m_up[xbin]=ybin + t1 - 1;
+ m_up[xbin]=ybin + t1;
}
}
max_x = m[xbin];//sum;
}
}
-
- if(max_x == 0)
- {
- //printf("\nmax_x == 0\n");
- return false;
- }
//printf("max_xbin %d max_x %d m_low %d m_up %d\n",max_xbin,max_x,m_low[max_xbin],m_up[max_xbin]);
- //printf("ylow %f yup %f\n",hist->GetYaxis()->GetBinCenter(m_low[max_xbin]),hist->GetYaxis()->GetBinCenter(m_up[max_xbin]));
+ //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[max_xbin]),hist->GetBinCenterY(m_up[max_xbin]));
//Determine a width in the x-direction
Int_t x_low=0,x_up=0;
for(Int_t xbin=max_xbin-1; xbin >= xmin; xbin--)
{
- if(m[xbin]*t2 < max_x)
+ if(m[xbin] < max_x*t2)
{
x_low = xbin+1;
break;
}
for(Int_t xbin = max_xbin+1; xbin <=xmax; xbin++)
{
- if(m[xbin]*t2 < max_x)
+ if(m[xbin] < max_x*t2)
{
x_up = xbin-1;
break;
}
}
-
+ printf("x_low %d x_up %d\n",x_low,x_up);
+
Double_t top=0,butt=0,value,x_peak;
+ if(x_up - x_low + 1 > t3)
+ {
+ t1 -= 1;
+ printf("\nxrange out if limit x_up %d x_low %d\n\n",x_low,x_up);
+ if(t1 > 1)
+ goto recompute;
+ else
+ {
+ x_peak = hist->GetBinCenterX(max_xbin);
+ goto moveon;
+ }
+ }
- // printf("xlow %f xup %f\n",hist->GetXaxis()->GetBinCenter(x_low),hist->GetXaxis()->GetBinCenter(x_up));
+ //printf("xlow %f xup %f\n",hist->GetBinCenterX(x_low),hist->GetBinCenterX(x_up));
//printf("Spread in x %d\n",x_up-x_low +1);
//Now, calculate the center of mass in x-direction
for(Int_t xbin=x_low; xbin <= x_up; xbin++)
{
- value = hist->GetXaxis()->GetBinCenter(xbin);
+ value = hist->GetBinCenterX(xbin);
top += value*m[xbin];
butt += m[xbin];
}
x_peak = top/butt;
- //printf("FOund x_peak %f\n",x_peak);
-
+ moveon:
+
//Find the peak in y direction:
- Int_t x_l = hist->GetXaxis()->FindBin(x_peak);
- if(hist->GetXaxis()->GetBinCenter(x_l) > x_peak)
+ Int_t x_l = hist->FindXbin(x_peak);
+ if(hist->GetBinCenterX(x_l) > x_peak)
x_l--;
Int_t x_u = x_l + 1;
- if(hist->GetXaxis()->GetBinCenter(x_l) > x_peak || hist->GetXaxis()->GetBinCenter(x_u) <= x_peak)
- printf("\nAliL3HoughMaxFinder::FindPeak : Wrong xrange %f %f %f\n\n",hist->GetXaxis()->GetBinCenter(x_l),x_peak,hist->GetXaxis()->GetBinCenter(x_u));
-
- //printf("\nxlow %f xup %f\n",hist->GetXaxis()->GetBinCenter(x_l),hist->GetXaxis()->GetBinCenter(x_u));
+ if(hist->GetBinCenterX(x_l) > x_peak || hist->GetBinCenterX(x_u) <= x_peak)
+ printf("\nAliL3HoughMaxFinder::FindPeak : Wrong xrange %f %f %f\n\n",hist->GetBinCenterX(x_l),x_peak,hist->GetBinCenterX(x_u));
+
+ //printf("\nxlow %f xup %f\n",hist->GetBinCenterX(x_l),hist->GetBinCenterX(x_u));
value=top=butt=0;
- //printf("ylow %f yup %f\n",hist->GetYaxis()->GetBinCenter(m_low[x_l]),hist->GetYaxis()->GetBinCenter(m_up[x_l]));
- //printf("ylow %f yup %f\n",hist->GetYaxis()->GetBinCenter(m_low[x_u]),hist->GetYaxis()->GetBinCenter(m_up[x_u]));
+ //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_l]),hist->GetBinCenterY(m_up[x_l]));
+ //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_u]),hist->GetBinCenterY(m_up[x_u]));
for(Int_t ybin=m_low[x_l]; ybin <= m_up[x_l]; ybin++)
{
- value = hist->GetYaxis()->GetBinCenter(ybin);
+ value = hist->GetBinCenterY(ybin);
bin = hist->GetBin(x_l,ybin);
top += value*hist->GetBinContent(bin);
butt += hist->GetBinContent(bin);
value=top=butt=0;
for(Int_t ybin=m_low[x_u]; ybin <= m_up[x_u]; ybin++)
{
- value = hist->GetYaxis()->GetBinCenter(ybin);
+ value = hist->GetBinCenterY(ybin);
bin = hist->GetBin(x_u,ybin);
top += value*hist->GetBinContent(bin);
butt += hist->GetBinContent(bin);
//printf("y_peak_up %f\n",y_peak_up);
- Double_t x_value_up = hist->GetXaxis()->GetBinCenter(x_u);
- Double_t x_value_low = hist->GetXaxis()->GetBinCenter(x_l);
+ Double_t x_value_up = hist->GetBinCenterX(x_u);
+ Double_t x_value_low = hist->GetBinCenterX(x_l);
Double_t y_peak = (y_peak_low*(x_value_up - x_peak) + y_peak_up*(x_peak - x_value_low))/(x_value_up - x_value_low);
+
//Find the weight:
bin = hist->FindBin(x_peak,y_peak);
Int_t weight = (Int_t)hist->GetBinContent(bin);
- if(weight==0) return false;
-
+ AliL3HoughTrack *track = new AliL3HoughTrack();
track->SetTrackParameters(x_peak,y_peak,weight);
+ //Reset area around peak
+ for(Int_t xbin=x_low; xbin<=x_up; xbin++)
+ {
+ for(Int_t ybin=m_low[xbin]; ybin<=m_up[xbin]; ybin++)
+ {
+ bin = hist->GetBin(xbin,ybin);
+ hist->SetBinContent(bin,0);
+ }
+ }
+
delete [] m;
- delete [] m_up;
delete [] m_low;
- return true;
+ delete [] m_up;
+
+ return track;
+
+
}
-AliL3TrackArray *AliL3HoughMaxFinder::FindPeak(TH2F *hist,Int_t t1,Double_t t2,Int_t t3)
+AliL3HoughTrack *AliL3HoughMaxFinder::FindPeak(Int_t t1,Double_t t2,Int_t t3)
{
//Attempt of a more sophisticated peak finder.
- AliL3TrackArray *tracks = new AliL3TrackArray("AliL3HoughTrack",1);
+ if(!fCurrentHisto)
+ {
+ printf("AliL3HoughMaxFinder::FindPeak : No histogram!!\n");
+ return 0;
+ }
+ AliL3Histogram *hist = fCurrentHisto;
- Int_t xmin = hist->GetXaxis()->GetFirst();
- Int_t xmax = hist->GetXaxis()->GetLast();
- Int_t ymin = hist->GetYaxis()->GetFirst();
- Int_t ymax = hist->GetYaxis()->GetLast();
- Int_t nbinsx = hist->GetXaxis()->GetNbins()+1;
+ Int_t xmin = hist->GetFirstXbin();
+ Int_t xmax = hist->GetLastXbin();
+ Int_t ymin = hist->GetFirstYbin();
+ Int_t ymax = hist->GetLastYbin();
+ Int_t nbinsx = hist->GetNbinsX()+1;
Int_t *m = new Int_t[nbinsx];
Int_t *m_low = new Int_t[nbinsx];
for(Int_t xbin=xmin; xbin<=xmax; xbin++)
{
- for(Int_t ybin=ymin; ybin < ymax - t1; ybin++)
+ for(Int_t ybin=ymin; ybin <= ymax - t1; ybin++)
{
sum = 0;
- for(Int_t y=ybin; y < ybin+t1; y++)
+ for(Int_t y=ybin; y <= ybin+t1; y++)
{
+ if(y>ymax) break;
//Inside window
bin = hist->GetBin(xbin,y);
sum += (Int_t)hist->GetBinContent(bin);
{
m[xbin]=sum;
m_low[xbin]=ybin;
- m_up[xbin]=ybin + t1 - 1;
+ m_up[xbin]=ybin + t1;
}
}
}
}
//printf("max_xbin %d max_x %d m_low %d m_up %d\n",max_xbin,max_x,m_low[max_xbin],m_up[max_xbin]);
- //printf("ylow %f yup %f\n",hist->GetYaxis()->GetBinCenter(m_low[max_xbin]),hist->GetYaxis()->GetBinCenter(m_up[max_xbin]));
+ //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[max_xbin]),hist->GetBinCenterY(m_up[max_xbin]));
//Determine a width in the x-direction
Int_t x_low=0,x_up=0;
+
for(Int_t xbin=max_xbin-1; xbin >= xmin; xbin--)
{
- if(m[xbin]*t2 < max_x)
+ if(m[xbin] < max_x*t2)
{
x_low = xbin+1;
break;
}
for(Int_t xbin = max_xbin+1; xbin <=xmax; xbin++)
{
- if(m[xbin]*t2 < max_x)
+ if(m[xbin] < max_x*t2)
{
x_up = xbin-1;
break;
if(x_up - x_low + 1 > t3)
{
t1 -= 1;
- printf("\nxrange out if limit x_up %d x_low %d\n\n",x_low,x_up);
+ printf("\nxrange out if limit x_up %d x_low %d t1 %d\n\n",x_low,x_up,t1);
if(t1 > 1)
goto recompute;
else
{
- x_peak = hist->GetXaxis()->GetBinCenter(max_xbin);
+ x_peak = hist->GetBinCenterX(max_xbin);
goto moveon;
}
}
- //printf("xlow %f xup %f\n",hist->GetXaxis()->GetBinCenter(x_low),hist->GetXaxis()->GetBinCenter(x_up));
+ //printf("xlow %f xup %f\n",hist->GetBinCenterX(x_low),hist->GetBinCenterX(x_up));
//printf("Spread in x %d\n",x_up-x_low +1);
//Now, calculate the center of mass in x-direction
for(Int_t xbin=x_low; xbin <= x_up; xbin++)
{
- value = hist->GetXaxis()->GetBinCenter(xbin);
+ value = hist->GetBinCenterX(xbin);
top += value*m[xbin];
butt += m[xbin];
}
moveon:
//Find the peak in y direction:
- Int_t x_l = hist->GetXaxis()->FindBin(x_peak);
- if(hist->GetXaxis()->GetBinCenter(x_l) > x_peak)
+ Int_t x_l = hist->FindXbin(x_peak);
+ if(hist->GetBinCenterX(x_l) > x_peak)
x_l--;
Int_t x_u = x_l + 1;
- if(hist->GetXaxis()->GetBinCenter(x_l) > x_peak || hist->GetXaxis()->GetBinCenter(x_u) <= x_peak)
- printf("\nAliL3HoughMaxFinder::FindPeak : Wrong xrange %f %f %f\n\n",hist->GetXaxis()->GetBinCenter(x_l),x_peak,hist->GetXaxis()->GetBinCenter(x_u));
+ if(hist->GetBinCenterX(x_l) > x_peak || hist->GetBinCenterX(x_u) <= x_peak)
+ printf("\nAliL3HoughMaxFinder::FindPeak : Wrong xrange %f %f %f\n\n",hist->GetBinCenterX(x_l),x_peak,hist->GetBinCenterX(x_u));
- //printf("\nxlow %f xup %f\n",hist->GetXaxis()->GetBinCenter(x_l),hist->GetXaxis()->GetBinCenter(x_u));
+ //printf("\nxlow %f xup %f\n",hist->GetBinCenterX(x_l),hist->GetBinCenterX(x_u));
value=top=butt=0;
- //printf("ylow %f yup %f\n",hist->GetYaxis()->GetBinCenter(m_low[x_l]),hist->GetYaxis()->GetBinCenter(m_up[x_l]));
- //printf("ylow %f yup %f\n",hist->GetYaxis()->GetBinCenter(m_low[x_u]),hist->GetYaxis()->GetBinCenter(m_up[x_u]));
+ //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_l]),hist->GetBinCenterY(m_up[x_l]));
+ //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_u]),hist->GetBinCenterY(m_up[x_u]));
for(Int_t ybin=m_low[x_l]; ybin <= m_up[x_l]; ybin++)
{
- value = hist->GetYaxis()->GetBinCenter(ybin);
+ value = hist->GetBinCenterY(ybin);
bin = hist->GetBin(x_l,ybin);
top += value*hist->GetBinContent(bin);
butt += hist->GetBinContent(bin);
value=top=butt=0;
for(Int_t ybin=m_low[x_u]; ybin <= m_up[x_u]; ybin++)
{
- value = hist->GetYaxis()->GetBinCenter(ybin);
+ value = hist->GetBinCenterY(ybin);
bin = hist->GetBin(x_u,ybin);
top += value*hist->GetBinContent(bin);
butt += hist->GetBinContent(bin);
//printf("y_peak_up %f\n",y_peak_up);
- Double_t x_value_up = hist->GetXaxis()->GetBinCenter(x_u);
- Double_t x_value_low = hist->GetXaxis()->GetBinCenter(x_l);
+ Double_t x_value_up = hist->GetBinCenterX(x_u);
+ Double_t x_value_low = hist->GetBinCenterX(x_l);
Double_t y_peak = (y_peak_low*(x_value_up - x_peak) + y_peak_up*(x_peak - x_value_low))/(x_value_up - x_value_low);
bin = hist->FindBin(x_peak,y_peak);
Int_t weight = (Int_t)hist->GetBinContent(bin);
- AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->NextTrack();
+ AliL3HoughTrack *track = new AliL3HoughTrack();
track->SetTrackParameters(x_peak,y_peak,weight);
delete [] m_low;
delete [] m_up;
- return tracks;
+ return track;
}
#include "AliL3RootTypes.h"
+class AliL3Histogram;
class AliL3TrackArray;
class AliL3HoughTrack;
+
class AliL3HoughMaxFinder : public TObject {
private:
Int_t fThreshold;
- //AliL3TrackArray *fTracks; //!
+ AliL3Histogram *fCurrentHisto;
Char_t fHistoType;
public:
AliL3HoughMaxFinder();
- AliL3HoughMaxFinder(Char_t *histotype);
+ AliL3HoughMaxFinder(Char_t *histotype,AliL3Histogram *hist=0);
virtual ~AliL3HoughMaxFinder();
- AliL3TrackArray *FindBigMaxima(TH2F *hist);
- AliL3TrackArray *FindMaxima(TH2F *hist,Int_t *rowrange=0,Int_t ref_row=0);
- AliL3TrackArray *LookForPeaks(TH2F *hist,Int_t nbins);
- AliL3TrackArray *LookInWindows(TH2F *hist,Int_t nbins,Int_t t1,Double_t t2,Int_t t3);
- Bool_t LocatePeak(TH2F *hist,AliL3HoughTrack *track,Int_t *xrange,Int_t *yrange,Int_t t1,Double_t t2,Int_t t3);
- AliL3TrackArray *FindPeak(TH2F *hist,Int_t t1,Double_t t2,Int_t t3);
+ Int_t *FindAbsMaxima();
+ AliL3TrackArray *FindBigMaxima(AliL3Histogram *hist);
+ AliL3TrackArray *FindMaxima(AliL3Histogram *hist,Int_t *rowrange=0,Int_t ref_row=0);
+ AliL3TrackArray *LookForPeaks(AliL3Histogram *hist,Int_t nbins);
+ AliL3TrackArray *LookInWindows(AliL3Histogram *hist,Int_t nbins,Int_t t1,Double_t t2,Int_t t3);
+ Bool_t LocatePeak(AliL3Histogram *hist,AliL3HoughTrack *track,Int_t *xrange,Int_t *yrange,Int_t t1,Double_t t2,Int_t t3);
+ AliL3HoughTrack *FindPeak(Int_t t1,Double_t t2,Int_t t3);
+ AliL3HoughTrack *CalculatePeakInWindow(Int_t *maxbin,Int_t t0,Int_t t1,Double_t t2,Int_t t3);
+
void SetThreshold(Int_t f) {fThreshold = f;}
-
+ void SetHistogram(AliL3Histogram *hist) {fCurrentHisto = hist;}
+
ClassDef(AliL3HoughMaxFinder,1)
};
+//Author: Anders Strand Vestbo
+//Last Modified: 28.6.01
+
#include <math.h>
+#include <TH2.h>
#include <TFile.h>
#include <TTree.h>
-#include <TH2.h>
#include "AliTPCParam.h"
#include "AliSimDigits.h"
+#include "AliL3Histogram.h"
+#include "AliL3Logging.h"
#include "AliL3Defs.h"
#include "AliL3Transform.h"
#include "AliL3HoughTransformer.h"
#include "AliL3HoughTrack.h"
#include "AliL3TrackArray.h"
+#include "AliL3DigitData.h"
ClassImp(AliL3HoughTransformer)
AliL3HoughTransformer::AliL3HoughTransformer()
{
//Default constructor
- fTransform = 0;
- fEtaMin = 0;
- fEtaMax = 0;
- fSlice = 0;
- fPatch = 0;
- fRowContainer = 0;
- fPhiRowContainer = 0;
- fVolume=0;
- fNumOfPadRows=0;
- fContainerBounds=0;
- fNDigits=0;
- fIndex = 0;
+
}
fSlice = slice;
fPatch = patch;
fNumOfPadRows = NRowsSlice;
+
+ fNRowsInPatch = NRows[fPatch][1]-NRows[fPatch][0] + 1;
+ fBinTableBounds = (fNRowsInPatch+1)*(MaxNPads+1);
+ fNDigitRowData = 0;
+ fDigitRowData = 0;
+ fHistoPt = 0;
}
AliL3HoughTransformer::~AliL3HoughTransformer()
{
//Destructor
- if(fRowContainer)
- delete [] fRowContainer;
+ if(fBinTable)
+ {
+ for(Int_t i=0; i<fBinTableBounds; i++)
+ delete [] fBinTable[i];
+ delete [] fBinTable;
+ }
+ if(fEtaIndex)
+ delete [] fEtaIndex;
+ if(fTrackTable)
+ {
+ for(Int_t i=0; i<fNumEtaSegments; i++)
+ delete [] fTrackTable[i];
+ delete [] fTrackTable;
+ }
if(fTransform)
delete fTransform;
- if(fPhiRowContainer)
- delete [] fPhiRowContainer;
- if(fVolume)
- delete [] fVolume;
- if(fIndex)
+
+
+}
+
+void AliL3HoughTransformer::SetInputData(UInt_t ndigits,AliL3DigitRowData *ptr)
+{
+ fNDigitRowData = ndigits;
+ fDigitRowData = ptr;
+}
+
+void AliL3HoughTransformer::InitTables()
+{
+ //Create LUT for the circle transform.
+ //the actual transformation is done in TransformTables.
+
+ AliL3Histogram *hist = fHistoPt;
+
+ Int_t nbinsy = hist->GetNbinsY();
+
+ fBinTable = new Int_t*[fBinTableBounds];
+ Int_t index;
+
+ Double_t etaslice = fEtaMax/fNumEtaSegments;
+
+ Int_t etabounds = (MaxNPads+1)*(fNRowsInPatch+1)*(MaxNTimeBins+1);
+
+ fEtaIndex = new Int_t[etabounds];
+ for(Int_t i=0; i<etabounds; i++)
+ fEtaIndex[i] = -1;
+
+ fTrackTable = new UChar_t*[fNumEtaSegments];
+ for(Int_t i=0; i<fNumEtaSegments; i++)
+ {
+ fTrackTable[i] = new UChar_t[fBinTableBounds];
+ for(Int_t j=0; j<fBinTableBounds; j++)
+ fTrackTable[i][j] = 0;
+ }
+
+ for(Int_t r=NRows[fPatch][0]; r<=NRows[fPatch][1]; r++)
+ {
+ Int_t prow = r - NRows[fPatch][0];
+
+ for(Int_t p=0; p<fTransform->GetNPads(r); p++)
+ {
+ Float_t xyz[3];
+ Int_t sector,row;
+ for(Int_t t=0; t<fTransform->GetNTimeBins(); t++)
+ {
+ fTransform->Slice2Sector(fSlice,r,sector,row);
+ fTransform->Raw2Local(xyz,sector,row,p,t);
+ Double_t eta = fTransform->GetEta(xyz);
+ if(eta < fEtaMin || eta > fEtaMax) continue;
+ Int_t ind = (prow<<17) + (p<<9) + t;
+ if(fEtaIndex[ind]>=0)
+ printf("AliL3HoughTransformer::InitTables : Overlapping indexes in eta!!\n");
+ Int_t eta_index = (Int_t)(eta/etaslice);
+ if(eta_index < 0 || eta_index > fNumEtaSegments)
+ continue;
+ fEtaIndex[ind] = eta_index;
+
+ }
+
+ Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
+ Double_t phi_pix = fTransform->GetPhi(xyz);
+ index = (prow<<8) + p;
+ fBinTable[index] = new Int_t[nbinsy+1];
+
+ for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
+ {
+ Double_t phi0 = hist->GetBinCenterY(b);
+ Double_t kappa = 2*sin(phi_pix-phi0)/r_pix;
+ Int_t bin = hist->FindBin(kappa,phi0);
+ if(fBinTable[index][b]!=0)
+ printf("AliL3HoughTransformer::InitTables : Overlapping indexes %d %d b %d\n",fBinTable[index][b],index,b);
+ fBinTable[index][b] = bin;
+ }
+ }
+
+ }
+
+}
+
+void AliL3HoughTransformer::TransformTables(AliL3Histogram **histos,AliL3Histogram **images)
+{
+ //Transform all the pixels while reading, and fill the corresponding histograms.
+ //Transform is done using LUT created in InitTables.
+ //fTrackTable : table telling whether a specific pixel is active (nonzero):
+ //fTrackTable = 0 -> no track
+ //fTrackindex = 1 -> track present
+ //fTrackindex = 2 -> track has been removed (already found)
+ //fEtaIndex : table storing the etaindex -> used to find correct histogram to fill
+ //fBinTable : table storing all the bins to fill for each nonzero pixel
+
+ Int_t eta_index;
+ AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fDigitRowData;
+ AliL3Histogram *hist;
+
+ if(!tempPt)
{
- for(Int_t i=0; i<fNDigits; i++)
- delete [] fIndex[i];
- delete [] fIndex;
+ LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformTables","data")<<
+ "Zero datapointer"<<ENDLOG;
+ return;
}
+
+ Int_t out_count=0,tot_count=0;
+ Int_t index,ind;
+
+ for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
+ {
+ Int_t prow = i - NRows[fPatch][0];
+ AliL3DigitData *bins = tempPt->fDigitData;
+ for(UInt_t dig=0; dig<tempPt->fNDigit; dig++)
+ {
+ index = (prow<<8) + bins[dig].fPad;
+ ind = (prow<<17) + (bins[dig].fPad<<9) + bins[dig].fTime;
+ eta_index = fEtaIndex[ind];
+ if(eta_index < 0) continue; //pixel out of etarange
+
+ if(fTrackTable[eta_index][index]==2) continue; //this pixel has already been removed.
+ fTrackTable[eta_index][index] = 1; //this pixel is now marked as active (it is nonzero)
+
+ tot_count++;
+ hist = histos[eta_index];
+
+ if(images)
+ {
+ //Display the transformed images.
+ AliL3Histogram *image = images[eta_index];
+ Float_t xyz_local[3];
+ Int_t sector,row;
+ fTransform->Slice2Sector(fSlice,i,sector,row);
+ fTransform->Raw2Local(xyz_local,sector,row,bins[dig].fPad,bins[dig].fTime);
+ image->Fill(xyz_local[0],xyz_local[1],bins[dig].fCharge);
+ }
+
+ if(!hist)
+ {
+ printf("Error getting histogram!\n");
+ continue;
+ }
+ for(Int_t p=hist->GetFirstYbin(); p<=hist->GetLastYbin(); p++)
+ hist->AddBinContent(fBinTable[index][p],bins[dig].fCharge);
+
+ }
+
+ Byte_t *tmp = (Byte_t*)tempPt;
+ Int_t size = sizeof(AliL3DigitRowData) + tempPt->fNDigit*sizeof(AliL3DigitData);
+ tmp += size;
+ tempPt = (AliL3DigitRowData*)tmp;
+ }
+
+}
+
+void AliL3HoughTransformer::WriteTables()
+{
+ //Write the tables to asci files.
+
+ AliL3Histogram *hist = fHistoPt;
+ Char_t name[100];
+ sprintf(name,"histogram_table_%d.txt",fPatch);
+ FILE *file = fopen(name,"w");
+
+ Int_t etabounds = (MaxNPads+1)*(fNRowsInPatch+1)*(MaxNTimeBins+1);
+ for(Int_t i=0; i<etabounds; i++)
+ {
+ if(fEtaIndex[i]<0) continue;
+ fprintf(file,"%d %d\n",i,fEtaIndex[i]);
+ }
+ fclose(file);
+
+ sprintf(name,"bin_table_%d.txt",fPatch);
+ FILE *file2 = fopen(name,"w");
+ for(Int_t i=0; i<fBinTableBounds; i++)
+ {
+ if(!fBinTable[i]) continue;
+ fprintf(file2,"%d ",i);
+ for(Int_t j=hist->GetFirstYbin(); j<=hist->GetLastYbin(); j++)
+ {
+ fprintf(file2,"%d ",fBinTable[i][j]);
+ }
+ fprintf(file2,"\n");
+ }
+ fclose(file2);
}
+/*
void AliL3HoughTransformer::InitTemplates(TH2F *hist)
{
}
-void AliL3HoughTransformer::CountBins()
-{
-
- Int_t middle_row = 87; //middle of the slice
-
- Double_t r_in_bundle = fTransform->Row2X(middle_row);
- // Double_t phi_min = (fSlice*20 - 10)*ToRad;
- //Double_t phi_max = (fSlice*20 + 10)*ToRad;
- Double_t phi_min = -15*ToRad;
- Double_t phi_max = 15*ToRad;
-
- Double_t phi_slice = (phi_max - phi_min)/fNPhiSegments;
- Double_t min_phi0 = 10000;
- Double_t max_phi0 = 0;
- Double_t min_kappa = 100000;
- Double_t max_kappa = 0;
-
- Int_t xbin = 60;
- Int_t ybin = 60;
- Float_t xrange[2] = {-0.006 , 0.006}; //Pt 0.2->
- Float_t yrange[2] = {-0.26 , 0.26}; //slice 2 0.55->0.88
-
- TH2F *histo = new TH2F("histo","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
-
- for(Int_t padrow=NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
- {
- for(Int_t pad=0; pad < fTransform->GetNPads(padrow); pad++)
- {
- for(Int_t time = 0; time < fTransform->GetNTimeBins(); time++)
- {
- Float_t xyz[3];
- Int_t sector,row;
- fTransform->Slice2Sector(fSlice,padrow,sector,row);
- fTransform->Raw2Global(xyz,sector,row,pad,time);
- Double_t eta = fTransform->GetEta(xyz);
- if(eta < fEtaMin || eta > fEtaMax) continue;
- fTransform->Global2Local(xyz,sector);
- Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
- Double_t phi_pix = fTransform->GetPhi(xyz);
-
- for(Int_t p=0; p<=fNPhiSegments; p++)
- {
- Double_t phi_in_bundle = phi_min + p*phi_slice;
-
- Double_t tanPhi0 = (r_pix*sin(phi_in_bundle)-r_in_bundle*sin(phi_pix))/(r_pix*cos(phi_in_bundle)-r_in_bundle*cos(phi_pix));
-
- Double_t phi0 = atan(tanPhi0);
- // if(phi0 < 0.55 || phi0 > 0.88) continue;
-
- //if(phi0 < 0) phi0 = phi0 +2*Pi;
- //Double_t kappa = sin(phi_in_bundle - phi0)*2/r_in_bundle;
-
- Double_t angle = phi_pix - phi0;
- Double_t kappa = 2*sin(angle)/r_pix;
- histo->Fill(kappa,phi0,1);
- if(phi0 < min_phi0)
- min_phi0 = phi0;
- if(phi0 > max_phi0)
- max_phi0 = phi0;
- if(kappa < min_kappa)
- min_kappa = kappa;
- if(kappa > max_kappa)
- max_kappa = kappa;
-
- }
-
- }
-
- }
-
- }
- Int_t count=0,bi=0;
-
- Int_t xmin = histo->GetXaxis()->GetFirst();
- Int_t xmax = histo->GetXaxis()->GetLast();
- Int_t ymin = histo->GetYaxis()->GetFirst();
- Int_t ymax = histo->GetYaxis()->GetLast();
-
-
- for(Int_t xbin=xmin+1; xbin<xmax; xbin++)
- {
- for(Int_t ybin=ymin+1; ybin<ymax; ybin++)
- {
- bi++;
- Int_t bin = histo->GetBin(xbin,ybin);
- if(histo->GetBinContent(bin)>0)
- count++;
- }
- }
-
-
- printf("Number of possible tracks in this region %d, bins %d\n",count,bi);
- printf("Phi, min %f max %f\n",min_phi0,max_phi0);
- printf("Kappa, min %f max %f\n",min_kappa,max_kappa);
- histo->Draw("box");
-}
-
-
void AliL3HoughTransformer::Transform2Circle(TH2F *hist,Int_t eta_index)
{
//Transformation is done with respect to local coordinates in slice.
//Transform every pixel into whole phirange, using parametrisation:
//kappa = 2*sin(phi-phi0)/R
- //Assumes you run InitTemplates firstly!!!!
+ //Assumes you run InitTemplates first!!!!
AliL3Digits *pix1;
//for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
for(pix1=(AliL3Digits*)fVolume[vol_index].first; pix1!=0; pix1=(AliL3Digits*)pix1->fNextVolumePixel)
{
-
+ Float_t xyz[3];
+ fTransform->Raw2Global(xyz,2,padrow,pix1->fPad,pix1->fTime);
+ Double_t eta = fTransform->GetEta(xyz);
+ if(eta < fEtaMin || eta > fEtaMax)
+ printf("\n Eta OUT OF RANGE\n");
+
Short_t signal = pix1->fCharge;
Int_t index = pix1->fIndex;
if(index < 0) continue; //This pixel has been removed.
printf("Total signal %d\n",totsignal);
}
-/*
- void AliL3HoughTransformer::Transform2Circle(TH2F *hist)
- {
- //Transformation is done with respect to local coordinates in slice.
- //Transform every pixel into whole phirange, using parametrisation:
- //kappa = 2*sin(phi-phi0)/R
-
- printf("Transforming 1 pixel only\n");
-
- AliL3Digits *pix1;
- Int_t sector,row;
-
- Int_t ymin = hist->GetYaxis()->GetFirst();
- Int_t ymax = hist->GetYaxis()->GetLast();
-
- for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
- {
-
- for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
- {
-
- Float_t xyz[3];
- fTransform->Slice2Sector(fSlice,padrow,sector,row);
- fTransform->Raw2Local(xyz,sector,row,pix1->fPad,pix1->fTime);
-
- Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
-
- Double_t phi_pix = fTransform->GetPhi(xyz);
- Short_t signal = pix1->fCharge;
-
- for(Int_t p=ymin+1; p<=ymax; p++)
- {
- Double_t phi0 = hist->GetYaxis()->GetBinCenter(p);
- Double_t kappa = 2*sin(phi_pix-phi0)/r_pix;
- hist->Fill(kappa,phi0,signal);
- }
-
- }
-
- }
-
+void AliL3HoughTransformer::Transform2Circle(TH2F **histos,Int_t n_eta_segments,UInt_t ndigits,AliL3DigitRowData *ptr)
+{
+ //Transform all the pixels while reading them, and fill the corresponding histograms.
+ //Everything is done in one go here.
+
+ Double_t etaslice = 0.9/n_eta_segments;
+ Int_t eta_index;
+ AliL3DigitRowData *tempPt = (AliL3DigitRowData*)ptr;
+ TH2F *hist;
+
+ Int_t out_count=0,tot_count=0;
+ for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
+ {
+ printf("doing row %d\n",i);
+ for(UInt_t dig=0; dig<tempPt->fNDigit; dig++)
+ {
+
+ AliL3DigitData *bins = tempPt->fDigitData;
+ Float_t xyz[3];
+ Int_t sector,row;
+ fTransform->Slice2Sector(fSlice,i,sector,row);
+ fTransform->Raw2Local(xyz,sector,row,bins[dig].fPad,bins[dig].fTime);
+ Double_t eta = fTransform->GetEta(xyz);
+ eta_index = (Int_t)(eta/etaslice);
+ if(eta_index < 0 || eta_index >= n_eta_segments)
+ {
+ //printf("Eta index out of range %d\n",eta_index);
+ out_count++;
+ continue;
+ }
+ tot_count++;
+ hist = histos[eta_index];
+ if(!hist)
+ {
+ printf("Error getting histogramm!\n");
+ return;
+ }
+
+ Int_t ymin = hist->GetYaxis()->GetFirst();
+ Int_t ymax = hist->GetYaxis()->GetLast();
+
+ Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
+ Double_t phi_pix = fTransform->GetPhi(xyz);
+ for(Int_t p=ymin; p<=ymax; p++)
+ {
+ Double_t phi0 = hist->GetYaxis()->GetBinCenter(p);
+ Double_t kappa = 2*sin(phi_pix-phi0)/r_pix;
+ //printf("kappa %f phi0 %f\n",kappa,phi0);
+ hist->Fill(kappa,phi0,bins[dig].fCharge);
+ }
+ }
+
+
+ Byte_t *tmp = (Byte_t*)tempPt;
+ Int_t size = sizeof(AliL3DigitRowData) + tempPt->fNDigit*sizeof(AliL3DigitData);
+ tmp += size;
+ tempPt = (AliL3DigitRowData*)tmp;
+
+ }
- }
-*/
+ printf("There were %d pixels out of range and %d inside\n", out_count,tot_count);
+
+}
+
void AliL3HoughTransformer::TransformLines2Circle(TH2F *hist,AliL3TrackArray *tracks)
{
}
+
void AliL3HoughTransformer::GetPixels(Char_t *rootfile,TH2F *hist)
{
+ //read data from rootfile. more or less obsolete code this.
+
TFile *file = new TFile(rootfile);
file->cd();
//fContainerBounds = (fNPhiSegments+1)*(fNumOfPadRows+1);
//printf("Allocating %d bytes to container of size %d\n",fContainerBounds*sizeof(AliL3HoughContainer),fContainerBounds);
- /*
- if(fPhiRowContainer)
- delete [] fPhiRowContainer;
- fPhiRowContainer = new AliL3HoughContainer[fContainerBounds];
- memset(fPhiRowContainer,0,fContainerBounds*sizeof(AliL3HoughContainer));
- */
+
+ // if(fPhiRowContainer)
+ // delete [] fPhiRowContainer;
+ // fPhiRowContainer = new AliL3HoughContainer[fContainerBounds];
+ // memset(fPhiRowContainer,0,fContainerBounds*sizeof(AliL3HoughContainer));
+
fContainerBounds = (fNumEtaSegments+1)*(fNumOfPadRows+1);
if(fVolume)
delete [] fVolume;
fVolume = new AliL3HoughContainer[fContainerBounds];
memset(fVolume,0,fContainerBounds*sizeof(AliL3HoughContainer));
- /*
- Double_t phi_min = -10*ToRad;
- Double_t phi_max = 10*ToRad;
- Double_t delta_phi = (phi_max-phi_min)/fNPhiSegments;
- */
+
+ // Double_t phi_min = -10*ToRad;
+ // Double_t phi_max = 10*ToRad;
+ // Double_t delta_phi = (phi_max-phi_min)/fNPhiSegments;
+
Double_t eta_slice = (fEtaMax-fEtaMin)/fNumEtaSegments;
Int_t index;
digit_counter=0;
- printf("\nLoading ALL pixels in slice\n\n");
+ //printf("\nLoading ALL pixels in slice\n\n");
for (Int_t i=0; i<num_of_entries; i++)
{
((AliL3Digits*)(fVolume[index].last))->fNextVolumePixel = dig;
fVolume[index].last = (void*)dig;
- /*
- Int_t phi_index = (Int_t)((phi-phi_min)/delta_phi);
- index = phi_index*fNumOfPadRows + padrow;
- if(phi_index > fContainerBounds || phi_index < 0)
- {
- printf("AliL3HoughTransform::GetPixels : index out of range %d\n",phi_index);
- continue;
- }
+
+ // Int_t phi_index = (Int_t)((phi-phi_min)/delta_phi);
+ // index = phi_index*fNumOfPadRows + padrow;
+ // if(phi_index > fContainerBounds || phi_index < 0)
+ // {
+ // printf("AliL3HoughTransform::GetPixels : index out of range %d\n",phi_index);
+ // continue;
+ // }
- if(fPhiRowContainer[index].first == NULL)
- fPhiRowContainer[index].first = (void*)dig;
- else
- ((AliL3Digits*)(fPhiRowContainer[index].last))->nextPhiRowPixel = dig;
- fPhiRowContainer[index].last=(void*)dig;
- */
+ // if(fPhiRowContainer[index].first == NULL)
+ // fPhiRowContainer[index].first = (void*)dig;
+ // else
+ // ((AliL3Digits*)(fPhiRowContainer[index].last))->nextPhiRowPixel = dig;
+ // fPhiRowContainer[index].last=(void*)dig;
+
}while (digarr->Next());
}
delete file;
}
+*/
void *last;
};
-class TH2F;
+
+
+class AliL3Histogram;
class AliL3Transform;
class AliL3HoughEvaluate;
class AliL3TrackArray;
+class AliL3DigitRowData;
class AliL3HoughTransformer : public TObject {
Float_t fEtaMax;
Int_t fNumEtaSegments;
Int_t fNumOfPadRows;
+ Int_t fNRowsInPatch;
+ Int_t fBinTableBounds;
+
+ UInt_t fNDigitRowData; //!
+ AliL3DigitRowData *fDigitRowData; //!
AliL3HoughContainer *fRowContainer; //!
AliL3HoughContainer *fPhiRowContainer; //!
AliL3HoughContainer *fVolume; //!
+
Int_t fContainerBounds;
Int_t fNDigits;
- Int_t **fIndex; //!
-
+ Int_t **fBinTable; //!
+ Int_t *fEtaIndex; //!
+ UChar_t **fTrackTable; //!
+ AliL3Histogram *fHistoPt;
Int_t fSlice;
Int_t fPatch;
AliL3HoughTransformer(Int_t slice,Int_t patch,Float_t *etarange);
AliL3HoughTransformer(Int_t slice,Int_t patch,Double_t *etarange=0,Int_t n_eta_segments=1);
virtual ~AliL3HoughTransformer();
-
- void Transform2Circle(TH2F *hist,Int_t eta_index);
- void Transform2Line(TH2F *hist,Int_t ref_row,Int_t *rowrange,Double_t *phirange,TH2F *raw=0);
- void TransformLines2Circle(TH2F *hist,AliL3TrackArray *tracks);
- void GetPixels(Char_t *rootfile,TH2F *hist=0);
- void InitTemplates(TH2F *hist);
- void CountBins();
+ void InitTables();
+ void TransformTables(AliL3Histogram **histos,AliL3Histogram **images=0);
+ void SetInputData(UInt_t ndigits,AliL3DigitRowData *ptr);
+ void WriteTables();
+ void SetHistogram(AliL3Histogram *hist) {fHistoPt = hist;}
+ /*
+ void Transform2Circle(TH2F *hist,Int_t eta_index);
+ void Transform2Circle(TH2F **histos,Int_t n_eta_segments,UInt_t ndigits,AliL3DigitRowData *ptr);
+ void Transform2Line(TH2F *hist,Int_t ref_row,Int_t *rowrange,Double_t *phirange,TH2F *raw=0);
+ void TransformLines2Circle(TH2F *hist,AliL3TrackArray *tracks);
+ void GetPixels(Char_t *rootfile,TH2F *hist=0);
+ void InitTemplates(TH2F *hist);
+ */
ClassDef(AliL3HoughTransformer,1)
};
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff