HLT/hough/AliL3HoughMaxFinder.cxx

   1 //$Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ASV
   5
   6 #include <math.h>
   7 #include <stdlib.h>
   8 #include <string.h>
   9
  10 #include <TNtuple.h>
  11 #include <TFile.h>
  12
  13 #include "AliL3Histogram.h"
  14 #include "AliL3TrackArray.h"
  15 #include "AliL3HoughTrack.h"
  16 #include "AliL3HoughMaxFinder.h"
  17
  18 //_____________________________________________________________
  19 // AliL3HoughMaxFinder
  20 //
  21 // Maximum finder
  22
  23 ClassImp(AliL3HoughMaxFinder)
  24
  25
  26 AliL3HoughMaxFinder::AliL3HoughMaxFinder()
  27 {
  28   //Default constructor
  29   fThreshold = 0;
  30   fHistoType=0;
  31   fXPeaks=0;
  32   fYPeaks=0;
  33   fNPeaks=0;
  34   fNMax=0;
  35   fNtuppel = 0;
  36 }
  37
  38
  39 AliL3HoughMaxFinder::AliL3HoughMaxFinder(Char_t *histotype,Int_t nmax,AliL3Histogram *hist)
  40 {
  41   //Constructor
  42
  43   //fTracks = new AliL3TrackArray("AliL3HoughTrack");
  44   if(strcmp(histotype,"KappaPhi")==0) fHistoType='c';
  45   if(strcmp(histotype,"DPsi")==0) fHistoType='l';
  46
  47   if(hist)
  48     fCurrentHisto = hist;
  49
  50   fNMax=nmax;
  51   fXPeaks = new Float_t[fNMax];
  52   fYPeaks = new Float_t[fNMax];
  53   fWeight = new Int_t[fNMax];
  54   fNtuppel = 0;
  55   fThreshold=0;
  56 }
  57
  58
  59 AliL3HoughMaxFinder::~AliL3HoughMaxFinder()
  60 {
  61   //Destructor
  62   if(fXPeaks)
  63     delete [] fXPeaks;
  64   if(fYPeaks)
  65     delete [] fYPeaks;
  66   if(fWeight)
  67     delete [] fWeight;
  68   if(fNtuppel)
  69     delete fNtuppel;
  70 }
  71
  72 void AliL3HoughMaxFinder::Reset()
  73 {
  74   for(Int_t i=0; i<fNMax; i++)
  75     {
  76       fXPeaks[i]=0;
  77       fYPeaks[i]=0;
  78       fWeight[i]=0;
  79     }
  80   fNPeaks=0;
  81 }
  82
  83 void AliL3HoughMaxFinder::CreateNtuppel()
  84 {
  85   fNtuppel = new TNtuple("ntuppel","Peak charateristics","kappa:phi0:weigth:content3:content5:content1:content7");
  86   //content#; neighbouring bins of the peak.
  87
  88 }
  89
  90 void AliL3HoughMaxFinder::WriteNtuppel(Char_t *filename)
  91 {
  92   TFile *file = TFile::Open(filename,"RECREATE");
  93   if(!file)
  94     {
  95       cerr<<"AliL3HoughMaxFinder::WriteNtuppel : Error opening file "<<filename<<endl;
  96       return;
  97     }
  98   fNtuppel->Write();
  99   file->Close();
 100 }
 101
 102 void AliL3HoughMaxFinder::FindAbsMaxima()
 103 {
 104
 105   if(!fCurrentHisto)
 106     {
 107       cerr<<"AliL3HoughMaxFinder::FindAbsMaxima : No histogram"<<endl;
 108       return;
 109     }
 110   AliL3Histogram *hist = fCurrentHisto;
 111
 112   Int_t xmin = hist->GetFirstXbin();
 113   Int_t xmax = hist->GetLastXbin();
 114   Int_t ymin = hist->GetFirstYbin();
 115   Int_t ymax = hist->GetLastYbin();
 116   Int_t bin;
 117   Double_t value,max_value=0;
 118
 119   Int_t max_xbin=0,max_ybin=0;
 120   for(Int_t xbin=xmin; xbin<=xmax; xbin++)
 121     {
 122       for(Int_t ybin=ymin; ybin<=ymax; ybin++)
 123         {
 124           bin = hist->GetBin(xbin,ybin);
 125           value = hist->GetBinContent(bin);
 126           if(value>max_value)
 127             {
 128               max_value = value;
 129               max_xbin = xbin;
 130               max_ybin = ybin;
 131             }
 132         }
 133     }
 134
 135   if(fNPeaks > fNMax)
 136     {
 137       cerr<<"AliL3HoughMaxFinder::FindAbsMaxima : Array out of range : "<<fNPeaks<<endl;
 138       return;
 139     }
 140
 141
 142   Double_t max_x = hist->GetBinCenterX(max_xbin);
 143   Double_t max_y = hist->GetBinCenterY(max_ybin);
 144   fXPeaks[fNPeaks] = max_x;
 145   fYPeaks[fNPeaks] = max_y;
 146   fWeight[fNPeaks] = (Int_t)max_value;
 147   fNPeaks++;
 148
 149   if(fNtuppel)
 150     {
 151       Int_t bin3 = hist->GetBin(max_xbin-1,max_ybin);
 152       Int_t bin5 = hist->GetBin(max_xbin+1,max_ybin);
 153       Int_t bin1 = hist->GetBin(max_xbin,max_ybin-1);
 154       Int_t bin7 = hist->GetBin(max_xbin,max_ybin+1);
 155
 156       fNtuppel->Fill(max_x,max_y,max_value,hist->GetBinContent(bin3),hist->GetBinContent(bin5),hist->GetBinContent(bin1),hist->GetBinContent(bin7));
 157     }
 158
 159 }
 160
 161 void AliL3HoughMaxFinder::FindBigMaxima()
 162 {
 163
 164   AliL3Histogram *hist = fCurrentHisto;
 165   Int_t xmin = hist->GetFirstXbin();
 166   Int_t xmax = hist->GetLastXbin();
 167   Int_t ymin = hist->GetFirstYbin();
 168   Int_t ymax = hist->GetLastYbin();
 169   Int_t bin[25],bin_index;
 170   Double_t value[25];
 171
 172   for(Int_t xbin=xmin+2; xbin<xmax-3; xbin++)
 173     {
 174       for(Int_t ybin=ymin+2; ybin<ymax-3; ybin++)
 175         {
 176           bin_index=0;
 177           for(Int_t xb=xbin-2; xb<=xbin+2; xb++)
 178             {
 179               for(Int_t yb=ybin-2; yb<=ybin+2; yb++)
 180                 {
 181                   bin[bin_index]=hist->GetBin(xb,yb);
 182                   value[bin_index]=hist->GetBinContent(bin[bin_index]);
 183                   bin_index++;
 184                 }
 185
 186             }
 187           if(value[12]==0) continue;
 188           Int_t b=0;
 189           while(1)
 190             {
 191               if(value[b] > value[12] || b==bin_index) break;
 192               b++;
 193               //printf("b %d\n",b);
 194             }
 195           if(b == bin_index)
 196             {
 197               //Found maxima
 198               if(fNPeaks > fNMax)
 199                 {
 200                   cerr<<"AliL3HoughMaxFinder::FindBigMaxima : Array out of range "<<fNPeaks<<endl;
 201                   return;
 202                 }
 203
 204               Double_t max_x = hist->GetBinCenterX(xbin);
 205               Double_t max_y = hist->GetBinCenterY(ybin);
 206               fXPeaks[fNPeaks] = max_x;
 207               fYPeaks[fNPeaks] = max_y;
 208               fNPeaks++;
 209             }
 210         }
 211     }
 212 }
 213
 214
 215 void AliL3HoughMaxFinder::FindMaxima(Double_t grad_x,Double_t grad_y)
 216 {
 217   //Locate all the maxima in input histogram.
 218   //Maxima is defined as bins with more entries than the
 219   //immediately neighbouring bins.
 220
 221   Int_t xmin = fCurrentHisto->GetFirstXbin();
 222   Int_t xmax = fCurrentHisto->GetLastXbin();
 223   Int_t ymin = fCurrentHisto->GetFirstYbin();
 224   Int_t ymax = fCurrentHisto->GetLastYbin();
 225   Int_t bin[9];
 226   Double_t value[9];
 227
 228   for(Int_t xbin=xmin+1; xbin<xmax-1; xbin++)
 229     {
 230       for(Int_t ybin=ymin+1; ybin<ymax-1; ybin++)
 231         {
 232           bin[0] = fCurrentHisto->GetBin(xbin-1,ybin-1);
 233           bin[1] = fCurrentHisto->GetBin(xbin,ybin-1);
 234           bin[2] = fCurrentHisto->GetBin(xbin+1,ybin-1);
 235           bin[3] = fCurrentHisto->GetBin(xbin-1,ybin);
 236           bin[4] = fCurrentHisto->GetBin(xbin,ybin);
 237           bin[5] = fCurrentHisto->GetBin(xbin+1,ybin);
 238           bin[6] = fCurrentHisto->GetBin(xbin-1,ybin+1);
 239           bin[7] = fCurrentHisto->GetBin(xbin,ybin+1);
 240           bin[8] = fCurrentHisto->GetBin(xbin+1,ybin+1);
 241           value[0] = fCurrentHisto->GetBinContent(bin[0]);
 242           value[1] = fCurrentHisto->GetBinContent(bin[1]);
 243           value[2] = fCurrentHisto->GetBinContent(bin[2]);
 244           value[3] = fCurrentHisto->GetBinContent(bin[3]);
 245           value[4] = fCurrentHisto->GetBinContent(bin[4]);
 246           value[5] = fCurrentHisto->GetBinContent(bin[5]);
 247           value[6] = fCurrentHisto->GetBinContent(bin[6]);
 248           value[7] = fCurrentHisto->GetBinContent(bin[7]);
 249           value[8] = fCurrentHisto->GetBinContent(bin[8]);
 250
 251
 252
 253           if(value[4]>value[0] && value[4]>value[1] && value[4]>value[2]
 254              && value[4]>value[3] && value[4]>value[5] && value[4]>value[6]
 255              && value[4]>value[7] && value[4]>value[8])
 256             {
 257               //Found a local maxima
 258               Float_t max_x = fCurrentHisto->GetBinCenterX(xbin);
 259               Float_t max_y = fCurrentHisto->GetBinCenterY(ybin);
 260
 261               cout<<"Checking for threshols "<<value[4]<<" "<<fThreshold<<endl;
 262               if((Int_t)value[4] <= fThreshold) continue;//central bin below threshold
 263
 264               if(fNPeaks > fNMax)
 265                 {
 266                   cerr<<"AliL3HoughMaxFinder::FindMaxima : Array out of range "<<fNPeaks<<endl;
 267                   return;
 268                 }
 269
 270               //Check the gradient:
 271               if(value[4]/value[3] < grad_x || value[4]/value[5] < grad_x ||
 272                  value[4]/value[1] < grad_y || value[4]/value[7] < grad_y)
 273                 continue;
 274
 275               fXPeaks[fNPeaks] = max_x;
 276               fYPeaks[fNPeaks] = max_y;
 277               fWeight[fNPeaks] = (Int_t)value[4];
 278               fNPeaks++;
 279
 280               /*
 281               //Check if the peak is overlapping with a previous:
 282               Bool_t bigger = kFALSE;
 283               for(Int_t p=0; p<entries; p++)
 284               {
 285                 if(fabs(max_x - xpeaks[p]) < kappa_overlap && fabs(max_y - ypeaks[p]) < phi_overlap)
 286               {
 287               bigger = kTRUE;
 288                       if(value[4] > weight[p]) //this peak is bigger.
 289                         {
 290                           xpeaks[p] = max_x;
 291                           ypeaks[p] = max_y;
 292                           weight[p] = (Int_t)value[4];
 293                         }
 294                       else
 295                         continue; //previous peak is bigger.
 296                     }
 297                 }
 298               if(!bigger) //there were no overlapping peaks.
 299                 {
 300                 xpeaks[entries] = max_x;
 301                   ypeaks[entries] = max_y;
 302                   weight[entries] = (Int_t)value[4];
 303                   entries++;
 304                 }
 305               */
 306             }
 307           else
 308             continue; //not a maxima
 309         }
 310     }
 311
 312 }
 313
 314
 315 AliL3TrackArray *AliL3HoughMaxFinder::LookForPeaks(AliL3Histogram *hist,Int_t nbins)
 316 {
 317
 318   AliL3TrackArray *tracks = new AliL3TrackArray("AliL3HoughTrack");
 319   AliL3HoughTrack *track;
 320   Int_t xmin = hist->GetFirstXbin();
 321   Int_t xmax = hist->GetLastXbin();
 322   Int_t ymin = hist->GetFirstYbin();
 323   Int_t ymax = hist->GetLastYbin();
 324
 325   Int_t weight_loc;
 326   for(Int_t xbin=xmin+nbins; xbin <= xmax - nbins; xbin++)
 327     {
 328       for(Int_t ybin=ymin+nbins; ybin <= ymax - nbins; ybin++)
 329         {
 330           weight_loc=0;
 331           for(Int_t xbin_loc = xbin-nbins; xbin_loc <= xbin+nbins; xbin_loc++)
 332             {
 333               for(Int_t ybin_loc = ybin-nbins; ybin_loc <= ybin+nbins; ybin_loc++)
 334                 {
 335                   Int_t bin_loc = hist->GetBin(xbin_loc,ybin_loc);
 336                   weight_loc += (Int_t)hist->GetBinContent(bin_loc);
 337                 }
 338             }
 339
 340           if(weight_loc > 0)
 341             {
 342               track = (AliL3HoughTrack*)tracks->NextTrack();
 343               track->SetTrackParameters(hist->GetBinCenterX(xbin),hist->GetBinCenterY(ybin),weight_loc);
 344             }
 345
 346         }
 347     }
 348   tracks->QSort();
 349
 350   AliL3HoughTrack *track1,*track2;
 351
 352   for(Int_t i=1; i<tracks->GetNTracks(); i++)
 353     {
 354       track1 = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
 355       if(!track1) continue;
 356       Int_t xbin1 = hist->FindXbin(track1->GetKappa());
 357       Int_t ybin1 = hist->FindXbin(track1->GetPhi0());
 358       for(Int_t j=0; j < i; j++)
 359         {
 360           track2 = (AliL3HoughTrack*)tracks->GetCheckedTrack(j);
 361           if(!track2) continue;
 362           Int_t xbin2 = hist->FindXbin(track2->GetKappa());
 363           Int_t ybin2 = hist->FindYbin(track2->GetPhi0());
 364           if(abs(xbin1-xbin2) < 10 && abs(ybin1-ybin2) < 10)
 365             {
 366               tracks->Remove(i);
 367               break;
 368             }
 369         }
 370
 371     }
 372   tracks->Compress();
 373   return tracks;
 374 }
 375
 376 AliL3HoughTrack *AliL3HoughMaxFinder::CalculatePeakInWindow(Int_t *maxbin,Int_t t0,Int_t t1,Double_t t2,Int_t t3)
 377 {
 378   //Try to expand the area around the maxbin +- t0
 379
 380   if(!fCurrentHisto)
 381     {
 382       printf("AliL3HoughMaxFinder::LocatePeak : No histogram\n");
 383       return 0;
 384     }
 385   AliL3Histogram *hist = fCurrentHisto;
 386   Int_t xmin = hist->GetFirstXbin();
 387   Int_t xmax = hist->GetLastXbin();
 388   Int_t ymin = hist->GetFirstYbin();
 389   Int_t ymax = hist->GetLastYbin();
 390
 391   Int_t xlow = maxbin[0]-t0;
 392   if(xlow < xmin)
 393     xlow = xmin;
 394   Int_t xup = maxbin[0]+t0;
 395   if(xup > xmax)
 396     xup = xmax;
 397   Int_t ylow = maxbin[1]-t0;
 398   if(ylow < ymin)
 399     ylow = ymin;
 400   Int_t yup = maxbin[1]+t0;
 401   if(yup > ymax)
 402     yup = ymax;
 403
 404   Int_t nbinsx = hist->GetNbinsX()+1;
 405
 406   Int_t *m = new Int_t[nbinsx];
 407   Int_t *m_low = new Int_t[nbinsx];
 408   Int_t *m_up = new Int_t[nbinsx];
 409
 410  recompute:
 411
 412   Int_t max_x=0,sum=0,max_xbin=0,bin;
 413
 414   for(Int_t xbin=xlow; xbin<=xup; xbin++)
 415     {
 416       for(Int_t ybin=ylow; ybin <= yup; ybin++)
 417         {
 418           sum = 0;
 419           for(Int_t y=ybin; y <= ybin+t1; y++)
 420             {
 421               if(y>yup) break; //reached the upper limit in y.
 422               //Inside window
 423               bin = hist->GetBin(xbin,y);
 424               sum += (Int_t)hist->GetBinContent(bin);
 425
 426             }
 427           if(sum > m[xbin]) //Max value locally in this xbin
 428             {
 429               m[xbin]=sum;
 430               m_low[xbin]=ybin;
 431               m_up[xbin]=ybin + t1;
 432             }
 433
 434         }
 435
 436       if(m[xbin] > max_x) //Max value globally in x-direction
 437         {
 438           max_xbin = xbin;
 439           max_x = m[xbin];//sum;
 440         }
 441     }
 442   //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]);
 443   //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[max_xbin]),hist->GetBinCenterY(m_up[max_xbin]));
 444
 445   //Determine a width in the x-direction
 446   Int_t x_low=0,x_up=0;
 447   for(Int_t xbin=max_xbin-1; xbin >= xmin; xbin--)
 448     {
 449       if(m[xbin] < max_x*t2)
 450         {
 451           x_low = xbin+1;
 452           break;
 453         }
 454     }
 455   for(Int_t xbin = max_xbin+1; xbin <=xmax; xbin++)
 456     {
 457       if(m[xbin] < max_x*t2)
 458         {
 459           x_up = xbin-1;
 460           break;
 461         }
 462     }
 463   printf("x_low %d x_up %d\n",x_low,x_up);
 464
 465   Double_t top=0,butt=0,value,x_peak;
 466   if(x_up - x_low + 1 > t3)
 467     {
 468       t1 -= 1;
 469       printf("\nxrange out if limit x_up %d x_low %d\n\n",x_low,x_up);
 470       if(t1 > 1)
 471         goto recompute;
 472       else
 473         {
 474           x_peak = hist->GetBinCenterX(max_xbin);
 475           goto moveon;
 476         }
 477     }
 478
 479   //printf("xlow %f xup %f\n",hist->GetBinCenterX(x_low),hist->GetBinCenterX(x_up));
 480   //printf("Spread in x %d\n",x_up-x_low +1);
 481
 482   //Now, calculate the center of mass in x-direction
 483   for(Int_t xbin=x_low; xbin <= x_up; xbin++)
 484     {
 485       value = hist->GetBinCenterX(xbin);
 486       top += value*m[xbin];
 487       butt += m[xbin];
 488     }
 489   x_peak = top/butt;
 490
 491  moveon:
 492
 493   //Find the peak in y direction:
 494   Int_t x_l = hist->FindXbin(x_peak);
 495   if(hist->GetBinCenterX(x_l) > x_peak)
 496     x_l--;
 497
 498   Int_t x_u = x_l + 1;
 499
 500   if(hist->GetBinCenterX(x_l) > x_peak || hist->GetBinCenterX(x_u) <= x_peak)
 501     printf("\nAliL3HoughMaxFinder::FindPeak : Wrong xrange %f %f %f\n\n",hist->GetBinCenterX(x_l),x_peak,hist->GetBinCenterX(x_u));
 502
 503     //printf("\nxlow %f xup %f\n",hist->GetBinCenterX(x_l),hist->GetBinCenterX(x_u));
 504
 505   value=top=butt=0;
 506
 507   //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_l]),hist->GetBinCenterY(m_up[x_l]));
 508   //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_u]),hist->GetBinCenterY(m_up[x_u]));
 509
 510   for(Int_t ybin=m_low[x_l]; ybin <= m_up[x_l]; ybin++)
 511     {
 512       value = hist->GetBinCenterY(ybin);
 513       bin = hist->GetBin(x_l,ybin);
 514       top += value*hist->GetBinContent(bin);
 515       butt += hist->GetBinContent(bin);
 516     }
 517   Double_t y_peak_low = top/butt;
 518
 519   //printf("y_peak_low %f\n",y_peak_low);
 520
 521   value=top=butt=0;
 522   for(Int_t ybin=m_low[x_u]; ybin <= m_up[x_u]; ybin++)
 523     {
 524       value = hist->GetBinCenterY(ybin);
 525       bin = hist->GetBin(x_u,ybin);
 526       top += value*hist->GetBinContent(bin);
 527       butt += hist->GetBinContent(bin);
 528     }
 529   Double_t y_peak_up = top/butt;
 530
 531   //printf("y_peak_up %f\n",y_peak_up);
 532
 533   Double_t x_value_up = hist->GetBinCenterX(x_u);
 534   Double_t x_value_low = hist->GetBinCenterX(x_l);
 535
 536   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);
 537
 538
 539   //Find the weight:
 540   bin = hist->FindBin(x_peak,y_peak);
 541   Int_t weight = (Int_t)hist->GetBinContent(bin);
 542
 543   AliL3HoughTrack *track = new AliL3HoughTrack();
 544   track->SetTrackParameters(x_peak,y_peak,weight);
 545
 546   //Reset area around peak
 547   for(Int_t xbin=x_low; xbin<=x_up; xbin++)
 548     {
 549       for(Int_t ybin=m_low[xbin]; ybin<=m_up[xbin]; ybin++)
 550         {
 551           bin = hist->GetBin(xbin,ybin);
 552           hist->SetBinContent(bin,0);
 553         }
 554     }
 555
 556   delete [] m;
 557   delete [] m_low;
 558   delete [] m_up;
 559
 560   return track;
 561
 562
 563 }
 564
 565 AliL3HoughTrack *AliL3HoughMaxFinder::FindPeakLine(Double_t rho,Double_t theta)
 566 {
 567   //Peak finder based on a second line transformation on kappa-phi space,
 568   //to use as a baseline.
 569
 570   if(!fCurrentHisto)
 571     {
 572       printf("AliL3HoughTransformer::FindPeakLine : No input histogram\n");
 573       return 0;
 574     }
 575
 576   //get the line parameters:
 577   Double_t a = -1./tan(theta);
 578   Double_t b = rho/sin(theta);
 579
 580   printf("rho %f theta %f\n",rho,theta);
 581   //now, start looking along the line.
 582   Int_t xmin = fCurrentHisto->GetFirstXbin();
 583   Int_t xmax = fCurrentHisto->GetLastXbin();
 584
 585   Int_t max_weight=0;
 586   Double_t max_bin[2];
 587   for(Int_t xbin=xmin; xbin<=xmax; xbin++)
 588     {
 589       Double_t x = fCurrentHisto->GetBinCenterX(xbin);
 590       Double_t y = a*x + b;
 591       Int_t bin = fCurrentHisto->FindBin(x,y);
 592       //printf("x %f y %f weight %d\n",x,y,fCurrentHisto->GetBinContent(bin));
 593       if(fCurrentHisto->GetBinContent(bin) > max_weight)
 594         {
 595           max_weight = (Int_t)fCurrentHisto->GetBinContent(bin);
 596           max_bin[0] = x;
 597           max_bin[1] = y;
 598         }
 599     }
 600
 601   AliL3HoughTrack *track = new AliL3HoughTrack();
 602   track->SetTrackParameters(max_bin[0],max_bin[1],max_weight);
 603   return track;
 604 }
 605
 606
 607
 608 void AliL3HoughMaxFinder::FindPeak1(Float_t *xpeaks,Float_t *ypeaks,Int_t *weight,Int_t &n,
 609                                     Int_t y_window,Int_t x_bin_sides)
 610 {
 611   //Testing mutliple peakfinding.
 612   //The algorithm searches the histogram for prepreaks by looking in windows
 613   //for each bin on the xaxis. The size of these windows is controlled by y_window.
 614   //Then the prepreaks are sorted according to their weight (sum inside window),
 615   //and the peak positions are calculated by taking the weighted mean in both
 616   //x and y direction. The size of the peak in x-direction is controlled by x_bin_sides.
 617
 618   Int_t max_entries = n;
 619   n=0;
 620   if(!fCurrentHisto)
 621     {
 622       printf("AliL3HoughMaxFinder::FindPeak1 : No input histogram\n");
 623       return;
 624     }
 625   //Int_t y_window=2;
 626   //Int_t x_bin_sides=1;
 627
 628   Float_t max_kappa = 0.001;
 629   Float_t max_phi0 = 0.05;
 630
 631   Int_t max_sum=0;
 632
 633   Int_t xmin = fCurrentHisto->GetFirstXbin();
 634   Int_t xmax = fCurrentHisto->GetLastXbin();
 635   Int_t ymin = fCurrentHisto->GetFirstYbin();
 636   Int_t ymax = fCurrentHisto->GetLastYbin();
 637   Int_t nbinsx = fCurrentHisto->GetNbinsX()+1;
 638
 639   AxisWindow **windowPt = new AxisWindow*[nbinsx];
 640   AxisWindow **anotherPt = new AxisWindow*[nbinsx];
 641
 642   for(Int_t i=0; i<nbinsx; i++)
 643     {
 644       windowPt[i] = new AxisWindow;
 645       bzero((void*)windowPt[i],sizeof(AxisWindow));
 646       anotherPt[i] = windowPt[i];
 647     }
 648
 649   for(Int_t xbin=xmin; xbin<=xmax; xbin++)
 650     {
 651       max_sum = 0;
 652       for(Int_t ybin=ymin; ybin<=ymax-y_window; ybin++)
 653         {
 654           Int_t sum_in_window=0;
 655           for(Int_t b=ybin; b<ybin+y_window; b++)
 656             {
 657               //inside window
 658               Int_t bin = fCurrentHisto->GetBin(xbin,b);
 659               sum_in_window += (Int_t)fCurrentHisto->GetBinContent(bin);
 660             }
 661
 662           if(sum_in_window > max_sum)
 663             {
 664               max_sum = sum_in_window;
 665               windowPt[xbin]->ymin = ybin;
 666               windowPt[xbin]->ymax = ybin + y_window;
 667               windowPt[xbin]->weight = sum_in_window;
 668               windowPt[xbin]->xbin = xbin;
 669             }
 670         }
 671     }
 672
 673   //Sort the windows according to the weight
 674   SortPeaks(windowPt,0,nbinsx);
 675
 676   Float_t top,butt;
 677   for(Int_t i=0; i<nbinsx; i++)
 678     {
 679       top=butt=0;
 680       Int_t xbin = windowPt[i]->xbin;
 681
 682       if(xbin<xmin || xbin > xmax-1) continue;
 683
 684       //Check if this is really a local maxima
 685       if(anotherPt[xbin-1]->weight > anotherPt[xbin]->weight ||
 686          anotherPt[xbin+1]->weight > anotherPt[xbin]->weight)
 687         continue;
 688
 689       for(Int_t j=windowPt[i]->ymin; j<windowPt[i]->ymax; j++)
 690         {
 691           //Calculate the mean in y direction:
 692           Int_t bin = fCurrentHisto->GetBin(windowPt[i]->xbin,j);
 693           top += (fCurrentHisto->GetBinCenterY(j))*(fCurrentHisto->GetBinContent(bin));
 694           butt += fCurrentHisto->GetBinContent(bin);
 695         }
 696       xpeaks[n] = fCurrentHisto->GetBinCenterX(windowPt[i]->xbin);
 697       ypeaks[n] = top/butt;
 698       weight[n] = (Int_t)butt;
 699       //cout<<"mean in y "<<ypeaks[n]<<" on x "<<windowPt[i]->xbin<<" content "<<butt<<endl;
 700       n++;
 701       if(n==max_entries) break;
 702     }
 703
 704
 705   //Improve the peaks by including the region around in x.
 706   Float_t ytop,ybutt;
 707   Int_t prev;
 708   Int_t w;
 709   for(Int_t i=0; i<n; i++)
 710     {
 711       Int_t xbin = fCurrentHisto->FindXbin(xpeaks[i]);
 712       if(xbin - x_bin_sides < xmin || xbin + x_bin_sides > xmax) continue;
 713       top=butt=0;
 714       ytop=0,ybutt=0;
 715       w=0;
 716       prev = xbin - x_bin_sides+1;
 717       for(Int_t j=xbin-x_bin_sides; j<=xbin+x_bin_sides; j++)
 718         {
 719           /*
 720           //Check if the windows are overlapping:
 721           if(anotherPt[j]->ymin > anotherPt[prev]->ymax) {prev=j; continue;}
 722           if(anotherPt[j]->ymax < anotherPt[prev]->ymin) {prev=j; continue;}
 723           prev = j;
 724           */
 725
 726           top += fCurrentHisto->GetBinCenterX(j)*anotherPt[j]->weight;
 727           butt += anotherPt[j]->weight;
 728
 729           for(Int_t k=anotherPt[j]->ymin; k<anotherPt[j]->ymax; k++)
 730             {
 731               Int_t bin = fCurrentHisto->GetBin(j,k);
 732               ytop += (fCurrentHisto->GetBinCenterY(k))*(fCurrentHisto->GetBinContent(bin));
 733               ybutt += fCurrentHisto->GetBinContent(bin);
 734               w+=(Int_t)fCurrentHisto->GetBinContent(bin);
 735             }
 736         }
 737
 738       xpeaks[i] = top/butt;
 739       ypeaks[i] = ytop/ybutt;
 740       weight[i] = w;
 741       //cout<<"Setting weight "<<w<<" kappa "<<xpeaks[i]<<" phi0 "<<ypeaks[i]<<endl;
 742
 743       //Check if this peak is overlapping with a previous:
 744       for(Int_t p=0; p<i-1; p++)
 745         {
 746           if(fabs(xpeaks[p] - xpeaks[i]) < max_kappa ||
 747              fabs(ypeaks[p] - ypeaks[i]) < max_phi0)
 748             {
 749               weight[i]=0;
 750               break;
 751             }
 752         }
 753
 754     }
 755
 756   for(Int_t i=0; i<nbinsx; i++)
 757     delete windowPt[i];
 758   delete [] windowPt;
 759   delete [] anotherPt;
 760
 761 }
 762
 763 void AliL3HoughMaxFinder::SortPeaks(struct AxisWindow **a,Int_t first,Int_t last)
 764 {
 765   //General sorting routine
 766   //Sort according to PeakCompare()
 767
 768   static struct AxisWindow *tmp;
 769   static int i;           // "static" to save stack space
 770   int j;
 771
 772   while (last - first > 1) {
 773     i = first;
 774     j = last;
 775     for (;;) {
 776       while (++i < last && PeakCompare(a[i], a[first]) < 0)
 777         ;
 778       while (--j > first && PeakCompare(a[j], a[first]) > 0)
 779         ;
 780       if (i >= j)
 781         break;
 782
 783       tmp  = a[i];
 784       a[i] = a[j];
 785       a[j] = tmp;
 786     }
 787     if (j == first) {
 788       ++first;
 789       continue;
 790     }
 791     tmp = a[first];
 792     a[first] = a[j];
 793     a[j] = tmp;
 794     if (j - first < last - (j + 1)) {
 795       SortPeaks(a, first, j);
 796       first = j + 1;   // QSort(j + 1, last);
 797     } else {
 798       SortPeaks(a, j + 1, last);
 799       last = j;        // QSort(first, j);
 800     }
 801   }
 802
 803 }
 804
 805 Int_t AliL3HoughMaxFinder::PeakCompare(struct AxisWindow *a,struct AxisWindow *b)
 806 {
 807   if(a->weight < b->weight) return 1;
 808   if(a->weight > b->weight) return -1;
 809   return 0;
 810
 811 }
 812
 813 void AliL3HoughMaxFinder::FindPeak(Int_t t1,Double_t t2,Int_t t3,Float_t &kappa,Float_t &phi0)
 814 {
 815   //Attempt of a more sophisticated peak finder.
 816   //Finds the best peak in the histogram, and returns the corresponding
 817   //track object.
 818
 819   if(!fCurrentHisto)
 820     {
 821       printf("AliL3HoughMaxFinder::FindPeak : No histogram!!\n");
 822       return;
 823     }
 824   AliL3Histogram *hist = fCurrentHisto;
 825
 826   Int_t xmin = hist->GetFirstXbin();
 827   Int_t xmax = hist->GetLastXbin();
 828   Int_t ymin = hist->GetFirstYbin();
 829   Int_t ymax = hist->GetLastYbin();
 830   Int_t nbinsx = hist->GetNbinsX()+1;
 831
 832   Int_t *m = new Int_t[nbinsx];
 833   Int_t *m_low = new Int_t[nbinsx];
 834   Int_t *m_up = new Int_t[nbinsx];
 835
 836
 837  recompute:  //this is a goto.
 838
 839   for(Int_t i=0; i<nbinsx; i++)
 840     {
 841       m[i]=0;
 842       m_low[i]=0;
 843       m_up[i]=0;
 844     }
 845
 846   Int_t max_x=0,sum=0,max_xbin=0,bin;
 847
 848   for(Int_t xbin=xmin; xbin<=xmax; xbin++)
 849     {
 850       for(Int_t ybin=ymin; ybin <= ymax - t1; ybin++)
 851         {
 852           sum = 0;
 853           for(Int_t y=ybin; y <= ybin+t1; y++)
 854             {
 855               if(y>ymax) break;
 856               //Inside window
 857               bin = hist->GetBin(xbin,y);
 858               sum += (Int_t)hist->GetBinContent(bin);
 859
 860             }
 861           if(sum > m[xbin]) //Max value locally in this xbin
 862             {
 863               m[xbin]=sum;
 864               m_low[xbin]=ybin;
 865               m_up[xbin]=ybin + t1;
 866             }
 867
 868         }
 869
 870       if(m[xbin] > max_x) //Max value globally in x-direction
 871         {
 872           max_xbin = xbin;
 873           max_x = m[xbin];//sum;
 874         }
 875     }
 876   //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]);
 877   //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[max_xbin]),hist->GetBinCenterY(m_up[max_xbin]));
 878
 879   //Determine a width in the x-direction
 880   Int_t x_low=0,x_up=0;
 881
 882   for(Int_t xbin=max_xbin-1; xbin >= xmin; xbin--)
 883     {
 884       if(m[xbin] < max_x*t2)
 885         {
 886           x_low = xbin+1;
 887           break;
 888         }
 889     }
 890   for(Int_t xbin = max_xbin+1; xbin <=xmax; xbin++)
 891     {
 892       if(m[xbin] < max_x*t2)
 893         {
 894           x_up = xbin-1;
 895           break;
 896         }
 897     }
 898
 899   Double_t top=0,butt=0,value,x_peak;
 900   if(x_up - x_low + 1 > t3)
 901     {
 902       t1 -= 1;
 903       printf("\nxrange out if limit x_up %d x_low %d t1 %d\n\n",x_low,x_up,t1);
 904       if(t1 > 1)
 905         goto recompute;
 906       else
 907         {
 908           x_peak = hist->GetBinCenterX(max_xbin);
 909           goto moveon;
 910         }
 911     }
 912
 913   //printf("xlow %f xup %f\n",hist->GetBinCenterX(x_low),hist->GetBinCenterX(x_up));
 914   //printf("Spread in x %d\n",x_up-x_low +1);
 915
 916   //Now, calculate the center of mass in x-direction
 917   for(Int_t xbin=x_low; xbin <= x_up; xbin++)
 918     {
 919       value = hist->GetBinCenterX(xbin);
 920       top += value*m[xbin];
 921       butt += m[xbin];
 922     }
 923   x_peak = top/butt;
 924
 925  moveon:
 926
 927   //Find the peak in y direction:
 928   Int_t x_l = hist->FindXbin(x_peak);
 929   if(hist->GetBinCenterX(x_l) > x_peak)
 930     x_l--;
 931
 932   Int_t x_u = x_l + 1;
 933
 934   if(hist->GetBinCenterX(x_l) > x_peak || hist->GetBinCenterX(x_u) <= x_peak)
 935     printf("\nAliL3HoughMaxFinder::FindPeak : Wrong xrange %f %f %f\n\n",hist->GetBinCenterX(x_l),x_peak,hist->GetBinCenterX(x_u));
 936
 937     //printf("\nxlow %f xup %f\n",hist->GetBinCenterX(x_l),hist->GetBinCenterX(x_u));
 938
 939   value=top=butt=0;
 940
 941   //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_l]),hist->GetBinCenterY(m_up[x_l]));
 942   //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_u]),hist->GetBinCenterY(m_up[x_u]));
 943
 944   for(Int_t ybin=m_low[x_l]; ybin <= m_up[x_l]; ybin++)
 945     {
 946       value = hist->GetBinCenterY(ybin);
 947       bin = hist->GetBin(x_l,ybin);
 948       top += value*hist->GetBinContent(bin);
 949       butt += hist->GetBinContent(bin);
 950     }
 951   Double_t y_peak_low = top/butt;
 952
 953   //printf("y_peak_low %f\n",y_peak_low);
 954
 955   value=top=butt=0;
 956   for(Int_t ybin=m_low[x_u]; ybin <= m_up[x_u]; ybin++)
 957     {
 958       value = hist->GetBinCenterY(ybin);
 959       bin = hist->GetBin(x_u,ybin);
 960       top += value*hist->GetBinContent(bin);
 961       butt += hist->GetBinContent(bin);
 962     }
 963   Double_t y_peak_up = top/butt;
 964
 965   //printf("y_peak_up %f\n",y_peak_up);
 966
 967   Double_t x_value_up = hist->GetBinCenterX(x_u);
 968   Double_t x_value_low = hist->GetBinCenterX(x_l);
 969
 970   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);
 971
 972
 973   //Find the weight:
 974   //bin = hist->FindBin(x_peak,y_peak);
 975   //Int_t weight = (Int_t)hist->GetBinContent(bin);
 976
 977   //AliL3HoughTrack *track = new AliL3HoughTrack();
 978   //track->SetTrackParameters(x_peak,y_peak,weight);
 979   kappa = x_peak;
 980   phi0 = y_peak;
 981
 982   delete [] m;
 983   delete [] m_low;
 984   delete [] m_up;
 985
 986   //return track;
 987
 988
 989 }
 990