HLT/hough/AliL3HoughMaxFinder.cxx

   1 //$Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ASV
   5
   6 #include "AliL3StandardIncludes.h"
   7
   8 #ifndef no_root
   9 #include <TNtuple.h>
  10 #include <TFile.h>
  11 #endif
  12
  13 #include "AliL3Logging.h"
  14 #include "AliL3HoughMaxFinder.h"
  15 #include "AliL3Histogram.h"
  16 #include "AliL3TrackArray.h"
  17 #include "AliL3HoughTrack.h"
  18
  19 #if GCCVERSION == 3
  20 using namespace std;
  21 #endif
  22
  23 //_____________________________________________________________
  24 // AliL3HoughMaxFinder
  25 //
  26 // Maximum finder
  27
  28 ClassImp(AliL3HoughMaxFinder)
  29
  30
  31 AliL3HoughMaxFinder::AliL3HoughMaxFinder()
  32 {
  33   //Default constructor
  34   fThreshold = 0;
  35   fHistoType=0;
  36   fXPeaks=0;
  37   fYPeaks=0;
  38   fNPeaks=0;
  39   fNMax=0;
  40 #ifndef no_root
  41   fNtuppel = 0;
  42 #endif
  43 }
  44
  45
  46 AliL3HoughMaxFinder::AliL3HoughMaxFinder(Char_t *histotype,Int_t nmax,AliL3Histogram *hist)
  47 {
  48   //Constructor
  49
  50   //fTracks = new AliL3TrackArray("AliL3HoughTrack");
  51   if(strcmp(histotype,"KappaPhi")==0) fHistoType='c';
  52   if(strcmp(histotype,"DPsi")==0) fHistoType='l';
  53
  54   if(hist)
  55     fCurrentHisto = hist;
  56
  57   fNMax=nmax;
  58   fXPeaks = new Float_t[fNMax];
  59   fYPeaks = new Float_t[fNMax];
  60   fWeight = new Int_t[fNMax];
  61 #ifndef no_root
  62   fNtuppel = 0;
  63 #endif
  64   fThreshold=0;
  65 }
  66
  67
  68 AliL3HoughMaxFinder::~AliL3HoughMaxFinder()
  69 {
  70   //Destructor
  71   if(fXPeaks)
  72     delete [] fXPeaks;
  73   if(fYPeaks)
  74     delete [] fYPeaks;
  75   if(fWeight)
  76     delete [] fWeight;
  77 #ifndef no_root
  78   if(fNtuppel)
  79     delete fNtuppel;
  80 #endif
  81 }
  82
  83 void AliL3HoughMaxFinder::Reset()
  84 {
  85   for(Int_t i=0; i<fNMax; i++)
  86     {
  87       fXPeaks[i]=0;
  88       fYPeaks[i]=0;
  89       fWeight[i]=0;
  90     }
  91   fNPeaks=0;
  92 }
  93
  94 void AliL3HoughMaxFinder::CreateNtuppel()
  95 {
  96 #ifndef no_root
  97   //content#; neighbouring bins of the peak.
  98   fNtuppel = new TNtuple("ntuppel","Peak charateristics","kappa:phi0:weigth:content3:content5:content1:content7");
  99 #endif
 100 }
 101
 102 void AliL3HoughMaxFinder::WriteNtuppel(Char_t *filename)
 103 {
 104 #ifndef no_root
 105   TFile *file = TFile::Open(filename,"RECREATE");
 106   if(!file)
 107     {
 108       cerr<<"AliL3HoughMaxFinder::WriteNtuppel : Error opening file "<<filename<<endl;
 109       return;
 110     }
 111   fNtuppel->Write();
 112   file->Close();
 113 #endif
 114 }
 115
 116 void AliL3HoughMaxFinder::FindAbsMaxima()
 117 {
 118
 119   if(!fCurrentHisto)
 120     {
 121       cerr<<"AliL3HoughMaxFinder::FindAbsMaxima : No histogram"<<endl;
 122       return;
 123     }
 124   AliL3Histogram *hist = fCurrentHisto;
 125
 126   Int_t xmin = hist->GetFirstXbin();
 127   Int_t xmax = hist->GetLastXbin();
 128   Int_t ymin = hist->GetFirstYbin();
 129   Int_t ymax = hist->GetLastYbin();
 130   Int_t bin;
 131   Double_t value,max_value=0;
 132
 133   Int_t max_xbin=0,max_ybin=0;
 134   for(Int_t xbin=xmin; xbin<=xmax; xbin++)
 135     {
 136       for(Int_t ybin=ymin; ybin<=ymax; ybin++)
 137         {
 138           bin = hist->GetBin(xbin,ybin);
 139           value = hist->GetBinContent(bin);
 140           if(value>max_value)
 141             {
 142               max_value = value;
 143               max_xbin = xbin;
 144               max_ybin = ybin;
 145             }
 146         }
 147     }
 148
 149   if(fNPeaks > fNMax)
 150     {
 151       cerr<<"AliL3HoughMaxFinder::FindAbsMaxima : Array out of range : "<<fNPeaks<<endl;
 152       return;
 153     }
 154
 155
 156   Double_t max_x = hist->GetBinCenterX(max_xbin);
 157   Double_t max_y = hist->GetBinCenterY(max_ybin);
 158   fXPeaks[fNPeaks] = max_x;
 159   fYPeaks[fNPeaks] = max_y;
 160   fWeight[fNPeaks] = (Int_t)max_value;
 161   fNPeaks++;
 162 #ifndef no_root
 163   if(fNtuppel)
 164     {
 165       Int_t bin3 = hist->GetBin(max_xbin-1,max_ybin);
 166       Int_t bin5 = hist->GetBin(max_xbin+1,max_ybin);
 167       Int_t bin1 = hist->GetBin(max_xbin,max_ybin-1);
 168       Int_t bin7 = hist->GetBin(max_xbin,max_ybin+1);
 169
 170       fNtuppel->Fill(max_x,max_y,max_value,hist->GetBinContent(bin3),hist->GetBinContent(bin5),hist->GetBinContent(bin1),hist->GetBinContent(bin7));
 171     }
 172 #endif
 173 }
 174
 175 void AliL3HoughMaxFinder::FindBigMaxima()
 176 {
 177
 178   AliL3Histogram *hist = fCurrentHisto;
 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[25],bin_index;
 184   Double_t value[25];
 185
 186   for(Int_t xbin=xmin+2; xbin<xmax-3; xbin++)
 187     {
 188       for(Int_t ybin=ymin+2; ybin<ymax-3; ybin++)
 189         {
 190           bin_index=0;
 191           for(Int_t xb=xbin-2; xb<=xbin+2; xb++)
 192             {
 193               for(Int_t yb=ybin-2; yb<=ybin+2; yb++)
 194                 {
 195                   bin[bin_index]=hist->GetBin(xb,yb);
 196                   value[bin_index]=hist->GetBinContent(bin[bin_index]);
 197                   bin_index++;
 198                 }
 199
 200             }
 201           if(value[12]==0) continue;
 202           Int_t b=0;
 203           while(1)
 204             {
 205               if(value[b] > value[12] || b==bin_index) break;
 206               b++;
 207               //printf("b %d\n",b);
 208             }
 209           if(b == bin_index)
 210             {
 211               //Found maxima
 212               if(fNPeaks > fNMax)
 213                 {
 214                   cerr<<"AliL3HoughMaxFinder::FindBigMaxima : Array out of range "<<fNPeaks<<endl;
 215                   return;
 216                 }
 217
 218               Double_t max_x = hist->GetBinCenterX(xbin);
 219               Double_t max_y = hist->GetBinCenterY(ybin);
 220               fXPeaks[fNPeaks] = max_x;
 221               fYPeaks[fNPeaks] = max_y;
 222               fNPeaks++;
 223             }
 224         }
 225     }
 226 }
 227
 228 void AliL3HoughMaxFinder::FindMaxima(Double_t grad_x,Double_t grad_y)
 229 {
 230   //Locate all the maxima in input histogram.
 231   //Maxima is defined as bins with more entries than the
 232   //immediately neighbouring bins.
 233
 234   Int_t xmin = fCurrentHisto->GetFirstXbin();
 235   Int_t xmax = fCurrentHisto->GetLastXbin();
 236   Int_t ymin = fCurrentHisto->GetFirstYbin();
 237   Int_t ymax = fCurrentHisto->GetLastYbin();
 238   Int_t bin[9];
 239   Double_t value[9];
 240
 241   for(Int_t xbin=xmin+1; xbin<xmax-1; xbin++)
 242     {
 243       for(Int_t ybin=ymin+1; ybin<ymax-1; ybin++)
 244         {
 245           bin[0] = fCurrentHisto->GetBin(xbin-1,ybin-1);
 246           bin[1] = fCurrentHisto->GetBin(xbin,ybin-1);
 247           bin[2] = fCurrentHisto->GetBin(xbin+1,ybin-1);
 248           bin[3] = fCurrentHisto->GetBin(xbin-1,ybin);
 249           bin[4] = fCurrentHisto->GetBin(xbin,ybin);
 250           bin[5] = fCurrentHisto->GetBin(xbin+1,ybin);
 251           bin[6] = fCurrentHisto->GetBin(xbin-1,ybin+1);
 252           bin[7] = fCurrentHisto->GetBin(xbin,ybin+1);
 253           bin[8] = fCurrentHisto->GetBin(xbin+1,ybin+1);
 254           value[0] = fCurrentHisto->GetBinContent(bin[0]);
 255           value[1] = fCurrentHisto->GetBinContent(bin[1]);
 256           value[2] = fCurrentHisto->GetBinContent(bin[2]);
 257           value[3] = fCurrentHisto->GetBinContent(bin[3]);
 258           value[4] = fCurrentHisto->GetBinContent(bin[4]);
 259           value[5] = fCurrentHisto->GetBinContent(bin[5]);
 260           value[6] = fCurrentHisto->GetBinContent(bin[6]);
 261           value[7] = fCurrentHisto->GetBinContent(bin[7]);
 262           value[8] = fCurrentHisto->GetBinContent(bin[8]);
 263
 264
 265
 266           if(value[4]>value[0] && value[4]>value[1] && value[4]>value[2]
 267              && value[4]>value[3] && value[4]>value[5] && value[4]>value[6]
 268              && value[4]>value[7] && value[4]>value[8])
 269             {
 270               //Found a local maxima
 271               Float_t max_x = fCurrentHisto->GetBinCenterX(xbin);
 272               Float_t max_y = fCurrentHisto->GetBinCenterY(ybin);
 273
 274               if((Int_t)value[4] <= fThreshold) continue;//central bin below threshold
 275               if(fNPeaks >= fNMax)
 276                 {
 277                   cout<<"AliL3HoughMaxFinder::FindMaxima : Array out of range "<<fNPeaks<<endl;
 278                   return;
 279                 }
 280
 281               /*
 282               //Check the gradient:
 283               if(value[3]/value[4] > 1./grad_x || value[5]/value[4] < 1./grad_x ||
 284                  value[1]/value[4] < 1./grad_y || value[7]/value[4] < 1./grad_y)
 285                 continue;
 286               */
 287
 288
 289               fXPeaks[fNPeaks] = max_x;
 290               fYPeaks[fNPeaks] = max_y;
 291               fWeight[fNPeaks] = (Int_t)value[4];
 292               fNPeaks++;
 293
 294               /*
 295               //Check if the peak is overlapping with a previous:
 296               Bool_t bigger = kFALSE;
 297               for(Int_t p=0; p<entries; p++)
 298               {
 299                 if(fabs(max_x - xpeaks[p]) < kappa_overlap && fabs(max_y - ypeaks[p]) < phi_overlap)
 300               {
 301               bigger = kTRUE;
 302                       if(value[4] > weight[p]) //this peak is bigger.
 303                         {
 304                           xpeaks[p] = max_x;
 305                           ypeaks[p] = max_y;
 306                           weight[p] = (Int_t)value[4];
 307                         }
 308                       else
 309                         continue; //previous peak is bigger.
 310                     }
 311                 }
 312               if(!bigger) //there were no overlapping peaks.
 313                 {
 314                 xpeaks[entries] = max_x;
 315                   ypeaks[entries] = max_y;
 316                   weight[entries] = (Int_t)value[4];
 317                   entries++;
 318                 }
 319               */
 320             }
 321           else
 322             continue; //not a maxima
 323         }
 324     }
 325
 326 }
 327
 328 AliL3TrackArray *AliL3HoughMaxFinder::LookForPeaks(AliL3Histogram *hist,Int_t nbins)
 329 {
 330
 331   AliL3TrackArray *tracks = new AliL3TrackArray("AliL3HoughTrack");
 332   AliL3HoughTrack *track;
 333   Int_t xmin = hist->GetFirstXbin();
 334   Int_t xmax = hist->GetLastXbin();
 335   Int_t ymin = hist->GetFirstYbin();
 336   Int_t ymax = hist->GetLastYbin();
 337
 338   Int_t weight_loc;
 339   for(Int_t xbin=xmin+nbins; xbin <= xmax - nbins; xbin++)
 340     {
 341       for(Int_t ybin=ymin+nbins; ybin <= ymax - nbins; ybin++)
 342         {
 343           weight_loc=0;
 344           for(Int_t xbin_loc = xbin-nbins; xbin_loc <= xbin+nbins; xbin_loc++)
 345             {
 346               for(Int_t ybin_loc = ybin-nbins; ybin_loc <= ybin+nbins; ybin_loc++)
 347                 {
 348                   Int_t bin_loc = hist->GetBin(xbin_loc,ybin_loc);
 349                   weight_loc += (Int_t)hist->GetBinContent(bin_loc);
 350                 }
 351             }
 352
 353           if(weight_loc > 0)
 354             {
 355               track = (AliL3HoughTrack*)tracks->NextTrack();
 356               track->SetTrackParameters(hist->GetBinCenterX(xbin),hist->GetBinCenterY(ybin),weight_loc);
 357             }
 358
 359         }
 360     }
 361   tracks->QSort();
 362
 363   AliL3HoughTrack *track1,*track2;
 364
 365   for(Int_t i=1; i<tracks->GetNTracks(); i++)
 366     {
 367       track1 = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
 368       if(!track1) continue;
 369       Int_t xbin1 = hist->FindXbin(track1->GetKappa());
 370       Int_t ybin1 = hist->FindXbin(track1->GetPhi0());
 371       for(Int_t j=0; j < i; j++)
 372         {
 373           track2 = (AliL3HoughTrack*)tracks->GetCheckedTrack(j);
 374           if(!track2) continue;
 375           Int_t xbin2 = hist->FindXbin(track2->GetKappa());
 376           Int_t ybin2 = hist->FindYbin(track2->GetPhi0());
 377           if(abs(xbin1-xbin2) < 10 && abs(ybin1-ybin2) < 10)
 378             {
 379               tracks->Remove(i);
 380               break;
 381             }
 382         }
 383
 384     }
 385   tracks->Compress();
 386   return tracks;
 387 }
 388
 389 AliL3HoughTrack *AliL3HoughMaxFinder::CalculatePeakInWindow(Int_t *maxbin,Int_t t0,Int_t t1,Double_t t2,Int_t t3)
 390 {
 391   //Try to expand the area around the maxbin +- t0
 392
 393   if(!fCurrentHisto)
 394     {
 395       printf("AliL3HoughMaxFinder::LocatePeak : No histogram\n");
 396       return 0;
 397     }
 398   AliL3Histogram *hist = fCurrentHisto;
 399   Int_t xmin = hist->GetFirstXbin();
 400   Int_t xmax = hist->GetLastXbin();
 401   Int_t ymin = hist->GetFirstYbin();
 402   Int_t ymax = hist->GetLastYbin();
 403
 404   Int_t xlow = maxbin[0]-t0;
 405   if(xlow < xmin)
 406     xlow = xmin;
 407   Int_t xup = maxbin[0]+t0;
 408   if(xup > xmax)
 409     xup = xmax;
 410   Int_t ylow = maxbin[1]-t0;
 411   if(ylow < ymin)
 412     ylow = ymin;
 413   Int_t yup = maxbin[1]+t0;
 414   if(yup > ymax)
 415     yup = ymax;
 416
 417   Int_t nbinsx = hist->GetNbinsX()+1;
 418
 419   Int_t *m = new Int_t[nbinsx];
 420   Int_t *m_low = new Int_t[nbinsx];
 421   Int_t *m_up = new Int_t[nbinsx];
 422
 423  recompute:
 424
 425   Int_t max_x=0,sum=0,max_xbin=0,bin;
 426
 427   for(Int_t xbin=xlow; xbin<=xup; xbin++)
 428     {
 429       for(Int_t ybin=ylow; ybin <= yup; ybin++)
 430         {
 431           sum = 0;
 432           for(Int_t y=ybin; y <= ybin+t1; y++)
 433             {
 434               if(y>yup) break; //reached the upper limit in y.
 435               //Inside window
 436               bin = hist->GetBin(xbin,y);
 437               sum += (Int_t)hist->GetBinContent(bin);
 438
 439             }
 440           if(sum > m[xbin]) //Max value locally in this xbin
 441             {
 442               m[xbin]=sum;
 443               m_low[xbin]=ybin;
 444               m_up[xbin]=ybin + t1;
 445             }
 446
 447         }
 448
 449       if(m[xbin] > max_x) //Max value globally in x-direction
 450         {
 451           max_xbin = xbin;
 452           max_x = m[xbin];//sum;
 453         }
 454     }
 455   //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]);
 456   //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[max_xbin]),hist->GetBinCenterY(m_up[max_xbin]));
 457
 458   //Determine a width in the x-direction
 459   Int_t x_low=0,x_up=0;
 460   for(Int_t xbin=max_xbin-1; xbin >= xmin; xbin--)
 461     {
 462       if(m[xbin] < max_x*t2)
 463         {
 464           x_low = xbin+1;
 465           break;
 466         }
 467     }
 468   for(Int_t xbin = max_xbin+1; xbin <=xmax; xbin++)
 469     {
 470       if(m[xbin] < max_x*t2)
 471         {
 472           x_up = xbin-1;
 473           break;
 474         }
 475     }
 476   printf("x_low %d x_up %d\n",x_low,x_up);
 477
 478   Double_t top=0,butt=0,value,x_peak;
 479   if(x_up - x_low + 1 > t3)
 480     {
 481       t1 -= 1;
 482       printf("\nxrange out if limit x_up %d x_low %d\n\n",x_low,x_up);
 483       if(t1 > 1)
 484         goto recompute;
 485       else
 486         {
 487           x_peak = hist->GetBinCenterX(max_xbin);
 488           goto moveon;
 489         }
 490     }
 491
 492   //printf("xlow %f xup %f\n",hist->GetBinCenterX(x_low),hist->GetBinCenterX(x_up));
 493   //printf("Spread in x %d\n",x_up-x_low +1);
 494
 495   //Now, calculate the center of mass in x-direction
 496   for(Int_t xbin=x_low; xbin <= x_up; xbin++)
 497     {
 498       value = hist->GetBinCenterX(xbin);
 499       top += value*m[xbin];
 500       butt += m[xbin];
 501     }
 502   x_peak = top/butt;
 503
 504  moveon:
 505
 506   //Find the peak in y direction:
 507   Int_t x_l = hist->FindXbin(x_peak);
 508   if(hist->GetBinCenterX(x_l) > x_peak)
 509     x_l--;
 510
 511   Int_t x_u = x_l + 1;
 512
 513   if(hist->GetBinCenterX(x_l) > x_peak || hist->GetBinCenterX(x_u) <= x_peak)
 514     printf("\nAliL3HoughMaxFinder::FindPeak : Wrong xrange %f %f %f\n\n",hist->GetBinCenterX(x_l),x_peak,hist->GetBinCenterX(x_u));
 515
 516     //printf("\nxlow %f xup %f\n",hist->GetBinCenterX(x_l),hist->GetBinCenterX(x_u));
 517
 518   value=top=butt=0;
 519
 520   //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_l]),hist->GetBinCenterY(m_up[x_l]));
 521   //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_u]),hist->GetBinCenterY(m_up[x_u]));
 522
 523   for(Int_t ybin=m_low[x_l]; ybin <= m_up[x_l]; ybin++)
 524     {
 525       value = hist->GetBinCenterY(ybin);
 526       bin = hist->GetBin(x_l,ybin);
 527       top += value*hist->GetBinContent(bin);
 528       butt += hist->GetBinContent(bin);
 529     }
 530   Double_t y_peak_low = top/butt;
 531
 532   //printf("y_peak_low %f\n",y_peak_low);
 533
 534   value=top=butt=0;
 535   for(Int_t ybin=m_low[x_u]; ybin <= m_up[x_u]; ybin++)
 536     {
 537       value = hist->GetBinCenterY(ybin);
 538       bin = hist->GetBin(x_u,ybin);
 539       top += value*hist->GetBinContent(bin);
 540       butt += hist->GetBinContent(bin);
 541     }
 542   Double_t y_peak_up = top/butt;
 543
 544   //printf("y_peak_up %f\n",y_peak_up);
 545
 546   Double_t x_value_up = hist->GetBinCenterX(x_u);
 547   Double_t x_value_low = hist->GetBinCenterX(x_l);
 548
 549   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);
 550
 551
 552   //Find the weight:
 553   bin = hist->FindBin(x_peak,y_peak);
 554   Int_t weight = (Int_t)hist->GetBinContent(bin);
 555
 556   AliL3HoughTrack *track = new AliL3HoughTrack();
 557   track->SetTrackParameters(x_peak,y_peak,weight);
 558
 559   //Reset area around peak
 560   for(Int_t xbin=x_low; xbin<=x_up; xbin++)
 561     {
 562       for(Int_t ybin=m_low[xbin]; ybin<=m_up[xbin]; ybin++)
 563         {
 564           bin = hist->GetBin(xbin,ybin);
 565           hist->SetBinContent(bin,0);
 566         }
 567     }
 568
 569   delete [] m;
 570   delete [] m_low;
 571   delete [] m_up;
 572
 573   return track;
 574
 575
 576 }
 577
 578 AliL3HoughTrack *AliL3HoughMaxFinder::FindPeakLine(Double_t rho,Double_t theta)
 579 {
 580   //Peak finder based on a second line transformation on kappa-phi space,
 581   //to use as a baseline.
 582
 583   if(!fCurrentHisto)
 584     {
 585       printf("AliL3HoughTransformer::FindPeakLine : No input histogram\n");
 586       return 0;
 587     }
 588
 589   //get the line parameters:
 590   Double_t a = -1./tan(theta);
 591   Double_t b = rho/sin(theta);
 592
 593   printf("rho %f theta %f\n",rho,theta);
 594   //now, start looking along the line.
 595   Int_t xmin = fCurrentHisto->GetFirstXbin();
 596   Int_t xmax = fCurrentHisto->GetLastXbin();
 597
 598   Int_t max_weight=0;
 599   Double_t max_bin[2];
 600   for(Int_t xbin=xmin; xbin<=xmax; xbin++)
 601     {
 602       Double_t x = fCurrentHisto->GetBinCenterX(xbin);
 603       Double_t y = a*x + b;
 604       Int_t bin = fCurrentHisto->FindBin(x,y);
 605       //printf("x %f y %f weight %d\n",x,y,fCurrentHisto->GetBinContent(bin));
 606       if(fCurrentHisto->GetBinContent(bin) > max_weight)
 607         {
 608           max_weight = (Int_t)fCurrentHisto->GetBinContent(bin);
 609           max_bin[0] = x;
 610           max_bin[1] = y;
 611         }
 612     }
 613
 614   AliL3HoughTrack *track = new AliL3HoughTrack();
 615   track->SetTrackParameters(max_bin[0],max_bin[1],max_weight);
 616   return track;
 617 }
 618
 619 void AliL3HoughMaxFinder::FindPeak1(Int_t y_window,Int_t x_bin_sides)
 620 {
 621   //Testing mutliple peakfinding.
 622   //The algorithm searches the histogram for prepreaks by looking in windows
 623   //for each bin on the xaxis. The size of these windows is controlled by y_window.
 624   //Then the prepreaks are sorted according to their weight (sum inside window),
 625   //and the peak positions are calculated by taking the weighted mean in both
 626   //x and y direction. The size of the peak in x-direction is controlled by x_bin_sides.
 627
 628   if(!fCurrentHisto)
 629     {
 630       printf("AliL3HoughMaxFinder::FindPeak1 : No input histogram\n");
 631       return;
 632     }
 633   //Int_t y_window=2;
 634   //Int_t x_bin_sides=1;
 635
 636   //Float_t max_kappa = 0.001;
 637   //Float_t max_phi0 = 0.05;
 638
 639   Int_t max_sum=0;
 640
 641   Int_t xmin = fCurrentHisto->GetFirstXbin();
 642   Int_t xmax = fCurrentHisto->GetLastXbin();
 643   Int_t ymin = fCurrentHisto->GetFirstYbin();
 644   Int_t ymax = fCurrentHisto->GetLastYbin();
 645   Int_t nbinsx = fCurrentHisto->GetNbinsX()+1;
 646
 647   AxisWindow **windowPt = new AxisWindow*[nbinsx];
 648   AxisWindow **anotherPt = new AxisWindow*[nbinsx];
 649
 650   for(Int_t i=0; i<nbinsx; i++)
 651     {
 652       windowPt[i] = new AxisWindow;
 653       bzero((void*)windowPt[i],sizeof(AxisWindow));
 654       anotherPt[i] = windowPt[i];
 655     }
 656
 657   for(Int_t xbin=xmin; xbin<=xmax; xbin++)
 658     {
 659       max_sum = 0;
 660       for(Int_t ybin=ymin; ybin<=ymax-y_window; ybin++)
 661         {
 662           Int_t sum_in_window=0;
 663           for(Int_t b=ybin; b<ybin+y_window; b++)
 664             {
 665               //inside window
 666               Int_t bin = fCurrentHisto->GetBin(xbin,b);
 667               sum_in_window += (Int_t)fCurrentHisto->GetBinContent(bin);
 668             }
 669
 670           if(sum_in_window > max_sum)
 671             {
 672               max_sum = sum_in_window;
 673               windowPt[xbin]->ymin = ybin;
 674               windowPt[xbin]->ymax = ybin + y_window;
 675               windowPt[xbin]->weight = sum_in_window;
 676               windowPt[xbin]->xbin = xbin;
 677             }
 678         }
 679     }
 680
 681   //Sort the windows according to the weight
 682   SortPeaks(windowPt,0,nbinsx);
 683
 684   Float_t top,butt;
 685   for(Int_t i=0; i<nbinsx; i++)
 686     {
 687       top=butt=0;
 688       Int_t xbin = windowPt[i]->xbin;
 689
 690       if(xbin<xmin || xbin > xmax-1) continue;
 691
 692       //Check if this is really a local maxima
 693       if(anotherPt[xbin-1]->weight > anotherPt[xbin]->weight ||
 694          anotherPt[xbin+1]->weight > anotherPt[xbin]->weight)
 695         continue;
 696
 697       for(Int_t j=windowPt[i]->ymin; j<windowPt[i]->ymax; j++)
 698         {
 699           //Calculate the mean in y direction:
 700           Int_t bin = fCurrentHisto->GetBin(windowPt[i]->xbin,j);
 701           top += (fCurrentHisto->GetBinCenterY(j))*(fCurrentHisto->GetBinContent(bin));
 702           butt += fCurrentHisto->GetBinContent(bin);
 703         }
 704
 705       if(butt < fThreshold)
 706         continue;
 707
 708       fXPeaks[fNPeaks] = fCurrentHisto->GetBinCenterX(windowPt[i]->xbin);
 709       fYPeaks[fNPeaks] = top/butt;
 710       fWeight[fNPeaks] = (Int_t)butt;
 711       //cout<<"mean in y "<<ypeaks[n]<<" on x "<<windowPt[i]->xbin<<" content "<<butt<<endl;
 712       fNPeaks++;
 713       if(fNPeaks==fNMax) break;
 714     }
 715
 716
 717   //Improve the peaks by including the region around in x.
 718   Float_t ytop,ybutt;
 719   Int_t prev;
 720   Int_t w;
 721   for(Int_t i=0; i<fNPeaks; i++)
 722     {
 723       Int_t xbin = fCurrentHisto->FindXbin(fXPeaks[i]);
 724       if(xbin - x_bin_sides < xmin || xbin + x_bin_sides > xmax) continue;
 725       top=butt=0;
 726       ytop=0,ybutt=0;
 727       w=0;
 728       prev = xbin - x_bin_sides+1;
 729       for(Int_t j=xbin-x_bin_sides; j<=xbin+x_bin_sides; j++)
 730         {
 731           /*
 732           //Check if the windows are overlapping:
 733           if(anotherPt[j]->ymin > anotherPt[prev]->ymax) {prev=j; continue;}
 734           if(anotherPt[j]->ymax < anotherPt[prev]->ymin) {prev=j; continue;}
 735           prev = j;
 736           */
 737
 738           top += fCurrentHisto->GetBinCenterX(j)*anotherPt[j]->weight;
 739           butt += anotherPt[j]->weight;
 740
 741           for(Int_t k=anotherPt[j]->ymin; k<anotherPt[j]->ymax; k++)
 742             {
 743               Int_t bin = fCurrentHisto->GetBin(j,k);
 744               ytop += (fCurrentHisto->GetBinCenterY(k))*(fCurrentHisto->GetBinContent(bin));
 745               ybutt += fCurrentHisto->GetBinContent(bin);
 746               w+=(Int_t)fCurrentHisto->GetBinContent(bin);
 747             }
 748         }
 749
 750       fXPeaks[i] = top/butt;
 751       fYPeaks[i] = ytop/ybutt;
 752       fWeight[i] = w;
 753       //cout<<"Setting weight "<<w<<" kappa "<<xpeaks[i]<<" phi0 "<<ypeaks[i]<<endl;
 754
 755       /*
 756       //Check if this peak is overlapping with a previous:
 757       for(Int_t p=0; p<i-1; p++)
 758         {
 759           if(fabs(fXPeaks[p] - fXPeaks[i]) < max_kappa ||
 760              fabs(fYPeaks[p] - fYPeaks[i]) < max_phi0)
 761             {
 762               fWeight[i]=0;
 763               break;
 764             }
 765         }
 766       */
 767     }
 768
 769   for(Int_t i=0; i<nbinsx; i++)
 770     delete windowPt[i];
 771   delete [] windowPt;
 772   delete [] anotherPt;
 773
 774 }
 775
 776 void AliL3HoughMaxFinder::SortPeaks(struct AxisWindow **a,Int_t first,Int_t last)
 777 {
 778   //General sorting routine
 779   //Sort according to PeakCompare()
 780
 781   static struct AxisWindow *tmp;
 782   static int i;           // "static" to save stack space
 783   int j;
 784
 785   while (last - first > 1) {
 786     i = first;
 787     j = last;
 788     for (;;) {
 789       while (++i < last && PeakCompare(a[i], a[first]) < 0)
 790         ;
 791       while (--j > first && PeakCompare(a[j], a[first]) > 0)
 792         ;
 793       if (i >= j)
 794         break;
 795
 796       tmp  = a[i];
 797       a[i] = a[j];
 798       a[j] = tmp;
 799     }
 800     if (j == first) {
 801       ++first;
 802       continue;
 803     }
 804     tmp = a[first];
 805     a[first] = a[j];
 806     a[j] = tmp;
 807     if (j - first < last - (j + 1)) {
 808       SortPeaks(a, first, j);
 809       first = j + 1;   // QSort(j + 1, last);
 810     } else {
 811       SortPeaks(a, j + 1, last);
 812       last = j;        // QSort(first, j);
 813     }
 814   }
 815
 816 }
 817
 818 Int_t AliL3HoughMaxFinder::PeakCompare(struct AxisWindow *a,struct AxisWindow *b)
 819 {
 820   if(a->weight < b->weight) return 1;
 821   if(a->weight > b->weight) return -1;
 822   return 0;
 823
 824 }
 825
 826 void AliL3HoughMaxFinder::FindPeak(Int_t t1,Double_t t2,Int_t t3)
 827 {
 828   //Attempt of a more sophisticated peak finder.
 829   //Finds the best peak in the histogram, and returns the corresponding
 830   //track object.
 831
 832   if(!fCurrentHisto)
 833     {
 834       printf("AliL3HoughMaxFinder::FindPeak : No histogram!!\n");
 835       return;
 836     }
 837   AliL3Histogram *hist = fCurrentHisto;
 838
 839   Int_t xmin = hist->GetFirstXbin();
 840   Int_t xmax = hist->GetLastXbin();
 841   Int_t ymin = hist->GetFirstYbin();
 842   Int_t ymax = hist->GetLastYbin();
 843   Int_t nbinsx = hist->GetNbinsX()+1;
 844
 845   Int_t *m = new Int_t[nbinsx];
 846   Int_t *m_low = new Int_t[nbinsx];
 847   Int_t *m_up = new Int_t[nbinsx];
 848
 849
 850  recompute:  //this is a goto.
 851
 852   for(Int_t i=0; i<nbinsx; i++)
 853     {
 854       m[i]=0;
 855       m_low[i]=0;
 856       m_up[i]=0;
 857     }
 858
 859   Int_t max_x=0,sum=0,max_xbin=0,bin=0;
 860
 861   for(Int_t xbin=xmin; xbin<=xmax; xbin++)
 862     {
 863       for(Int_t ybin=ymin; ybin <= ymax - t1; ybin++)
 864         {
 865           sum = 0;
 866           for(Int_t y=ybin; y <= ybin+t1; y++)
 867             {
 868               if(y>ymax) break;
 869               //Inside window
 870               bin = hist->GetBin(xbin,y);
 871               sum += (Int_t)hist->GetBinContent(bin);
 872
 873             }
 874           if(sum > m[xbin]) //Max value locally in this xbin
 875             {
 876               m[xbin]=sum;
 877               m_low[xbin]=ybin;
 878               m_up[xbin]=ybin + t1;
 879             }
 880
 881         }
 882
 883       if(m[xbin] > max_x) //Max value globally in x-direction
 884         {
 885           max_xbin = xbin;
 886           max_x = m[xbin];//sum;
 887         }
 888     }
 889   //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]);
 890   //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[max_xbin]),hist->GetBinCenterY(m_up[max_xbin]));
 891
 892   //Determine a width in the x-direction
 893   Int_t x_low=0,x_up=0;
 894
 895   for(Int_t xbin=max_xbin-1; xbin >= xmin; xbin--)
 896     {
 897       if(m[xbin] < max_x*t2)
 898         {
 899           x_low = xbin+1;
 900           break;
 901         }
 902     }
 903   for(Int_t xbin = max_xbin+1; xbin <=xmax; xbin++)
 904     {
 905       if(m[xbin] < max_x*t2)
 906         {
 907           x_up = xbin-1;
 908           break;
 909         }
 910     }
 911
 912   Double_t top=0,butt=0,value,x_peak;
 913   if(x_up - x_low + 1 > t3)
 914     {
 915       t1 -= 1;
 916       printf("\nxrange out if limit x_up %d x_low %d t1 %d\n\n",x_low,x_up,t1);
 917       if(t1 > 1)
 918         goto recompute;
 919       else
 920         {
 921           x_peak = hist->GetBinCenterX(max_xbin);
 922           goto moveon;
 923         }
 924     }
 925
 926   //printf("xlow %f xup %f\n",hist->GetBinCenterX(x_low),hist->GetBinCenterX(x_up));
 927   //printf("Spread in x %d\n",x_up-x_low +1);
 928
 929   //Now, calculate the center of mass in x-direction
 930   for(Int_t xbin=x_low; xbin <= x_up; xbin++)
 931     {
 932       value = hist->GetBinCenterX(xbin);
 933       top += value*m[xbin];
 934       butt += m[xbin];
 935     }
 936   x_peak = top/butt;
 937
 938  moveon:
 939
 940   //Find the peak in y direction:
 941   Int_t x_l = hist->FindXbin(x_peak);
 942   if(hist->GetBinCenterX(x_l) > x_peak)
 943     x_l--;
 944
 945   Int_t x_u = x_l + 1;
 946
 947   if(hist->GetBinCenterX(x_l) > x_peak || hist->GetBinCenterX(x_u) <= x_peak)
 948     printf("\nAliL3HoughMaxFinder::FindPeak : Wrong xrange %f %f %f\n\n",hist->GetBinCenterX(x_l),x_peak,hist->GetBinCenterX(x_u));
 949
 950     //printf("\nxlow %f xup %f\n",hist->GetBinCenterX(x_l),hist->GetBinCenterX(x_u));
 951
 952   value=top=butt=0;
 953
 954   //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_l]),hist->GetBinCenterY(m_up[x_l]));
 955   //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_u]),hist->GetBinCenterY(m_up[x_u]));
 956
 957   for(Int_t ybin=m_low[x_l]; ybin <= m_up[x_l]; ybin++)
 958     {
 959       value = hist->GetBinCenterY(ybin);
 960       bin = hist->GetBin(x_l,ybin);
 961       top += value*hist->GetBinContent(bin);
 962       butt += hist->GetBinContent(bin);
 963     }
 964   Double_t y_peak_low = top/butt;
 965
 966   //printf("y_peak_low %f\n",y_peak_low);
 967
 968   value=top=butt=0;
 969   for(Int_t ybin=m_low[x_u]; ybin <= m_up[x_u]; ybin++)
 970     {
 971       value = hist->GetBinCenterY(ybin);
 972       bin = hist->GetBin(x_u,ybin);
 973       top += value*hist->GetBinContent(bin);
 974       butt += hist->GetBinContent(bin);
 975     }
 976   Double_t y_peak_up = top/butt;
 977
 978   //printf("y_peak_up %f\n",y_peak_up);
 979
 980   Double_t x_value_up = hist->GetBinCenterX(x_u);
 981   Double_t x_value_low = hist->GetBinCenterX(x_l);
 982
 983   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);
 984
 985
 986   //Find the weight:
 987   //bin = hist->FindBin(x_peak,y_peak);
 988   //Int_t weight = (Int_t)hist->GetBinContent(bin);
 989
 990   //AliL3HoughTrack *track = new AliL3HoughTrack();
 991   //track->SetTrackParameters(x_peak,y_peak,weight);
 992   fXPeaks[fNPeaks]=x_peak;
 993   fYPeaks[fNPeaks]=y_peak;
 994   fWeight[fNPeaks]=(Int_t)hist->GetBinContent(bin);
 995   fNPeaks++;
 996
 997   delete [] m;
 998   delete [] m_low;
 999   delete [] m_up;
1000
1001   //return track;
1002
1003
1004 }
1005