HLT/hough/AliL3HoughMaxFinder.cxx

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