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