HLT/hough/AliL3HoughMaxFinder.cxx

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