Added new peakfinder

diff --git a/HLT/hough/AliL3HoughMaxFinder.cxx b/HLT/hough/AliL3HoughMaxFinder.cxx
index d36b28c6e36009a2b38c2f211fae7b231b101229..bb6b8764aa09b13358d39f49c022aa9b531c5a18 100644 (file)
--- a/HLT/hough/AliL3HoughMaxFinder.cxx
+++ b/HLT/hough/AliL3HoughMaxFinder.cxx
@@ -111,3 +111,140 @@ AliL3TrackArray *AliL3HoughMaxFinder::FindMaxima(TH2F *hist,Int_t *rowrange,Int_
   tracks->QSort();
   return tracks;
 }
+
+void AliL3HoughMaxFinder::FindPeak(TH2F *hist,Double_t *peak)
+{
+  //Attempt of a more sophisticated peak finder.
+  
+  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();
+  
+  Int_t *m = new Int_t[nbinsx];
+  Int_t *m_low = new Int_t[nbinsx];
+  Int_t *m_up = new Int_t[nbinsx];
+  for(Int_t i=0; i<nbinsx; i++)
+    {
+      m[i]=0;
+      m_low[i]=0;
+      m_up[i]=0;
+    }
+  
+
+  Int_t t1 = 2;
+  Double_t t2 = 0.95;
+
+  Int_t max_x=0,sum=0,max_xbin=0,bin;
+
+  for(Int_t xbin=xmin; xbin<=xmax; xbin++)
+    {
+      for(Int_t ybin=ymin; ybin <= ymax - t1; ybin += t1)
+       {
+         sum = 0;
+         for(Int_t y=ybin; y<=ybin+t1; y++)
+           {
+             //Inside window
+             bin = hist->GetBin(xbin,y);
+             sum += (Int_t)hist->GetBinContent(bin);
+             
+           }
+         if(sum > m[xbin]) //Max value locally in this xbin
+           {
+             m[xbin]=sum;
+             m_low[xbin]=ybin;
+             m_up[xbin]=ybin+t1;
+           }
+         
+       }
+      
+      if(m[xbin] > max_x) //Max value globally in x-direction
+       {
+         max_xbin = xbin;
+         max_x = m[xbin];//sum;
+       }
+    }
+  //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]));
+
+  //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)
+       {
+         x_low = xbin;
+         break;
+       }
+    }
+  for(Int_t xbin = max_xbin+1; xbin <=xmax; xbin++)
+    {
+      if(m[xbin]*t2 < max_x)
+       {
+         x_up = xbin;
+         break;
+       }
+    }
+  
+  //  printf("xlow %f xup %f\n",hist->GetXaxis()->GetBinCenter(x_low),hist->GetXaxis()->GetBinCenter(x_up));
+
+  //Now, calculate the center of mass in x-direction
+  Double_t top=0,butt=0,x_peak,value;
+  for(Int_t xbin=x_low; xbin <= x_up; xbin++)
+    {
+      value = hist->GetXaxis()->GetBinCenter(xbin);
+      top += value*m[xbin];
+      butt += m[xbin];
+    }
+  x_peak = top/butt;
+  
+  
+  //Find the peak in y direction:
+  Int_t x_l = hist->GetXaxis()->FindBin(x_peak);
+  Int_t x_u = x_l + 1;
+    
+  printf("xlow %f xup %f\n",hist->GetXaxis()->GetBinCenter(x_l),hist->GetXaxis()->GetBinCenter(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]));
+  
+  for(Int_t ybin=m_low[x_l]; ybin <= m_up[x_l]; ybin++)
+    {
+      value = hist->GetYaxis()->GetBinCenter(ybin);
+      bin = hist->GetBin(x_l,ybin);
+      top += value*hist->GetBinContent(bin);
+      butt += hist->GetBinContent(bin);
+    }
+  Double_t y_peak_low = top/butt;
+  
+  printf("y_peak_low %f\n",y_peak_low);
+
+  value=top=butt=0;
+  for(Int_t ybin=m_low[x_u]; ybin <= m_up[x_u]; ybin++)
+    {
+      value = hist->GetYaxis()->GetBinCenter(ybin);
+      bin = hist->GetBin(x_u,ybin);
+      top += value*hist->GetBinContent(bin);
+      butt += hist->GetBinContent(bin);
+    }
+  Double_t y_peak_up = top/butt;
+  
+  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 y_peak = (y_peak_low*(x_value_up - x_peak) + y_peak_up*(x_peak - x_value_low))/(x_value_up - x_value_low);
+
+
+  peak[0] = x_peak;
+  peak[1] = y_peak;
+
+  //delete [] m;
+  //delete [] m_low;
+  //delete [] m_up;
+  
+}

