Major changes
authorvestbo <vestbo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 29 Jun 2001 14:43:32 +0000 (14:43 +0000)
committervestbo <vestbo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 29 Jun 2001 14:43:32 +0000 (14:43 +0000)
HLT/hough/AliL3Histogram.cxx
HLT/hough/AliL3Histogram.h
HLT/hough/AliL3HoughMaxFinder.cxx
HLT/hough/AliL3HoughMaxFinder.h
HLT/hough/AliL3HoughTransformer.cxx
HLT/hough/AliL3HoughTransformer.h

index 382481b..16c3365 100644 (file)
@@ -17,12 +17,16 @@ AliL3Histogram::AliL3Histogram()
   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);
@@ -35,7 +39,12 @@ AliL3Histogram::AliL3Histogram(Char_t *name,Char_t *id,Int_t nxbin,Double_t xmin
   fXmax = xmax;
   fYmax = ymax;
   fEntries = 0;
-  
+  fFirstXbin = 1;
+  fFirstYbin = 1;
+  fLastXbin = nxbin;
+  fLastYbin = nybin;
+  fRootHisto = 0;
+
   fContent = new Double_t[fNcells];
   Reset();
 }
@@ -45,6 +54,8 @@ AliL3Histogram::~AliL3Histogram()
   //Destructor
   if(fContent)
     delete [] fContent;
+  if(fRootHisto)
+    delete fRootHisto;
 }
 
 
@@ -65,29 +76,91 @@ void AliL3Histogram::Fill(Double_t x,Double_t y,Int_t weight)
 
 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::FindBin","array")<<AliL3Log::kDec<<
+      LOG(AliL3Log::kError,"AliL3Histogram::FindXbin","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::FindBin","array")<<AliL3Log::kDec<<
+      LOG(AliL3Log::kError,"AliL3Histogram::FindYbin","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);
 
 }
@@ -104,7 +177,7 @@ void AliL3Histogram::AddBinContent(Int_t bin,Int_t weight)
   fContent[bin] += weight;
 }
 
-Double_t AliL3Histogram::GetXBinCenter(Int_t xbin)
+Double_t AliL3Histogram::GetBinCenterX(Int_t xbin)
 {
   
   Double_t binwidth = (fXmax - fXmin) / fNxbins;
@@ -112,7 +185,7 @@ Double_t AliL3Histogram::GetXBinCenter(Int_t xbin)
   
 }
 
-Double_t AliL3Histogram::GetYBinCenter(Int_t ybin)
+Double_t AliL3Histogram::GetBinCenterY(Int_t ybin)
 {
   
   Double_t binwidth = (fYmax - fYmin) / fNybins;
@@ -121,16 +194,15 @@ Double_t AliL3Histogram::GetYBinCenter(Int_t ybin)
 }
 
 
-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);
+  
 }
index 3558b27..f5d14ff 100644 (file)
@@ -4,7 +4,8 @@
 #include "AliL3RootTypes.h"
 #include <TH2.h>
 
-class AliL3Histogram : public TH2F {
+
+class AliL3Histogram : public TObject {
   
  private:
   
@@ -14,12 +15,17 @@ class AliL3Histogram : public TH2F {
   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();
@@ -29,20 +35,29 @@ class AliL3Histogram : public TH2F {
   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)
index 21c9874..441ff92 100644 (file)
@@ -1,7 +1,9 @@
-#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"
@@ -18,7 +20,7 @@ AliL3HoughMaxFinder::AliL3HoughMaxFinder()
 }
 
 
-AliL3HoughMaxFinder::AliL3HoughMaxFinder(Char_t *histotype)
+AliL3HoughMaxFinder::AliL3HoughMaxFinder(Char_t *histotype,AliL3Histogram *hist)
 {
   //Constructor
 
@@ -26,6 +28,8 @@ AliL3HoughMaxFinder::AliL3HoughMaxFinder(Char_t *histotype)
   if(strcmp(histotype,"KappaPhi")==0) fHistoType='c';
   if(strcmp(histotype,"DPsi")==0) fHistoType='l';
   
+  if(hist)
+    fCurrentHisto = hist;
 }
 
 
@@ -35,12 +39,47 @@ AliL3HoughMaxFinder::~AliL3HoughMaxFinder()
 
 }
 
-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];
   
@@ -52,9 +91,9 @@ AliL3TrackArray *AliL3HoughMaxFinder::FindBigMaxima(TH2F *hist)
       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]);
@@ -73,8 +112,8 @@ AliL3TrackArray *AliL3HoughMaxFinder::FindBigMaxima(TH2F *hist)
          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]);
            }
@@ -86,16 +125,16 @@ AliL3TrackArray *AliL3HoughMaxFinder::FindBigMaxima(TH2F *hist)
 }
 
 
-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];
   
@@ -132,8 +171,8 @@ AliL3TrackArray *AliL3HoughMaxFinder::FindMaxima(TH2F *hist,Int_t *rowrange,Int_
             && 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();
              
@@ -163,15 +202,15 @@ AliL3TrackArray *AliL3HoughMaxFinder::FindMaxima(TH2F *hist,Int_t *rowrange,Int_
 }
 
 
-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;
@@ -192,7 +231,7 @@ AliL3TrackArray *AliL3HoughMaxFinder::LookForPeaks(TH2F *hist,Int_t nbins)
          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);
            }
          
        }
@@ -205,14 +244,14 @@ AliL3TrackArray *AliL3HoughMaxFinder::LookForPeaks(TH2F *hist,Int_t nbins)
     {
       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);
