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