Updating changes before the weekend
authorvestbo <vestbo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 30 Mar 2001 14:55:23 +0000 (14:55 +0000)
committervestbo <vestbo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 30 Mar 2001 14:55:23 +0000 (14:55 +0000)
12 files changed:
HLT/hough/AliL3Defs.h
HLT/hough/AliL3HoughEval.cxx
HLT/hough/AliL3HoughEval.h
HLT/hough/AliL3HoughLinkDef.h
HLT/hough/AliL3HoughMaxFinder.cxx
HLT/hough/AliL3HoughMaxFinder.h
HLT/hough/AliL3HoughMerge.cxx
HLT/hough/AliL3HoughTransformer.cxx
HLT/hough/AliL3HoughTransformer.h
HLT/hough/Makefile
HLT/hough/hough_merge.C
HLT/hough/testPF.C

index c9bfdf052e5414dc684dea2c52d197116d792901..a59c4c5525c9bfc295dfc2686081a58625db2f48 100644 (file)
@@ -5,6 +5,7 @@
 
 const Int_t NRowsSlice = 174;
 const Int_t NRows[5][2] = {{ 0, 45},{46,77},{78,109},{110,141},{142,173}};
+//const Int_t NRows[5][2] = {{ 0, 31},{46,173},{78,109},{110,141},{142,173}};
 const Double_t Pi = 3.14159265358979323846;
 const Double_t ToRad = Pi/180.;
 
index b7c608f252256e448d083e53ac6228364f6f8fe9..5bdfcc92d345aa0e6b37a43d3c4967df810a6cb0 100644 (file)
@@ -40,7 +40,7 @@ AliL3HoughEval::~AliL3HoughEval()
 }
 
 