@@ -225,16 +264,16 @@ AliL3TrackArray *AliL3HoughMaxFinder::LookForPeaks(TH2F *hist,Int_t nbins)
   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)
     {
@@ -257,57 +296,52 @@ AliL3TrackArray *AliL3HoughMaxFinder::LookInWindows(TH2F *hist,Int_t nbins,Int_t
   
 }
 
-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);
@@ -317,7 +351,7 @@ Bool_t AliL3HoughMaxFinder::LocatePeak(TH2F *hist,AliL3HoughTrack *track,Int_t *
            {
              m[xbin]=sum;
              m_low[xbin]=ybin;
-             m_up[xbin]=ybin + t1 - 1;
+             m_up[xbin]=ybin + t1;
            }
          
        }
@@ -328,20 +362,14 @@ Bool_t AliL3HoughMaxFinder::LocatePeak(TH2F *hist,AliL3HoughTrack *track,Int_t *
          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;
@@ -349,49 +377,62 @@ Bool_t AliL3HoughMaxFinder::LocatePeak(TH2F *hist,AliL3HoughTrack *track,Int_t *
     }
   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);
@@ -403,7 +444,7 @@ Bool_t AliL3HoughMaxFinder::LocatePeak(TH2F *hist,AliL3HoughTrack *track,Int_t *
   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);
@@ -412,37 +453,55 @@ Bool_t AliL3HoughMaxFinder::LocatePeak(TH2F *hist,AliL3HoughTrack *track,Int_t *
   
   //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];
@@ -462,11 +521,12 @@ AliL3TrackArray *AliL3HoughMaxFinder::FindPeak(TH2F *hist,Int_t t1,Double_t t2,I
 
   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);
@@ -476,7 +536,7 @@ AliL3TrackArray *AliL3HoughMaxFinder::FindPeak(TH2F *hist,Int_t t1,Double_t t2,I
            {
              m[xbin]=sum;
              m_low[xbin]=ybin;
-             m_up[xbin]=ybin + t1 - 1;
+             m_up[xbin]=ybin + t1;
            }
          
        }
@@ -488,13 +548,14 @@ AliL3TrackArray *AliL3HoughMaxFinder::FindPeak(TH2F *hist,Int_t t1,Double_t t2,I
        }
     }
   //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;
@@ -502,7 +563,7 @@ AliL3TrackArray *AliL3HoughMaxFinder::FindPeak(TH2F *hist,Int_t t1,Double_t t2,I
     }
   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;
@@ -513,23 +574,23 @@ AliL3TrackArray *AliL3HoughMaxFinder::FindPeak(TH2F *hist,Int_t t1,Double_t t2,I
   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];
     }
@@ -538,25 +599,25 @@ AliL3TrackArray *AliL3HoughMaxFinder::FindPeak(TH2F *hist,Int_t t1,Double_t t2,I
  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);
@@ -568,7 +629,7 @@ AliL3TrackArray *AliL3HoughMaxFinder::FindPeak(TH2F *hist,Int_t t1,Double_t t2,I
   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);
@@ -577,8 +638,8 @@ AliL3TrackArray *AliL3HoughMaxFinder::FindPeak(TH2F *hist,Int_t t1,Double_t t2,I
   
   //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);
 
@@ -587,7 +648,7 @@ AliL3TrackArray *AliL3HoughMaxFinder::FindPeak(TH2F *hist,Int_t t1,Double_t t2,I
   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);
   
   
@@ -595,7 +656,7 @@ AliL3TrackArray *AliL3HoughMaxFinder::FindPeak(TH2F *hist,Int_t t1,Double_t t2,I
   delete [] m_low;
   delete [] m_up;
   
-  return tracks;
+  return track;
 
     
 }
index d778813..2776e4f 100644 (file)
@@ -3,32 +3,38 @@
 
 #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)
 
 };
index 8da50d1..6e61b53 100644 (file)
@@ -1,35 +1,30 @@
+//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;
+  
 }
 
 
@@ -66,28 +61,220 @@ AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Double_t *e
   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)
 {
 
@@ -135,110 +322,12 @@ 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;
     
@@ -254,7 +343,12 @@ void AliL3HoughTransformer::Transform2Circle(TH2F *hist,Int_t eta_index)
       //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.
@@ -269,50 +363,70 @@ void AliL3HoughTransformer::Transform2Circle(TH2F *hist,Int_t eta_index)
   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)
 {
@@ -419,9 +533,12 @@ void AliL3HoughTransformer::Transform2Line(TH2F *hist,Int_t ref_row,Int_t *rowra
     
 }
 
+
 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();
 
@@ -448,27 +565,27 @@ void AliL3HoughTransformer::GetPixels(Char_t *rootfile,TH2F *hist)
   //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++) 
     { 
@@ -532,21 +649,21 @@ void AliL3HoughTransformer::GetPixels(Char_t *rootfile,TH2F *hist)
          ((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());
       
     }
@@ -558,3 +675,4 @@ void AliL3HoughTransformer::GetPixels(Char_t *rootfile,TH2F *hist)
   delete file;
   
 }
+*/
index 3e6b790..6e702f3 100644 (file)
@@ -20,10 +20,13 @@ struct AliL3HoughContainer
   void *last;
 };
 
-class TH2F;
+
+
+class AliL3Histogram;
 class AliL3Transform;
 class AliL3HoughEvaluate;
 class AliL3TrackArray;
+class AliL3DigitRowData;
 
 class AliL3HoughTransformer : public TObject {
   
@@ -38,14 +41,22 @@ 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;
@@ -55,14 +66,20 @@ class AliL3HoughTransformer : public TObject {
   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)
 
 };