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