diff --git a/HLT/hough/AliL3HoughMaxFinder.h b/HLT/hough/AliL3HoughMaxFinder.h
index e937dda1cc6638bcb9befc6a319d57088d2fd281..e4d5f924cea1eb40af5fac26401079f3c31a8fbc 100644 (file)
--- a/HLT/hough/AliL3HoughMaxFinder.h
+++ b/HLT/hough/AliL3HoughMaxFinder.h
@@ -20,9 +20,9 @@ class AliL3HoughMaxFinder : public TObject {
   virtual ~AliL3HoughMaxFinder();
 
   AliL3TrackArray *FindMaxima(TH2F *hist,Int_t *rowrange=0,Int_t ref_row=0);
+  void FindPeak(TH2F *hist,Double_t *peak);
   void SetThreshold(Int_t f) {fThreshold = f;}
-
-  //AliL3TrackArray *GetTracks() {return fTracks;}
+  
 
   ClassDef(AliL3HoughMaxFinder,1)
 

diff --git a/HLT/hough/AliL3HoughTransformer.cxx b/HLT/hough/AliL3HoughTransformer.cxx
index 5dfeccdea15a7edd05bdc85a8c337b76533f85cb..39bf6d71691ddf42f90ea0547a349ba7c7b30cbe 100644 (file)
--- a/HLT/hough/AliL3HoughTransformer.cxx
+++ b/HLT/hough/AliL3HoughTransformer.cxx
@@ -92,7 +92,6 @@ void AliL3HoughTransformer::InitTemplates(TH2F *hist)
          Double_t r_pix = sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
          
          Double_t phi_pix = fTransform->GetPhi(xyz);
-         Short_t signal = pixel->fCharge;
          Int_t index = pixel->fIndex;
          if(index >= fNDigits)
            printf("AliL3HoughTransformer::InitTemplates : Index error! %d\n",index);
@@ -290,14 +289,10 @@ 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");
+  //Assumes you run InitTemplates firstly!!!!
 
   AliL3Digits *pix1;
-  Int_t sector,row;
-  
-  Int_t ymin = hist->GetYaxis()->GetFirst();
-  Int_t ymax = hist->GetYaxis()->GetLast();
+    
   Int_t nbinsy = hist->GetNbinsY();
 
   for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)

diff --git a/HLT/hough/hough.C b/HLT/hough/hough.C
index 8b1332d4e8297b6c2a1423fcd6f4d7f57e4773cc..e58f4b57b6c245b2e6b9e428b261f8c59d7e4474 100644 (file)
--- a/HLT/hough/hough.C
+++ b/HLT/hough/hough.C
@@ -24,8 +24,8 @@ void hough(char *rootfile,int patch,Bool_t MC=false)
   peaks->SetMarkerColor(2);
   
   int slice = 2,n_phi_segments=30;
-  //float eta[2] = {0.3,0.4};
-  float eta[2] = {0.03,0.04};
+  float eta[2] = {0.3,0.4};
+  //float eta[2] = {0.03,0.04};
   
   AliL3HoughTransformer *a = new AliL3HoughTransformer(slice,patch,eta,n_phi_segments);
   a->GetPixels(rootfile,raw);
@@ -39,8 +39,18 @@ void hough(char *rootfile,int patch,Bool_t MC=false)
   
   AliL3HoughMaxFinder *b = new AliL3HoughMaxFinder("KappaPhi");
   //b->SetThreshold(10000);
-  AliL3TrackArray *tracks = b->FindMaxima(hist);
-    
+  //AliL3TrackArray *tracks = b->FindMaxima(hist);
+  
+  Double_t peak[2];
+  b->FindPeak(hist,peak);
+  printf("Found peak %f %f\n",peak[0],peak[1]);
+  peaks->Fill(peak[0],peak[1],1);
+  peaks->SetMarkerStyle(3);
+  peaks->SetMarkerColor(3);
+  hist->Draw("box");
+  peaks->Draw("same");
+  return;
+
   AliL3HoughEval *c = new AliL3HoughEval(a);
   printf("number of peaks before %d\n",tracks->GetNTracks());