]>
Commit | Line | Data |
---|---|---|
1 | // @(#) $Id$ | |
2 | ||
3 | // Author: Anders Vestbo <mailto:vestbo@fi.uib.no> | |
4 | //*-- Copyright © ALICE HLT Group | |
5 | ||
6 | #include <strings.h> | |
7 | #include "AliL3StandardIncludes.h" | |
8 | ||
9 | #ifndef no_root | |
10 | #include <TNtuple.h> | |
11 | #include <TFile.h> | |
12 | #endif | |
13 | ||
14 | #include "AliL3Logging.h" | |
15 | #include "AliL3HoughMaxFinder.h" | |
16 | #include "AliL3Histogram.h" | |
17 | #include "AliL3TrackArray.h" | |
18 | #include "AliL3HoughTrack.h" | |
19 | ||
20 | #if __GNUC__ == 3 | |
21 | using namespace std; | |
22 | #endif | |
23 | ||
24 | /** \class AliL3HoughMaxFinder | |
25 | <pre> | |
26 | //_____________________________________________________________ | |
27 | // AliL3HoughMaxFinder | |
28 | // | |
29 | // Maximum finder | |
30 | // | |
31 | </pre> | |
32 | */ | |
33 | ||
34 | ClassImp(AliL3HoughMaxFinder) | |
35 | ||
36 | ||
37 | AliL3HoughMaxFinder::AliL3HoughMaxFinder() | |
38 | { | |
39 | //Default constructor | |
40 | fThreshold = 0; | |
41 | fHistoType=0; | |
42 | fXPeaks=0; | |
43 | fYPeaks=0; | |
44 | fNPeaks=0; | |
45 | fNMax=0; | |
46 | fGradX=1; | |
47 | fGradY=1; | |
48 | #ifndef no_root | |
49 | fNtuppel = 0; | |
50 | #endif | |
51 | } | |
52 | ||
53 | AliL3HoughMaxFinder::AliL3HoughMaxFinder(Char_t *histotype,Int_t nmax,AliL3Histogram *hist) | |
54 | { | |
55 | //Constructor | |
56 | ||
57 | //fTracks = new AliL3TrackArray("AliL3HoughTrack"); | |
58 | if(strcmp(histotype,"KappaPhi")==0) fHistoType='c'; | |
59 | if(strcmp(histotype,"DPsi")==0) fHistoType='l'; | |
60 | ||
61 | if(hist) | |
62 | fCurrentHisto = hist; | |
63 | ||
64 | fGradX=1; | |
65 | fGradY=1; | |
66 | fNMax=nmax; | |
67 | fXPeaks = new Float_t[fNMax]; | |
68 | fYPeaks = new Float_t[fNMax]; | |
69 | fWeight = new Int_t[fNMax]; | |
70 | #ifndef no_root | |
71 | fNtuppel = 0; | |
72 | #endif | |
73 | fThreshold=0; | |
74 | } | |
75 | ||
76 | AliL3HoughMaxFinder::~AliL3HoughMaxFinder() | |
77 | { | |
78 | //Destructor | |
79 | if(fXPeaks) | |
80 | delete [] fXPeaks; | |
81 | if(fYPeaks) | |
82 | delete [] fYPeaks; | |
83 | if(fWeight) | |
84 | delete [] fWeight; | |
85 | #ifndef no_root | |
86 | if(fNtuppel) | |
87 | delete fNtuppel; | |
88 | #endif | |
89 | } | |
90 | ||
91 | void AliL3HoughMaxFinder::Reset() | |
92 | { | |
93 | for(Int_t i=0; i<fNMax; i++) | |
94 | { | |
95 | fXPeaks[i]=0; | |
96 | fYPeaks[i]=0; | |
97 | fWeight[i]=0; | |
98 | } | |
99 | fNPeaks=0; | |
100 | } | |
101 | ||
102 | void AliL3HoughMaxFinder::CreateNtuppel() | |
103 | { | |
104 | #ifndef no_root | |
105 | //content#; neighbouring bins of the peak. | |
106 | fNtuppel = new TNtuple("ntuppel","Peak charateristics","kappa:phi0:weigth:content3:content5:content1:content7"); | |
107 | fNtuppel->SetDirectory(0); | |
108 | #endif | |
109 | } | |
110 | ||
111 | void AliL3HoughMaxFinder::WriteNtuppel(Char_t *filename) | |
112 | { | |
113 | #ifndef no_root | |
114 | TFile *file = TFile::Open(filename,"RECREATE"); | |
115 | if(!file) | |
116 | { | |
117 | cerr<<"AliL3HoughMaxFinder::WriteNtuppel : Error opening file "<<filename<<endl; | |
118 | return; | |
119 | } | |
120 | fNtuppel->Write(); | |
121 | file->Close(); | |
122 | #endif | |
123 | } | |
124 | ||
125 | void AliL3HoughMaxFinder::FindAbsMaxima() | |
126 | { | |
127 | ||
128 | if(!fCurrentHisto) | |
129 | { | |
130 | cerr<<"AliL3HoughMaxFinder::FindAbsMaxima : No histogram"<<endl; | |
131 | return; | |
132 | } | |
133 | AliL3Histogram *hist = fCurrentHisto; | |
134 | ||
135 | if(hist->GetNEntries() == 0) | |
136 | return; | |
137 | ||
138 | Int_t xmin = hist->GetFirstXbin(); | |
139 | Int_t xmax = hist->GetLastXbin(); | |
140 | Int_t ymin = hist->GetFirstYbin(); | |
141 | Int_t ymax = hist->GetLastYbin(); | |
142 | Int_t bin; | |
143 | Double_t value,max_value=0; | |
144 | ||
145 | Int_t max_xbin=0,max_ybin=0; | |
146 | for(Int_t xbin=xmin; xbin<=xmax; xbin++) | |
147 | { | |
148 | for(Int_t ybin=ymin; ybin<=ymax; ybin++) | |
149 | { | |
150 | bin = hist->GetBin(xbin,ybin); | |
151 | value = hist->GetBinContent(bin); | |
152 | if(value>max_value) | |
153 | { | |
154 | max_value = value; | |
155 | max_xbin = xbin; | |
156 | max_ybin = ybin; | |
157 | } | |
158 | } | |
159 | } | |
160 | ||
161 | if(max_value == 0) | |
162 | return; | |
163 | ||
164 | if(fNPeaks > fNMax) | |
165 | { | |
166 | cerr<<"AliL3HoughMaxFinder::FindAbsMaxima : Array out of range : "<<fNPeaks<<endl; | |
167 | return; | |
168 | } | |
169 | ||
170 | Double_t max_x = hist->GetBinCenterX(max_xbin); | |
171 | Double_t max_y = hist->GetBinCenterY(max_ybin); | |
172 | fXPeaks[fNPeaks] = max_x; | |
173 | fYPeaks[fNPeaks] = max_y; | |
174 | fWeight[fNPeaks] = (Int_t)max_value; | |
175 | ||
176 | fNPeaks++; | |
177 | #ifndef no_root | |
178 | if(fNtuppel) | |
179 | { | |
180 | Int_t bin3 = hist->GetBin(max_xbin-1,max_ybin); | |
181 | Int_t bin5 = hist->GetBin(max_xbin+1,max_ybin); | |
182 | Int_t bin1 = hist->GetBin(max_xbin,max_ybin-1); | |
183 | Int_t bin7 = hist->GetBin(max_xbin,max_ybin+1); | |
184 | ||
185 | fNtuppel->Fill(max_x,max_y,max_value,hist->GetBinContent(bin3),hist->GetBinContent(bin5),hist->GetBinContent(bin1),hist->GetBinContent(bin7)); | |
186 | } | |
187 | #endif | |
188 | } | |
189 | ||
190 | void AliL3HoughMaxFinder::FindBigMaxima() | |
191 | { | |
192 | ||
193 | AliL3Histogram *hist = fCurrentHisto; | |
194 | ||
195 | if(hist->GetNEntries() == 0) | |
196 | return; | |
197 | ||
198 | Int_t xmin = hist->GetFirstXbin(); | |
199 | Int_t xmax = hist->GetLastXbin(); | |
200 | Int_t ymin = hist->GetFirstYbin(); | |
201 | Int_t ymax = hist->GetLastYbin(); | |
202 | Int_t bin[25],bin_index; | |
203 | Double_t value[25]; | |
204 | ||
205 | for(Int_t xbin=xmin+2; xbin<xmax-3; xbin++) | |
206 | { | |
207 | for(Int_t ybin=ymin+2; ybin<ymax-3; ybin++) | |
208 | { | |
209 | bin_index=0; | |
210 | for(Int_t xb=xbin-2; xb<=xbin+2; xb++) | |
211 | { | |
212 | for(Int_t yb=ybin-2; yb<=ybin+2; yb++) | |
213 | { | |
214 | bin[bin_index]=hist->GetBin(xb,yb); | |
215 | value[bin_index]=hist->GetBinContent(bin[bin_index]); | |
216 | bin_index++; | |
217 | } | |
218 | } | |
219 | if(value[12]==0) continue; | |
220 | Int_t b=0; | |
221 | while(1) | |
222 | { | |
223 | if(value[b] > value[12] || b==bin_index) break; | |
224 | b++; | |
225 | //printf("b %d\n",b); | |
226 | } | |
227 | if(b == bin_index) | |
228 | { | |
229 | //Found maxima | |
230 | if(fNPeaks > fNMax) | |
231 | { | |
232 | cerr<<"AliL3HoughMaxFinder::FindBigMaxima : Array out of range "<<fNPeaks<<endl; | |
233 | return; | |
234 | } | |
235 | ||
236 | Double_t max_x = hist->GetBinCenterX(xbin); | |
237 | Double_t max_y = hist->GetBinCenterY(ybin); | |
238 | fXPeaks[fNPeaks] = max_x; | |
239 | fYPeaks[fNPeaks] = max_y; | |
240 | fNPeaks++; | |
241 | } | |
242 | } | |
243 | } | |
244 | } | |
245 | ||
246 | void AliL3HoughMaxFinder::FindMaxima(Int_t threshold) | |
247 | { | |
248 | //Locate all the maxima in input histogram. | |
249 | //Maxima is defined as bins with more entries than the | |
250 | //immediately neighbouring bins. | |
251 | ||
252 | if(fCurrentHisto->GetNEntries() == 0) | |
253 | return; | |
254 | ||
255 | Int_t xmin = fCurrentHisto->GetFirstXbin(); | |
256 | Int_t xmax = fCurrentHisto->GetLastXbin(); | |
257 | Int_t ymin = fCurrentHisto->GetFirstYbin(); | |
258 | Int_t ymax = fCurrentHisto->GetLastYbin(); | |
259 | Int_t bin[9]; | |
260 | Double_t value[9]; | |
261 | ||
262 | //Float_t max_kappa = 0.001; | |
263 | //Float_t max_phi0 = 0.08; | |
264 | ||
265 | for(Int_t xbin=xmin+1; xbin<=xmax-1; xbin++) | |
266 | { | |
267 | for(Int_t ybin=ymin+1; ybin<=ymax-1; ybin++) | |
268 | { | |
269 | bin[0] = fCurrentHisto->GetBin(xbin-1,ybin-1); | |
270 | bin[1] = fCurrentHisto->GetBin(xbin,ybin-1); | |
271 | bin[2] = fCurrentHisto->GetBin(xbin+1,ybin-1); | |
272 | bin[3] = fCurrentHisto->GetBin(xbin-1,ybin); | |
273 | bin[4] = fCurrentHisto->GetBin(xbin,ybin); | |
274 | bin[5] = fCurrentHisto->GetBin(xbin+1,ybin); | |
275 | bin[6] = fCurrentHisto->GetBin(xbin-1,ybin+1); | |
276 | bin[7] = fCurrentHisto->GetBin(xbin,ybin+1); | |
277 | bin[8] = fCurrentHisto->GetBin(xbin+1,ybin+1); | |
278 | value[0] = fCurrentHisto->GetBinContent(bin[0]); | |
279 | value[1] = fCurrentHisto->GetBinContent(bin[1]); | |
280 | value[2] = fCurrentHisto->GetBinContent(bin[2]); | |
281 | value[3] = fCurrentHisto->GetBinContent(bin[3]); | |
282 | value[4] = fCurrentHisto->GetBinContent(bin[4]); | |
283 | value[5] = fCurrentHisto->GetBinContent(bin[5]); | |
284 | value[6] = fCurrentHisto->GetBinContent(bin[6]); | |
285 | value[7] = fCurrentHisto->GetBinContent(bin[7]); | |
286 | value[8] = fCurrentHisto->GetBinContent(bin[8]); | |
287 | ||
288 | ||
289 | if(value[4]>value[0] && value[4]>value[1] && value[4]>value[2] | |
290 | && value[4]>value[3] && value[4]>value[5] && value[4]>value[6] | |
291 | && value[4]>value[7] && value[4]>value[8]) | |
292 | { | |
293 | //Found a local maxima | |
294 | Float_t max_x = fCurrentHisto->GetBinCenterX(xbin); | |
295 | Float_t max_y = fCurrentHisto->GetBinCenterY(ybin); | |
296 | ||
297 | if((Int_t)value[4] <= threshold) continue;//central bin below threshold | |
298 | if(fNPeaks >= fNMax) | |
299 | { | |
300 | cout<<"AliL3HoughMaxFinder::FindMaxima : Array out of range "<<fNPeaks<<endl; | |
301 | return; | |
302 | } | |
303 | ||
304 | //Check the gradient: | |
305 | if(value[3]/value[4] > fGradX && value[5]/value[4] > fGradX) | |
306 | continue; | |
307 | ||
308 | if(value[1]/value[4] > fGradY && value[7]/value[4] > fGradY) | |
309 | continue; | |
310 | ||
311 | fXPeaks[fNPeaks] = max_x; | |
312 | fYPeaks[fNPeaks] = max_y; | |
313 | fWeight[fNPeaks] = (Int_t)value[4]; | |
314 | fNPeaks++; | |
315 | ||
316 | /* | |
317 | //Check if the peak is overlapping with a previous: | |
318 | Bool_t bigger = kFALSE; | |
319 | for(Int_t p=0; p<fNPeaks; p++) | |
320 | { | |
321 | if(fabs(max_x - fXPeaks[p]) < max_kappa && fabs(max_y - fYPeaks[p]) < max_phi0) | |
322 | { | |
323 | bigger = kTRUE; | |
324 | if(value[4] > fWeight[p]) //this peak is bigger. | |
325 | { | |
326 | fXPeaks[p] = max_x; | |
327 | fYPeaks[p] = max_y; | |
328 | fWeight[p] = (Int_t)value[4]; | |
329 | } | |
330 | else | |
331 | continue; //previous peak is bigger. | |
332 | } | |
333 | } | |
334 | if(!bigger) //there were no overlapping peaks. | |
335 | { | |
336 | fXPeaks[fNPeaks] = max_x; | |
337 | fYPeaks[fNPeaks] = max_y; | |
338 | fWeight[fNPeaks] = (Int_t)value[4]; | |
339 | fNPeaks++; | |
340 | } | |
341 | */ | |
342 | } | |
343 | } | |
344 | } | |
345 | ||
346 | } | |
347 | ||
348 | struct Window | |
349 | { | |
350 | Int_t start; | |
351 | Int_t sum; | |
352 | }; | |
353 | ||
354 | void AliL3HoughMaxFinder::FindAdaptedPeaks(Int_t kappawindow,Float_t cut_ratio) | |
355 | { | |
356 | //Peak finder which looks for peaks with a certain shape. | |
357 | //The first step involves a pre-peak finder, which looks for peaks | |
358 | //in windows (size controlled by kappawindow) summing over each psi-bin. | |
359 | //These pre-preaks are then matched between neighbouring kappa-bins to | |
360 | //look for real 2D peaks exhbiting the typical cross-shape in the Hough circle transform. | |
361 | //The maximum bin within this region is marked as the peak itself, and | |
362 | //a few checks is performed to avoid the clear fake peaks (asymmetry check etc.) | |
363 | ||
364 | ||
365 | AliL3Histogram *hist = fCurrentHisto; | |
366 | ||
367 | if(!hist) | |
368 | { | |
369 | cerr<<"AliL3HoughMaxFinder : No histogram!"<<endl; | |
370 | return; | |
371 | } | |
372 | ||
373 | if(hist->GetNEntries() == 0) | |
374 | return; | |
375 | ||
376 | Int_t xmin = hist->GetFirstXbin(); | |
377 | Int_t xmax = hist->GetLastXbin(); | |
378 | Int_t ymin = hist->GetFirstYbin(); | |
379 | Int_t ymax = hist->GetLastYbin(); | |
380 | ||
381 | ||
382 | //Start by looking for pre-peaks: | |
383 | ||
384 | Window **local_maxima = new Window*[hist->GetNbinsY()]; | |
385 | ||
386 | Short_t *nmaxs = new Short_t[hist->GetNbinsY()]; | |
387 | Int_t n,last_sum,sum; | |
388 | Bool_t sum_was_rising; | |
389 | for(Int_t ybin=ymin; ybin<=ymax; ybin++) | |
390 | { | |
391 | local_maxima[ybin-ymin] = new Window[hist->GetNbinsX()]; | |
392 | nmaxs[ybin-ymin] = 0; | |
393 | sum_was_rising=0; | |
394 | last_sum=0; | |
395 | n=0; | |
396 | for(Int_t xbin=xmin; xbin<=xmax-kappawindow; xbin++) | |
397 | { | |
398 | sum=0; | |
399 | for(Int_t lbin=xbin; lbin<xbin+kappawindow; lbin++) | |
400 | sum += hist->GetBinContent(hist->GetBin(lbin,ybin)); | |
401 | ||
402 | if(sum < last_sum) | |
403 | { | |
404 | if(sum > fThreshold) | |
405 | if(sum_was_rising)//Previous sum was a local maxima | |
406 | { | |
407 | local_maxima[ybin-ymin][nmaxs[ybin-ymin]].start = xbin-1; | |
408 | local_maxima[ybin-ymin][nmaxs[ybin-ymin]].sum = last_sum; | |
409 | nmaxs[ybin-ymin]++; | |
410 | } | |
411 | ||
412 | sum_was_rising=0; | |
413 | } | |
414 | else if(sum > 0) | |
415 | sum_was_rising=1; | |
416 | last_sum=sum; | |
417 | } | |
418 | } | |
419 | ||
420 | Int_t match=0; | |
421 | Int_t *starts = new Int_t[hist->GetNbinsY()+1]; | |
422 | Int_t *maxs = new Int_t[hist->GetNbinsY()+1]; | |
423 | ||
424 | for(Int_t ybin=ymax; ybin >= ymin+1; ybin--) | |
425 | { | |
426 | for(Int_t i=0; i<nmaxs[ybin-ymin]; i++) | |
427 | { | |
428 | Int_t lw = local_maxima[ybin-ymin][i].sum; | |
429 | ||
430 | if(lw<0) | |
431 | continue; //already used | |
432 | ||
433 | Int_t maxvalue=0,maxybin=0,maxxbin=0,maxwindow=0; | |
434 | for(Int_t k=local_maxima[ybin-ymin][i].start; k<local_maxima[ybin-ymin][i].start + kappawindow; k++) | |
435 | if(hist->GetBinContent(hist->GetBin(k,ybin)) > maxvalue) | |
436 | { | |
437 | maxvalue = hist->GetBinContent(hist->GetBin(k,ybin)); | |
438 | maxybin = ybin; | |
439 | maxxbin = k; | |
440 | } | |
441 | ||
442 | //start expanding in the psi-direction: | |
443 | ||
444 | Int_t lb = local_maxima[ybin-ymin][i].start; | |
445 | //Int_t ystart=ybin; | |
446 | starts[ybin] = local_maxima[ybin-ymin][i].start; | |
447 | maxs[ybin] = maxxbin; | |
448 | Int_t yl=ybin-1,nybins=1; | |
449 | ||
450 | //cout<<"Starting search at ybin "<<ybin<<" start "<<lb<<" with sum "<<local_maxima[ybin-ymin][i].sum<<endl; | |
451 | while(yl >= ymin) | |
452 | { | |
453 | Bool_t found=0; | |
454 | for(Int_t j=0; j<nmaxs[yl-ymin]; j++) | |
455 | { | |
456 | if( local_maxima[yl-ymin][j].start - lb < 0) continue; | |
457 | if( local_maxima[yl-ymin][j].start < lb + kappawindow + match && | |
458 | local_maxima[yl-ymin][j].start >= lb && local_maxima[yl-ymin][j].sum > 0) | |
459 | { | |
460 | ||
461 | //cout<<"match at ybin "<<yl<<" yvalue "<<hist->GetBinCenterY(yl)<<" start "<<local_maxima[yl-ymin][j].start<<" sum "<<local_maxima[yl-ymin][j].sum<<endl; | |
462 | ||
463 | Int_t lmaxvalue=0,lmaxxbin=0; | |
464 | for(Int_t k=local_maxima[yl-ymin][j].start; k<local_maxima[yl-ymin][j].start + kappawindow; k++) | |
465 | { | |
466 | if(hist->GetBinContent(hist->GetBin(k,yl)) > maxvalue) | |
467 | { | |
468 | maxvalue = hist->GetBinContent(hist->GetBin(k,yl)); | |
469 | maxxbin = k; | |
470 | maxybin = yl; | |
471 | maxwindow = j; | |
472 | } | |
473 | if(hist->GetBinContent(hist->GetBin(k,yl)) > lmaxvalue)//local maxima value | |
474 | { | |
475 | lmaxvalue=hist->GetBinContent(hist->GetBin(k,yl)); | |
476 | lmaxxbin=k; | |
477 | } | |
478 | } | |
479 | nybins++; | |
480 | starts[yl] = local_maxima[yl-ymin][j].start; | |
481 | maxs[yl] = lmaxxbin; | |
482 | local_maxima[yl-ymin][j].sum=-1; //Mark as used | |
483 | found=1; | |
484 | lb = local_maxima[yl-ymin][j].start; | |
485 | break;//Since we found a match in this bin, we dont have to search it anymore, goto next bin. | |
486 | } | |
487 | } | |
488 | if(!found || yl == ymin)//no more local maximas to be matched, so write the final peak and break the expansion: | |
489 | { | |
490 | if(nybins > 4) | |
491 | { | |
492 | //cout<<"ystart "<<ystart<<" and nybins "<<nybins<<endl; | |
493 | ||
494 | Bool_t truepeak=kTRUE; | |
495 | ||
496 | //cout<<"Maxima found at xbin "<<maxxbin<<" ybin "<<maxybin<<" value "<<maxvalue<<endl; | |
497 | //cout<<"Starting to sum at xbin "<<starts[maxybin-ymin]<<endl; | |
498 | ||
499 | ||
500 | //Look in a window on both sides to probe the asymmetry | |
501 | Float_t right=0,left=0; | |
502 | for(Int_t w=maxxbin+1; w<=maxxbin+3; w++) | |
503 | { | |
504 | for(Int_t r=maxybin+1; r<=maxybin+3; r++) | |
505 | { | |
506 | right += (Float_t)hist->GetBinContent(hist->GetBin(w,r)); | |
507 | } | |
508 | } | |
509 | ||
510 | for(Int_t w=maxxbin-1; w>=maxxbin-3; w--) | |
511 | { | |
512 | for(Int_t r=maxybin+1; r<=maxybin+3; r++) | |
513 | { | |
514 | left += (Float_t)hist->GetBinContent(hist->GetBin(w,r)); | |
515 | } | |
516 | } | |
517 | ||
518 | //cout<<"ratio "<<right/left<<endl; | |
519 | ||
520 | Float_t upper_ratio=1,lower_ratio=1; | |
521 | if(left) | |
522 | upper_ratio = right/left; | |
523 | ||
524 | right=left=0; | |
525 | for(Int_t w=maxxbin+1; w<=maxxbin+3; w++) | |
526 | { | |
527 | for(Int_t r=maxybin-1; r>=maxybin-3; r--) | |
528 | { | |
529 | right += (Float_t)hist->GetBinContent(hist->GetBin(w,r)); | |
530 | } | |
531 | } | |
532 | ||
533 | for(Int_t w=maxxbin-1; w>=maxxbin-3; w--) | |
534 | { | |
535 | for(Int_t r=maxybin-1; r>=maxybin-3; r--) | |
536 | { | |
537 | left += (Float_t)hist->GetBinContent(hist->GetBin(w,r)); | |
538 | } | |
539 | } | |
540 | ||
541 | //cout<<"ratio "<<left/right<<endl; | |
542 | ||
543 | if(right) | |
544 | lower_ratio = left/right; | |
545 | ||
546 | if(upper_ratio > cut_ratio || lower_ratio > cut_ratio) | |
547 | truepeak=kFALSE; | |
548 | ||
549 | if(truepeak) | |
550 | { | |
551 | ||
552 | fXPeaks[fNPeaks] = hist->GetBinCenterX(maxxbin); | |
553 | fYPeaks[fNPeaks] = hist->GetBinCenterY(maxybin); | |
554 | fWeight[fNPeaks] = maxvalue; | |
555 | fNPeaks++; | |
556 | ||
557 | /* | |
558 | //Calculate the peak using weigthed means: | |
559 | Float_t sum=0; | |
560 | fYPeaks[fNPeaks]=0; | |
561 | for(Int_t k=maxybin-1; k<=maxybin+1; k++) | |
562 | { | |
563 | Float_t lsum = 0; | |
564 | for(Int_t l=starts[k]; l<starts[k]+kappawindow; l++) | |
565 | { | |
566 | lsum += (Float_t)hist->GetBinContent(hist->GetBin(l,k)); | |
567 | sum += (Float_t)hist->GetBinContent(hist->GetBin(l,k)); | |
568 | } | |
569 | fYPeaks[fNPeaks] += lsum*hist->GetBinCenterY(k); | |
570 | } | |
571 | fYPeaks[fNPeaks] /= sum; | |
572 | Int_t ybin1,ybin2; | |
573 | if(fYPeaks[fNPeaks] < hist->GetBinCenterY(hist->FindYbin(fYPeaks[fNPeaks]))) | |
574 | { | |
575 | ybin1 = hist->FindYbin(fYPeaks[fNPeaks])-1; | |
576 | ybin2 = ybin1+1; | |
577 | } | |
578 | else | |
579 | { | |
580 | ybin1 = hist->FindYbin(fYPeaks[fNPeaks]); | |
581 | ybin2 = ybin1+1; | |
582 | } | |
583 | ||
584 | Float_t kappa1=0,kappa2=0; | |
585 | sum=0; | |
586 | for(Int_t k=starts[ybin1]; k<starts[ybin1] + kappawindow; k++) | |
587 | { | |
588 | kappa1 += hist->GetBinCenterX(k)*hist->GetBinContent(hist->GetBin(k,ybin1)); | |
589 | sum += (Float_t)hist->GetBinContent(hist->GetBin(k,ybin1)); | |
590 | } | |
591 | kappa1 /= sum; | |
592 | sum=0; | |
593 | for(Int_t k=starts[ybin2]; k<starts[ybin2] + kappawindow; k++) | |
594 | { | |
595 | kappa2 += hist->GetBinCenterX(k)*hist->GetBinContent(hist->GetBin(k,ybin2)); | |
596 | sum += (Float_t)hist->GetBinContent(hist->GetBin(k,ybin2)); | |
597 | } | |
598 | kappa2 /= sum; | |
599 | ||
600 | fXPeaks[fNPeaks] = ( kappa1*( hist->GetBinCenterY(ybin2) - fYPeaks[fNPeaks] ) + | |
601 | kappa2*( fYPeaks[fNPeaks] - hist->GetBinCenterY(ybin1) ) ) / | |
602 | (hist->GetBinCenterY(ybin2) - hist->GetBinCenterY(ybin1)); | |
603 | ||
604 | fNPeaks++; | |
605 | */ | |
606 | } | |
607 | } | |
608 | break; | |
609 | } | |
610 | else | |
611 | yl--;//Search continues... | |
612 | } | |
613 | } | |
614 | } | |
615 | ||
616 | for(Int_t i=0; i<hist->GetNbinsY(); i++) | |
617 | delete local_maxima[i]; | |
618 | ||
619 | delete [] local_maxima; | |
620 | delete [] nmaxs; | |
621 | delete [] starts; | |
622 | delete [] maxs; | |
623 | } | |
624 | ||
625 | struct PreYPeak | |
626 | { | |
627 | Int_t start_position; | |
628 | Int_t end_position; | |
629 | Int_t value; | |
630 | Int_t prev_value; | |
631 | }; | |
632 | ||
633 | struct Pre2DPeak | |
634 | { | |
635 | Int_t start_x_position; | |
636 | Int_t end_x_position; | |
637 | Int_t start_y_position; | |
638 | Int_t end_y_position; | |
639 | Int_t value; | |
640 | }; | |
641 | ||
642 | void AliL3HoughMaxFinder::FindAdaptedRowPeaks(Int_t kappawindow,Int_t /*xsize*/,Int_t ysize) | |
643 | { | |
644 | ||
645 | AliL3Histogram *hist = fCurrentHisto; | |
646 | ||
647 | if(!hist) | |
648 | { | |
649 | cerr<<"AliL3HoughMaxFinder : No histogram!"<<endl; | |
650 | return; | |
651 | } | |
652 | ||
653 | if(hist->GetNEntries() == 0) | |
654 | return; | |
655 | ||
656 | Int_t xmin = hist->GetFirstXbin(); | |
657 | Int_t xmax = hist->GetLastXbin(); | |
658 | Int_t ymin = hist->GetFirstYbin(); | |
659 | Int_t ymax = hist->GetLastYbin(); | |
660 | ||
661 | //Start by looking for pre-peaks: | |
662 | ||
663 | PreYPeak **local_maxima = new PreYPeak*[hist->GetNbinsY()]; | |
664 | ||
665 | Short_t *nmaxs = new Short_t[hist->GetNbinsY()]; | |
666 | Int_t n2dmax = 0; | |
667 | Int_t last_value,value; | |
668 | for(Int_t ybin=ymin; ybin<=ymax; ybin++) | |
669 | { | |
670 | local_maxima[ybin-ymin] = new PreYPeak[hist->GetNbinsX()]; | |
671 | nmaxs[ybin-ymin] = 0; | |
672 | last_value = 0; | |
673 | Bool_t found = 0; | |
674 | for(Int_t xbin=xmin; xbin<=xmax; xbin++) | |
675 | { | |
676 | value = hist->GetBinContent(hist->GetBin(xbin,ybin)); | |
677 | if(value >= fThreshold) | |
678 | { | |
679 | if(value > last_value) | |
680 | // if(value > (last_value + 1)) | |
681 | { | |
682 | local_maxima[ybin-ymin][nmaxs[ybin-ymin]].start_position = xbin; | |
683 | local_maxima[ybin-ymin][nmaxs[ybin-ymin]].end_position = xbin; | |
684 | local_maxima[ybin-ymin][nmaxs[ybin-ymin]].value = value; | |
685 | local_maxima[ybin-ymin][nmaxs[ybin-ymin]].prev_value = 0; | |
686 | found = 1; | |
687 | } | |
688 | if(value == last_value) | |
689 | // if(abs(value - last_value) <= 1) | |
690 | if(found) | |
691 | local_maxima[ybin-ymin][nmaxs[ybin-ymin]].end_position = xbin; | |
692 | } | |
693 | ||
694 | if((value < fThreshold) || (value < last_value)) | |
695 | // if((value < fThreshold) || (value < (last_value - 1))) | |
696 | { | |
697 | if(found) | |
698 | { | |
699 | nmaxs[ybin-ymin]++; | |
700 | found = 0; | |
701 | } | |
702 | } | |
703 | last_value = value; | |
704 | ||
705 | } | |
706 | n2dmax += nmaxs[ybin-ymin]; | |
707 | } | |
708 | ||
709 | // Pre2DPeak *maxima = new Pre2DPeak[n2dmax]; | |
710 | for(Int_t ybin=ymax; ybin >= ymin; ybin--) | |
711 | { | |
712 | for(Int_t i=0; i<nmaxs[ybin-ymin]; i++) | |
713 | { | |
714 | Int_t local_value = local_maxima[ybin-ymin][i].value; | |
715 | Int_t local_prev_value = local_maxima[ybin-ymin][i].prev_value; | |
716 | Int_t local_next_value = 0; | |
717 | ||
718 | if(local_value<0) | |
719 | continue; //already used | |
720 | ||
721 | //start expanding in the psi-direction: | |
722 | ||
723 | Int_t local_x_start = local_maxima[ybin-ymin][i].start_position; | |
724 | Int_t local_x_end = local_maxima[ybin-ymin][i].end_position; | |
725 | Int_t temp_x_start = local_maxima[ybin-ymin][i].start_position; | |
726 | Int_t temp_x_end = local_maxima[ybin-ymin][i].end_position; | |
727 | ||
728 | Int_t local_y=ybin-1,nybins=1; | |
729 | ||
730 | while(local_y >= ymin) | |
731 | { | |
732 | Bool_t found=0; | |
733 | for(Int_t j=0; j<nmaxs[local_y-ymin]; j++) | |
734 | { | |
735 | Int_t adapted_kappawindow; | |
736 | Float_t adapted_x,adapted_y; | |
737 | adapted_x = ((Float_t)local_x_start+(Float_t)local_x_end)/2.0; | |
738 | adapted_y = ((Float_t)ybin+(Float_t)(local_y+1))/2.0; | |
739 | // if(adapted_y<(-0.46*adapted_x+86) && adapted_y>(-0.46*adapted_x+56) && adapted_y<(-0.215*adapted_x+68) && adapted_y>(-0.215*adapted_x+40)) | |
740 | // if(adapted_y<(-0.215*adapted_x+68) && adapted_y>(-0.215*adapted_x+40)) | |
741 | // if((adapted_y<(-0.46*adapted_x+86) && adapted_y>(-0.46*adapted_x+56)) && !(adapted_y<(-0.215*adapted_x+68) && adapted_y>(-0.215*adapted_x+40))) | |
742 | if(((adapted_y<(-0.46*adapted_x+86) && adapted_y>(-0.46*adapted_x+56)) && !(adapted_y<(-0.215*adapted_x+68) && adapted_y>(-0.215*adapted_x+40))) || | |
743 | ((adapted_y<(-0.215*adapted_x+68) && adapted_y>(-0.215*adapted_x+40))&& !(adapted_y<(-0.46*adapted_x+86) && adapted_y>(-0.46*adapted_x+56)))) | |
744 | adapted_kappawindow = 1; | |
745 | else | |
746 | adapted_kappawindow = 0; | |
747 | //temprorary here | |
748 | // adapted_kappawindow = 0; | |
749 | ||
750 | if( (local_maxima[local_y-ymin][j].start_position <= (temp_x_end + kappawindow + adapted_kappawindow)) && (local_maxima[local_y-ymin][j].end_position >= temp_x_start)) | |
751 | { | |
752 | if( local_maxima[local_y-ymin][j].value == local_value ) | |
753 | { | |
754 | local_x_end = local_maxima[local_y-ymin][j].end_position; | |
755 | temp_x_start = local_maxima[local_y-ymin][j].start_position; | |
756 | temp_x_end = local_maxima[local_y-ymin][j].end_position; | |
757 | local_maxima[local_y-ymin][j].value = -1; | |
758 | found = 1; | |
759 | nybins++; | |
760 | break; | |
761 | } | |
762 | else | |
763 | { | |
764 | local_maxima[local_y-ymin][j].prev_value = local_value; | |
765 | local_next_value = local_maxima[local_y-ymin][j].value; | |
766 | } | |
767 | } | |
768 | } | |
769 | if(!found || local_y == ymin)//no more local maximas to be matched, so write the final peak and break the expansion: | |
770 | { | |
771 | Int_t adapted_xsize; | |
772 | Float_t adapted_x,adapted_y; | |
773 | adapted_x = ((Float_t)local_x_start+(Float_t)local_x_end)/2.0; | |
774 | adapted_y = ((Float_t)ybin+(Float_t)(local_y+1))/2.0; | |
775 | // if(adapted_y<(-0.46*adapted_x+86) && adapted_y>(-0.46*adapted_x+56) && adapted_y<(-0.215*adapted_x+68) && adapted_y>(-0.215*adapted_x+40)) | |
776 | // if(adapted_y<(-0.215*adapted_x+68) && adapted_y>(-0.215*adapted_x+40)) | |
777 | if((adapted_y<(-0.46*adapted_x+86) && adapted_y>(-0.46*adapted_x+56)) && !(adapted_y<(-0.215*adapted_x+68) && adapted_y>(-0.215*adapted_x+40))) | |
778 | adapted_xsize = 1; | |
779 | else | |
780 | adapted_xsize = 2; | |
781 | //temprorary here | |
782 | // adapted_xsize = 1; | |
783 | if((nybins > ysize) && ((local_x_end-local_x_start+1) > adapted_xsize) && (local_value > local_prev_value) && (local_value > local_next_value)) | |
784 | { | |
785 | fXPeaks[fNPeaks] = hist->GetPreciseBinCenterX(adapted_x); | |
786 | fYPeaks[fNPeaks] = hist->GetPreciseBinCenterY(adapted_y); | |
787 | fWeight[fNPeaks] = local_value; | |
788 | #ifdef do_mc | |
789 | cout<<"Peak found at: "<<((Float_t)local_x_start+(Float_t)local_x_end)/2.0<<" "<<((Float_t)ybin+(Float_t)(local_y+1))/2.0<<" "<<fXPeaks[fNPeaks]<<" "<<fYPeaks[fNPeaks]<<" with weight "<<fWeight[fNPeaks]<<" and size "<<local_x_end-local_x_start+1<<" by "<<nybins<<endl; | |
790 | #endif | |
791 | fNPeaks++; | |
792 | } | |
793 | break; | |
794 | } | |
795 | else | |
796 | local_y--;//Search continues... | |
797 | } | |
798 | } | |
799 | } | |
800 | ||
801 | for(Int_t i=0; i<hist->GetNbinsY(); i++) | |
802 | delete local_maxima[i]; | |
803 | ||
804 | delete [] local_maxima; | |
805 | delete [] nmaxs; | |
806 | //delete [] maxima; | |
807 | } | |
808 | ||
809 | void AliL3HoughMaxFinder::FindPeak1(Int_t y_window,Int_t x_bin_sides) | |
810 | { | |
811 | //Testing mutliple peakfinding. | |
812 | //The algorithm searches the histogram for prepreaks by looking in windows | |
813 | //for each bin on the xaxis. The size of these windows is controlled by y_window. | |
814 | //Then the prepreaks are sorted according to their weight (sum inside window), | |
815 | //and the peak positions are calculated by taking the weighted mean in both | |
816 | //x and y direction. The size of the peak in x-direction is controlled by x_bin_sides. | |
817 | ||
818 | if(!fCurrentHisto) | |
819 | { | |
820 | printf("AliL3HoughMaxFinder::FindPeak1 : No input histogram\n"); | |
821 | return; | |
822 | } | |
823 | if(fCurrentHisto->GetNEntries()==0) | |
824 | return; | |
825 | ||
826 | //Int_t y_window=2; | |
827 | //Int_t x_bin_sides=1; | |
828 | ||
829 | //Float_t max_kappa = 0.001; | |
830 | //Float_t max_phi0 = 0.08; | |
831 | ||
832 | Int_t max_sum=0; | |
833 | ||
834 | Int_t xmin = fCurrentHisto->GetFirstXbin(); | |
835 | Int_t xmax = fCurrentHisto->GetLastXbin(); | |
836 | Int_t ymin = fCurrentHisto->GetFirstYbin(); | |
837 | Int_t ymax = fCurrentHisto->GetLastYbin(); | |
838 | Int_t nbinsx = fCurrentHisto->GetNbinsX()+1; | |
839 | ||
840 | AxisWindow **windowPt = new AxisWindow*[nbinsx]; | |
841 | AxisWindow **anotherPt = new AxisWindow*[nbinsx]; | |
842 | ||
843 | for(Int_t i=0; i<nbinsx; i++) | |
844 | { | |
845 | windowPt[i] = new AxisWindow; | |
846 | #if defined(__DECCXX) | |
847 | bzero((char *)windowPt[i],sizeof(AxisWindow)); | |
848 | #else | |
849 | bzero((void*)windowPt[i],sizeof(AxisWindow)); | |
850 | #endif | |
851 | anotherPt[i] = windowPt[i]; | |
852 | } | |
853 | ||
854 | for(Int_t xbin=xmin; xbin<=xmax; xbin++) | |
855 | { | |
856 | max_sum = 0; | |
857 | for(Int_t ybin=ymin; ybin<=ymax-y_window; ybin++) | |
858 | { | |
859 | Int_t sum_in_window=0; | |
860 | for(Int_t b=ybin; b<ybin+y_window; b++) | |
861 | { | |
862 | //inside window | |
863 | Int_t bin = fCurrentHisto->GetBin(xbin,b); | |
864 | sum_in_window += (Int_t)fCurrentHisto->GetBinContent(bin); | |
865 | } | |
866 | ||
867 | if(sum_in_window > max_sum) | |
868 | { | |
869 | max_sum = sum_in_window; | |
870 | windowPt[xbin]->ymin = ybin; | |
871 | windowPt[xbin]->ymax = ybin + y_window; | |
872 | windowPt[xbin]->weight = sum_in_window; | |
873 | windowPt[xbin]->xbin = xbin; | |
874 | } | |
875 | } | |
876 | } | |
877 | ||
878 | //Sort the windows according to the weight | |
879 | SortPeaks(windowPt,0,nbinsx); | |
880 | ||
881 | Float_t top,butt; | |
882 | for(Int_t i=0; i<nbinsx; i++) | |
883 | { | |
884 | top=butt=0; | |
885 | Int_t xbin = windowPt[i]->xbin; | |
886 | ||
887 | if(xbin<xmin || xbin > xmax-1) continue; | |
888 | ||
889 | //Check if this is really a local maxima | |
890 | if(anotherPt[xbin-1]->weight > anotherPt[xbin]->weight || | |
891 | anotherPt[xbin+1]->weight > anotherPt[xbin]->weight) | |
892 | continue; | |
893 | ||
894 | for(Int_t j=windowPt[i]->ymin; j<windowPt[i]->ymax; j++) | |
895 | { | |
896 | //Calculate the mean in y direction: | |
897 | Int_t bin = fCurrentHisto->GetBin(windowPt[i]->xbin,j); | |
898 | top += (fCurrentHisto->GetBinCenterY(j))*(fCurrentHisto->GetBinContent(bin)); | |
899 | butt += fCurrentHisto->GetBinContent(bin); | |
900 | } | |
901 | ||
902 | if(butt < fThreshold) | |
903 | continue; | |
904 | ||
905 | fXPeaks[fNPeaks] = fCurrentHisto->GetBinCenterX(windowPt[i]->xbin); | |
906 | fYPeaks[fNPeaks] = top/butt; | |
907 | fWeight[fNPeaks] = (Int_t)butt; | |
908 | //cout<<"mean in y "<<ypeaks[n]<<" on x "<<windowPt[i]->xbin<<" content "<<butt<<endl; | |
909 | fNPeaks++; | |
910 | if(fNPeaks==fNMax) | |
911 | { | |
912 | cerr<<"AliL3HoughMaxFinder::FindPeak1 : Peak array out of range!!!"<<endl; | |
913 | break; | |
914 | } | |
915 | } | |
916 | ||
917 | ||
918 | //Improve the peaks by including the region around in x. | |
919 | Float_t ytop,ybutt; | |
920 | Int_t prev; | |
921 | Int_t w; | |
922 | for(Int_t i=0; i<fNPeaks; i++) | |
923 | { | |
924 | Int_t xbin = fCurrentHisto->FindXbin(fXPeaks[i]); | |
925 | if(xbin - x_bin_sides < xmin || xbin + x_bin_sides > xmax) continue; | |
926 | top=butt=0; | |
927 | ytop=0,ybutt=0; | |
928 | w=0; | |
929 | prev = xbin - x_bin_sides+1; | |
930 | for(Int_t j=xbin-x_bin_sides; j<=xbin+x_bin_sides; j++) | |
931 | { | |
932 | /* | |
933 | //Check if the windows are overlapping: | |
934 | if(anotherPt[j]->ymin > anotherPt[prev]->ymax) {prev=j; continue;} | |
935 | if(anotherPt[j]->ymax < anotherPt[prev]->ymin) {prev=j; continue;} | |
936 | prev = j; | |
937 | */ | |
938 | ||
939 | top += fCurrentHisto->GetBinCenterX(j)*anotherPt[j]->weight; | |
940 | butt += anotherPt[j]->weight; | |
941 | ||
942 | for(Int_t k=anotherPt[j]->ymin; k<anotherPt[j]->ymax; k++) | |
943 | { | |
944 | Int_t bin = fCurrentHisto->GetBin(j,k); | |
945 | ytop += (fCurrentHisto->GetBinCenterY(k))*(fCurrentHisto->GetBinContent(bin)); | |
946 | ybutt += fCurrentHisto->GetBinContent(bin); | |
947 | w+=(Int_t)fCurrentHisto->GetBinContent(bin); | |
948 | } | |
949 | } | |
950 | ||
951 | fXPeaks[i] = top/butt; | |
952 | fYPeaks[i] = ytop/ybutt; | |
953 | fWeight[i] = w; | |
954 | //cout<<"Setting weight "<<w<<" kappa "<<fXPeaks[i]<<" phi0 "<<fYPeaks[i]<<endl; | |
955 | ||
956 | /* | |
957 | //Check if this peak is overlapping with a previous: | |
958 | for(Int_t p=0; p<i; p++) | |
959 | { | |
960 | //cout<<fabs(fXPeaks[p] - fXPeaks[i])<<" "<<fabs(fYPeaks[p] - fYPeaks[i])<<endl; | |
961 | if(fabs(fXPeaks[p] - fXPeaks[i]) < max_kappa && | |
962 | fabs(fYPeaks[p] - fYPeaks[i]) < max_phi0) | |
963 | { | |
964 | fWeight[i]=0; | |
965 | //break; | |
966 | } | |
967 | } | |
968 | */ | |
969 | } | |
970 | ||
971 | for(Int_t i=0; i<nbinsx; i++) | |
972 | delete windowPt[i]; | |
973 | delete [] windowPt; | |
974 | delete [] anotherPt; | |
975 | } | |
976 | ||
977 | void AliL3HoughMaxFinder::SortPeaks(struct AxisWindow **a,Int_t first,Int_t last) | |
978 | { | |
979 | //General sorting routine | |
980 | //Sort according to PeakCompare() | |
981 | ||
982 | static struct AxisWindow *tmp; | |
983 | static int i; // "static" to save stack space | |
984 | int j; | |
985 | ||
986 | while (last - first > 1) { | |
987 | i = first; | |
988 | j = last; | |
989 | for (;;) { | |
990 | while (++i < last && PeakCompare(a[i], a[first]) < 0) | |
991 | ; | |
992 | while (--j > first && PeakCompare(a[j], a[first]) > 0) | |
993 | ; | |
994 | if (i >= j) | |
995 | break; | |
996 | ||
997 | tmp = a[i]; | |
998 | a[i] = a[j]; | |
999 | a[j] = tmp; | |
1000 | } | |
1001 | if (j == first) { | |
1002 | ++first; | |
1003 | continue; | |
1004 | } | |
1005 | tmp = a[first]; | |
1006 | a[first] = a[j]; | |
1007 | a[j] = tmp; | |
1008 | if (j - first < last - (j + 1)) { | |
1009 | SortPeaks(a, first, j); | |
1010 | first = j + 1; // QSort(j + 1, last); | |
1011 | } else { | |
1012 | SortPeaks(a, j + 1, last); | |
1013 | last = j; // QSort(first, j); | |
1014 | } | |
1015 | } | |
1016 | ||
1017 | } | |
1018 | ||
1019 | Int_t AliL3HoughMaxFinder::PeakCompare(struct AxisWindow *a,struct AxisWindow *b) | |
1020 | { | |
1021 | if(a->weight < b->weight) return 1; | |
1022 | if(a->weight > b->weight) return -1; | |
1023 | return 0; | |
1024 | } | |
1025 | ||
1026 | void AliL3HoughMaxFinder::FindPeak(Int_t t1,Double_t t2,Int_t t3) | |
1027 | { | |
1028 | //Attempt of a more sophisticated peak finder. | |
1029 | //Finds the best peak in the histogram, and returns the corresponding | |
1030 | //track object. | |
1031 | ||
1032 | if(!fCurrentHisto) | |
1033 | { | |
1034 | printf("AliL3HoughMaxFinder::FindPeak : No histogram!!\n"); | |
1035 | return; | |
1036 | } | |
1037 | AliL3Histogram *hist = fCurrentHisto; | |
1038 | if(hist->GetNEntries()==0) | |
1039 | return; | |
1040 | ||
1041 | Int_t xmin = hist->GetFirstXbin(); | |
1042 | Int_t xmax = hist->GetLastXbin(); | |
1043 | Int_t ymin = hist->GetFirstYbin(); | |
1044 | Int_t ymax = hist->GetLastYbin(); | |
1045 | Int_t nbinsx = hist->GetNbinsX()+1; | |
1046 | ||
1047 | Int_t *m = new Int_t[nbinsx]; | |
1048 | Int_t *m_low = new Int_t[nbinsx]; | |
1049 | Int_t *m_up = new Int_t[nbinsx]; | |
1050 | ||
1051 | ||
1052 | recompute: //this is a goto. | |
1053 | ||
1054 | for(Int_t i=0; i<nbinsx; i++) | |
1055 | { | |
1056 | m[i]=0; | |
1057 | m_low[i]=0; | |
1058 | m_up[i]=0; | |
1059 | } | |
1060 | ||
1061 | Int_t max_x=0,sum=0,max_xbin=0,bin=0; | |
1062 | ||
1063 | for(Int_t xbin=xmin; xbin<=xmax; xbin++) | |
1064 | { | |
1065 | for(Int_t ybin=ymin; ybin <= ymax - t1; ybin++) | |
1066 | { | |
1067 | sum = 0; | |
1068 | for(Int_t y=ybin; y <= ybin+t1; y++) | |
1069 | { | |
1070 | if(y>ymax) break; | |
1071 | //Inside window | |
1072 | bin = hist->GetBin(xbin,y); | |
1073 | sum += (Int_t)hist->GetBinContent(bin); | |
1074 | ||
1075 | } | |
1076 | if(sum > m[xbin]) //Max value locally in this xbin | |
1077 | { | |
1078 | m[xbin]=sum; | |
1079 | m_low[xbin]=ybin; | |
1080 | m_up[xbin]=ybin + t1; | |
1081 | } | |
1082 | ||
1083 | } | |
1084 | ||
1085 | if(m[xbin] > max_x) //Max value globally in x-direction | |
1086 | { | |
1087 | max_xbin = xbin; | |
1088 | max_x = m[xbin];//sum; | |
1089 | } | |
1090 | } | |
1091 | //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]); | |
1092 | //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[max_xbin]),hist->GetBinCenterY(m_up[max_xbin])); | |
1093 | ||
1094 | //Determine a width in the x-direction | |
1095 | Int_t x_low=0,x_up=0; | |
1096 | ||
1097 | for(Int_t xbin=max_xbin-1; xbin >= xmin; xbin--) | |
1098 | { | |
1099 | if(m[xbin] < max_x*t2) | |
1100 | { | |
1101 | x_low = xbin+1; | |
1102 | break; | |
1103 | } | |
1104 | } | |
1105 | for(Int_t xbin = max_xbin+1; xbin <=xmax; xbin++) | |
1106 | { | |
1107 | if(m[xbin] < max_x*t2) | |
1108 | { | |
1109 | x_up = xbin-1; | |
1110 | break; | |
1111 | } | |
1112 | } | |
1113 | ||
1114 | Double_t top=0,butt=0,value,x_peak; | |
1115 | if(x_up - x_low + 1 > t3) | |
1116 | { | |
1117 | t1 -= 1; | |
1118 | printf("\nxrange out if limit x_up %d x_low %d t1 %d\n\n",x_low,x_up,t1); | |
1119 | if(t1 > 1) | |
1120 | goto recompute; | |
1121 | else | |
1122 | { | |
1123 | x_peak = hist->GetBinCenterX(max_xbin); | |
1124 | goto moveon; | |
1125 | } | |
1126 | } | |
1127 | ||
1128 | //printf("xlow %f xup %f\n",hist->GetBinCenterX(x_low),hist->GetBinCenterX(x_up)); | |
1129 | //printf("Spread in x %d\n",x_up-x_low +1); | |
1130 | ||
1131 | //Now, calculate the center of mass in x-direction | |
1132 | for(Int_t xbin=x_low; xbin <= x_up; xbin++) | |
1133 | { | |
1134 | value = hist->GetBinCenterX(xbin); | |
1135 | top += value*m[xbin]; | |
1136 | butt += m[xbin]; | |
1137 | } | |
1138 | x_peak = top/butt; | |
1139 | ||
1140 | moveon: | |
1141 | ||
1142 | //Find the peak in y direction: | |
1143 | Int_t x_l = hist->FindXbin(x_peak); | |
1144 | if(hist->GetBinCenterX(x_l) > x_peak) | |
1145 | x_l--; | |
1146 | ||
1147 | Int_t x_u = x_l + 1; | |
1148 | ||
1149 | if(hist->GetBinCenterX(x_l) > x_peak || hist->GetBinCenterX(x_u) <= x_peak) | |
1150 | printf("\nAliL3HoughMaxFinder::FindPeak : Wrong xrange %f %f %f\n\n",hist->GetBinCenterX(x_l),x_peak,hist->GetBinCenterX(x_u)); | |
1151 | ||
1152 | //printf("\nxlow %f xup %f\n",hist->GetBinCenterX(x_l),hist->GetBinCenterX(x_u)); | |
1153 | ||
1154 | value=top=butt=0; | |
1155 | ||
1156 | //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_l]),hist->GetBinCenterY(m_up[x_l])); | |
1157 | //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_u]),hist->GetBinCenterY(m_up[x_u])); | |
1158 | ||
1159 | for(Int_t ybin=m_low[x_l]; ybin <= m_up[x_l]; ybin++) | |
1160 | { | |
1161 | value = hist->GetBinCenterY(ybin); | |
1162 | bin = hist->GetBin(x_l,ybin); | |
1163 | top += value*hist->GetBinContent(bin); | |
1164 | butt += hist->GetBinContent(bin); | |
1165 | } | |
1166 | Double_t y_peak_low = top/butt; | |
1167 | ||
1168 | //printf("y_peak_low %f\n",y_peak_low); | |
1169 | ||
1170 | value=top=butt=0; | |
1171 | for(Int_t ybin=m_low[x_u]; ybin <= m_up[x_u]; ybin++) | |
1172 | { | |
1173 | value = hist->GetBinCenterY(ybin); | |
1174 | bin = hist->GetBin(x_u,ybin); | |
1175 | top += value*hist->GetBinContent(bin); | |
1176 | butt += hist->GetBinContent(bin); | |
1177 | } | |
1178 | Double_t y_peak_up = top/butt; | |
1179 | ||
1180 | //printf("y_peak_up %f\n",y_peak_up); | |
1181 | ||
1182 | Double_t x_value_up = hist->GetBinCenterX(x_u); | |
1183 | Double_t x_value_low = hist->GetBinCenterX(x_l); | |
1184 | ||
1185 | 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); | |
1186 | ||
1187 | ||
1188 | //Find the weight: | |
1189 | //bin = hist->FindBin(x_peak,y_peak); | |
1190 | //Int_t weight = (Int_t)hist->GetBinContent(bin); | |
1191 | ||
1192 | //AliL3HoughTrack *track = new AliL3HoughTrack(); | |
1193 | //track->SetTrackParameters(x_peak,y_peak,weight); | |
1194 | fXPeaks[fNPeaks]=x_peak; | |
1195 | fYPeaks[fNPeaks]=y_peak; | |
1196 | fWeight[fNPeaks]=(Int_t)hist->GetBinContent(bin); | |
1197 | fNPeaks++; | |
1198 | ||
1199 | delete [] m; | |
1200 | delete [] m_low; | |
1201 | delete [] m_up; | |
1202 | ||
1203 | //return track; | |
1204 | } | |
1205 |