HLT/TPCLib/tracking/AliHLTTPCHoughMaxFinder.cxx

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