HLT/hough/AliL3HoughMaxFinder.cxx

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