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