HLT/hough/AliL3HoughMaxFinder.cxx

   1 // @(#) $Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ALICE HLT Group
   5
   6 #include <strings.h>
   7 #include "AliL3StandardIncludes.h"
   8
   9 #ifndef no_root
  10 #include <TNtuple.h>
  11 #include <TFile.h>
  12 #endif
  13
  14 #include "AliL3Logging.h"
  15 #include "AliL3HoughMaxFinder.h"
  16 #include "AliL3Histogram.h"
  17 #include "AliL3TrackArray.h"
  18 #include "AliL3HoughTrack.h"
  19 #include "AliL3Transform.h"
  20 #include "AliL3HoughTransformerRow.h"
  21
  22 #if __GNUC__ == 3
  23 using namespace std;
  24 #endif
  25
  26 /** \class AliL3HoughMaxFinder
  27 <pre>
  28 //_____________________________________________________________
  29 // AliL3HoughMaxFinder
  30 //
  31 // Maximum finder
  32 //
  33 </pre>
  34 */
  35
  36 ClassImp(AliL3HoughMaxFinder)
  37
  38
  39 AliL3HoughMaxFinder::AliL3HoughMaxFinder()
  40 {
  41   //Default constructor
  42   fThreshold = 0;
  43   fCurrentEtaSlice = -1;
  44   fHistoType=0;
  45   fXPeaks=0;
  46   fYPeaks=0;
  47   fWeight=0;
  48   fNPeaks=0;
  49   fN1PeaksPrevEtaSlice=0;
  50   fN2PeaksPrevEtaSlice=0;
  51   fSTARTXPeaks=0;
  52   fSTARTYPeaks=0;
  53   fENDXPeaks=0;
  54   fENDYPeaks=0;
  55   fSTARTETAPeaks=0;
  56   fENDETAPeaks=0;
  57   fNMax=0;
  58   fGradX=1;
  59   fGradY=1;
  60 #ifndef no_root
  61   fNtuppel = 0;
  62 #endif
  63 }
  64
  65 AliL3HoughMaxFinder::AliL3HoughMaxFinder(Char_t *histotype,Int_t nmax,AliL3Histogram *hist)
  66 {
  67   //Constructor
  68
  69   //fTracks = new AliL3TrackArray("AliL3HoughTrack");
  70   if(strcmp(histotype,"KappaPhi")==0) fHistoType='c';
  71   if(strcmp(histotype,"DPsi")==0) fHistoType='l';
  72
  73   fCurrentEtaSlice = -1;
  74
  75   if(hist)
  76     fCurrentHisto = hist;
  77
  78   fGradX=1;
  79   fGradY=1;
  80   fNMax=nmax;
  81   fXPeaks = new Float_t[fNMax];
  82   fYPeaks = new Float_t[fNMax];
  83   fWeight = new Int_t[fNMax];
  84   fSTARTXPeaks = new Int_t[fNMax];
  85   fSTARTYPeaks = new Int_t[fNMax];
  86   fENDXPeaks = new Int_t[fNMax];
  87   fENDYPeaks = new Int_t[fNMax];
  88   fSTARTETAPeaks = new Int_t[fNMax];
  89   fENDETAPeaks = new Int_t[fNMax];
  90 #ifndef no_root
  91   fNtuppel = 0;
  92 #endif
  93   fThreshold=0;
  94 }
  95
  96 AliL3HoughMaxFinder::~AliL3HoughMaxFinder()
  97 {
  98   //Destructor
  99   if(fXPeaks)
 100     delete [] fXPeaks;
 101   if(fYPeaks)
 102     delete [] fYPeaks;
 103   if(fWeight)
 104     delete [] fWeight;
 105   if(fSTARTXPeaks)
 106     delete [] fSTARTXPeaks;
 107   if(fSTARTYPeaks)
 108     delete [] fSTARTYPeaks;
 109   if(fENDXPeaks)
 110     delete [] fENDXPeaks;
 111   if(fENDYPeaks)
 112     delete [] fENDYPeaks;
 113   if(fSTARTETAPeaks)
 114     delete [] fSTARTETAPeaks;
 115   if(fENDETAPeaks)
 116     delete [] fENDETAPeaks;
 117 #ifndef no_root
 118   if(fNtuppel)
 119     delete fNtuppel;
 120 #endif
 121 }
 122
 123 void AliL3HoughMaxFinder::Reset()
 124 {
 125   for(Int_t i=0; i<fNMax; i++)
 126     {
 127       fXPeaks[i]=0;
 128       fYPeaks[i]=0;
 129       fWeight[i]=0;
 130       fSTARTXPeaks[i]=0;
 131       fSTARTYPeaks[i]=0;
 132       fENDXPeaks[i]=0;
 133       fENDYPeaks[i]=0;
 134       fSTARTETAPeaks[i]=0;
 135       fENDETAPeaks[i]=0;
 136     }
 137   fNPeaks=0;
 138   fN1PeaksPrevEtaSlice=0;
 139   fN2PeaksPrevEtaSlice=0;
 140 }
 141
 142 void AliL3HoughMaxFinder::CreateNtuppel()
 143 {
 144 #ifndef no_root
 145   //content#; neighbouring bins of the peak.
 146   fNtuppel = new TNtuple("ntuppel","Peak charateristics","kappa:phi0:weigth:content3:content5:content1:content7");
 147   fNtuppel->SetDirectory(0);
 148 #endif
 149 }
 150
 151 void AliL3HoughMaxFinder::WriteNtuppel(Char_t *filename)
 152 {
 153 #ifndef no_root
 154   TFile *file = TFile::Open(filename,"RECREATE");
 155   if(!file)
 156     {
 157       cerr<<"AliL3HoughMaxFinder::WriteNtuppel : Error opening file "<<filename<<endl;
 158       return;
 159     }
 160   fNtuppel->Write();
 161   file->Close();
 162 #endif
 163 }
 164
 165 void AliL3HoughMaxFinder::FindAbsMaxima()
 166 {
 167
 168   if(!fCurrentHisto)
 169     {
 170       cerr<<"AliL3HoughMaxFinder::FindAbsMaxima : No histogram"<<endl;
 171       return;
 172     }
 173   AliL3Histogram *hist = fCurrentHisto;
 174
 175   if(hist->GetNEntries() == 0)
 176     return;
 177
 178   Int_t xmin = hist->GetFirstXbin();
 179   Int_t xmax = hist->GetLastXbin();
 180   Int_t ymin = hist->GetFirstYbin();
 181   Int_t ymax = hist->GetLastYbin();
 182   Int_t bin;
 183   Double_t value,max_value=0;
 184
 185   Int_t max_xbin=0,max_ybin=0;
 186   for(Int_t xbin=xmin; xbin<=xmax; xbin++)
 187     {
 188       for(Int_t ybin=ymin; ybin<=ymax; ybin++)
 189         {
 190           bin = hist->GetBin(xbin,ybin);
 191           value = hist->GetBinContent(bin);
 192           if(value>max_value)
 193             {
 194               max_value = value;
 195               max_xbin = xbin;
 196               max_ybin = ybin;
 197             }
 198         }
 199     }
 200
 201   if(max_value == 0)
 202     return;
 203
 204   if(fNPeaks > fNMax)
 205     {
 206       cerr<<"AliL3HoughMaxFinder::FindAbsMaxima : Array out of range : "<<fNPeaks<<endl;
 207       return;
 208     }
 209
 210   Double_t max_x = hist->GetBinCenterX(max_xbin);
 211   Double_t max_y = hist->GetBinCenterY(max_ybin);
 212   fXPeaks[fNPeaks] = max_x;
 213   fYPeaks[fNPeaks] = max_y;
 214   fWeight[fNPeaks] = (Int_t)max_value;
 215
 216   fNPeaks++;
 217 #ifndef no_root
 218   if(fNtuppel)
 219     {
 220       Int_t bin3 = hist->GetBin(max_xbin-1,max_ybin);
 221       Int_t bin5 = hist->GetBin(max_xbin+1,max_ybin);
 222       Int_t bin1 = hist->GetBin(max_xbin,max_ybin-1);
 223       Int_t bin7 = hist->GetBin(max_xbin,max_ybin+1);
 224
 225       fNtuppel->Fill(max_x,max_y,max_value,hist->GetBinContent(bin3),hist->GetBinContent(bin5),hist->GetBinContent(bin1),hist->GetBinContent(bin7));
 226     }
 227 #endif
 228 }
 229
 230 void AliL3HoughMaxFinder::FindBigMaxima()
 231 {
 232
 233   AliL3Histogram *hist = fCurrentHisto;
 234
 235   if(hist->GetNEntries() == 0)
 236     return;
 237
 238   Int_t xmin = hist->GetFirstXbin();
 239   Int_t xmax = hist->GetLastXbin();
 240   Int_t ymin = hist->GetFirstYbin();
 241   Int_t ymax = hist->GetLastYbin();
 242   Int_t bin[25],bin_index;
 243   Double_t value[25];
 244
 245   for(Int_t xbin=xmin+2; xbin<xmax-3; xbin++)
 246     {
 247       for(Int_t ybin=ymin+2; ybin<ymax-3; ybin++)
 248         {
 249           bin_index=0;
 250           for(Int_t xb=xbin-2; xb<=xbin+2; xb++)
 251             {
 252               for(Int_t yb=ybin-2; yb<=ybin+2; yb++)
 253                 {
 254                   bin[bin_index]=hist->GetBin(xb,yb);
 255                   value[bin_index]=hist->GetBinContent(bin[bin_index]);
 256                   bin_index++;
 257                 }
 258             }
 259           if(value[12]==0) continue;
 260           Int_t b=0;
 261           while(1)
 262             {
 263               if(value[b] > value[12] || b==bin_index) break;
 264               b++;
 265               //printf("b %d\n",b);
 266             }
 267           if(b == bin_index)
 268             {
 269               //Found maxima
 270               if(fNPeaks > fNMax)
 271                 {
 272                   cerr<<"AliL3HoughMaxFinder::FindBigMaxima : Array out of range "<<fNPeaks<<endl;
 273                   return;
 274                 }
 275
 276               Double_t max_x = hist->GetBinCenterX(xbin);
 277               Double_t max_y = hist->GetBinCenterY(ybin);
 278               fXPeaks[fNPeaks] = max_x;
 279               fYPeaks[fNPeaks] = max_y;
 280               fNPeaks++;
 281             }
 282         }
 283     }
 284 }
 285
 286 void AliL3HoughMaxFinder::FindMaxima(Int_t threshold)
 287 {
 288   //Locate all the maxima in input histogram.
 289   //Maxima is defined as bins with more entries than the
 290   //immediately neighbouring bins.
 291
 292   if(fCurrentHisto->GetNEntries() == 0)
 293     return;
 294
 295   Int_t xmin = fCurrentHisto->GetFirstXbin();
 296   Int_t xmax = fCurrentHisto->GetLastXbin();
 297   Int_t ymin = fCurrentHisto->GetFirstYbin();
 298   Int_t ymax = fCurrentHisto->GetLastYbin();
 299   Int_t bin[9];
 300   Double_t value[9];
 301
 302   //Float_t max_kappa = 0.001;
 303   //Float_t max_phi0 = 0.08;
 304
 305   for(Int_t xbin=xmin+1; xbin<=xmax-1; xbin++)
 306     {
 307       for(Int_t ybin=ymin+1; ybin<=ymax-1; ybin++)
 308         {
 309           bin[0] = fCurrentHisto->GetBin(xbin-1,ybin-1);
 310           bin[1] = fCurrentHisto->GetBin(xbin,ybin-1);
 311           bin[2] = fCurrentHisto->GetBin(xbin+1,ybin-1);
 312           bin[3] = fCurrentHisto->GetBin(xbin-1,ybin);
 313           bin[4] = fCurrentHisto->GetBin(xbin,ybin);
 314           bin[5] = fCurrentHisto->GetBin(xbin+1,ybin);
 315           bin[6] = fCurrentHisto->GetBin(xbin-1,ybin+1);
 316           bin[7] = fCurrentHisto->GetBin(xbin,ybin+1);
 317           bin[8] = fCurrentHisto->GetBin(xbin+1,ybin+1);
 318           value[0] = fCurrentHisto->GetBinContent(bin[0]);
 319           value[1] = fCurrentHisto->GetBinContent(bin[1]);
 320           value[2] = fCurrentHisto->GetBinContent(bin[2]);
 321           value[3] = fCurrentHisto->GetBinContent(bin[3]);
 322           value[4] = fCurrentHisto->GetBinContent(bin[4]);
 323           value[5] = fCurrentHisto->GetBinContent(bin[5]);
 324           value[6] = fCurrentHisto->GetBinContent(bin[6]);
 325           value[7] = fCurrentHisto->GetBinContent(bin[7]);
 326           value[8] = fCurrentHisto->GetBinContent(bin[8]);
 327
 328
 329           if(value[4]>value[0] && value[4]>value[1] && value[4]>value[2]
 330              && value[4]>value[3] && value[4]>value[5] && value[4]>value[6]
 331              && value[4]>value[7] && value[4]>value[8])
 332             {
 333               //Found a local maxima
 334               Float_t max_x = fCurrentHisto->GetBinCenterX(xbin);
 335               Float_t max_y = fCurrentHisto->GetBinCenterY(ybin);
 336
 337               if((Int_t)value[4] <= threshold) continue;//central bin below threshold
 338               if(fNPeaks >= fNMax)
 339                 {
 340                   cout<<"AliL3HoughMaxFinder::FindMaxima : Array out of range "<<fNPeaks<<endl;
 341                   return;
 342                 }
 343
 344               //Check the gradient:
 345               if(value[3]/value[4] > fGradX && value[5]/value[4] > fGradX)
 346                 continue;
 347
 348               if(value[1]/value[4] > fGradY && value[7]/value[4] > fGradY)
 349                 continue;
 350
 351               fXPeaks[fNPeaks] = max_x;
 352               fYPeaks[fNPeaks] = max_y;
 353               fWeight[fNPeaks] = (Int_t)value[4];
 354               fNPeaks++;
 355
 356               /*
 357               //Check if the peak is overlapping with a previous:
 358               Bool_t bigger = kFALSE;
 359               for(Int_t p=0; p<fNPeaks; p++)
 360                 {
 361                   if(fabs(max_x - fXPeaks[p]) < max_kappa && fabs(max_y - fYPeaks[p]) < max_phi0)
 362                     {
 363                       bigger = kTRUE;
 364                       if(value[4] > fWeight[p]) //this peak is bigger.
 365                         {
 366                           fXPeaks[p] = max_x;
 367                           fYPeaks[p] = max_y;
 368                           fWeight[p] = (Int_t)value[4];
 369                         }
 370                       else
 371                         continue; //previous peak is bigger.
 372                     }
 373                 }
 374               if(!bigger) //there were no overlapping peaks.
 375                 {
 376                   fXPeaks[fNPeaks] = max_x;
 377                   fYPeaks[fNPeaks] = max_y;
 378                   fWeight[fNPeaks] = (Int_t)value[4];
 379                   fNPeaks++;
 380                 }
 381               */
 382             }
 383         }
 384     }
 385
 386 }
 387
 388 struct Window
 389 {
 390   Int_t start;
 391   Int_t sum;
 392 };
 393
 394 void AliL3HoughMaxFinder::FindAdaptedPeaks(Int_t kappawindow,Float_t cut_ratio)
 395 {
 396   //Peak finder which looks for peaks with a certain shape.
 397   //The first step involves a pre-peak finder, which looks for peaks
 398   //in windows (size controlled by kappawindow) summing over each psi-bin.
 399   //These pre-preaks are then matched between neighbouring kappa-bins to
 400   //look for real 2D peaks exhbiting the typical cross-shape in the Hough circle transform.
 401   //The maximum bin within this region is marked as the peak itself, and
 402   //a few checks is performed to avoid the clear fake peaks (asymmetry check etc.)
 403
 404
 405   AliL3Histogram *hist = fCurrentHisto;
 406
 407   if(!hist)
 408     {
 409       cerr<<"AliL3HoughMaxFinder : No histogram!"<<endl;
 410       return;
 411     }
 412
 413   if(hist->GetNEntries() == 0)
 414     return;
 415
 416   Int_t xmin = hist->GetFirstXbin();
 417   Int_t xmax = hist->GetLastXbin();
 418   Int_t ymin = hist->GetFirstYbin();
 419   Int_t ymax = hist->GetLastYbin();
 420
 421
 422   //Start by looking for pre-peaks:
 423
 424   Window **local_maxima = new Window*[hist->GetNbinsY()];
 425
 426   Short_t *nmaxs = new Short_t[hist->GetNbinsY()];
 427   Int_t n,last_sum,sum;
 428   Bool_t sum_was_rising;
 429   for(Int_t ybin=ymin; ybin<=ymax; ybin++)
 430     {
 431       local_maxima[ybin-ymin] = new Window[hist->GetNbinsX()];
 432       nmaxs[ybin-ymin] = 0;
 433       sum_was_rising=0;
 434       last_sum=0;
 435       n=0;
 436       for(Int_t xbin=xmin; xbin<=xmax-kappawindow; xbin++)
 437         {
 438           sum=0;
 439           for(Int_t lbin=xbin; lbin<xbin+kappawindow; lbin++)
 440             sum += hist->GetBinContent(hist->GetBin(lbin,ybin));
 441
 442           if(sum < last_sum)
 443             {
 444               if(sum > fThreshold)
 445                 if(sum_was_rising)//Previous sum was a local maxima
 446                   {
 447                     local_maxima[ybin-ymin][nmaxs[ybin-ymin]].start = xbin-1;
 448                     local_maxima[ybin-ymin][nmaxs[ybin-ymin]].sum = last_sum;
 449                     nmaxs[ybin-ymin]++;
 450                   }
 451
 452               sum_was_rising=0;
 453             }
 454           else if(sum > 0)
 455             sum_was_rising=1;
 456           last_sum=sum;
 457         }
 458     }
 459
 460   Int_t match=0;
 461   Int_t *starts = new Int_t[hist->GetNbinsY()+1];
 462   Int_t *maxs = new Int_t[hist->GetNbinsY()+1];
 463
 464   for(Int_t ybin=ymax; ybin >= ymin+1; ybin--)
 465     {
 466       for(Int_t i=0; i<nmaxs[ybin-ymin]; i++)
 467         {
 468           Int_t lw = local_maxima[ybin-ymin][i].sum;
 469
 470           if(lw<0)
 471             continue; //already used
 472
 473           Int_t maxvalue=0,maxybin=0,maxxbin=0,maxwindow=0;
 474           for(Int_t k=local_maxima[ybin-ymin][i].start; k<local_maxima[ybin-ymin][i].start + kappawindow; k++)
 475             if(hist->GetBinContent(hist->GetBin(k,ybin)) > maxvalue)
 476               {
 477                 maxvalue = hist->GetBinContent(hist->GetBin(k,ybin));
 478                 maxybin = ybin;
 479                 maxxbin = k;
 480               }
 481
 482           //start expanding in the psi-direction:
 483
 484           Int_t lb = local_maxima[ybin-ymin][i].start;
 485           //Int_t ystart=ybin;
 486           starts[ybin] = local_maxima[ybin-ymin][i].start;
 487           maxs[ybin] = maxxbin;
 488           Int_t yl=ybin-1,nybins=1;
 489
 490           //cout<<"Starting search at ybin "<<ybin<<" start "<<lb<<" with sum "<<local_maxima[ybin-ymin][i].sum<<endl;
 491           while(yl >= ymin)
 492             {
 493               Bool_t found=0;
 494               for(Int_t j=0; j<nmaxs[yl-ymin]; j++)
 495                 {
 496                   if( local_maxima[yl-ymin][j].start - lb < 0) continue;
 497                   if( local_maxima[yl-ymin][j].start < lb + kappawindow + match &&
 498                       local_maxima[yl-ymin][j].start >= lb && local_maxima[yl-ymin][j].sum > 0)
 499                     {
 500
 501                       //cout<<"match at ybin "<<yl<<" yvalue "<<hist->GetBinCenterY(yl)<<" start "<<local_maxima[yl-ymin][j].start<<" sum "<<local_maxima[yl-ymin][j].sum<<endl;
 502
 503                       Int_t lmaxvalue=0,lmaxxbin=0;
 504                       for(Int_t k=local_maxima[yl-ymin][j].start; k<local_maxima[yl-ymin][j].start + kappawindow; k++)
 505                         {
 506                           if(hist->GetBinContent(hist->GetBin(k,yl)) > maxvalue)
 507                             {
 508                               maxvalue = hist->GetBinContent(hist->GetBin(k,yl));
 509                               maxxbin = k;
 510                               maxybin = yl;
 511                               maxwindow = j;
 512                             }
 513                           if(hist->GetBinContent(hist->GetBin(k,yl)) > lmaxvalue)//local maxima value
 514                             {
 515                               lmaxvalue=hist->GetBinContent(hist->GetBin(k,yl));
 516                               lmaxxbin=k;
 517                             }
 518                         }
 519                       nybins++;
 520                       starts[yl] = local_maxima[yl-ymin][j].start;
 521                       maxs[yl] = lmaxxbin;
 522                       local_maxima[yl-ymin][j].sum=-1; //Mark as used
 523                       found=1;
 524                       lb = local_maxima[yl-ymin][j].start;
 525                       break;//Since we found a match in this bin, we dont have to search it anymore, goto next bin.
 526                     }
 527                 }
 528               if(!found || yl == ymin)//no more local maximas to be matched, so write the final peak and break the expansion:
 529                 {
 530                   if(nybins > 4)
 531                     {
 532                       //cout<<"ystart "<<ystart<<" and nybins "<<nybins<<endl;
 533
 534                       Bool_t truepeak=kTRUE;
 535
 536                       //cout<<"Maxima found at xbin "<<maxxbin<<" ybin "<<maxybin<<" value "<<maxvalue<<endl;
 537                       //cout<<"Starting to sum at xbin "<<starts[maxybin-ymin]<<endl;
 538
 539
 540                       //Look in a window on both sides to probe the asymmetry
 541                       Float_t right=0,left=0;
 542                       for(Int_t w=maxxbin+1; w<=maxxbin+3; w++)
 543                         {
 544                           for(Int_t r=maxybin+1; r<=maxybin+3; r++)
 545                             {
 546                               right += (Float_t)hist->GetBinContent(hist->GetBin(w,r));
 547                             }
 548                         }
 549
 550                       for(Int_t w=maxxbin-1; w>=maxxbin-3; w--)
 551                         {
 552                           for(Int_t r=maxybin+1; r<=maxybin+3; r++)
 553                             {
 554                               left += (Float_t)hist->GetBinContent(hist->GetBin(w,r));
 555                             }
 556                         }
 557
 558                       //cout<<"ratio "<<right/left<<endl;
 559
 560                       Float_t upper_ratio=1,lower_ratio=1;
 561                       if(left)
 562                         upper_ratio = right/left;
 563
 564                       right=left=0;
 565                       for(Int_t w=maxxbin+1; w<=maxxbin+3; w++)
 566                         {
 567                           for(Int_t r=maxybin-1; r>=maxybin-3; r--)
 568                             {
 569                               right += (Float_t)hist->GetBinContent(hist->GetBin(w,r));
 570                             }
 571                         }
 572
 573                       for(Int_t w=maxxbin-1; w>=maxxbin-3; w--)
 574                         {
 575                           for(Int_t r=maxybin-1; r>=maxybin-3; r--)
 576                             {
 577                               left += (Float_t)hist->GetBinContent(hist->GetBin(w,r));
 578                             }
 579                         }
 580
 581                       //cout<<"ratio "<<left/right<<endl;
 582
 583                       if(right)
 584                         lower_ratio = left/right;
 585
 586                       if(upper_ratio > cut_ratio || lower_ratio > cut_ratio)
 587                         truepeak=kFALSE;
 588
 589                       if(truepeak)
 590                         {
 591
 592                           fXPeaks[fNPeaks] = hist->GetBinCenterX(maxxbin);
 593                           fYPeaks[fNPeaks] = hist->GetBinCenterY(maxybin);
 594                           fWeight[fNPeaks] = maxvalue;
 595                           fNPeaks++;
 596
 597                           /*
 598                           //Calculate the peak using weigthed means:
 599                           Float_t sum=0;
 600                           fYPeaks[fNPeaks]=0;
 601                           for(Int_t k=maxybin-1; k<=maxybin+1; k++)
 602                             {
 603                               Float_t lsum = 0;
 604                               for(Int_t l=starts[k]; l<starts[k]+kappawindow; l++)
 605                                 {
 606                                   lsum += (Float_t)hist->GetBinContent(hist->GetBin(l,k));
 607                                   sum += (Float_t)hist->GetBinContent(hist->GetBin(l,k));
 608                                 }
 609                               fYPeaks[fNPeaks] += lsum*hist->GetBinCenterY(k);
 610                             }
 611                           fYPeaks[fNPeaks] /= sum;
 612                           Int_t ybin1,ybin2;
 613                           if(fYPeaks[fNPeaks] < hist->GetBinCenterY(hist->FindYbin(fYPeaks[fNPeaks])))
 614                             {
 615                               ybin1 = hist->FindYbin(fYPeaks[fNPeaks])-1;
 616                               ybin2 = ybin1+1;
 617                             }
 618                           else
 619                             {
 620                               ybin1 = hist->FindYbin(fYPeaks[fNPeaks]);
 621                               ybin2 = ybin1+1;
 622                             }
 623
 624                           Float_t kappa1=0,kappa2=0;
 625                           sum=0;
 626                           for(Int_t k=starts[ybin1]; k<starts[ybin1] + kappawindow; k++)
 627                             {
 628                               kappa1 += hist->GetBinCenterX(k)*hist->GetBinContent(hist->GetBin(k,ybin1));
 629                               sum += (Float_t)hist->GetBinContent(hist->GetBin(k,ybin1));
 630                             }
 631                           kappa1 /= sum;
 632                           sum=0;
 633                           for(Int_t k=starts[ybin2]; k<starts[ybin2] + kappawindow; k++)
 634                             {
 635                               kappa2 += hist->GetBinCenterX(k)*hist->GetBinContent(hist->GetBin(k,ybin2));
 636                               sum += (Float_t)hist->GetBinContent(hist->GetBin(k,ybin2));
 637                             }
 638                           kappa2 /= sum;
 639
 640                           fXPeaks[fNPeaks] = ( kappa1*( hist->GetBinCenterY(ybin2) - fYPeaks[fNPeaks] ) +
 641                                                kappa2*( fYPeaks[fNPeaks] - hist->GetBinCenterY(ybin1) ) )  /
 642                             (hist->GetBinCenterY(ybin2) - hist->GetBinCenterY(ybin1));
 643
 644                           fNPeaks++;
 645                           */
 646                         }
 647                     }
 648                   break;
 649                 }
 650               else
 651                 yl--;//Search continues...
 652             }
 653         }
 654     }
 655
 656   for(Int_t i=0; i<hist->GetNbinsY(); i++)
 657     delete local_maxima[i];
 658
 659   delete [] local_maxima;
 660   delete [] nmaxs;
 661   delete [] starts;
 662   delete [] maxs;
 663 }
 664
 665 struct PreYPeak
 666 {
 667   Int_t start_position;
 668   Int_t end_position;
 669   Int_t min_value;
 670   Int_t max_value;
 671   Int_t prev_value;
 672   Int_t left_value;
 673   Int_t right_value;
 674 };
 675
 676 struct Pre2DPeak
 677 {
 678   Float_t x;
 679   Float_t y;
 680   Float_t size_x;
 681   Float_t size_y;
 682   Int_t start_x;
 683   Int_t start_y;
 684   Int_t end_x;
 685   Int_t end_y;
 686   Float_t weight;
 687 };
 688
 689 void AliL3HoughMaxFinder::FindAdaptedRowPeaks(Int_t kappawindow,Int_t xsize,Int_t ysize)
 690 {
 691
 692   AliL3Histogram *hist = fCurrentHisto;
 693
 694   if(!hist)
 695     {
 696       cerr<<"AliL3HoughMaxFinder : No histogram!"<<endl;
 697       return;
 698     }
 699
 700   if(hist->GetNEntries() == 0)
 701     return;
 702
 703   Int_t xmin = hist->GetFirstXbin();
 704   Int_t xmax = hist->GetLastXbin();
 705   Int_t ymin = hist->GetFirstYbin();
 706   Int_t ymax = hist->GetLastYbin();
 707
 708   //Start by looking for pre-peaks:
 709
 710   PreYPeak **local_maxima = new PreYPeak*[hist->GetNbinsY()];
 711
 712   Short_t *nmaxs = new Short_t[hist->GetNbinsY()];
 713   Int_t last_value=0,value=0;
 714   for(Int_t ybin=ymin; ybin<=ymax; ybin++)
 715     {
 716       local_maxima[ybin-ymin] = new PreYPeak[hist->GetNbinsX()];
 717       nmaxs[ybin-ymin] = 0;
 718       last_value = 0;
 719       Bool_t found = 0;
 720       for(Int_t xbin=xmin; xbin<=xmax; xbin++)
 721         {
 722           value = hist->GetBinContent(hist->GetBin(xbin,ybin));
 723           if(value > 0)
 724             {
 725               if(abs(value - last_value) > 1)
 726                 {
 727                   if(found) {
 728                     local_maxima[ybin-ymin][nmaxs[ybin-ymin]].right_value = value;
 729                     nmaxs[ybin-ymin]++;
 730                   }
 731                   local_maxima[ybin-ymin][nmaxs[ybin-ymin]].start_position = xbin;
 732                   local_maxima[ybin-ymin][nmaxs[ybin-ymin]].end_position = xbin;
 733                   local_maxima[ybin-ymin][nmaxs[ybin-ymin]].min_value = value;
 734                   local_maxima[ybin-ymin][nmaxs[ybin-ymin]].max_value = value;
 735                   local_maxima[ybin-ymin][nmaxs[ybin-ymin]].prev_value = 0;
 736                   local_maxima[ybin-ymin][nmaxs[ybin-ymin]].left_value = last_value;
 737                   found = 1;
 738                 }
 739               if(abs(value - last_value) <= 1)
 740                 {
 741                     local_maxima[ybin-ymin][nmaxs[ybin-ymin]].end_position = xbin;
 742                     if(value>local_maxima[ybin-ymin][nmaxs[ybin-ymin]].max_value)
 743                       local_maxima[ybin-ymin][nmaxs[ybin-ymin]].max_value = value;
 744                     if(value<local_maxima[ybin-ymin][nmaxs[ybin-ymin]].min_value)
 745                       local_maxima[ybin-ymin][nmaxs[ybin-ymin]].min_value = value;
 746                 }
 747             }
 748           last_value = value;
 749
 750         }
 751       if(found) {
 752         local_maxima[ybin-ymin][nmaxs[ybin-ymin]].right_value = value;
 753         nmaxs[ybin-ymin]++;
 754       }
 755     }
 756
 757   Pre2DPeak maxima[500];
 758   Int_t nmaxima = 0;
 759
 760   for(Int_t ybin=ymax; ybin >= ymin; ybin--)
 761     {
 762       for(Int_t i=0; i<nmaxs[ybin-ymin]; i++)
 763         {
 764           Int_t local_min_value = local_maxima[ybin-ymin][i].min_value;
 765           Int_t local_max_value = local_maxima[ybin-ymin][i].max_value;
 766           Int_t local_prev_value = local_maxima[ybin-ymin][i].prev_value;
 767           Int_t local_next_value = 0;
 768           Int_t local_left_value = local_maxima[ybin-ymin][i].left_value;
 769           Int_t local_right_value = local_maxima[ybin-ymin][i].right_value;
 770
 771           if(local_min_value<0)
 772             continue; //already used
 773
 774           //start expanding in the psi-direction:
 775
 776           Int_t local_x_start = local_maxima[ybin-ymin][i].start_position;
 777           Int_t local_x_end = local_maxima[ybin-ymin][i].end_position;
 778           Int_t temp_x_start = local_maxima[ybin-ymin][i].start_position;
 779           Int_t temp_x_end = local_maxima[ybin-ymin][i].end_position;
 780
 781           Int_t local_y=ybin-1,nybins=1;
 782
 783           while(local_y >= ymin)
 784             {
 785               Bool_t found=0;
 786               for(Int_t j=0; j<nmaxs[local_y-ymin]; j++)
 787                 {
 788                   if( (local_maxima[local_y-ymin][j].start_position <= (temp_x_end + kappawindow)) && (local_maxima[local_y-ymin][j].end_position >= (temp_x_start - kappawindow)))
 789                     {
 790                       if(((local_maxima[local_y-ymin][j].min_value <= local_max_value) && (local_maxima[local_y-ymin][j].min_value >= local_min_value)) ||
 791                          ((local_maxima[local_y-ymin][j].max_value >= local_min_value) && (local_maxima[local_y-ymin][j].max_value <= local_max_value)))
 792                         {
 793                           if(local_maxima[local_y-ymin][j].end_position > local_x_end)
 794                             local_x_end = local_maxima[local_y-ymin][j].end_position;
 795                           if(local_maxima[local_y-ymin][j].start_position < local_x_start)
 796                             local_x_start = local_maxima[local_y-ymin][j].start_position;
 797                           temp_x_start = local_maxima[local_y-ymin][j].start_position;
 798                           temp_x_end = local_maxima[local_y-ymin][j].end_position;
 799                           if(local_maxima[local_y-ymin][j].min_value < local_min_value)
 800                             local_min_value = local_maxima[local_y-ymin][j].min_value;
 801                           if(local_maxima[local_y-ymin][j].max_value > local_max_value)
 802                             local_max_value = local_maxima[local_y-ymin][j].max_value;
 803                           if(local_maxima[local_y-ymin][j].right_value > local_right_value)
 804                             local_right_value = local_maxima[local_y-ymin][j].right_value;
 805                           if(local_maxima[local_y-ymin][j].left_value > local_left_value)
 806                             local_left_value = local_maxima[local_y-ymin][j].left_value;
 807                           local_maxima[local_y-ymin][j].min_value = -1;
 808                           found = 1;
 809                           nybins++;
 810                           break;
 811                         }
 812                       else
 813                         {
 814                           if(local_max_value > local_maxima[local_y-ymin][j].prev_value)
 815                             local_maxima[local_y-ymin][j].prev_value = local_max_value;
 816                           if(local_maxima[local_y-ymin][j].max_value > local_next_value)
 817                             local_next_value = local_maxima[local_y-ymin][j].max_value;
 818                         }
 819                     }
 820                 }
 821               if(!found || local_y == ymin)//no more local maximas to be matched, so write the final peak and break the expansion:
 822                 {
 823                   if((nybins > ysize) && ((local_x_end-local_x_start+1) > xsize) && (local_max_value > local_prev_value) && (local_max_value > local_next_value) && (local_max_value > local_left_value) && (local_max_value > local_right_value))
 824                   //              if((nybins > ysize) && ((local_x_end-local_x_start+1) > xsize))
 825                     {
 826                       maxima[nmaxima].x = ((Float_t)local_x_start+(Float_t)local_x_end)/2.0;
 827                       maxima[nmaxima].y = ((Float_t)ybin+(Float_t)(local_y+1))/2.0;
 828                       maxima[nmaxima].size_x = local_x_end-local_x_start+1;
 829                       maxima[nmaxima].size_y = nybins;
 830                       maxima[nmaxima].weight = (local_min_value+local_max_value)/2;
 831                       maxima[nmaxima].start_x = local_x_start;
 832                       maxima[nmaxima].end_x = local_x_end;
 833                       maxima[nmaxima].start_y = local_y +1;
 834                       maxima[nmaxima].end_y = ybin;
 835 #ifdef do_mc
 836                       //                      cout<<"Peak found at: "<<((Float_t)local_x_start+(Float_t)local_x_end)/2.0<<" "<<((Float_t)ybin+(Float_t)(local_y+1))/2.0<<" "<<local_min_value<<" "<<local_max_value<<" "<<" with weight "<<(local_min_value+local_max_value)/2<<" and size "<<local_x_end-local_x_start+1<<" by "<<nybins<<endl;
 837 #endif
 838                       nmaxima++;
 839                     }
 840                   break;
 841                 }
 842               else
 843                 local_y--;//Search continues...
 844             }
 845         }
 846     }
 847
 848   //remove fake tracks
 849   Int_t nxbins = hist->GetNbinsX()+2;
 850
 851   for(Int_t i = 0; i < (nmaxima - 1); i++)
 852     {
 853       if(maxima[i].weight < 0) continue;
 854       for(Int_t j = i + 1; j < nmaxima; j++)
 855         {
 856           if(maxima[j].weight < 0) continue;
 857           Int_t xtrack1=0,xtrack2=0,ytrack1=0,ytrack2=0;
 858           Int_t deltax = 9999;
 859           for(Int_t ix1 = maxima[i].start_x; ix1 <= maxima[i].end_x; ix1++) {
 860             for(Int_t ix2 = maxima[j].start_x; ix2 <= maxima[j].end_x; ix2++) {
 861               if(abs(ix1 - ix2) < deltax) {
 862                 deltax = abs(ix1 - ix2);
 863                 xtrack1 = ix1;
 864                 xtrack2 = ix2;
 865               }
 866             }
 867           }
 868           Int_t deltay = 9999;
 869           for(Int_t iy1 = maxima[i].start_y; iy1 <= maxima[i].end_y; iy1++) {
 870             for(Int_t iy2 = maxima[j].start_y; iy2 <= maxima[j].end_y; iy2++) {
 871               if(abs(iy1 - iy2) < deltay) {
 872                 deltay = abs(iy1 - iy2);
 873                 ytrack1 = iy1;
 874                 ytrack2 = iy2;
 875               }
 876             }
 877           }
 878           Int_t firstrow1 = AliL3HoughTransformerRow::GetTrackFirstRow()[xtrack1 + nxbins*ytrack1];
 879           Int_t lastrow1 = AliL3HoughTransformerRow::GetTrackLastRow()[xtrack1 + nxbins*ytrack1];
 880           Int_t firstrow2 = AliL3HoughTransformerRow::GetTrackFirstRow()[xtrack1 + nxbins*ytrack1];
 881           Int_t lastrow2 = AliL3HoughTransformerRow::GetTrackLastRow()[xtrack1 + nxbins*ytrack1];
 882           Int_t firstrow,lastrow;
 883           if(firstrow1 < firstrow2)
 884             firstrow = firstrow2;
 885           else
 886             firstrow = firstrow1;
 887
 888           if(lastrow1 > lastrow2)
 889             lastrow = lastrow2;
 890           else
 891             lastrow = lastrow1;
 892
 893           AliL3HoughTrack track1;
 894           Float_t x1 = hist->GetPreciseBinCenterX(xtrack1);
 895           Float_t y1 = hist->GetPreciseBinCenterY(ytrack1);
 896           Float_t psi1 = atan((x1-y1)/(AliL3HoughTransformerRow::GetBeta1()-AliL3HoughTransformerRow::GetBeta2()));
 897           Float_t kappa1 = 2.0*(x1*cos(psi1)-AliL3HoughTransformerRow::GetBeta1()*sin(psi1));
 898           track1.SetTrackParameters(kappa1,psi1,1);
 899           Float_t firsthit1[3];
 900           if(!track1.GetCrossingPoint(firstrow,firsthit1)) continue;
 901           Float_t lasthit1[3];
 902           if(!track1.GetCrossingPoint(lastrow,lasthit1)) continue;
 903
 904           AliL3HoughTrack track2;
 905           Float_t x2 = hist->GetPreciseBinCenterX(xtrack2);
 906           Float_t y2 = hist->GetPreciseBinCenterY(ytrack2);
 907           Float_t psi2 = atan((x2-y2)/(AliL3HoughTransformerRow::GetBeta1()-AliL3HoughTransformerRow::GetBeta2()));
 908           Float_t kappa2 = 2.0*(x2*cos(psi2)-AliL3HoughTransformerRow::GetBeta1()*sin(psi2));
 909           track2.SetTrackParameters(kappa2,psi2,1);
 910           Float_t firsthit2[3];
 911           if(!track2.GetCrossingPoint(firstrow,firsthit2)) continue;
 912           Float_t lasthit2[3];
 913           if(!track2.GetCrossingPoint(lastrow,lasthit2)) continue;
 914
 915           Float_t padpitchlow = AliL3Transform::GetPadPitchWidth(AliL3Transform::GetPatch(firstrow));
 916           Float_t padpitchup = AliL3Transform::GetPadPitchWidth(AliL3Transform::GetPatch(lastrow));
 917           // check the distance between tracks at the edges
 918 #ifdef do_mc
 919           //      cout<<"DEBUG Merge peaks "<<i<<" "<<j<<" "<<firsthit1[1]<<" "<<firsthit2[1]<<"     "<<lasthit1[1]<<" "<<lasthit2[1]<<endl;
 920 #endif
 921           //cvetan please check!!! I added a cast to Int_t
 922           if((fabs(firsthit1[1]-firsthit2[1]) < 3.0*padpitchlow) && (fabs(lasthit1[1]-lasthit2[1]) < 3.0*padpitchup)) {
 923             if(maxima[i].size_x*maxima[i].size_y > maxima[j].size_x*maxima[j].size_y)
 924               maxima[j].weight = -maxima[j].weight;
 925             if(maxima[i].size_x*maxima[i].size_y < maxima[j].size_x*maxima[j].size_y)
 926               maxima[i].weight = -maxima[i].weight;
 927 #ifdef do_mc
 928             //      cout<<"Merge peaks "<<i<<" "<<j<<" "<<maxima[i].weight<<" "<<maxima[j].weight<<endl;
 929 #endif
 930           }
 931         }
 932     }
 933
 934   //merge tracks in neighbour eta slices
 935   /*
 936   for(Int_t i = 0; i < nmaxima; i++) {
 937     if(maxima[i].weight > 0) {
 938       fXPeaks[fNPeaks] = hist->GetPreciseBinCenterX(maxima[i].x);
 939       fYPeaks[fNPeaks] = hist->GetPreciseBinCenterY(maxima[i].y);
 940       fWeight[fNPeaks] = (Int_t)maxima[i].weight;
 941 #ifdef do_mc
 942       cout<<"Final Peak found at: "<<maxima[i].x<<" "<<maxima[i].y<<" "<<" "<<fXPeaks[fNPeaks]<<" "<<fYPeaks[fNPeaks]<<" with weight "<<fWeight[fNPeaks]<<" and size "<<maxima[i].size_x<<" by "<<maxima[i].size_y<<endl;
 943 #endif
 944       fNPeaks++;
 945     }
 946   }
 947   */
 948
 949   Int_t currentnpeaks = fNPeaks;
 950   for(Int_t i = 0; i < nmaxima; i++) {
 951     if(maxima[i].weight < 0) continue;
 952     Bool_t merged = kFALSE;
 953     for(Int_t j = fN1PeaksPrevEtaSlice; j < fN2PeaksPrevEtaSlice; j++) {
 954       if(fWeight[j] < 0) continue;
 955       if((fENDETAPeaks[j]-fSTARTETAPeaks[j]) >= 1) continue;
 956       if((maxima[i].start_x <= fENDXPeaks[j]+1) && (maxima[i].end_x >= fSTARTXPeaks[j]-1)) {
 957         if((maxima[i].start_y <= fENDYPeaks[j]+1) && (maxima[i].end_y >= fSTARTYPeaks[j]-1)) {
 958           //merge
 959           merged = kTRUE;
 960           fXPeaks[fNPeaks] = (hist->GetPreciseBinCenterX(maxima[i].x)+(fENDETAPeaks[j]-fSTARTETAPeaks[j]+1)*fXPeaks[j])/(fENDETAPeaks[j]-fSTARTETAPeaks[j]+2);
 961           fYPeaks[fNPeaks] = (hist->GetPreciseBinCenterY(maxima[i].y)+(fENDETAPeaks[j]-fSTARTETAPeaks[j]+1)*fYPeaks[j])/(fENDETAPeaks[j]-fSTARTETAPeaks[j]+2);
 962           fWeight[fNPeaks] = (Int_t)maxima[i].weight + fWeight[j];
 963           fSTARTXPeaks[fNPeaks] = maxima[i].start_x;
 964           fSTARTYPeaks[fNPeaks] = maxima[i].start_y;
 965           fENDXPeaks[fNPeaks] = maxima[i].end_x;
 966           fENDYPeaks[fNPeaks] = maxima[i].end_y;
 967           fSTARTETAPeaks[fNPeaks] = fSTARTETAPeaks[j];
 968           fENDETAPeaks[fNPeaks] = fCurrentEtaSlice;
 969           fNPeaks++;
 970           fWeight[j] = -fWeight[j];
 971         }
 972       }
 973     }
 974     fXPeaks[fNPeaks] = hist->GetPreciseBinCenterX(maxima[i].x);
 975     fYPeaks[fNPeaks] = hist->GetPreciseBinCenterY(maxima[i].y);
 976     if(!merged)
 977       fWeight[fNPeaks] = (Int_t)maxima[i].weight;
 978     else
 979       fWeight[fNPeaks] = -(Int_t)maxima[i].weight;
 980     fSTARTXPeaks[fNPeaks] = maxima[i].start_x;
 981     fSTARTYPeaks[fNPeaks] = maxima[i].start_y;
 982     fENDXPeaks[fNPeaks] = maxima[i].end_x;
 983     fENDYPeaks[fNPeaks] = maxima[i].end_y;
 984     fSTARTETAPeaks[fNPeaks] = fCurrentEtaSlice;
 985     fENDETAPeaks[fNPeaks] = fCurrentEtaSlice;
 986     fNPeaks++;
 987   }
 988   fN1PeaksPrevEtaSlice = currentnpeaks;
 989   fN2PeaksPrevEtaSlice = fNPeaks;
 990
 991   for(Int_t i=0; i<hist->GetNbinsY(); i++)
 992     delete local_maxima[i];
 993
 994   delete [] local_maxima;
 995   delete [] nmaxs;
 996 }
 997
 998 void AliL3HoughMaxFinder::FindPeak1(Int_t y_window,Int_t x_bin_sides)
 999 {
1000   //Testing mutliple peakfinding.
1001   //The algorithm searches the histogram for prepreaks by looking in windows
1002   //for each bin on the xaxis. The size of these windows is controlled by y_window.
1003   //Then the prepreaks are sorted according to their weight (sum inside window),
1004   //and the peak positions are calculated by taking the weighted mean in both
1005   //x and y direction. The size of the peak in x-direction is controlled by x_bin_sides.
1006
1007   if(!fCurrentHisto)
1008     {
1009       printf("AliL3HoughMaxFinder::FindPeak1 : No input histogram\n");
1010       return;
1011     }
1012   if(fCurrentHisto->GetNEntries()==0)
1013     return;
1014
1015   //Int_t y_window=2;
1016   //Int_t x_bin_sides=1;
1017
1018   //Float_t max_kappa = 0.001;
1019   //Float_t max_phi0 = 0.08;
1020
1021   Int_t max_sum=0;
1022
1023   Int_t xmin = fCurrentHisto->GetFirstXbin();
1024   Int_t xmax = fCurrentHisto->GetLastXbin();
1025   Int_t ymin = fCurrentHisto->GetFirstYbin();
1026   Int_t ymax = fCurrentHisto->GetLastYbin();
1027   Int_t nbinsx = fCurrentHisto->GetNbinsX()+1;
1028
1029   AxisWindow **windowPt = new AxisWindow*[nbinsx];
1030   AxisWindow **anotherPt = new AxisWindow*[nbinsx];
1031
1032   for(Int_t i=0; i<nbinsx; i++)
1033     {
1034       windowPt[i] = new AxisWindow;
1035 #if defined(__DECCXX)
1036       bzero((char *)windowPt[i],sizeof(AxisWindow));
1037 #else
1038       bzero((void*)windowPt[i],sizeof(AxisWindow));
1039 #endif
1040       anotherPt[i] = windowPt[i];
1041     }
1042
1043   for(Int_t xbin=xmin; xbin<=xmax; xbin++)
1044     {
1045       max_sum = 0;
1046       for(Int_t ybin=ymin; ybin<=ymax-y_window; ybin++)
1047         {
1048           Int_t sum_in_window=0;
1049           for(Int_t b=ybin; b<ybin+y_window; b++)
1050             {
1051               //inside window
1052               Int_t bin = fCurrentHisto->GetBin(xbin,b);
1053               sum_in_window += (Int_t)fCurrentHisto->GetBinContent(bin);
1054             }
1055
1056           if(sum_in_window > max_sum)
1057             {
1058               max_sum = sum_in_window;
1059               windowPt[xbin]->ymin = ybin;
1060               windowPt[xbin]->ymax = ybin + y_window;
1061               windowPt[xbin]->weight = sum_in_window;
1062               windowPt[xbin]->xbin = xbin;
1063             }
1064         }
1065     }
1066
1067   //Sort the windows according to the weight
1068   SortPeaks(windowPt,0,nbinsx);
1069
1070   Float_t top,butt;
1071   for(Int_t i=0; i<nbinsx; i++)
1072     {
1073       top=butt=0;
1074       Int_t xbin = windowPt[i]->xbin;
1075
1076       if(xbin<xmin || xbin > xmax-1) continue;
1077
1078       //Check if this is really a local maxima
1079       if(anotherPt[xbin-1]->weight > anotherPt[xbin]->weight ||
1080          anotherPt[xbin+1]->weight > anotherPt[xbin]->weight)
1081         continue;
1082
1083       for(Int_t j=windowPt[i]->ymin; j<windowPt[i]->ymax; j++)
1084         {
1085           //Calculate the mean in y direction:
1086           Int_t bin = fCurrentHisto->GetBin(windowPt[i]->xbin,j);
1087           top += (fCurrentHisto->GetBinCenterY(j))*(fCurrentHisto->GetBinContent(bin));
1088           butt += fCurrentHisto->GetBinContent(bin);
1089         }
1090
1091       if(butt < fThreshold)
1092         continue;
1093
1094       fXPeaks[fNPeaks] = fCurrentHisto->GetBinCenterX(windowPt[i]->xbin);
1095       fYPeaks[fNPeaks] = top/butt;
1096       fWeight[fNPeaks] = (Int_t)butt;
1097       //cout<<"mean in y "<<ypeaks[n]<<" on x "<<windowPt[i]->xbin<<" content "<<butt<<endl;
1098       fNPeaks++;
1099       if(fNPeaks==fNMax)
1100         {
1101           cerr<<"AliL3HoughMaxFinder::FindPeak1 : Peak array out of range!!!"<<endl;
1102           break;
1103         }
1104     }
1105
1106
1107   //Improve the peaks by including the region around in x.
1108   Float_t ytop,ybutt;
1109   Int_t prev;
1110   Int_t w;
1111   for(Int_t i=0; i<fNPeaks; i++)
1112     {
1113       Int_t xbin = fCurrentHisto->FindXbin(fXPeaks[i]);
1114       if(xbin - x_bin_sides < xmin || xbin + x_bin_sides > xmax) continue;
1115       top=butt=0;
1116       ytop=0,ybutt=0;
1117       w=0;
1118       prev = xbin - x_bin_sides+1;
1119       for(Int_t j=xbin-x_bin_sides; j<=xbin+x_bin_sides; j++)
1120         {
1121           /*
1122           //Check if the windows are overlapping:
1123           if(anotherPt[j]->ymin > anotherPt[prev]->ymax) {prev=j; continue;}
1124           if(anotherPt[j]->ymax < anotherPt[prev]->ymin) {prev=j; continue;}
1125           prev = j;
1126           */
1127
1128           top += fCurrentHisto->GetBinCenterX(j)*anotherPt[j]->weight;
1129           butt += anotherPt[j]->weight;
1130
1131           for(Int_t k=anotherPt[j]->ymin; k<anotherPt[j]->ymax; k++)
1132             {
1133               Int_t bin = fCurrentHisto->GetBin(j,k);
1134               ytop += (fCurrentHisto->GetBinCenterY(k))*(fCurrentHisto->GetBinContent(bin));
1135               ybutt += fCurrentHisto->GetBinContent(bin);
1136               w+=(Int_t)fCurrentHisto->GetBinContent(bin);
1137             }
1138         }
1139
1140       fXPeaks[i] = top/butt;
1141       fYPeaks[i] = ytop/ybutt;
1142       fWeight[i] = w;
1143       //cout<<"Setting weight "<<w<<" kappa "<<fXPeaks[i]<<" phi0 "<<fYPeaks[i]<<endl;
1144
1145       /*
1146       //Check if this peak is overlapping with a previous:
1147       for(Int_t p=0; p<i; p++)
1148         {
1149           //cout<<fabs(fXPeaks[p] - fXPeaks[i])<<" "<<fabs(fYPeaks[p] - fYPeaks[i])<<endl;
1150           if(fabs(fXPeaks[p] - fXPeaks[i]) < max_kappa &&
1151              fabs(fYPeaks[p] - fYPeaks[i]) < max_phi0)
1152             {
1153               fWeight[i]=0;
1154               //break;
1155             }
1156         }
1157       */
1158     }
1159
1160   for(Int_t i=0; i<nbinsx; i++)
1161     delete windowPt[i];
1162   delete [] windowPt;
1163   delete [] anotherPt;
1164 }
1165
1166 void AliL3HoughMaxFinder::SortPeaks(struct AxisWindow **a,Int_t first,Int_t last)
1167 {
1168   //General sorting routine
1169   //Sort according to PeakCompare()
1170
1171   static struct AxisWindow *tmp;
1172   static int i;           // "static" to save stack space
1173   int j;
1174
1175   while (last - first > 1) {
1176     i = first;
1177     j = last;
1178     for (;;) {
1179       while (++i < last && PeakCompare(a[i], a[first]) < 0)
1180         ;
1181       while (--j > first && PeakCompare(a[j], a[first]) > 0)
1182         ;
1183       if (i >= j)
1184         break;
1185
1186       tmp  = a[i];
1187       a[i] = a[j];
1188       a[j] = tmp;
1189     }
1190     if (j == first) {
1191       ++first;
1192       continue;
1193     }
1194     tmp = a[first];
1195     a[first] = a[j];
1196     a[j] = tmp;
1197     if (j - first < last - (j + 1)) {
1198       SortPeaks(a, first, j);
1199       first = j + 1;   // QSort(j + 1, last);
1200     } else {
1201       SortPeaks(a, j + 1, last);
1202       last = j;        // QSort(first, j);
1203     }
1204   }
1205
1206 }
1207
1208 Int_t AliL3HoughMaxFinder::PeakCompare(struct AxisWindow *a,struct AxisWindow *b)
1209 {
1210   if(a->weight < b->weight) return 1;
1211   if(a->weight > b->weight) return -1;
1212   return 0;
1213 }
1214
1215 void AliL3HoughMaxFinder::FindPeak(Int_t t1,Double_t t2,Int_t t3)
1216 {
1217   //Attempt of a more sophisticated peak finder.
1218   //Finds the best peak in the histogram, and returns the corresponding
1219   //track object.
1220
1221   if(!fCurrentHisto)
1222     {
1223       printf("AliL3HoughMaxFinder::FindPeak : No histogram!!\n");
1224       return;
1225     }
1226   AliL3Histogram *hist = fCurrentHisto;
1227   if(hist->GetNEntries()==0)
1228     return;
1229
1230   Int_t xmin = hist->GetFirstXbin();
1231   Int_t xmax = hist->GetLastXbin();
1232   Int_t ymin = hist->GetFirstYbin();
1233   Int_t ymax = hist->GetLastYbin();
1234   Int_t nbinsx = hist->GetNbinsX()+1;
1235
1236   Int_t *m = new Int_t[nbinsx];
1237   Int_t *m_low = new Int_t[nbinsx];
1238   Int_t *m_up = new Int_t[nbinsx];
1239
1240
1241  recompute:  //this is a goto.
1242
1243   for(Int_t i=0; i<nbinsx; i++)
1244     {
1245       m[i]=0;
1246       m_low[i]=0;
1247       m_up[i]=0;
1248     }
1249
1250   Int_t max_x=0,sum=0,max_xbin=0,bin=0;
1251
1252   for(Int_t xbin=xmin; xbin<=xmax; xbin++)
1253     {
1254       for(Int_t ybin=ymin; ybin <= ymax - t1; ybin++)
1255         {
1256           sum = 0;
1257           for(Int_t y=ybin; y <= ybin+t1; y++)
1258             {
1259               if(y>ymax) break;
1260               //Inside window
1261               bin = hist->GetBin(xbin,y);
1262               sum += (Int_t)hist->GetBinContent(bin);
1263
1264             }
1265           if(sum > m[xbin]) //Max value locally in this xbin
1266             {
1267               m[xbin]=sum;
1268               m_low[xbin]=ybin;
1269               m_up[xbin]=ybin + t1;
1270             }
1271
1272         }
1273
1274       if(m[xbin] > max_x) //Max value globally in x-direction
1275         {
1276           max_xbin = xbin;
1277           max_x = m[xbin];//sum;
1278         }
1279     }
1280   //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]);
1281   //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[max_xbin]),hist->GetBinCenterY(m_up[max_xbin]));
1282
1283   //Determine a width in the x-direction
1284   Int_t x_low=0,x_up=0;
1285
1286   for(Int_t xbin=max_xbin-1; xbin >= xmin; xbin--)
1287     {
1288       if(m[xbin] < max_x*t2)
1289         {
1290           x_low = xbin+1;
1291           break;
1292         }
1293     }
1294   for(Int_t xbin = max_xbin+1; xbin <=xmax; xbin++)
1295     {
1296       if(m[xbin] < max_x*t2)
1297         {
1298           x_up = xbin-1;
1299           break;
1300         }
1301     }
1302
1303   Double_t top=0,butt=0,value,x_peak;
1304   if(x_up - x_low + 1 > t3)
1305     {
1306       t1 -= 1;
1307       printf("\nxrange out if limit x_up %d x_low %d t1 %d\n\n",x_low,x_up,t1);
1308       if(t1 > 1)
1309         goto recompute;
1310       else
1311         {
1312           x_peak = hist->GetBinCenterX(max_xbin);
1313           goto moveon;
1314         }
1315     }
1316
1317   //printf("xlow %f xup %f\n",hist->GetBinCenterX(x_low),hist->GetBinCenterX(x_up));
1318   //printf("Spread in x %d\n",x_up-x_low +1);
1319
1320   //Now, calculate the center of mass in x-direction
1321   for(Int_t xbin=x_low; xbin <= x_up; xbin++)
1322     {
1323       value = hist->GetBinCenterX(xbin);
1324       top += value*m[xbin];
1325       butt += m[xbin];
1326     }
1327   x_peak = top/butt;
1328
1329  moveon:
1330
1331   //Find the peak in y direction:
1332   Int_t x_l = hist->FindXbin(x_peak);
1333   if(hist->GetBinCenterX(x_l) > x_peak)
1334     x_l--;
1335
1336   Int_t x_u = x_l + 1;
1337
1338   if(hist->GetBinCenterX(x_l) > x_peak || hist->GetBinCenterX(x_u) <= x_peak)
1339     printf("\nAliL3HoughMaxFinder::FindPeak : Wrong xrange %f %f %f\n\n",hist->GetBinCenterX(x_l),x_peak,hist->GetBinCenterX(x_u));
1340
1341     //printf("\nxlow %f xup %f\n",hist->GetBinCenterX(x_l),hist->GetBinCenterX(x_u));
1342
1343   value=top=butt=0;
1344
1345   //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_l]),hist->GetBinCenterY(m_up[x_l]));
1346   //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_u]),hist->GetBinCenterY(m_up[x_u]));
1347
1348   for(Int_t ybin=m_low[x_l]; ybin <= m_up[x_l]; ybin++)
1349     {
1350       value = hist->GetBinCenterY(ybin);
1351       bin = hist->GetBin(x_l,ybin);
1352       top += value*hist->GetBinContent(bin);
1353       butt += hist->GetBinContent(bin);
1354     }
1355   Double_t y_peak_low = top/butt;
1356
1357   //printf("y_peak_low %f\n",y_peak_low);
1358
1359   value=top=butt=0;
1360   for(Int_t ybin=m_low[x_u]; ybin <= m_up[x_u]; ybin++)
1361     {
1362       value = hist->GetBinCenterY(ybin);
1363       bin = hist->GetBin(x_u,ybin);
1364       top += value*hist->GetBinContent(bin);
1365       butt += hist->GetBinContent(bin);
1366     }
1367   Double_t y_peak_up = top/butt;
1368
1369   //printf("y_peak_up %f\n",y_peak_up);
1370
1371   Double_t x_value_up = hist->GetBinCenterX(x_u);
1372   Double_t x_value_low = hist->GetBinCenterX(x_l);
1373
1374   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);
1375
1376
1377   //Find the weight:
1378   //bin = hist->FindBin(x_peak,y_peak);
1379   //Int_t weight = (Int_t)hist->GetBinContent(bin);
1380
1381   //AliL3HoughTrack *track = new AliL3HoughTrack();
1382   //track->SetTrackParameters(x_peak,y_peak,weight);
1383   fXPeaks[fNPeaks]=x_peak;
1384   fYPeaks[fNPeaks]=y_peak;
1385   fWeight[fNPeaks]=(Int_t)hist->GetBinContent(bin);
1386   fNPeaks++;
1387
1388   delete [] m;
1389   delete [] m_low;
1390   delete [] m_up;
1391
1392   //return track;
1393 }
1394
1395 Float_t AliL3HoughMaxFinder::GetXPeakSize(Int_t i)
1396 {
1397   if(i<0 || i>fNMax)
1398     {
1399       STDCERR<<"AliL3HoughMaxFinder::GetXPeakSize : Invalid index "<<i<<STDENDL;
1400       return 0;
1401     }
1402   Float_t BinWidth = fCurrentHisto->GetBinWidthX();
1403   return BinWidth*(fENDXPeaks[i]-fSTARTXPeaks[i]+1);
1404 }
1405
1406 Float_t AliL3HoughMaxFinder::GetYPeakSize(Int_t i)
1407 {
1408   if(i<0 || i>fNMax)
1409     {
1410       STDCERR<<"AliL3HoughMaxFinder::GetYPeak : Invalid index "<<i<<STDENDL;
1411       return 0;
1412     }
1413   Float_t BinWidth = fCurrentHisto->GetBinWidthY();
1414   return BinWidth*(fENDYPeaks[i]-fSTARTYPeaks[i]+1);
1415 }