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