-void AliL3HoughEval::LookInsideRoad(AliL3TrackArray *tracks,TH2F *hist,TH2F *fake)
+void AliL3HoughEval::LookInsideRoad(AliL3TrackArray *tracks,Bool_t remove,TH2F *hist)
 {
   //Look at rawdata along the road specified by the track candidates.
 
@@ -75,8 +75,6 @@ void AliL3HoughEval::LookInsideRoad(AliL3TrackArray *tracks,TH2F *hist,TH2F *fak
              continue;
            }
          
-         if(fake)
-           fake->Fill(xyz[0],xyz[1],1);
          Int_t npixs = 0;
          
          //Get the pixels along the track candidate
@@ -93,6 +91,9 @@ void AliL3HoughEval::LookInsideRoad(AliL3TrackArray *tracks,TH2F *hist,TH2F *fak
              if(fabs(xyz_pix[1] - xyz[1]) > num_of_pads_to_look*fTransform->GetPadPitchWidthLow()) continue; 
              npixs++;
              
+             if(remove) //Remove this pixel from image
+               pixel->fIndex = -1;
+
              if(hist)
                {
                  //fTransform->Local2Global(xyz_pix,slice);
index 1d64d91f74c901e1d43af2f83fae780bbe071184..eb56a8acefc329c91ec6ee27d38ebd22faf2e27d 100644 (file)
@@ -25,7 +25,7 @@ class AliL3HoughEval : public TObject {
   virtual ~AliL3HoughEval();
 
   TClonesArray *GetParticles(Char_t *rootfile);
-  void LookInsideRoad(AliL3TrackArray *tracks,TH2F *hist=0,TH2F *fake=0);
+  void LookInsideRoad(AliL3TrackArray *tracks,Bool_t remove=(Bool_t)false,TH2F *hist=0);
   void DisplaySlice(TH2F *hist);
   void CompareMC(Char_t *rootfile,AliL3TrackArray *merged_tracks,Float_t *eta);
   
index 534862b547fd3907f84267813b8883902590391c..18a9a234a4ec6c3dd16e34d48d800468de1ad9c8 100644 (file)
@@ -4,9 +4,8 @@
 #pragma link off all classes;
 #pragma link off all functions;
 
+#pragma link C++ class AliL3Hough; 
 #pragma link C++ class AliL3HoughTransformer; 
-#pragma link C++ class AliL3HoughPixel; 
-//#pragma link C++ class AliL3HoughTrack; 
 #pragma link C++ class AliL3HoughMaxFinder;
 #pragma link C++ class AliL3HoughEval;
 #pragma link C++ class AliL3HoughMerge;
index 9a502f0a03be841a39dde48ef9432778087ec92c..ee8473a4e3b24040d27ec9f3ddb873439a223337 100644 (file)
@@ -35,6 +35,56 @@ AliL3HoughMaxFinder::~AliL3HoughMaxFinder()
 
 }
 
+AliL3TrackArray *AliL3HoughMaxFinder::FindBigMaxima(TH2F *hist)
+{
+  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 bin[25],bin_index;
+  Stat_t value[25];
+  
+  AliL3TrackArray *tracks = new AliL3TrackArray("AliL3HoughTrack");
+  AliL3HoughTrack *track;
+  
+  for(Int_t xbin=xmin+2; xbin<xmax-3; xbin++)
+    {
+      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 yb=ybin-2; yb<ybin+3; yb++)
+               {
+                 bin[bin_index]=hist->GetBin(xb,yb);
+                 value[bin_index]=hist->GetBinContent(bin[bin_index]);
+                 bin_index++;
+               }
+             
+           }
+         if(value[12]==0) continue;
+         Int_t b=0;
+         while(1)
+           {
+             if(value[b] > value[12] || b==bin_index) break;
+             b++;
+             printf("b %d\n",b);
+           }
+         if(b == bin_index)
+           {
+             //Found maxima
+             Double_t max_x = hist->GetXaxis()->GetBinCenter(xbin);
+             Double_t max_y = hist->GetYaxis()->GetBinCenter(ybin);
+             track = (AliL3HoughTrack*)tracks->NextTrack();
+             track->SetTrackParameters(max_x,max_y,value[12]);
+           }
+       }
+    }
+  
+  tracks->QSort();
+  return tracks;
+}
+
 
 AliL3TrackArray *AliL3HoughMaxFinder::FindMaxima(TH2F *hist,Int_t *rowrange,Int_t ref_row)
 {
@@ -112,17 +162,283 @@ AliL3TrackArray *AliL3HoughMaxFinder::FindMaxima(TH2F *hist,Int_t *rowrange,Int_
   return tracks;
 }
 
+
+AliL3TrackArray *AliL3HoughMaxFinder::LookForPeaks(TH2F *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 weight_loc;
+  Int_t candidate=0;
+  for(Int_t xbin=xmin+nbins; xbin <= xmax - nbins; xbin++)
+    {
+      for(Int_t ybin=ymin+nbins; ybin <= ymax - nbins; ybin++)
+       {
+         weight_loc=0;
+         for(Int_t xbin_loc = xbin-nbins; xbin_loc <= xbin+nbins; xbin_loc++)
+           {
+             for(Int_t ybin_loc = ybin-nbins; ybin_loc <= ybin+nbins; ybin_loc++)
+               {
+                 Int_t bin_loc = hist->GetBin(xbin_loc,ybin_loc);
+                 weight_loc += (Int_t)hist->GetBinContent(bin_loc);
+               }
+           }
+         
+         if(weight_loc > 0)
+           {
+             track = (AliL3HoughTrack*)tracks->NextTrack();
+             track->SetTrackParameters(hist->GetXaxis()->GetBinCenter(xbin),hist->GetYaxis()->GetBinCenter(ybin),weight_loc);
+           }
+         
+       }
+    }
+  tracks->QSort();
+  
+  AliL3HoughTrack *track1,*track2;
+
+  for(Int_t i=1; i<tracks->GetNTracks(); i++)
+    {
+      track1 = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
+      if(!track1) continue;
+      Int_t xbin1 = hist->GetXaxis()->FindBin(track1->GetKappa());
+      Int_t ybin1 = hist->GetYaxis()->FindBin(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());
+         if(abs(xbin1-xbin2) < 10 && abs(ybin1-ybin2) < 10)
+           {
+             tracks->Remove(i);
+             break;
+           }
+       }
+
+    }
+  tracks->Compress();
+  return tracks;
+}
+
+AliL3TrackArray *AliL3HoughMaxFinder::LookInWindows(TH2F *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();
+  
+  for(Int_t xbin=xmin+nbins; xbin < xmax - nbins; xbin += nbins)
+    {
+      for(Int_t ybin=ymin+nbins; ybin<ymax-nbins; ybin += nbins)
+       {
+         //Int_t bin = hist->GetBin(xbin,ybin);
+         //if((Int_t)hist->GetBinContent(bin)==0) continue;
+
+         Int_t xrange[2] = {xbin-nbins,xbin+nbins};
+         Int_t yrange[2] = {ybin-nbins,ybin+nbins};
+         if(LocatePeak(hist,track,xrange,yrange,t1,t2,t3))
+           track = (AliL3HoughTrack*)tracks->NextTrack();
+         else
+           continue;
+       }
+    }
+  tracks->QSort();
+  return tracks;
+  
+}
+
+Bool_t AliL3HoughMaxFinder::LocatePeak(TH2F *hist,AliL3HoughTrack *track,Int_t *xrange,Int_t *yrange,Int_t t1,Double_t t2,Int_t t3)
+{
+
+  /*  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++)
+    {
+      m[i]=0;
+      m_low[i]=0;
+      m_up[i]=0;
+    }
+
+  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++)
+       {
+         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 - 1;
+           }
+         
+       }
+      
+      if(m[xbin] > max_x) //Max value globally in x-direction
+       {
+         max_xbin = xbin;
+         max_x = m[xbin];//sum;
+       }
+    }
+  
+  if(max_x == 0) 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]));
+
+  //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+1;
+         break;
+       }
+    }
+  for(Int_t xbin = max_xbin+1; xbin <=xmax; xbin++)
+    {
+      if(m[xbin]*t2 < max_x)
+       {
+         x_up = xbin-1;
+         break;
+       }
+    }
+  
+  Double_t top=0,butt=0,value,x_peak;
+  
+  // printf("xlow %f xup %f\n",hist->GetXaxis()->GetBinCenter(x_low),hist->GetXaxis()->GetBinCenter(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);
+      top += value*m[xbin];
+      butt += m[xbin];
+    }
+  x_peak = top/butt;
+  
+  //printf("FOund x_peak %f\n",x_peak);
+
+  //Find the peak in y direction:
+  Int_t x_l = hist->GetXaxis()->FindBin(x_peak);
+  if(hist->GetXaxis()->GetBinCenter(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));
+
+  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);
+
+  //Find the weight:
+  bin = hist->FindBin(x_peak,y_peak);
+  Int_t weight = (Int_t)hist->GetBinContent(bin);
+
+  if(weight==0) return false;
+  
+  track->SetTrackParameters(x_peak,y_peak,weight);
+  
+  delete [] m;
+  delete [] m_up;
+  delete [] m_low;
+  return true;
+}
+
+
 AliL3TrackArray *AliL3HoughMaxFinder::FindPeak(TH2F *hist,Int_t t1,Double_t t2,Int_t t3)
 {
   //Attempt of a more sophisticated peak finder.
   
-  AliL3TrackArray *tracks = new AliL3TrackArray("AliL3HoughTrack");
+  AliL3TrackArray *tracks = new AliL3TrackArray("AliL3HoughTrack",1);
 
   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 nbinsx = hist->GetXaxis()->GetNbins()+1;
   
   Int_t *m = new Int_t[nbinsx];
   Int_t *m_low = new Int_t[nbinsx];
@@ -270,11 +586,13 @@ AliL3TrackArray *AliL3HoughMaxFinder::FindPeak(TH2F *hist,Int_t t1,Double_t t2,I
   AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->NextTrack();
   track->SetTrackParameters(x_peak,y_peak,weight);
   
+  
+  delete [] m;
+  delete [] m_low;
+  delete [] m_up;
+  
   return tracks;
 
-  //delete [] m;
-  //delete [] m_low;
-  //delete [] m_up;
-  
+    
 }
 
index 4eca14eb5d99b428bb592ff0ed005eefdd6d49d4..d7788136aeb1a36dedd37aff71e7c7e6f0a4c44c 100644 (file)
@@ -4,6 +4,7 @@
 #include "AliL3RootTypes.h"
 
 class AliL3TrackArray;
+class AliL3HoughTrack;
 
 class AliL3HoughMaxFinder : public TObject {
   
@@ -19,7 +20,11 @@ class AliL3HoughMaxFinder : public TObject {
   AliL3HoughMaxFinder(Char_t *histotype);
   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);
   void SetThreshold(Int_t f) {fThreshold = f;}
   
index 534e329c8c3d7a4466d24c7840ec5079c54e651e..81604dae3fc47f0525cfdad8e8675428ae4c8d7c 100644 (file)
@@ -61,11 +61,3 @@ void AliL3HoughMerge::FillHisto(TH2F *merge_hist)
     }
 }
 
-/*
-void AliL3HoughMerge::MergeLines()
-{
-
-  
-
-}
-*/
index 39bf6d71691ddf42f90ea0547a349ba7c7b30cbe..8208830515edab69aaed1f7c060eacd22a57533f 100644 (file)
@@ -26,6 +26,7 @@ AliL3HoughTransformer::AliL3HoughTransformer()
   fPatch = 0;
   fRowContainer = 0;
   fPhiRowContainer = 0;
+  fVolume=0;
   fNumOfPadRows=0;
   fContainerBounds=0;
   fNDigits=0;
@@ -33,21 +34,41 @@ AliL3HoughTransformer::AliL3HoughTransformer()
 }
 
 
-AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Float_t *etarange,Int_t phi_segments)
+AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Float_t *etarange)
 {
   //Constructor
   
   fTransform = new AliL3Transform();
-  
 
   fEtaMin = etarange[0];
   fEtaMax = etarange[1];
   fSlice = slice; 
   fPatch = patch;
-  fNPhiSegments = phi_segments;
   fNumOfPadRows=NRowsSlice;
 }
 
+AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Double_t *etarange,Int_t n_eta_segments)
+{
+  
+  fTransform = new AliL3Transform();
+  if(etarange)
+    {
+      //assumes you want to look at a given etaslice only
+      fEtaMin = etarange[0];
+      fEtaMax = etarange[1];
+    }
+  else
+    {
+      //transform in all etaslices
+      fEtaMin = 0;
+      fEtaMax = slice < 18 ? 0.9 : -0.9;
+    }
+  fNumEtaSegments = n_eta_segments;
+  fSlice = slice;
+  fPatch = patch;
+  fNumOfPadRows = NRowsSlice;
+}
+
 
 AliL3HoughTransformer::~AliL3HoughTransformer()
 {
@@ -58,6 +79,8 @@ AliL3HoughTransformer::~AliL3HoughTransformer()
     delete fTransform;
   if(fPhiRowContainer)
     delete [] fPhiRowContainer;
+  if(fVolume)
+    delete [] fVolume;
   if(fIndex)
     {
       for(Int_t i=0; i<fNDigits; i++)
@@ -80,8 +103,9 @@ void AliL3HoughTransformer::InitTemplates(TH2F *hist)
     fIndex[i] = new Int_t[nbinsy+1];
     
   Int_t sector,row;
+  
   for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
-    {
+    {        
       
       for(pixel=(AliL3Digits*)fRowContainer[padrow].first; pixel!=0; pixel=(AliL3Digits*)pixel->nextRowPixel)
        {
@@ -100,7 +124,7 @@ void AliL3HoughTransformer::InitTemplates(TH2F *hist)
              
              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);
              Int_t bin = hist->FindBin(kappa,phi0);
              if(fIndex[index][p]!=0)
                printf("AliL3HoughTransformer::InitTemplates : Overlapping indexes\n");
@@ -210,81 +234,7 @@ void AliL3HoughTransformer::CountBins()
 }
 
 
-void AliL3HoughTransformer::Transform2Circle(TH2F *hist,Int_t middle_row)
-{
-  //Transformation is done with respect to local coordinates in slice.
-
-
-  AliL3Digits *pix1;
-  Int_t sector,row;
-  
-  //Define a common point
-  /*
-    Double_t rowdist1 = fTransform->Row2X(middle_row-1);
-    Double_t rowdist2 = fTransform->Row2X(middle_row);
-    Double_t r_in_bundle = rowdist1 + (rowdist1-rowdist2)/2;
-  */
-  
-  Double_t r_in_bundle = fTransform->Row2X(middle_row);
-  //Make overlap between slices
-  Double_t phi_min = -15*ToRad;
-  Double_t phi_max = 15*ToRad;
-  
-  Double_t phi_slice = (phi_max - phi_min)/fNPhiSegments;
-  
-  for(Int_t p=0; p <= fNPhiSegments; p++)
-    {
-      Double_t phi_in_bundle = phi_min + p*phi_slice;
-      //printf("phi %f in slice %d patch %d middle row %f\n",phi_in_bundle/ToRad,fSlice,fPatch,r_in_bundle);
-      
-      for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
-       //for(Int_t padrow = middle_row; padrow <= 173; padrow++)
-       //for(Int_t padrow = 0; padrow <= middle_row; 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);
-             //fTransform->Raw2Global(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);
-                             
-             //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 tanPhi0 = (r_in_bundle*sin(phi_pix)-r_pix*sin(phi_in_bundle))/(r_in_bundle*cos(phi_pix)-r_pix*cos(phi_in_bundle));
-             
-             Double_t phi0 = atan(tanPhi0);
-             //if(padrow > middle_row)
-             //phi0 = -phi0;
-             //if(phi0 < 0.55 || phi0 > 0.88) continue;
-             
-             Double_t angle = phi_pix - phi0;
-             Double_t kappa = 2*sin(angle)/r_pix;
-             
-             //Double_t angle = phi_in_bundle - phi0;
-             //Double_t kappa = 2*sin(angle)/r_in_bundle;
-
-             //if(kappa < -0.006 || kappa > 0.006) continue;
-             
-             Short_t signal = pix1->fCharge;
-             
-             hist->Fill(kappa,phi0,signal);
-             
-             
-           }
-         
-       }
-    }
-  
-}
-
-
-void AliL3HoughTransformer::Transform2Circle(TH2F *hist)
+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:
@@ -295,15 +245,21 @@ void AliL3HoughTransformer::Transform2Circle(TH2F *hist)
     
   Int_t nbinsy = hist->GetNbinsY();
 
+  Int_t totsignal=0,vol_index;
+  if(fNumEtaSegments==1)
+    eta_index=0; //only looking in one etaslice.
+
   for(Int_t padrow = NRows[fPatch][0]; padrow <= NRows[fPatch][1]; padrow++)
     {
-      
-      for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
+      vol_index = eta_index*fNumOfPadRows + padrow;
+      //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)
        {
          
          Short_t signal = pix1->fCharge;
          Int_t index = pix1->fIndex;
-         
+         if(index < 0) continue; //This pixel has been removed.
+         totsignal += signal;
          for(Int_t p=0; p <= nbinsy; p++)
            hist->AddBinContent(fIndex[index][p],signal);
                          
@@ -311,7 +267,7 @@ void AliL3HoughTransformer::Transform2Circle(TH2F *hist)
       
     }
     
-  
+  printf("Total signal %d\n",totsignal);
 }
 
 /*
@@ -431,8 +387,8 @@ void AliL3HoughTransformer::Transform2Line(TH2F *hist,Int_t ref_row,Int_t *rowra
              printf("AliL3HoughTransformer::Transform2Line : index %d out of range \n",index);
              return;
            }
-         //for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
-         for(pix1=(AliL3Digits*)fPhiRowContainer[index].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextPhiRowPixel)
+         for(pix1=(AliL3Digits*)fRowContainer[padrow].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextRowPixel)
+         //for(pix1=(AliL3Digits*)fPhiRowContainer[index].first; pix1!=0; pix1=(AliL3Digits*)pix1->nextPhiRowPixel)
            {
              //printf("Transforming pixel in index %d pad %d time %d padrow %d\n",index,pix1->fPad,pix1->fTime,padrow);
              Float_t xyz[3];
@@ -479,7 +435,7 @@ void AliL3HoughTransformer::GetPixels(Char_t *rootfile,TH2F *hist)
 
   Int_t digit_counter=0;
   Float_t xyz[3];
-  Double_t eta,phi;
+  Double_t eta;
   
   Int_t nrows = NRows[fPatch][1] - NRows[fPatch][0] + 1;
   printf("nrows %d slice %d patch %d\n",nrows,fSlice,fPatch);
@@ -490,19 +446,31 @@ void AliL3HoughTransformer::GetPixels(Char_t *rootfile,TH2F *hist)
   memset(fRowContainer,0,fNumOfPadRows*sizeof(AliL3HoughContainer));
 
 
-  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));
+  //fContainerBounds = (fNPhiSegments+1)*(fNumOfPadRows+1);
+  //printf("Allocating %d bytes to container of size %d\n",fContainerBounds*sizeof(AliL3HoughContainer),fContainerBounds);
 
-  Double_t phi_min = -10*ToRad;
-  Double_t phi_max = 10*ToRad;
-  Double_t delta_phi = (phi_max-phi_min)/fNPhiSegments;
+  /*
+    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 eta_slice = (fEtaMax-fEtaMin)/fNumEtaSegments;
   Int_t index;
   digit_counter=0;
 
+  printf("\nLoading ALL pixels in slice\n\n");
+
   for (Int_t i=0; i<num_of_entries; i++) 
     { 
       t->GetEvent(i);
@@ -527,11 +495,12 @@ void AliL3HoughTransformer::GetPixels(Char_t *rootfile,TH2F *hist)
        
        fTransform->Raw2Global(xyz,sector,row,pad,time);
        eta = fTransform->GetEta(xyz);
+       
        if(eta < fEtaMin || eta > fEtaMax) continue;
        fTransform->Global2Local(xyz,sector);
        
        
-       phi = fTransform->GetPhi(xyz);
+       //phi = fTransform->GetPhi(xyz);
        if(hist)
          hist->Fill(xyz[0],xyz[1],signal);
        
@@ -548,20 +517,37 @@ void AliL3HoughTransformer::GetPixels(Char_t *rootfile,TH2F *hist)
          ((AliL3Digits*)(fRowContainer[padrow].last))->nextRowPixel=dig;
        fRowContainer[padrow].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)
+       //thisHit->etaIndex=(Int_t)((thisHit->GetEta()-fEtaMin)/etaSlice + 1);
+       Int_t eta_index = (Int_t)((eta-fEtaMin)/eta_slice);
+       index = eta_index*fNumOfPadRows + padrow;
+       if(index > fContainerBounds || index < 0)
          {
-           printf("AliL3HoughTransform::GetPixels : index out of range %d\n",phi_index);
+           
+           //printf("AliL3HoughTransformer::GetPixels : index out of range %d %d eta_index %d padrow %d\n",index,fContainerBounds,eta_index,padrow);
            continue;
          }
-               
-       if(fPhiRowContainer[index].first == NULL)
-         fPhiRowContainer[index].first = (void*)dig;
+       
+       if(fVolume[index].first == NULL)
+         fVolume[index].first = (void*)dig;
        else
-         ((AliL3Digits*)(fPhiRowContainer[index].last))->nextPhiRowPixel = dig;
-       fPhiRowContainer[index].last=(void*)dig;
+         ((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;
+         }
+         
+         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());
       
     }
index 2833e6a40237969f48d42df32a70531c954afe47..3e6b79045ce86d2efb3a1e355b9e8e028ec8c6a7 100644 (file)
@@ -9,7 +9,8 @@ struct AliL3Digits
   UChar_t fPad;
   UShort_t fTime;
   Int_t fIndex;
-  AliL3Digits *nextPhiRowPixel;
+  AliL3Digits *fNextVolumePixel;
+  //AliL3Digits *nextPhiRowPixel;
   AliL3Digits *nextRowPixel;
 };
 
@@ -40,6 +41,7 @@ class AliL3HoughTransformer : public TObject {
   
   AliL3HoughContainer *fRowContainer; //!
   AliL3HoughContainer *fPhiRowContainer; //!
+  AliL3HoughContainer *fVolume; //!
   Int_t fContainerBounds;
   Int_t fNDigits;
   Int_t **fIndex; //!
@@ -50,11 +52,11 @@ class AliL3HoughTransformer : public TObject {
 
  public:
   AliL3HoughTransformer(); 
-  AliL3HoughTransformer(Int_t slice,Int_t patch,Float_t *etarange,Int_t phi_segments);
+  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 middle_row);
-  void Transform2Circle(TH2F *hist);
+  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);
index a4b97f7efb66c7a32fc95425dacef72c068d2531..6268a76fefea01056fd7aca4d0f026dfb759abc6 100644 (file)
@@ -10,7 +10,7 @@ PACKAGE = AliL3Hough
 # C++ sources
 
 
-SRCS          = AliL3HoughTransformer.cxx AliL3HoughPixel.cxx  \
+SRCS          = AliL3HoughTransformer.cxx AliL3Hough.cxx \
                 AliL3HoughMaxFinder.cxx AliL3HoughEval.cxx AliL3HoughMerge.cxx 
                 
 
index 74c7c79b8770602d8c1889c0232f60a164ffed44..db1f017f7e818f4231194f36c607803fbb1925d0 100644 (file)
@@ -2,8 +2,8 @@ void hough_merge(char *rootfile,Bool_t MC=false)
 {
   gStyle->SetOptStat(0);
   
-  Int_t xbin = 60;
-  Int_t ybin = 60;
+  Int_t xbin = 70;
+  Int_t ybin = 70;
   Float_t xrange[2] = {-0.006 , 0.006}; //Pt 0.1->4GeV 0.006-0.006
   //Float_t yrange[2] = {-0.17 , 0.17}; //slice 2 0.55->0.88
   Float_t yrange[2] = {-0.26 , 0.26}; //slice 2 0.55->0.88
@@ -11,91 +11,61 @@ void hough_merge(char *rootfile,Bool_t MC=false)
 
   Int_t xr[2] = {0,250};
   Int_t yr[2] = {-125,125};
-  TH1F *ntracks = new TH1F("ntracks","",100,0,200);
-  TH2F *raw = new TH2F("raw","",250,xr[0],xr[1],250,yr[0],yr[1]);
-  TH2F *road = new TH2F("road","",250,0,250,250,yr[0],yr[1]);
-  TH2F *fake = new TH2F("fake","",300,0,300,250,-125,125);
-  TH2F *peaks = new TH2F("peaks","",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
+  ntracks = new TH1F("ntracks","",100,0,200);
+  raw = new TH2F("raw","",250,xr[0],xr[1],250,yr[0],yr[1]);
+  road = new TH2F("road","",250,0,250,250,yr[0],yr[1]);
+  fake = new TH2F("fake","",300,0,300,250,-125,125);
+  peaks = new TH2F("peaks","",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
   peaks->SetMarkerStyle(3);
   peaks->SetMarkerColor(2);
+  road->SetMarkerStyle(6);
+
+  int slice = 2;
+  //float eta[2] = {0.3,0.4};
+  float eta[2] = {0.04,0.05};
   
-  int slice = 2,patch=0,n_phi_segments=40;
-  float eta[2] = {0.3,0.4};
-  // float eta[2] = {0.04,0.05};
-  
-  AliL3HoughTransformer *a;
-  AliL3HoughMaxFinder *b;
-  AliL3HoughEval *c = new AliL3HoughEval(slice);
-  AliL3HoughMerge *d = new AliL3HoughMerge(slice);
-  TH2F *hist;
+  b = new AliL3HoughMaxFinder("KappaPhi");
+  //hist = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
+  TH2F **histos = new TH2F*[5];
   for(int pat=0; pat<5; pat++)
     {
-      a = new AliL3HoughTransformer(slice,pat,eta,n_phi_segments);
-      hist = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
+      a = new AliL3HoughTransformer(slice,pat,eta);
+      histos[pat] = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
       a->GetPixels(rootfile,raw);
-      a->Transform2Circle(hist,87);
-      
-      b = new AliL3HoughMaxFinder("KappaPhi");
-      AliL3TrackArray *tracks = b->FindMaxima(hist);
-      c->SetTransformer(a);
-      c->LookInsideRoad(tracks);
-      d->FillTracks(tracks,pat);
+      a->InitTemplates(hist);
+      a->Transform2Circle(hist);
       delete a;
-      delete b;
     }
-  return;
-  
-
-  TH2F *merge_hist = new TH2F("merge_hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
-  
-  d->FillHisto(merge_hist);
-  
   
-  b = new AliL3HoughMaxFinder(slice,pat);
-  b->FindMaxima(merge_hist);
-  
-  AliL3TrackArray *mtrs = b->GetTracks();
+  hist = new TH2F("hist","Parameter space",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
+  for(Int_t i=0; i<5; i++)
+    {
+      hist->Add(histos[i],1);
+    }
+  tracks = (AliL3TrackArray*)b->FindPeak(hist,4,0.95,5);
+  track = (AliL3HoughTrack*)tracks->GetCheckedTrack(0);
+  peaks->Fill(track->GetKappa(),track->GetPhi0(),1);
   
-  c->CompareMC(rootfile,mtrs,eta);
   
-  AliL3Transform *transform = new AliL3Transform();
-  for(int tr=0; tr<mtrs->GetNTracks(); tr++)
+  for(Int_t padrow=0; padrow<174; padrow++)
     {
-      AliL3HoughTrack *t = (AliL3HoughTrack*)mtrs->GetCheckedTrack(tr);
-      if(!t) {printf("NO TRACK11\n"); continue;}
-      if(t->GetMCid() < 0) {printf("Track %d was fake, weigth %d\n",tr,t->GetNHits()); continue;}
-      printf("Found pt %f weigth %d\n",t->GetPt(),t->GetNHits());
-      //printf("charge %f\n",t->GetCharge());
-      for(Int_t j=0; j<174; j++)
-       {
-         Float_t xyz[3];
-         if(!t->GetCrossingPoint(j,xyz)) 
-           {
-             printf("Track does not cross line\n");
-             continue;
-           }
-         road->Fill(xyz[0],xyz[1],1);
-         
-       }
-       
-      
+      Float_t xyz[3];
+      track->GetCrossingPoint(padrow,xyz);
+      road->Fill(xyz[0],xyz[1],1);
     }
-  
 
-  TCanvas *c1 = new TCanvas("c1","",800,800);
-  SetCanvasOptions(c1);
+  c1 = new TCanvas("c1","",1000,1000);
   c1->Divide(2,2);
   c1->cd(1);
-  merge_hist->Draw("box");
-  
-  //TCanvas *c2 = new TCanvas("c2","",2);
-  c1->cd(3);
+  hist->Draw("box");
+  peaks->Draw("same");
+  c1->cd(2);
   raw->Draw();
-
-  //TCanvas *c3 = new TCanvas("c3","",2);
-  c1->cd(4);
-  road->SetMarkerStyle(7);
+  c1->cd(3);
   road->Draw();
+  printf("Pt %f phi0 %f\n",track->GetPt(),track->GetPhi0());
+  return;
+  
 
   delete b;
   delete d;
index a00161c453b7b8bba43f20cc38cb67d54f4f11dd..cdbbd1abf45e6c49c27ac02cdd460efd2be54bfc 100644 (file)
@@ -22,26 +22,59 @@ void testPF(char *rootfile,int patch)
   TH2F *peaks = new TH2F("peaks","",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
   peaks->SetMarkerStyle(3);
   peaks->SetMarkerColor(2);
+  road->SetMarkerStyle(6);
+  road->SetMarkerColor(2);
   
   real_peaks = new TH2F("real_peaks","",xbin,xrange[0],xrange[1],ybin,yrange[0],yrange[1]);
   real_peaks->SetMarkerStyle(3);
   real_peaks->SetMarkerColor(4);
 
   int slice = 2,n_phi_segments=30;
-  float eta[2] = {0.3,0.4};
-  //float eta[2] = {0.03,0.04};
+  //double eta[2] = {0.3,0.4};
+  double eta[2] = {0.04,0.05};
   
-  a = new AliL3HoughTransformer(slice,patch,eta,n_phi_segments);
+  a = new AliL3HoughTransformer(slice,0,eta,1);
   a->GetPixels(rootfile,raw);
   a->InitTemplates(hist);
-  a->Transform2Circle(hist);
-    
+  a->Transform2Circle(hist,3);
+  
   b = new AliL3HoughMaxFinder("KappaPhi");
-  tracks = (AliL3TrackArray*)b->FindPeak(hist,3,0.95,5);
-
   c = new AliL3HoughEval(a);
-  c->LookInsideRoad(tracks);
-
+  //  tracks = (AliL3TrackArray*)b->FindMaxima(hist);
+  /*
+  int n_iter=1;
+  AliL3TrackArray **track_cand = new AliL3TrackArray*[n_iter];
+  for(int k=0; k<n_iter; k++)
+    {
+      hist->Reset();
+      a->Transform2Circle(hist);
+      track_cand[k] = (AliL3TrackArray*)b->FindPeak(hist,3,0.95,5);
+      c->LookInsideRoad(track_cand[k],true);
+    }
+  
+  //tracks = (AliL3TrackArray*)b->LookInWindows(hist,4,4,0.95,5);
+  
+  //tracks = (AliL3TrackArray*)b->FindMaxima(hist);
+  //cout << "Found "<<tracks->GetNTracks()<<" peaks"<<endl;
+  
+  for(int j=0; j<n_iter; j++)
+    {
+      track = (AliL3HoughTrack*)track_cand[j]->GetCheckedTrack(0);
+      if(!track) continue;
+      printf("Pt %f phi0 %f weight %d\n",track->GetPt(),track->GetPhi0(),track->GetWeight());
+      for(int padrow=0; padrow<174; padrow++)
+       {
+         Float_t xyz[3];
+         track->GetCrossingPoint(padrow,xyz);
+         road->Fill(xyz[0],xyz[1],1);
+       }
+      peaks->Fill(track->GetKappa(),track->GetPhi0(),1);
+    }
+  */
+  
+  printf("Looking for big maxima\n");
+  tracks = (AliL3TrackArray*)b->FindPeak(hist,3,0.95,5);
+  cout<<"NTracks after checking "<<tracks->GetNTracks()<<endl;
   for(int i=0; i<tracks->GetNTracks(); i++)
     {
       track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
@@ -50,37 +83,21 @@ void testPF(char *rootfile,int patch)
       printf("Pt %f phi0 %f weight %d\n",track->GetPt(),track->GetPhi0(),track->GetWeight());
     }  
   
-  c2 = new TCanvas("c2","",800,800);
-  //  c2->Divide(2,2);
-  //c2->cd(1);
-  hist->Draw("box");
-  peaks->Draw("same");
-  real_peaks->Draw("same");
-  
   
-  /*
-  AliL3TrackArray *tracks = b->FindMaxima(hist);
-  
-  track = (AliL3HoughTrack*)tracks->GetCheckedTrack(0);
-  peaks->Fill(track->GetKappa(),track->GetPhi0(),track->GetWeight());
-  int xbin = hist->GetXaxis()->FindBin(track->GetKappa());
-  int ybin = hist->GetYaxis()->FindBin(track->GetPhi0());
-  
-  TH1D *proj_kappa = hist->ProjectionX("proj_kappa",ybin,ybin);
-  TH1D *proj_phi = hist->ProjectionY("proj_phi",xbin,xbin);
-  
-  TCanvas *c1 = new TCanvas("c1","",800,800);
-  c1->Divide(2,2);
-  c1->cd(1);
+  c2 = new TCanvas("c2","",900,900);
+  c2->Divide(2,2);
+  c2->cd(1);
+  raw->Draw("");
+  road->Draw("same");
+  c2->cd(2);
+  raw->DrawCopy();
+  c2->cd(3);
   hist->Draw("box");
   peaks->Draw("same");
-  
-  c1->cd(3);
-  proj_kappa->Draw();
-  c1->cd(4);
-  proj_phi->Draw();
-  return;
-  */
+  //real_peaks->Draw("same");
+  c2->cd(4);
+  hist->DrawCopy("lego1");
+
   delete a;
   delete b;
   delete c;