HLT/hough/AliL3HoughMaxFinder.cxx

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