HLT/hough/AliL3HoughMaxFinder.cxx

   1 #include <string.h>
   2
   3 #include <TH2.h>
   4
   5 #include "AliL3TrackArray.h"
   6 #include "AliL3HoughTrack.h"
   7 #include "AliL3HoughMaxFinder.h"
   8
   9 ClassImp(AliL3HoughMaxFinder)
  10
  11
  12 AliL3HoughMaxFinder::AliL3HoughMaxFinder()
  13 {
  14   //Default constructor
  15   fThreshold = 0;
  16   //fTracks = 0;
  17   fHistoType=0;
  18 }
  19
  20
  21 AliL3HoughMaxFinder::AliL3HoughMaxFinder(Char_t *histotype)
  22 {
  23   //Constructor
  24
  25   //fTracks = new AliL3TrackArray("AliL3HoughTrack");
  26   if(strcmp(histotype,"KappaPhi")==0) fHistoType='c';
  27   if(strcmp(histotype,"DPsi")==0) fHistoType='l';
  28
  29 }
  30
  31
  32 AliL3HoughMaxFinder::~AliL3HoughMaxFinder()
  33 {
  34   //Destructor
  35
  36 }
  37
  38
  39 AliL3TrackArray *AliL3HoughMaxFinder::FindMaxima(TH2F *hist,Int_t *rowrange,Int_t ref_row)
  40 {
  41   //Locate all the maxima in input histogram.
  42   //Maxima is defined as bins with more entries than the
  43   //immediately neighbouring bins.
  44
  45   Int_t xmin = hist->GetXaxis()->GetFirst();
  46   Int_t xmax = hist->GetXaxis()->GetLast();
  47   Int_t ymin = hist->GetYaxis()->GetFirst();
  48   Int_t ymax = hist->GetYaxis()->GetLast();
  49   Int_t bin[9],track_counter=0;
  50   Stat_t value[9];
  51
  52   AliL3TrackArray *tracks = new AliL3TrackArray("AliL3HoughTrack");
  53   AliL3HoughTrack *track;
  54
  55   for(Int_t xbin=xmin+1; xbin<xmax-1; xbin++)
  56     {
  57       for(Int_t ybin=ymin+1; ybin<ymax-1; ybin++)
  58         {
  59           bin[0] = hist->GetBin(xbin-1,ybin-1);
  60           bin[1] = hist->GetBin(xbin,ybin-1);
  61           bin[2] = hist->GetBin(xbin+1,ybin-1);
  62           bin[3] = hist->GetBin(xbin-1,ybin);
  63           bin[4] = hist->GetBin(xbin,ybin);
  64           bin[5] = hist->GetBin(xbin+1,ybin);
  65           bin[6] = hist->GetBin(xbin-1,ybin+1);
  66           bin[7] = hist->GetBin(xbin,ybin+1);
  67           bin[8] = hist->GetBin(xbin+1,ybin+1);
  68           value[0] = hist->GetBinContent(bin[0]);
  69           value[1] = hist->GetBinContent(bin[1]);
  70           value[2] = hist->GetBinContent(bin[2]);
  71           value[3] = hist->GetBinContent(bin[3]);
  72           value[4] = hist->GetBinContent(bin[4]);
  73           value[5] = hist->GetBinContent(bin[5]);
  74           value[6] = hist->GetBinContent(bin[6]);
  75           value[7] = hist->GetBinContent(bin[7]);
  76           value[8] = hist->GetBinContent(bin[8]);
  77
  78           if(value[4] <= fThreshold) continue;//central bin below threshold
  79
  80           if(value[4]>value[0] && value[4]>value[1] && value[4]>value[2]
  81              && value[4]>value[3] && value[4]>value[5] && value[4]>value[6]
  82              && value[4]>value[7] && value[4]>value[8])
  83             {
  84               //Found a local maxima
  85               Float_t max_x = hist->GetXaxis()->GetBinCenter(xbin);
  86               Float_t max_y = hist->GetYaxis()->GetBinCenter(ybin);
  87
  88               track = (AliL3HoughTrack*)tracks->NextTrack();
  89
  90
  91               switch(fHistoType)
  92                 {
  93                 case 'c':
  94                   track->SetTrackParameters(max_x,max_y,(Int_t)value[4]);
  95                   break;
  96                 case 'l':
  97                   track->SetLineParameters(max_x,max_y,(Int_t)value[4],rowrange,ref_row);
  98                   break;
  99                 default:
 100                   printf("AliL3HoughMaxFinder: Error in tracktype\n");
 101                 }
 102
 103               track_counter++;
 104
 105
 106             }
 107           else
 108             continue; //not a maxima
 109         }
 110     }
 111   tracks->QSort();
 112   return tracks;
 113 }
 114
 115 AliL3TrackArray *AliL3HoughMaxFinder::FindPeak(TH2F *hist,Int_t t1,Double_t t2,Int_t t3)
 116 {
 117   //Attempt of a more sophisticated peak finder.
 118
 119   AliL3TrackArray *tracks = new AliL3TrackArray("AliL3HoughTrack");
 120
 121   Int_t xmin = hist->GetXaxis()->GetFirst();
 122   Int_t xmax = hist->GetXaxis()->GetLast();
 123   Int_t ymin = hist->GetYaxis()->GetFirst();
 124   Int_t ymax = hist->GetYaxis()->GetLast();
 125   Int_t nbinsx = hist->GetXaxis()->GetNbins();
 126
 127   Int_t *m = new Int_t[nbinsx];
 128   Int_t *m_low = new Int_t[nbinsx];
 129   Int_t *m_up = new Int_t[nbinsx];
 130
 131
 132  recompute:  //this is a goto.
 133
 134   for(Int_t i=0; i<nbinsx; i++)
 135     {
 136       m[i]=0;
 137       m_low[i]=0;
 138       m_up[i]=0;
 139     }
 140
 141   Int_t max_x=0,sum=0,max_xbin=0,bin;
 142
 143   for(Int_t xbin=xmin; xbin<=xmax; xbin++)
 144     {
 145       for(Int_t ybin=ymin; ybin < ymax - t1; ybin++)
 146         {
 147           sum = 0;
 148           for(Int_t y=ybin; y < ybin+t1; y++)
 149             {
 150               //Inside window
 151               bin = hist->GetBin(xbin,y);
 152               sum += (Int_t)hist->GetBinContent(bin);
 153
 154             }
 155           if(sum > m[xbin]) //Max value locally in this xbin
 156             {
 157               m[xbin]=sum;
 158               m_low[xbin]=ybin;
 159               m_up[xbin]=ybin + t1 - 1;
 160             }
 161
 162         }
 163
 164       if(m[xbin] > max_x) //Max value globally in x-direction
 165         {
 166           max_xbin = xbin;
 167           max_x = m[xbin];//sum;
 168         }
 169     }
 170   //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]);
 171   //printf("ylow %f yup %f\n",hist->GetYaxis()->GetBinCenter(m_low[max_xbin]),hist->GetYaxis()->GetBinCenter(m_up[max_xbin]));
 172
 173   //Determine a width in the x-direction
 174   Int_t x_low=0,x_up=0;
 175   for(Int_t xbin=max_xbin-1; xbin >= xmin; xbin--)
 176     {
 177       if(m[xbin]*t2 < max_x)
 178         {
 179           x_low = xbin+1;
 180           break;
 181         }
 182     }
 183   for(Int_t xbin = max_xbin+1; xbin <=xmax; xbin++)
 184     {
 185       if(m[xbin]*t2 < max_x)
 186         {
 187           x_up = xbin-1;
 188           break;
 189         }
 190     }
 191
 192   Double_t top=0,butt=0,value,x_peak;
 193   if(x_up - x_low + 1 > t3)
 194     {
 195       t1 -= 1;
 196       printf("\nxrange out if limit x_up %d x_low %d\n\n",x_low,x_up);
 197       if(t1 > 1)
 198         goto recompute;
 199       else
 200         {
 201           x_peak = hist->GetXaxis()->GetBinCenter(max_xbin);
 202           goto moveon;
 203         }
 204     }
 205
 206   //printf("xlow %f xup %f\n",hist->GetXaxis()->GetBinCenter(x_low),hist->GetXaxis()->GetBinCenter(x_up));
 207   //printf("Spread in x %d\n",x_up-x_low +1);
 208
 209   //Now, calculate the center of mass in x-direction
 210   for(Int_t xbin=x_low; xbin <= x_up; xbin++)
 211     {
 212       value = hist->GetXaxis()->GetBinCenter(xbin);
 213       top += value*m[xbin];
 214       butt += m[xbin];
 215     }
 216   x_peak = top/butt;
 217
 218  moveon:
 219
 220   //Find the peak in y direction:
 221   Int_t x_l = hist->GetXaxis()->FindBin(x_peak);
 222   if(hist->GetXaxis()->GetBinCenter(x_l) > x_peak)
 223     x_l--;
 224
 225   Int_t x_u = x_l + 1;
 226
 227   if(hist->GetXaxis()->GetBinCenter(x_l) > x_peak || hist->GetXaxis()->GetBinCenter(x_u) <= x_peak)
 228     printf("\nAliL3HoughMaxFinder::FindPeak : Wrong xrange %f %f %f\n\n",hist->GetXaxis()->GetBinCenter(x_l),x_peak,hist->GetXaxis()->GetBinCenter(x_u));
 229
 230     //printf("\nxlow %f xup %f\n",hist->GetXaxis()->GetBinCenter(x_l),hist->GetXaxis()->GetBinCenter(x_u));
 231
 232   value=top=butt=0;
 233
 234   //printf("ylow %f yup %f\n",hist->GetYaxis()->GetBinCenter(m_low[x_l]),hist->GetYaxis()->GetBinCenter(m_up[x_l]));
 235   //printf("ylow %f yup %f\n",hist->GetYaxis()->GetBinCenter(m_low[x_u]),hist->GetYaxis()->GetBinCenter(m_up[x_u]));
 236
 237   for(Int_t ybin=m_low[x_l]; ybin <= m_up[x_l]; ybin++)
 238     {
 239       value = hist->GetYaxis()->GetBinCenter(ybin);
 240       bin = hist->GetBin(x_l,ybin);
 241       top += value*hist->GetBinContent(bin);
 242       butt += hist->GetBinContent(bin);
 243     }
 244   Double_t y_peak_low = top/butt;
 245
 246   //printf("y_peak_low %f\n",y_peak_low);
 247
 248   value=top=butt=0;
 249   for(Int_t ybin=m_low[x_u]; ybin <= m_up[x_u]; ybin++)
 250     {
 251       value = hist->GetYaxis()->GetBinCenter(ybin);
 252       bin = hist->GetBin(x_u,ybin);
 253       top += value*hist->GetBinContent(bin);
 254       butt += hist->GetBinContent(bin);
 255     }
 256   Double_t y_peak_up = top/butt;
 257
 258   //printf("y_peak_up %f\n",y_peak_up);
 259
 260   Double_t x_value_up = hist->GetXaxis()->GetBinCenter(x_u);
 261   Double_t x_value_low = hist->GetXaxis()->GetBinCenter(x_l);
 262
 263   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);
 264
 265
 266   //Find the weight:
 267   bin = hist->FindBin(x_peak,y_peak);
 268   Int_t weight = (Int_t)hist->GetBinContent(bin);
 269
 270   AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->NextTrack();
 271   track->SetTrackParameters(x_peak,y_peak,weight);
 272
 273   return tracks;
 274
 275   //delete [] m;
 276   //delete [] m_low;
 277   //delete [] m_up;
 278
 279 }
 280