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