Compilation on Sun
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughMaxFinder.cxx
CommitLineData
3e87ef69 1// @(#) $Id$
b5a207b4 2
b1886074 3// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
3e87ef69 4//*-- Copyright &copy ALICE HLT Group
4de874d1 5
62bb4b3d 6#include <strings.h>
e06900d5 7#include "AliL3StandardIncludes.h"
8
2061190a 9#ifndef no_root
208b54c5 10#include <TNtuple.h>
11#include <TFile.h>
2061190a 12#endif
208b54c5 13
e06900d5 14#include "AliL3Logging.h"
15#include "AliL3HoughMaxFinder.h"
4cafa5fc 16#include "AliL3Histogram.h"
4de874d1 17#include "AliL3TrackArray.h"
18#include "AliL3HoughTrack.h"
e06900d5 19
0bd0c1ef 20#if __GNUC__ == 3
e06900d5 21using namespace std;
22#endif
4de874d1 23
0bd0c1ef 24/** \class AliL3HoughMaxFinder
25<pre>
b1886074 26//_____________________________________________________________
27// AliL3HoughMaxFinder
28//
29// Maximum finder
0bd0c1ef 30//
31</pre>
32*/
b1886074 33
4de874d1 34ClassImp(AliL3HoughMaxFinder)
35
36
37AliL3HoughMaxFinder::AliL3HoughMaxFinder()
38{
39 //Default constructor
40 fThreshold = 0;
4de874d1 41 fHistoType=0;
3fe49b5b 42 fXPeaks=0;
43 fYPeaks=0;
44 fNPeaks=0;
45 fNMax=0;
b2a02bce 46 fGradX=1;
47 fGradY=1;
2061190a 48#ifndef no_root
208b54c5 49 fNtuppel = 0;
2061190a 50#endif
4de874d1 51}
52
3fe49b5b 53AliL3HoughMaxFinder::AliL3HoughMaxFinder(Char_t *histotype,Int_t nmax,AliL3Histogram *hist)
4de874d1 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
4cafa5fc 61 if(hist)
62 fCurrentHisto = hist;
3fe49b5b 63
b2a02bce 64 fGradX=1;
65 fGradY=1;
3fe49b5b 66 fNMax=nmax;
67 fXPeaks = new Float_t[fNMax];
68 fYPeaks = new Float_t[fNMax];
69 fWeight = new Int_t[fNMax];
2061190a 70#ifndef no_root
208b54c5 71 fNtuppel = 0;
2061190a 72#endif
208b54c5 73 fThreshold=0;
4de874d1 74}
75
4de874d1 76AliL3HoughMaxFinder::~AliL3HoughMaxFinder()
77{
78 //Destructor
3fe49b5b 79 if(fXPeaks)
80 delete [] fXPeaks;
81 if(fYPeaks)
82 delete [] fYPeaks;
83 if(fWeight)
84 delete [] fWeight;
2061190a 85#ifndef no_root
208b54c5 86 if(fNtuppel)
87 delete fNtuppel;
2061190a 88#endif
3fe49b5b 89}
4de874d1 90
3fe49b5b 91void 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;
4de874d1 100}
101
208b54c5 102void AliL3HoughMaxFinder::CreateNtuppel()
103{
2061190a 104#ifndef no_root
208b54c5 105 //content#; neighbouring bins of the peak.
2061190a 106 fNtuppel = new TNtuple("ntuppel","Peak charateristics","kappa:phi0:weigth:content3:content5:content1:content7");
3e87ef69 107 fNtuppel->SetDirectory(0);
2061190a 108#endif
208b54c5 109}
110
111void AliL3HoughMaxFinder::WriteNtuppel(Char_t *filename)
112{
2061190a 113#ifndef no_root
208b54c5 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();
2061190a 122#endif
208b54c5 123}
124
3fe49b5b 125void AliL3HoughMaxFinder::FindAbsMaxima()
52a2a604 126{
208b54c5 127
4cafa5fc 128 if(!fCurrentHisto)
129 {
208b54c5 130 cerr<<"AliL3HoughMaxFinder::FindAbsMaxima : No histogram"<<endl;
4fc9a6a4 131 return;
4cafa5fc 132 }
133 AliL3Histogram *hist = fCurrentHisto;
4fc9a6a4 134
b2a02bce 135 if(hist->GetNEntries() == 0)
136 return;
137
4cafa5fc 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;
95a00d93 143 Double_t value,max_value=0;
3fe49b5b 144
145 Int_t max_xbin=0,max_ybin=0;
4cafa5fc 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;
4fc9a6a4 155 max_xbin = xbin;
156 max_ybin = ybin;
4cafa5fc 157 }
158 }
159 }
4fc9a6a4 160
b2a02bce 161 if(max_value == 0)
162 return;
163
3fe49b5b 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;
b2a02bce 175
3fe49b5b 176 fNPeaks++;
2061190a 177#ifndef no_root
208b54c5 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 }
2061190a 187#endif
4cafa5fc 188}
189
3fe49b5b 190void AliL3HoughMaxFinder::FindBigMaxima()
4cafa5fc 191{
192
3fe49b5b 193 AliL3Histogram *hist = fCurrentHisto;
b2a02bce 194
195 if(hist->GetNEntries() == 0)
196 return;
197
4cafa5fc 198 Int_t xmin = hist->GetFirstXbin();
199 Int_t xmax = hist->GetLastXbin();
200 Int_t ymin = hist->GetFirstYbin();
201 Int_t ymax = hist->GetLastYbin();
52a2a604 202 Int_t bin[25],bin_index;
95a00d93 203 Double_t value[25];
52a2a604 204
52a2a604 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;
4cafa5fc 210 for(Int_t xb=xbin-2; xb<=xbin+2; xb++)
52a2a604 211 {
4cafa5fc 212 for(Int_t yb=ybin-2; yb<=ybin+2; yb++)
52a2a604 213 {
214 bin[bin_index]=hist->GetBin(xb,yb);
215 value[bin_index]=hist->GetBinContent(bin[bin_index]);
216 bin_index++;
217 }
52a2a604 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++;
3fe49b5b 225 //printf("b %d\n",b);
52a2a604 226 }
227 if(b == bin_index)
228 {
229 //Found maxima
3fe49b5b 230 if(fNPeaks > fNMax)
231 {
232 cerr<<"AliL3HoughMaxFinder::FindBigMaxima : Array out of range "<<fNPeaks<<endl;
233 return;
234 }
235
4cafa5fc 236 Double_t max_x = hist->GetBinCenterX(xbin);
237 Double_t max_y = hist->GetBinCenterY(ybin);
3fe49b5b 238 fXPeaks[fNPeaks] = max_x;
239 fYPeaks[fNPeaks] = max_y;
240 fNPeaks++;
52a2a604 241 }
242 }
243 }
52a2a604 244}
245
b2a02bce 246void AliL3HoughMaxFinder::FindMaxima(Int_t threshold)
4de874d1 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
b2a02bce 252 if(fCurrentHisto->GetNEntries() == 0)
253 return;
254
4ab9f8f0 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];
95a00d93 260 Double_t value[9];
4de874d1 261
3e87ef69 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++)
4de874d1 266 {
3e87ef69 267 for(Int_t ybin=ymin+1; ybin<=ymax-1; ybin++)
4de874d1 268 {
4ab9f8f0 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
4de874d1 288
4de874d1 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
4ab9f8f0 294 Float_t max_x = fCurrentHisto->GetBinCenterX(xbin);
295 Float_t max_y = fCurrentHisto->GetBinCenterY(ybin);
4de874d1 296
b2a02bce 297 if((Int_t)value[4] <= threshold) continue;//central bin below threshold
0309a5ee 298 if(fNPeaks >= fNMax)
4de874d1 299 {
3060f008 300 cout<<"AliL3HoughMaxFinder::FindMaxima : Array out of range "<<fNPeaks<<endl;
4ab9f8f0 301 return;
4de874d1 302 }
3e87ef69 303
3fe49b5b 304 //Check the gradient:
b2a02bce 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
3fe49b5b 311 fXPeaks[fNPeaks] = max_x;
312 fYPeaks[fNPeaks] = max_y;
313 fWeight[fNPeaks] = (Int_t)value[4];
314 fNPeaks++;
3e87ef69 315
3fe49b5b 316 /*
4ab9f8f0 317 //Check if the peak is overlapping with a previous:
318 Bool_t bigger = kFALSE;
3e87ef69 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.
4ab9f8f0 325 {
3e87ef69 326 fXPeaks[p] = max_x;
327 fYPeaks[p] = max_y;
328 fWeight[p] = (Int_t)value[4];
4ab9f8f0 329 }
330 else
331 continue; //previous peak is bigger.
332 }
333 }
334 if(!bigger) //there were no overlapping peaks.
335 {
3e87ef69 336 fXPeaks[fNPeaks] = max_x;
337 fYPeaks[fNPeaks] = max_y;
338 fWeight[fNPeaks] = (Int_t)value[4];
339 fNPeaks++;
4ab9f8f0 340 }
3fe49b5b 341 */
4de874d1 342 }
4de874d1 343 }
344 }
4ab9f8f0 345
4de874d1 346}
e4c21048 347
b2a02bce 348struct Window
349{
350 Int_t start;
351 Int_t sum;
352};
353
045549b7 354void AliL3HoughMaxFinder::FindAdaptedPeaks(Int_t kappawindow,Float_t cut_ratio)
52a2a604 355{
b2a02bce 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.)
52a2a604 363
52a2a604 364
b2a02bce 365 AliL3Histogram *hist = fCurrentHisto;
52a2a604 366
b2a02bce 367 if(!hist)
52a2a604 368 {
b2a02bce 369 cerr<<"AliL3HoughMaxFinder : No histogram!"<<endl;
370 return;
52a2a604 371 }
b2a02bce 372
373 if(hist->GetNEntries() == 0)
374 return;
0bd0c1ef 375
4cafa5fc 376 Int_t xmin = hist->GetFirstXbin();
377 Int_t xmax = hist->GetLastXbin();
378 Int_t ymin = hist->GetFirstYbin();
379 Int_t ymax = hist->GetLastYbin();
0bd0c1ef 380
381
b2a02bce 382 //Start by looking for pre-peaks:
4cafa5fc 383
b2a02bce 384 Window **local_maxima = new Window*[hist->GetNbinsY()];
4cafa5fc 385
b2a02bce 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++)
52a2a604 390 {
b2a02bce 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++)
52a2a604 397 {
b2a02bce 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)
52a2a604 403 {
b2a02bce 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 }
52a2a604 411
b2a02bce 412 sum_was_rising=0;
52a2a604 413 }
b2a02bce 414 else if(sum > 0)
415 sum_was_rising=1;
416 last_sum=sum;
4cafa5fc 417 }
418 }
52a2a604 419
b2a02bce 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];
52a2a604 423
b2a02bce 424 for(Int_t ybin=ymax; ybin >= ymin+1; ybin--)
4cafa5fc 425 {
b2a02bce 426 for(Int_t i=0; i<nmaxs[ybin-ymin]; i++)
4cafa5fc 427 {
b2a02bce 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 }
4cafa5fc 613 }
614 }
4cafa5fc 615
b2a02bce 616 for(Int_t i=0; i<hist->GetNbinsY(); i++)
617 delete local_maxima[i];
52a2a604 618
b2a02bce 619 delete [] local_maxima;
620 delete [] nmaxs;
621 delete [] starts;
622 delete [] maxs;
4fc9a6a4 623}
624
0bd0c1ef 625struct PreYPeak
626{
627 Int_t start_position;
628 Int_t end_position;
629 Int_t value;
630 Int_t prev_value;
631};
632
633struct 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
642void 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}
b2a02bce 808
afd8fed4 809void AliL3HoughMaxFinder::FindPeak1(Int_t y_window,Int_t x_bin_sides)
4fc9a6a4 810{
b1886074 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
4fc9a6a4 818 if(!fCurrentHisto)
819 {
820 printf("AliL3HoughMaxFinder::FindPeak1 : No input histogram\n");
821 return;
822 }
b2a02bce 823 if(fCurrentHisto->GetNEntries()==0)
824 return;
825
b1886074 826 //Int_t y_window=2;
827 //Int_t x_bin_sides=1;
828
afd8fed4 829 //Float_t max_kappa = 0.001;
3e87ef69 830 //Float_t max_phi0 = 0.08;
b1886074 831
4fc9a6a4 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 bzero((void*)windowPt[i],sizeof(AxisWindow));
847 anotherPt[i] = windowPt[i];
848 }
849
850 for(Int_t xbin=xmin; xbin<=xmax; xbin++)
851 {
852 max_sum = 0;
853 for(Int_t ybin=ymin; ybin<=ymax-y_window; ybin++)
854 {
855 Int_t sum_in_window=0;
856 for(Int_t b=ybin; b<ybin+y_window; b++)
857 {
858 //inside window
859 Int_t bin = fCurrentHisto->GetBin(xbin,b);
860 sum_in_window += (Int_t)fCurrentHisto->GetBinContent(bin);
861 }
862
863 if(sum_in_window > max_sum)
864 {
865 max_sum = sum_in_window;
866 windowPt[xbin]->ymin = ybin;
867 windowPt[xbin]->ymax = ybin + y_window;
868 windowPt[xbin]->weight = sum_in_window;
869 windowPt[xbin]->xbin = xbin;
870 }
871 }
872 }
873
874 //Sort the windows according to the weight
875 SortPeaks(windowPt,0,nbinsx);
876
4fc9a6a4 877 Float_t top,butt;
878 for(Int_t i=0; i<nbinsx; i++)
879 {
880 top=butt=0;
881 Int_t xbin = windowPt[i]->xbin;
882
883 if(xbin<xmin || xbin > xmax-1) continue;
884
b1886074 885 //Check if this is really a local maxima
886 if(anotherPt[xbin-1]->weight > anotherPt[xbin]->weight ||
887 anotherPt[xbin+1]->weight > anotherPt[xbin]->weight)
4fc9a6a4 888 continue;
889
890 for(Int_t j=windowPt[i]->ymin; j<windowPt[i]->ymax; j++)
891 {
892 //Calculate the mean in y direction:
893 Int_t bin = fCurrentHisto->GetBin(windowPt[i]->xbin,j);
894 top += (fCurrentHisto->GetBinCenterY(j))*(fCurrentHisto->GetBinContent(bin));
895 butt += fCurrentHisto->GetBinContent(bin);
896 }
afd8fed4 897
b2a02bce 898 if(butt < fThreshold)
899 continue;
afd8fed4 900
901 fXPeaks[fNPeaks] = fCurrentHisto->GetBinCenterX(windowPt[i]->xbin);
902 fYPeaks[fNPeaks] = top/butt;
903 fWeight[fNPeaks] = (Int_t)butt;
4ab9f8f0 904 //cout<<"mean in y "<<ypeaks[n]<<" on x "<<windowPt[i]->xbin<<" content "<<butt<<endl;
afd8fed4 905 fNPeaks++;
3e87ef69 906 if(fNPeaks==fNMax)
907 {
908 cerr<<"AliL3HoughMaxFinder::FindPeak1 : Peak array out of range!!!"<<endl;
909 break;
910 }
4fc9a6a4 911 }
912
913
914 //Improve the peaks by including the region around in x.
915 Float_t ytop,ybutt;
b1886074 916 Int_t prev;
917 Int_t w;
afd8fed4 918 for(Int_t i=0; i<fNPeaks; i++)
4fc9a6a4 919 {
afd8fed4 920 Int_t xbin = fCurrentHisto->FindXbin(fXPeaks[i]);
4fc9a6a4 921 if(xbin - x_bin_sides < xmin || xbin + x_bin_sides > xmax) continue;
922 top=butt=0;
923 ytop=0,ybutt=0;
b1886074 924 w=0;
925 prev = xbin - x_bin_sides+1;
4fc9a6a4 926 for(Int_t j=xbin-x_bin_sides; j<=xbin+x_bin_sides; j++)
927 {
4ab9f8f0 928 /*
b1886074 929 //Check if the windows are overlapping:
930 if(anotherPt[j]->ymin > anotherPt[prev]->ymax) {prev=j; continue;}
931 if(anotherPt[j]->ymax < anotherPt[prev]->ymin) {prev=j; continue;}
932 prev = j;
4ab9f8f0 933 */
934
4fc9a6a4 935 top += fCurrentHisto->GetBinCenterX(j)*anotherPt[j]->weight;
936 butt += anotherPt[j]->weight;
b1886074 937
4fc9a6a4 938 for(Int_t k=anotherPt[j]->ymin; k<anotherPt[j]->ymax; k++)
939 {
940 Int_t bin = fCurrentHisto->GetBin(j,k);
941 ytop += (fCurrentHisto->GetBinCenterY(k))*(fCurrentHisto->GetBinContent(bin));
942 ybutt += fCurrentHisto->GetBinContent(bin);
95a00d93 943 w+=(Int_t)fCurrentHisto->GetBinContent(bin);
4fc9a6a4 944 }
945 }
b1886074 946
afd8fed4 947 fXPeaks[i] = top/butt;
948 fYPeaks[i] = ytop/ybutt;
949 fWeight[i] = w;
3e87ef69 950 //cout<<"Setting weight "<<w<<" kappa "<<fXPeaks[i]<<" phi0 "<<fYPeaks[i]<<endl;
b5a207b4 951
afd8fed4 952 /*
b1886074 953 //Check if this peak is overlapping with a previous:
3e87ef69 954 for(Int_t p=0; p<i; p++)
b1886074 955 {
3e87ef69 956 //cout<<fabs(fXPeaks[p] - fXPeaks[i])<<" "<<fabs(fYPeaks[p] - fYPeaks[i])<<endl;
957 if(fabs(fXPeaks[p] - fXPeaks[i]) < max_kappa &&
afd8fed4 958 fabs(fYPeaks[p] - fYPeaks[i]) < max_phi0)
b1886074 959 {
afd8fed4 960 fWeight[i]=0;
3e87ef69 961 //break;
b1886074 962 }
963 }
afd8fed4 964 */
4fc9a6a4 965 }
966
967 for(Int_t i=0; i<nbinsx; i++)
968 delete windowPt[i];
969 delete [] windowPt;
970 delete [] anotherPt;
4fc9a6a4 971}
972
973void AliL3HoughMaxFinder::SortPeaks(struct AxisWindow **a,Int_t first,Int_t last)
974{
975 //General sorting routine
976 //Sort according to PeakCompare()
977
978 static struct AxisWindow *tmp;
979 static int i; // "static" to save stack space
980 int j;
981
982 while (last - first > 1) {
983 i = first;
984 j = last;
985 for (;;) {
986 while (++i < last && PeakCompare(a[i], a[first]) < 0)
987 ;
988 while (--j > first && PeakCompare(a[j], a[first]) > 0)
989 ;
990 if (i >= j)
991 break;
992
993 tmp = a[i];
994 a[i] = a[j];
995 a[j] = tmp;
996 }
997 if (j == first) {
998 ++first;
999 continue;
1000 }
1001 tmp = a[first];
1002 a[first] = a[j];
1003 a[j] = tmp;
1004 if (j - first < last - (j + 1)) {
1005 SortPeaks(a, first, j);
1006 first = j + 1; // QSort(j + 1, last);
1007 } else {
1008 SortPeaks(a, j + 1, last);
1009 last = j; // QSort(first, j);
1010 }
1011 }
1012
1013}
1014
1015Int_t AliL3HoughMaxFinder::PeakCompare(struct AxisWindow *a,struct AxisWindow *b)
1016{
1017 if(a->weight < b->weight) return 1;
1018 if(a->weight > b->weight) return -1;
1019 return 0;
4fc9a6a4 1020}
52a2a604 1021
7f66cda9 1022void AliL3HoughMaxFinder::FindPeak(Int_t t1,Double_t t2,Int_t t3)
e4c21048 1023{
1024 //Attempt of a more sophisticated peak finder.
4fc9a6a4 1025 //Finds the best peak in the histogram, and returns the corresponding
1026 //track object.
1027
4cafa5fc 1028 if(!fCurrentHisto)
1029 {
1030 printf("AliL3HoughMaxFinder::FindPeak : No histogram!!\n");
b1886074 1031 return;
4cafa5fc 1032 }
1033 AliL3Histogram *hist = fCurrentHisto;
b2a02bce 1034 if(hist->GetNEntries()==0)
1035 return;
81516384 1036
4cafa5fc 1037 Int_t xmin = hist->GetFirstXbin();
1038 Int_t xmax = hist->GetLastXbin();
1039 Int_t ymin = hist->GetFirstYbin();
1040 Int_t ymax = hist->GetLastYbin();
1041 Int_t nbinsx = hist->GetNbinsX()+1;
e4c21048 1042
1043 Int_t *m = new Int_t[nbinsx];
1044 Int_t *m_low = new Int_t[nbinsx];
1045 Int_t *m_up = new Int_t[nbinsx];
81516384 1046
1047
1048 recompute: //this is a goto.
1049
e4c21048 1050 for(Int_t i=0; i<nbinsx; i++)
1051 {
1052 m[i]=0;
1053 m_low[i]=0;
1054 m_up[i]=0;
1055 }
e4c21048 1056
e06900d5 1057 Int_t max_x=0,sum=0,max_xbin=0,bin=0;
e4c21048 1058
1059 for(Int_t xbin=xmin; xbin<=xmax; xbin++)
1060 {
4cafa5fc 1061 for(Int_t ybin=ymin; ybin <= ymax - t1; ybin++)
e4c21048 1062 {
1063 sum = 0;
4cafa5fc 1064 for(Int_t y=ybin; y <= ybin+t1; y++)
e4c21048 1065 {
4cafa5fc 1066 if(y>ymax) break;
e4c21048 1067 //Inside window
1068 bin = hist->GetBin(xbin,y);
1069 sum += (Int_t)hist->GetBinContent(bin);
1070
1071 }
1072 if(sum > m[xbin]) //Max value locally in this xbin
1073 {
1074 m[xbin]=sum;
1075 m_low[xbin]=ybin;
4cafa5fc 1076 m_up[xbin]=ybin + t1;
e4c21048 1077 }
1078
1079 }
1080
1081 if(m[xbin] > max_x) //Max value globally in x-direction
1082 {
1083 max_xbin = xbin;
1084 max_x = m[xbin];//sum;
1085 }
1086 }
1087 //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]);
4cafa5fc 1088 //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[max_xbin]),hist->GetBinCenterY(m_up[max_xbin]));
e4c21048 1089
1090 //Determine a width in the x-direction
1091 Int_t x_low=0,x_up=0;
4cafa5fc 1092
e4c21048 1093 for(Int_t xbin=max_xbin-1; xbin >= xmin; xbin--)
1094 {
4cafa5fc 1095 if(m[xbin] < max_x*t2)
e4c21048 1096 {
81516384 1097 x_low = xbin+1;
e4c21048 1098 break;
1099 }
1100 }
1101 for(Int_t xbin = max_xbin+1; xbin <=xmax; xbin++)
1102 {
4cafa5fc 1103 if(m[xbin] < max_x*t2)
e4c21048 1104 {
81516384 1105 x_up = xbin-1;
e4c21048 1106 break;
1107 }
1108 }
1109
81516384 1110 Double_t top=0,butt=0,value,x_peak;
1111 if(x_up - x_low + 1 > t3)
1112 {
1113 t1 -= 1;
4cafa5fc 1114 printf("\nxrange out if limit x_up %d x_low %d t1 %d\n\n",x_low,x_up,t1);
81516384 1115 if(t1 > 1)
1116 goto recompute;
1117 else
1118 {
4cafa5fc 1119 x_peak = hist->GetBinCenterX(max_xbin);
81516384 1120 goto moveon;
1121 }
1122 }
1123
4cafa5fc 1124 //printf("xlow %f xup %f\n",hist->GetBinCenterX(x_low),hist->GetBinCenterX(x_up));
81516384 1125 //printf("Spread in x %d\n",x_up-x_low +1);
e4c21048 1126
1127 //Now, calculate the center of mass in x-direction
e4c21048 1128 for(Int_t xbin=x_low; xbin <= x_up; xbin++)
1129 {
4cafa5fc 1130 value = hist->GetBinCenterX(xbin);
e4c21048 1131 top += value*m[xbin];
1132 butt += m[xbin];
1133 }
1134 x_peak = top/butt;
1135
81516384 1136 moveon:
e4c21048 1137
1138 //Find the peak in y direction:
4cafa5fc 1139 Int_t x_l = hist->FindXbin(x_peak);
1140 if(hist->GetBinCenterX(x_l) > x_peak)
81516384 1141 x_l--;
1142
e4c21048 1143 Int_t x_u = x_l + 1;
81516384 1144
4cafa5fc 1145 if(hist->GetBinCenterX(x_l) > x_peak || hist->GetBinCenterX(x_u) <= x_peak)
1146 printf("\nAliL3HoughMaxFinder::FindPeak : Wrong xrange %f %f %f\n\n",hist->GetBinCenterX(x_l),x_peak,hist->GetBinCenterX(x_u));
e4c21048 1147
4cafa5fc 1148 //printf("\nxlow %f xup %f\n",hist->GetBinCenterX(x_l),hist->GetBinCenterX(x_u));
e4c21048 1149
1150 value=top=butt=0;
1151
4cafa5fc 1152 //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_l]),hist->GetBinCenterY(m_up[x_l]));
1153 //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_u]),hist->GetBinCenterY(m_up[x_u]));
e4c21048 1154
1155 for(Int_t ybin=m_low[x_l]; ybin <= m_up[x_l]; ybin++)
1156 {
4cafa5fc 1157 value = hist->GetBinCenterY(ybin);
e4c21048 1158 bin = hist->GetBin(x_l,ybin);
1159 top += value*hist->GetBinContent(bin);
1160 butt += hist->GetBinContent(bin);
1161 }
1162 Double_t y_peak_low = top/butt;
1163
81516384 1164 //printf("y_peak_low %f\n",y_peak_low);
e4c21048 1165
1166 value=top=butt=0;
1167 for(Int_t ybin=m_low[x_u]; ybin <= m_up[x_u]; ybin++)
1168 {
4cafa5fc 1169 value = hist->GetBinCenterY(ybin);
e4c21048 1170 bin = hist->GetBin(x_u,ybin);
1171 top += value*hist->GetBinContent(bin);
1172 butt += hist->GetBinContent(bin);
1173 }
1174 Double_t y_peak_up = top/butt;
1175
81516384 1176 //printf("y_peak_up %f\n",y_peak_up);
e4c21048 1177
4cafa5fc 1178 Double_t x_value_up = hist->GetBinCenterX(x_u);
1179 Double_t x_value_low = hist->GetBinCenterX(x_l);
e4c21048 1180
1181 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);
1182
1183
81516384 1184 //Find the weight:
95a00d93 1185 //bin = hist->FindBin(x_peak,y_peak);
1186 //Int_t weight = (Int_t)hist->GetBinContent(bin);
81516384 1187
b1886074 1188 //AliL3HoughTrack *track = new AliL3HoughTrack();
1189 //track->SetTrackParameters(x_peak,y_peak,weight);
7f66cda9 1190 fXPeaks[fNPeaks]=x_peak;
1191 fYPeaks[fNPeaks]=y_peak;
1192 fWeight[fNPeaks]=(Int_t)hist->GetBinContent(bin);
1193 fNPeaks++;
52a2a604 1194
1195 delete [] m;
1196 delete [] m_low;
1197 delete [] m_up;
1198
b1886074 1199 //return track;
e4c21048 1200}
81516384 1201