]>
Commit | Line | Data |
---|---|---|
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 | ||
13 | ClassImp(AliL3HoughMaxFinder) | |
14 | ||
15 | ||
16 | AliL3HoughMaxFinder::AliL3HoughMaxFinder() | |
17 | { | |
18 | //Default constructor | |
19 | fThreshold = 0; | |
20 | //fTracks = 0; | |
21 | fHistoType=0; | |
22 | } | |
23 | ||
24 | ||
4cafa5fc | 25 | AliL3HoughMaxFinder::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 | ||
38 | AliL3HoughMaxFinder::~AliL3HoughMaxFinder() | |
39 | { | |
40 | //Destructor | |
41 | ||
42 | } | |
43 | ||
4fc9a6a4 | 44 | void 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 | ||
77 | AliL3TrackArray *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 | 129 | AliL3TrackArray *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 | 206 | AliL3TrackArray *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 | 267 | AliL3HoughTrack *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 | 456 | AliL3HoughTrack *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 | ||
499 | void 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 | ||
619 | void 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 | ||
661 | Int_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 | 669 | AliL3HoughTrack *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 |