]>
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 | |
62bb4b3d | 6 | #include <strings.h> |
e06900d5 | 7 | #include "AliL3StandardIncludes.h" |
8 | ||
2061190a | 9 | #ifndef no_root |
208b54c5 | 10 | #include <TNtuple.h> |
11 | #include <TFile.h> | |
2061190a | 12 | #endif |
208b54c5 | 13 | |
e06900d5 | 14 | #include "AliL3Logging.h" |
15 | #include "AliL3HoughMaxFinder.h" | |
4cafa5fc | 16 | #include "AliL3Histogram.h" |
4de874d1 | 17 | #include "AliL3TrackArray.h" |
18 | #include "AliL3HoughTrack.h" | |
de3c3890 | 19 | #include "AliL3Transform.h" |
20 | #include "AliL3HoughTransformerRow.h" | |
e06900d5 | 21 | |
0bd0c1ef | 22 | #if __GNUC__ == 3 |
e06900d5 | 23 | using namespace std; |
24 | #endif | |
4de874d1 | 25 | |
0bd0c1ef | 26 | /** \class AliL3HoughMaxFinder |
27 | <pre> | |
b1886074 | 28 | //_____________________________________________________________ |
29 | // AliL3HoughMaxFinder | |
30 | // | |
31 | // Maximum finder | |
0bd0c1ef | 32 | // |
33 | </pre> | |
34 | */ | |
b1886074 | 35 | |
4de874d1 | 36 | ClassImp(AliL3HoughMaxFinder) |
37 | ||
38 | ||
39 | AliL3HoughMaxFinder::AliL3HoughMaxFinder() | |
40 | { | |
41 | //Default constructor | |
42 | fThreshold = 0; | |
de3c3890 | 43 | fCurrentEtaSlice = -1; |
4de874d1 | 44 | fHistoType=0; |
3fe49b5b | 45 | fXPeaks=0; |
46 | fYPeaks=0; | |
de3c3890 | 47 | fWeight=0; |
3fe49b5b | 48 | fNPeaks=0; |
de3c3890 | 49 | fN1PeaksPrevEtaSlice=0; |
50 | fN2PeaksPrevEtaSlice=0; | |
51 | fSTARTXPeaks=0; | |
52 | fSTARTYPeaks=0; | |
53 | fENDXPeaks=0; | |
54 | fENDYPeaks=0; | |
55 | fSTARTETAPeaks=0; | |
56 | fENDETAPeaks=0; | |
3fe49b5b | 57 | fNMax=0; |
b2a02bce | 58 | fGradX=1; |
59 | fGradY=1; | |
2061190a | 60 | #ifndef no_root |
208b54c5 | 61 | fNtuppel = 0; |
2061190a | 62 | #endif |
4de874d1 | 63 | } |
64 | ||
3fe49b5b | 65 | AliL3HoughMaxFinder::AliL3HoughMaxFinder(Char_t *histotype,Int_t nmax,AliL3Histogram *hist) |
4de874d1 | 66 | { |
67 | //Constructor | |
68 | ||
69 | //fTracks = new AliL3TrackArray("AliL3HoughTrack"); | |
70 | if(strcmp(histotype,"KappaPhi")==0) fHistoType='c'; | |
71 | if(strcmp(histotype,"DPsi")==0) fHistoType='l'; | |
de3c3890 | 72 | |
73 | fCurrentEtaSlice = -1; | |
4de874d1 | 74 | |
4cafa5fc | 75 | if(hist) |
76 | fCurrentHisto = hist; | |
3fe49b5b | 77 | |
b2a02bce | 78 | fGradX=1; |
79 | fGradY=1; | |
3fe49b5b | 80 | fNMax=nmax; |
81 | fXPeaks = new Float_t[fNMax]; | |
82 | fYPeaks = new Float_t[fNMax]; | |
83 | fWeight = new Int_t[fNMax]; | |
de3c3890 | 84 | fSTARTXPeaks = new Int_t[fNMax]; |
85 | fSTARTYPeaks = new Int_t[fNMax]; | |
86 | fENDXPeaks = new Int_t[fNMax]; | |
87 | fENDYPeaks = new Int_t[fNMax]; | |
88 | fSTARTETAPeaks = new Int_t[fNMax]; | |
89 | fENDETAPeaks = new Int_t[fNMax]; | |
2061190a | 90 | #ifndef no_root |
208b54c5 | 91 | fNtuppel = 0; |
2061190a | 92 | #endif |
208b54c5 | 93 | fThreshold=0; |
4de874d1 | 94 | } |
95 | ||
4de874d1 | 96 | AliL3HoughMaxFinder::~AliL3HoughMaxFinder() |
97 | { | |
98 | //Destructor | |
3fe49b5b | 99 | if(fXPeaks) |
100 | delete [] fXPeaks; | |
101 | if(fYPeaks) | |
102 | delete [] fYPeaks; | |
103 | if(fWeight) | |
104 | delete [] fWeight; | |
de3c3890 | 105 | if(fSTARTXPeaks) |
106 | delete [] fSTARTXPeaks; | |
107 | if(fSTARTYPeaks) | |
108 | delete [] fSTARTYPeaks; | |
109 | if(fENDXPeaks) | |
110 | delete [] fENDXPeaks; | |
111 | if(fENDYPeaks) | |
112 | delete [] fENDYPeaks; | |
113 | if(fSTARTETAPeaks) | |
114 | delete [] fSTARTETAPeaks; | |
115 | if(fENDETAPeaks) | |
116 | delete [] fENDETAPeaks; | |
2061190a | 117 | #ifndef no_root |
208b54c5 | 118 | if(fNtuppel) |
119 | delete fNtuppel; | |
2061190a | 120 | #endif |
3fe49b5b | 121 | } |
4de874d1 | 122 | |
3fe49b5b | 123 | void AliL3HoughMaxFinder::Reset() |
124 | { | |
917e711b | 125 | // Method to reinit the Peak Finder |
3fe49b5b | 126 | for(Int_t i=0; i<fNMax; i++) |
127 | { | |
128 | fXPeaks[i]=0; | |
129 | fYPeaks[i]=0; | |
130 | fWeight[i]=0; | |
de3c3890 | 131 | fSTARTXPeaks[i]=0; |
132 | fSTARTYPeaks[i]=0; | |
133 | fENDXPeaks[i]=0; | |
134 | fENDYPeaks[i]=0; | |
135 | fSTARTETAPeaks[i]=0; | |
136 | fENDETAPeaks[i]=0; | |
3fe49b5b | 137 | } |
138 | fNPeaks=0; | |
de3c3890 | 139 | fN1PeaksPrevEtaSlice=0; |
140 | fN2PeaksPrevEtaSlice=0; | |
4de874d1 | 141 | } |
142 | ||
208b54c5 | 143 | void AliL3HoughMaxFinder::CreateNtuppel() |
144 | { | |
917e711b | 145 | // Fill a NTuple with the peak parameters |
2061190a | 146 | #ifndef no_root |
208b54c5 | 147 | //content#; neighbouring bins of the peak. |
2061190a | 148 | fNtuppel = new TNtuple("ntuppel","Peak charateristics","kappa:phi0:weigth:content3:content5:content1:content7"); |
3e87ef69 | 149 | fNtuppel->SetDirectory(0); |
2061190a | 150 | #endif |
208b54c5 | 151 | } |
152 | ||
153 | void AliL3HoughMaxFinder::WriteNtuppel(Char_t *filename) | |
154 | { | |
917e711b | 155 | // Write the NTuple with the peak parameters |
2061190a | 156 | #ifndef no_root |
208b54c5 | 157 | TFile *file = TFile::Open(filename,"RECREATE"); |
158 | if(!file) | |
159 | { | |
160 | cerr<<"AliL3HoughMaxFinder::WriteNtuppel : Error opening file "<<filename<<endl; | |
161 | return; | |
162 | } | |
163 | fNtuppel->Write(); | |
164 | file->Close(); | |
2061190a | 165 | #endif |
208b54c5 | 166 | } |
167 | ||
3fe49b5b | 168 | void AliL3HoughMaxFinder::FindAbsMaxima() |
52a2a604 | 169 | { |
917e711b | 170 | // Simple Peak Finder in the Hough space |
4cafa5fc | 171 | if(!fCurrentHisto) |
172 | { | |
208b54c5 | 173 | cerr<<"AliL3HoughMaxFinder::FindAbsMaxima : No histogram"<<endl; |
4fc9a6a4 | 174 | return; |
4cafa5fc | 175 | } |
176 | AliL3Histogram *hist = fCurrentHisto; | |
4fc9a6a4 | 177 | |
b2a02bce | 178 | if(hist->GetNEntries() == 0) |
179 | return; | |
180 | ||
4cafa5fc | 181 | Int_t xmin = hist->GetFirstXbin(); |
182 | Int_t xmax = hist->GetLastXbin(); | |
183 | Int_t ymin = hist->GetFirstYbin(); | |
184 | Int_t ymax = hist->GetLastYbin(); | |
185 | Int_t bin; | |
917e711b | 186 | Double_t value,maxvalue=0; |
3fe49b5b | 187 | |
917e711b | 188 | Int_t maxxbin=0,maxybin=0; |
4cafa5fc | 189 | for(Int_t xbin=xmin; xbin<=xmax; xbin++) |
190 | { | |
191 | for(Int_t ybin=ymin; ybin<=ymax; ybin++) | |
192 | { | |
193 | bin = hist->GetBin(xbin,ybin); | |
194 | value = hist->GetBinContent(bin); | |
917e711b | 195 | if(value>maxvalue) |
4cafa5fc | 196 | { |
917e711b | 197 | maxvalue = value; |
198 | maxxbin = xbin; | |
199 | maxybin = ybin; | |
4cafa5fc | 200 | } |
201 | } | |
202 | } | |
4fc9a6a4 | 203 | |
917e711b | 204 | if(maxvalue == 0) |
b2a02bce | 205 | return; |
206 | ||
3fe49b5b | 207 | if(fNPeaks > fNMax) |
208 | { | |
209 | cerr<<"AliL3HoughMaxFinder::FindAbsMaxima : Array out of range : "<<fNPeaks<<endl; | |
210 | return; | |
211 | } | |
212 | ||
917e711b | 213 | Double_t maxx = hist->GetBinCenterX(maxxbin); |
214 | Double_t maxy = hist->GetBinCenterY(maxybin); | |
215 | fXPeaks[fNPeaks] = maxx; | |
216 | fYPeaks[fNPeaks] = maxy; | |
217 | fWeight[fNPeaks] = (Int_t)maxvalue; | |
b2a02bce | 218 | |
3fe49b5b | 219 | fNPeaks++; |
2061190a | 220 | #ifndef no_root |
208b54c5 | 221 | if(fNtuppel) |
222 | { | |
917e711b | 223 | Int_t bin3 = hist->GetBin(maxxbin-1,maxybin); |
224 | Int_t bin5 = hist->GetBin(maxxbin+1,maxybin); | |
225 | Int_t bin1 = hist->GetBin(maxxbin,maxybin-1); | |
226 | Int_t bin7 = hist->GetBin(maxxbin,maxybin+1); | |
208b54c5 | 227 | |
917e711b | 228 | fNtuppel->Fill(maxx,maxy,maxvalue,hist->GetBinContent(bin3),hist->GetBinContent(bin5),hist->GetBinContent(bin1),hist->GetBinContent(bin7)); |
208b54c5 | 229 | } |
2061190a | 230 | #endif |
4cafa5fc | 231 | } |
232 | ||
3fe49b5b | 233 | void AliL3HoughMaxFinder::FindBigMaxima() |
4cafa5fc | 234 | { |
917e711b | 235 | // Another Peak finder |
3fe49b5b | 236 | AliL3Histogram *hist = fCurrentHisto; |
b2a02bce | 237 | |
238 | if(hist->GetNEntries() == 0) | |
239 | return; | |
240 | ||
4cafa5fc | 241 | Int_t xmin = hist->GetFirstXbin(); |
242 | Int_t xmax = hist->GetLastXbin(); | |
243 | Int_t ymin = hist->GetFirstYbin(); | |
244 | Int_t ymax = hist->GetLastYbin(); | |
917e711b | 245 | Int_t bin[25],binindex; |
95a00d93 | 246 | Double_t value[25]; |
52a2a604 | 247 | |
52a2a604 | 248 | for(Int_t xbin=xmin+2; xbin<xmax-3; xbin++) |
249 | { | |
250 | for(Int_t ybin=ymin+2; ybin<ymax-3; ybin++) | |
251 | { | |
917e711b | 252 | binindex=0; |
4cafa5fc | 253 | for(Int_t xb=xbin-2; xb<=xbin+2; xb++) |
52a2a604 | 254 | { |
4cafa5fc | 255 | for(Int_t yb=ybin-2; yb<=ybin+2; yb++) |
52a2a604 | 256 | { |
917e711b | 257 | bin[binindex]=hist->GetBin(xb,yb); |
258 | value[binindex]=hist->GetBinContent(bin[binindex]); | |
259 | binindex++; | |
52a2a604 | 260 | } |
52a2a604 | 261 | } |
262 | if(value[12]==0) continue; | |
263 | Int_t b=0; | |
264 | while(1) | |
265 | { | |
917e711b | 266 | if(value[b] > value[12] || b==binindex) break; |
52a2a604 | 267 | b++; |
3fe49b5b | 268 | //printf("b %d\n",b); |
52a2a604 | 269 | } |
917e711b | 270 | if(b == binindex) |
52a2a604 | 271 | { |
272 | //Found maxima | |
3fe49b5b | 273 | if(fNPeaks > fNMax) |
274 | { | |
275 | cerr<<"AliL3HoughMaxFinder::FindBigMaxima : Array out of range "<<fNPeaks<<endl; | |
276 | return; | |
277 | } | |
278 | ||
917e711b | 279 | Double_t maxx = hist->GetBinCenterX(xbin); |
280 | Double_t maxy = hist->GetBinCenterY(ybin); | |
281 | fXPeaks[fNPeaks] = maxx; | |
282 | fYPeaks[fNPeaks] = maxy; | |
3fe49b5b | 283 | fNPeaks++; |
52a2a604 | 284 | } |
285 | } | |
286 | } | |
52a2a604 | 287 | } |
288 | ||
b2a02bce | 289 | void AliL3HoughMaxFinder::FindMaxima(Int_t threshold) |
4de874d1 | 290 | { |
291 | //Locate all the maxima in input histogram. | |
292 | //Maxima is defined as bins with more entries than the | |
293 | //immediately neighbouring bins. | |
294 | ||
b2a02bce | 295 | if(fCurrentHisto->GetNEntries() == 0) |
296 | return; | |
297 | ||
4ab9f8f0 | 298 | Int_t xmin = fCurrentHisto->GetFirstXbin(); |
299 | Int_t xmax = fCurrentHisto->GetLastXbin(); | |
300 | Int_t ymin = fCurrentHisto->GetFirstYbin(); | |
301 | Int_t ymax = fCurrentHisto->GetLastYbin(); | |
302 | Int_t bin[9]; | |
95a00d93 | 303 | Double_t value[9]; |
4de874d1 | 304 | |
3e87ef69 | 305 | //Float_t max_kappa = 0.001; |
306 | //Float_t max_phi0 = 0.08; | |
307 | ||
308 | for(Int_t xbin=xmin+1; xbin<=xmax-1; xbin++) | |
4de874d1 | 309 | { |
3e87ef69 | 310 | for(Int_t ybin=ymin+1; ybin<=ymax-1; ybin++) |
4de874d1 | 311 | { |
4ab9f8f0 | 312 | bin[0] = fCurrentHisto->GetBin(xbin-1,ybin-1); |
313 | bin[1] = fCurrentHisto->GetBin(xbin,ybin-1); | |
314 | bin[2] = fCurrentHisto->GetBin(xbin+1,ybin-1); | |
315 | bin[3] = fCurrentHisto->GetBin(xbin-1,ybin); | |
316 | bin[4] = fCurrentHisto->GetBin(xbin,ybin); | |
317 | bin[5] = fCurrentHisto->GetBin(xbin+1,ybin); | |
318 | bin[6] = fCurrentHisto->GetBin(xbin-1,ybin+1); | |
319 | bin[7] = fCurrentHisto->GetBin(xbin,ybin+1); | |
320 | bin[8] = fCurrentHisto->GetBin(xbin+1,ybin+1); | |
321 | value[0] = fCurrentHisto->GetBinContent(bin[0]); | |
322 | value[1] = fCurrentHisto->GetBinContent(bin[1]); | |
323 | value[2] = fCurrentHisto->GetBinContent(bin[2]); | |
324 | value[3] = fCurrentHisto->GetBinContent(bin[3]); | |
325 | value[4] = fCurrentHisto->GetBinContent(bin[4]); | |
326 | value[5] = fCurrentHisto->GetBinContent(bin[5]); | |
327 | value[6] = fCurrentHisto->GetBinContent(bin[6]); | |
328 | value[7] = fCurrentHisto->GetBinContent(bin[7]); | |
329 | value[8] = fCurrentHisto->GetBinContent(bin[8]); | |
330 | ||
4de874d1 | 331 | |
4de874d1 | 332 | if(value[4]>value[0] && value[4]>value[1] && value[4]>value[2] |
333 | && value[4]>value[3] && value[4]>value[5] && value[4]>value[6] | |
334 | && value[4]>value[7] && value[4]>value[8]) | |
335 | { | |
336 | //Found a local maxima | |
917e711b | 337 | Float_t maxx = fCurrentHisto->GetBinCenterX(xbin); |
338 | Float_t maxy = fCurrentHisto->GetBinCenterY(ybin); | |
4de874d1 | 339 | |
b2a02bce | 340 | if((Int_t)value[4] <= threshold) continue;//central bin below threshold |
0309a5ee | 341 | if(fNPeaks >= fNMax) |
4de874d1 | 342 | { |
3060f008 | 343 | cout<<"AliL3HoughMaxFinder::FindMaxima : Array out of range "<<fNPeaks<<endl; |
4ab9f8f0 | 344 | return; |
4de874d1 | 345 | } |
3e87ef69 | 346 | |
3fe49b5b | 347 | //Check the gradient: |
b2a02bce | 348 | if(value[3]/value[4] > fGradX && value[5]/value[4] > fGradX) |
349 | continue; | |
350 | ||
351 | if(value[1]/value[4] > fGradY && value[7]/value[4] > fGradY) | |
352 | continue; | |
353 | ||
917e711b | 354 | fXPeaks[fNPeaks] = maxx; |
355 | fYPeaks[fNPeaks] = maxy; | |
3fe49b5b | 356 | fWeight[fNPeaks] = (Int_t)value[4]; |
357 | fNPeaks++; | |
3e87ef69 | 358 | |
3fe49b5b | 359 | /* |
4ab9f8f0 | 360 | //Check if the peak is overlapping with a previous: |
361 | Bool_t bigger = kFALSE; | |
3e87ef69 | 362 | for(Int_t p=0; p<fNPeaks; p++) |
363 | { | |
917e711b | 364 | if(fabs(maxx - fXPeaks[p]) < max_kappa && fabs(maxy - fYPeaks[p]) < max_phi0) |
3e87ef69 | 365 | { |
366 | bigger = kTRUE; | |
367 | if(value[4] > fWeight[p]) //this peak is bigger. | |
4ab9f8f0 | 368 | { |
917e711b | 369 | fXPeaks[p] = maxx; |
370 | fYPeaks[p] = maxy; | |
3e87ef69 | 371 | fWeight[p] = (Int_t)value[4]; |
4ab9f8f0 | 372 | } |
373 | else | |
374 | continue; //previous peak is bigger. | |
375 | } | |
376 | } | |
377 | if(!bigger) //there were no overlapping peaks. | |
378 | { | |
917e711b | 379 | fXPeaks[fNPeaks] = maxx; |
380 | fYPeaks[fNPeaks] = maxy; | |
3e87ef69 | 381 | fWeight[fNPeaks] = (Int_t)value[4]; |
382 | fNPeaks++; | |
4ab9f8f0 | 383 | } |
3fe49b5b | 384 | */ |
4de874d1 | 385 | } |
4de874d1 | 386 | } |
387 | } | |
4ab9f8f0 | 388 | |
4de874d1 | 389 | } |
e4c21048 | 390 | |
917e711b | 391 | struct AliL3Window |
b2a02bce | 392 | { |
917e711b | 393 | Int_t fStart; // Start |
394 | Int_t fSum; // Sum | |
b2a02bce | 395 | }; |
396 | ||
917e711b | 397 | void AliL3HoughMaxFinder::FindAdaptedPeaks(Int_t kappawindow,Float_t cutratio) |
52a2a604 | 398 | { |
b2a02bce | 399 | //Peak finder which looks for peaks with a certain shape. |
400 | //The first step involves a pre-peak finder, which looks for peaks | |
401 | //in windows (size controlled by kappawindow) summing over each psi-bin. | |
402 | //These pre-preaks are then matched between neighbouring kappa-bins to | |
403 | //look for real 2D peaks exhbiting the typical cross-shape in the Hough circle transform. | |
404 | //The maximum bin within this region is marked as the peak itself, and | |
405 | //a few checks is performed to avoid the clear fake peaks (asymmetry check etc.) | |
52a2a604 | 406 | |
52a2a604 | 407 | |
b2a02bce | 408 | AliL3Histogram *hist = fCurrentHisto; |
52a2a604 | 409 | |
b2a02bce | 410 | if(!hist) |
52a2a604 | 411 | { |
b2a02bce | 412 | cerr<<"AliL3HoughMaxFinder : No histogram!"<<endl; |
413 | return; | |
52a2a604 | 414 | } |
b2a02bce | 415 | |
416 | if(hist->GetNEntries() == 0) | |
417 | return; | |
0bd0c1ef | 418 | |
4cafa5fc | 419 | Int_t xmin = hist->GetFirstXbin(); |
420 | Int_t xmax = hist->GetLastXbin(); | |
421 | Int_t ymin = hist->GetFirstYbin(); | |
422 | Int_t ymax = hist->GetLastYbin(); | |
0bd0c1ef | 423 | |
424 | ||
b2a02bce | 425 | //Start by looking for pre-peaks: |
4cafa5fc | 426 | |
917e711b | 427 | AliL3Window **localmaxima = new AliL3Window*[hist->GetNbinsY()]; |
4cafa5fc | 428 | |
b2a02bce | 429 | Short_t *nmaxs = new Short_t[hist->GetNbinsY()]; |
917e711b | 430 | Int_t n,lastsum,sum; |
431 | Bool_t sumwasrising; | |
b2a02bce | 432 | for(Int_t ybin=ymin; ybin<=ymax; ybin++) |
52a2a604 | 433 | { |
917e711b | 434 | localmaxima[ybin-ymin] = new AliL3Window[hist->GetNbinsX()]; |
b2a02bce | 435 | nmaxs[ybin-ymin] = 0; |
917e711b | 436 | sumwasrising=0; |
437 | lastsum=0; | |
b2a02bce | 438 | n=0; |
439 | for(Int_t xbin=xmin; xbin<=xmax-kappawindow; xbin++) | |
52a2a604 | 440 | { |
b2a02bce | 441 | sum=0; |
442 | for(Int_t lbin=xbin; lbin<xbin+kappawindow; lbin++) | |
443 | sum += hist->GetBinContent(hist->GetBin(lbin,ybin)); | |
444 | ||
917e711b | 445 | if(sum < lastsum) |
52a2a604 | 446 | { |
b2a02bce | 447 | if(sum > fThreshold) |
917e711b | 448 | if(sumwasrising)//Previous sum was a local maxima |
b2a02bce | 449 | { |
917e711b | 450 | localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fStart = xbin-1; |
451 | localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fSum = lastsum; | |
b2a02bce | 452 | nmaxs[ybin-ymin]++; |
453 | } | |
52a2a604 | 454 | |
917e711b | 455 | sumwasrising=0; |
52a2a604 | 456 | } |
b2a02bce | 457 | else if(sum > 0) |
917e711b | 458 | sumwasrising=1; |
459 | lastsum=sum; | |
4cafa5fc | 460 | } |
461 | } | |
52a2a604 | 462 | |
b2a02bce | 463 | Int_t match=0; |
464 | Int_t *starts = new Int_t[hist->GetNbinsY()+1]; | |
465 | Int_t *maxs = new Int_t[hist->GetNbinsY()+1]; | |
52a2a604 | 466 | |
b2a02bce | 467 | for(Int_t ybin=ymax; ybin >= ymin+1; ybin--) |
4cafa5fc | 468 | { |
b2a02bce | 469 | for(Int_t i=0; i<nmaxs[ybin-ymin]; i++) |
4cafa5fc | 470 | { |
917e711b | 471 | Int_t lw = localmaxima[ybin-ymin][i].fSum; |
b2a02bce | 472 | |
473 | if(lw<0) | |
474 | continue; //already used | |
475 | ||
476 | Int_t maxvalue=0,maxybin=0,maxxbin=0,maxwindow=0; | |
917e711b | 477 | for(Int_t k=localmaxima[ybin-ymin][i].fStart; k<localmaxima[ybin-ymin][i].fStart + kappawindow; k++) |
b2a02bce | 478 | if(hist->GetBinContent(hist->GetBin(k,ybin)) > maxvalue) |
479 | { | |
480 | maxvalue = hist->GetBinContent(hist->GetBin(k,ybin)); | |
481 | maxybin = ybin; | |
482 | maxxbin = k; | |
483 | } | |
484 | ||
485 | //start expanding in the psi-direction: | |
486 | ||
917e711b | 487 | Int_t lb = localmaxima[ybin-ymin][i].fStart; |
b2a02bce | 488 | //Int_t ystart=ybin; |
917e711b | 489 | starts[ybin] = localmaxima[ybin-ymin][i].fStart; |
b2a02bce | 490 | maxs[ybin] = maxxbin; |
491 | Int_t yl=ybin-1,nybins=1; | |
492 | ||
917e711b | 493 | //cout<<"Starting search at ybin "<<ybin<<" start "<<lb<<" with sum "<<localmaxima[ybin-ymin][i].sum<<endl; |
b2a02bce | 494 | while(yl >= ymin) |
495 | { | |
496 | Bool_t found=0; | |
497 | for(Int_t j=0; j<nmaxs[yl-ymin]; j++) | |
498 | { | |
917e711b | 499 | if( localmaxima[yl-ymin][j].fStart - lb < 0) continue; |
500 | if( localmaxima[yl-ymin][j].fStart < lb + kappawindow + match && | |
501 | localmaxima[yl-ymin][j].fStart >= lb && localmaxima[yl-ymin][j].fSum > 0) | |
b2a02bce | 502 | { |
503 | ||
917e711b | 504 | //cout<<"match at ybin "<<yl<<" yvalue "<<hist->GetBinCenterY(yl)<<" start "<<localmaxima[yl-ymin][j].start<<" sum "<<localmaxima[yl-ymin][j].sum<<endl; |
b2a02bce | 505 | |
506 | Int_t lmaxvalue=0,lmaxxbin=0; | |
917e711b | 507 | for(Int_t k=localmaxima[yl-ymin][j].fStart; k<localmaxima[yl-ymin][j].fStart + kappawindow; k++) |
b2a02bce | 508 | { |
509 | if(hist->GetBinContent(hist->GetBin(k,yl)) > maxvalue) | |
510 | { | |
511 | maxvalue = hist->GetBinContent(hist->GetBin(k,yl)); | |
512 | maxxbin = k; | |
513 | maxybin = yl; | |
514 | maxwindow = j; | |
515 | } | |
516 | if(hist->GetBinContent(hist->GetBin(k,yl)) > lmaxvalue)//local maxima value | |
517 | { | |
518 | lmaxvalue=hist->GetBinContent(hist->GetBin(k,yl)); | |
519 | lmaxxbin=k; | |
520 | } | |
521 | } | |
522 | nybins++; | |
917e711b | 523 | starts[yl] = localmaxima[yl-ymin][j].fStart; |
b2a02bce | 524 | maxs[yl] = lmaxxbin; |
917e711b | 525 | localmaxima[yl-ymin][j].fSum=-1; //Mark as used |
b2a02bce | 526 | found=1; |
917e711b | 527 | lb = localmaxima[yl-ymin][j].fStart; |
b2a02bce | 528 | break;//Since we found a match in this bin, we dont have to search it anymore, goto next bin. |
529 | } | |
530 | } | |
531 | if(!found || yl == ymin)//no more local maximas to be matched, so write the final peak and break the expansion: | |
532 | { | |
533 | if(nybins > 4) | |
534 | { | |
535 | //cout<<"ystart "<<ystart<<" and nybins "<<nybins<<endl; | |
536 | ||
537 | Bool_t truepeak=kTRUE; | |
538 | ||
539 | //cout<<"Maxima found at xbin "<<maxxbin<<" ybin "<<maxybin<<" value "<<maxvalue<<endl; | |
540 | //cout<<"Starting to sum at xbin "<<starts[maxybin-ymin]<<endl; | |
541 | ||
542 | ||
543 | //Look in a window on both sides to probe the asymmetry | |
544 | Float_t right=0,left=0; | |
545 | for(Int_t w=maxxbin+1; w<=maxxbin+3; w++) | |
546 | { | |
547 | for(Int_t r=maxybin+1; r<=maxybin+3; r++) | |
548 | { | |
549 | right += (Float_t)hist->GetBinContent(hist->GetBin(w,r)); | |
550 | } | |
551 | } | |
552 | ||
553 | for(Int_t w=maxxbin-1; w>=maxxbin-3; w--) | |
554 | { | |
555 | for(Int_t r=maxybin+1; r<=maxybin+3; r++) | |
556 | { | |
557 | left += (Float_t)hist->GetBinContent(hist->GetBin(w,r)); | |
558 | } | |
559 | } | |
560 | ||
561 | //cout<<"ratio "<<right/left<<endl; | |
562 | ||
917e711b | 563 | Float_t upperratio=1,lowerratio=1; |
b2a02bce | 564 | if(left) |
917e711b | 565 | upperratio = right/left; |
b2a02bce | 566 | |
567 | right=left=0; | |
568 | for(Int_t w=maxxbin+1; w<=maxxbin+3; w++) | |
569 | { | |
570 | for(Int_t r=maxybin-1; r>=maxybin-3; r--) | |
571 | { | |
572 | right += (Float_t)hist->GetBinContent(hist->GetBin(w,r)); | |
573 | } | |
574 | } | |
575 | ||
576 | for(Int_t w=maxxbin-1; w>=maxxbin-3; w--) | |
577 | { | |
578 | for(Int_t r=maxybin-1; r>=maxybin-3; r--) | |
579 | { | |
580 | left += (Float_t)hist->GetBinContent(hist->GetBin(w,r)); | |
581 | } | |
582 | } | |
583 | ||
584 | //cout<<"ratio "<<left/right<<endl; | |
585 | ||
586 | if(right) | |
917e711b | 587 | lowerratio = left/right; |
b2a02bce | 588 | |
917e711b | 589 | if(upperratio > cutratio || lowerratio > cutratio) |
b2a02bce | 590 | truepeak=kFALSE; |
591 | ||
592 | if(truepeak) | |
593 | { | |
594 | ||
595 | fXPeaks[fNPeaks] = hist->GetBinCenterX(maxxbin); | |
596 | fYPeaks[fNPeaks] = hist->GetBinCenterY(maxybin); | |
597 | fWeight[fNPeaks] = maxvalue; | |
598 | fNPeaks++; | |
599 | ||
600 | /* | |
601 | //Calculate the peak using weigthed means: | |
602 | Float_t sum=0; | |
603 | fYPeaks[fNPeaks]=0; | |
604 | for(Int_t k=maxybin-1; k<=maxybin+1; k++) | |
605 | { | |
606 | Float_t lsum = 0; | |
607 | for(Int_t l=starts[k]; l<starts[k]+kappawindow; l++) | |
608 | { | |
609 | lsum += (Float_t)hist->GetBinContent(hist->GetBin(l,k)); | |
610 | sum += (Float_t)hist->GetBinContent(hist->GetBin(l,k)); | |
611 | } | |
612 | fYPeaks[fNPeaks] += lsum*hist->GetBinCenterY(k); | |
613 | } | |
614 | fYPeaks[fNPeaks] /= sum; | |
615 | Int_t ybin1,ybin2; | |
616 | if(fYPeaks[fNPeaks] < hist->GetBinCenterY(hist->FindYbin(fYPeaks[fNPeaks]))) | |
617 | { | |
618 | ybin1 = hist->FindYbin(fYPeaks[fNPeaks])-1; | |
619 | ybin2 = ybin1+1; | |
620 | } | |
621 | else | |
622 | { | |
623 | ybin1 = hist->FindYbin(fYPeaks[fNPeaks]); | |
624 | ybin2 = ybin1+1; | |
625 | } | |
626 | ||
627 | Float_t kappa1=0,kappa2=0; | |
628 | sum=0; | |
629 | for(Int_t k=starts[ybin1]; k<starts[ybin1] + kappawindow; k++) | |
630 | { | |
631 | kappa1 += hist->GetBinCenterX(k)*hist->GetBinContent(hist->GetBin(k,ybin1)); | |
632 | sum += (Float_t)hist->GetBinContent(hist->GetBin(k,ybin1)); | |
633 | } | |
634 | kappa1 /= sum; | |
635 | sum=0; | |
636 | for(Int_t k=starts[ybin2]; k<starts[ybin2] + kappawindow; k++) | |
637 | { | |
638 | kappa2 += hist->GetBinCenterX(k)*hist->GetBinContent(hist->GetBin(k,ybin2)); | |
639 | sum += (Float_t)hist->GetBinContent(hist->GetBin(k,ybin2)); | |
640 | } | |
641 | kappa2 /= sum; | |
642 | ||
643 | fXPeaks[fNPeaks] = ( kappa1*( hist->GetBinCenterY(ybin2) - fYPeaks[fNPeaks] ) + | |
644 | kappa2*( fYPeaks[fNPeaks] - hist->GetBinCenterY(ybin1) ) ) / | |
645 | (hist->GetBinCenterY(ybin2) - hist->GetBinCenterY(ybin1)); | |
646 | ||
647 | fNPeaks++; | |
648 | */ | |
649 | } | |
650 | } | |
651 | break; | |
652 | } | |
653 | else | |
654 | yl--;//Search continues... | |
655 | } | |
4cafa5fc | 656 | } |
657 | } | |
4cafa5fc | 658 | |
b2a02bce | 659 | for(Int_t i=0; i<hist->GetNbinsY(); i++) |
917e711b | 660 | delete localmaxima[i]; |
52a2a604 | 661 | |
917e711b | 662 | delete [] localmaxima; |
b2a02bce | 663 | delete [] nmaxs; |
664 | delete [] starts; | |
665 | delete [] maxs; | |
4fc9a6a4 | 666 | } |
667 | ||
917e711b | 668 | struct AliL3PreYPeak |
0bd0c1ef | 669 | { |
917e711b | 670 | Int_t fStartPosition; // Start position in X |
671 | Int_t fEndPosition; // End position in X | |
672 | Int_t fMinValue; // Minimum value inside the prepeak | |
673 | Int_t fMaxValue; // Maximum value inside the prepeak | |
674 | Int_t fPrevValue; // Neighbour values | |
675 | Int_t fLeftValue; // Neighbour values | |
676 | Int_t fRightValue; // Neighbour values | |
0bd0c1ef | 677 | }; |
678 | ||
de3c3890 | 679 | void AliL3HoughMaxFinder::FindAdaptedRowPeaks(Int_t kappawindow,Int_t xsize,Int_t ysize) |
0bd0c1ef | 680 | { |
917e711b | 681 | // Peak finder which is working over the Hough Space provided by the AliL3HoughTransformerRow class |
0bd0c1ef | 682 | AliL3Histogram *hist = fCurrentHisto; |
683 | ||
684 | if(!hist) | |
685 | { | |
686 | cerr<<"AliL3HoughMaxFinder : No histogram!"<<endl; | |
687 | return; | |
688 | } | |
689 | ||
690 | if(hist->GetNEntries() == 0) | |
691 | return; | |
692 | ||
693 | Int_t xmin = hist->GetFirstXbin(); | |
694 | Int_t xmax = hist->GetLastXbin(); | |
695 | Int_t ymin = hist->GetFirstYbin(); | |
696 | Int_t ymax = hist->GetLastYbin(); | |
a8ffd46b | 697 | Int_t nxbins = hist->GetNbinsX()+2; |
698 | Int_t *content = hist->GetContentArray(); | |
699 | ||
0bd0c1ef | 700 | //Start by looking for pre-peaks: |
701 | ||
917e711b | 702 | AliL3PreYPeak **localmaxima = new AliL3PreYPeak*[hist->GetNbinsY()]; |
0bd0c1ef | 703 | |
704 | Short_t *nmaxs = new Short_t[hist->GetNbinsY()]; | |
1e75562a | 705 | memset(nmaxs,0,hist->GetNbinsY()*sizeof(Short_t)); |
917e711b | 706 | Int_t lastvalue=0,value=0; |
5929c18d | 707 | for(Int_t ybin=fNextRow[ymin]; ybin<=ymax; ybin = fNextRow[ybin+1]) |
0bd0c1ef | 708 | { |
a8ffd46b | 709 | localmaxima[ybin-ymin] = new AliL3PreYPeak[nxbins-2]; |
917e711b | 710 | lastvalue = 0; |
0bd0c1ef | 711 | Bool_t found = 0; |
712 | for(Int_t xbin=xmin; xbin<=xmax; xbin++) | |
713 | { | |
a8ffd46b | 714 | value = content[xbin + nxbins*ybin]; //value = hist->GetBinContent(xbin + nxbins*ybin); //value = hist->GetBinContent(hist->GetBin(xbin,ybin)); |
de3c3890 | 715 | if(value > 0) |
0bd0c1ef | 716 | { |
cdb3bb24 | 717 | if((value - lastvalue) > 1) |
0bd0c1ef | 718 | { |
917e711b | 719 | localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fStartPosition = xbin; |
720 | localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fEndPosition = xbin; | |
721 | localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMinValue = value; | |
722 | localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMaxValue = value; | |
723 | localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fPrevValue = 0; | |
724 | localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fLeftValue = lastvalue; | |
0bd0c1ef | 725 | found = 1; |
726 | } | |
917e711b | 727 | if(abs(value - lastvalue) <= 1) |
0bd0c1ef | 728 | { |
cdb3bb24 | 729 | if(found) { |
917e711b | 730 | localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fEndPosition = xbin; |
731 | if(value>localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMaxValue) | |
732 | localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMaxValue = value; | |
733 | if(value<localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMinValue) | |
734 | localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fMinValue = value; | |
cdb3bb24 | 735 | } |
0bd0c1ef | 736 | } |
cdb3bb24 | 737 | if((value - lastvalue) < -1) |
738 | { | |
739 | if(found) { | |
740 | localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fRightValue = value; | |
741 | nmaxs[ybin-ymin]++; | |
742 | found = 0; | |
743 | } | |
744 | } | |
745 | } | |
746 | else | |
747 | { | |
748 | if(found) { | |
749 | localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fRightValue = value; | |
750 | nmaxs[ybin-ymin]++; | |
751 | found = 0; | |
752 | } | |
0bd0c1ef | 753 | } |
917e711b | 754 | lastvalue = value; |
0bd0c1ef | 755 | |
756 | } | |
de3c3890 | 757 | if(found) { |
cdb3bb24 | 758 | localmaxima[ybin-ymin][nmaxs[ybin-ymin]].fRightValue = 0; |
de3c3890 | 759 | nmaxs[ybin-ymin]++; |
760 | } | |
0bd0c1ef | 761 | } |
762 | ||
917e711b | 763 | AliL3Pre2DPeak maxima[500]; |
de3c3890 | 764 | Int_t nmaxima = 0; |
765 | ||
0bd0c1ef | 766 | for(Int_t ybin=ymax; ybin >= ymin; ybin--) |
767 | { | |
768 | for(Int_t i=0; i<nmaxs[ybin-ymin]; i++) | |
769 | { | |
917e711b | 770 | Int_t localminvalue = localmaxima[ybin-ymin][i].fMinValue; |
771 | Int_t localmaxvalue = localmaxima[ybin-ymin][i].fMaxValue; | |
772 | Int_t localprevvalue = localmaxima[ybin-ymin][i].fPrevValue; | |
773 | Int_t localnextvalue = 0; | |
774 | Int_t localleftvalue = localmaxima[ybin-ymin][i].fLeftValue; | |
775 | Int_t localrightvalue = localmaxima[ybin-ymin][i].fRightValue; | |
776 | ||
777 | if(localminvalue<0) | |
0bd0c1ef | 778 | continue; //already used |
779 | ||
780 | //start expanding in the psi-direction: | |
781 | ||
917e711b | 782 | Int_t localxstart = localmaxima[ybin-ymin][i].fStartPosition; |
783 | Int_t localxend = localmaxima[ybin-ymin][i].fEndPosition; | |
784 | Int_t tempxstart = localmaxima[ybin-ymin][i].fStartPosition; | |
785 | Int_t tempxend = localmaxima[ybin-ymin][i].fEndPosition; | |
0bd0c1ef | 786 | |
917e711b | 787 | Int_t localy=ybin-1,nybins=1; |
de3c3890 | 788 | |
917e711b | 789 | while(localy >= ymin) |
0bd0c1ef | 790 | { |
791 | Bool_t found=0; | |
917e711b | 792 | for(Int_t j=0; j<nmaxs[localy-ymin]; j++) |
0bd0c1ef | 793 | { |
917e711b | 794 | if( (localmaxima[localy-ymin][j].fStartPosition <= (tempxend + kappawindow)) && (localmaxima[localy-ymin][j].fEndPosition >= (tempxstart - kappawindow))) |
0bd0c1ef | 795 | { |
cdb3bb24 | 796 | if((localmaxima[localy-ymin][j].fMinValue <= localmaxvalue) && (localmaxima[localy-ymin][j].fMaxValue >= localminvalue)) |
0bd0c1ef | 797 | { |
917e711b | 798 | if(localmaxima[localy-ymin][j].fEndPosition > localxend) |
799 | localxend = localmaxima[localy-ymin][j].fEndPosition; | |
800 | if(localmaxima[localy-ymin][j].fStartPosition < localxstart) | |
801 | localxstart = localmaxima[localy-ymin][j].fStartPosition; | |
802 | tempxstart = localmaxima[localy-ymin][j].fStartPosition; | |
803 | tempxend = localmaxima[localy-ymin][j].fEndPosition; | |
804 | if(localmaxima[localy-ymin][j].fMinValue < localminvalue) | |
805 | localminvalue = localmaxima[localy-ymin][j].fMinValue; | |
806 | if(localmaxima[localy-ymin][j].fMaxValue > localmaxvalue) | |
807 | localmaxvalue = localmaxima[localy-ymin][j].fMaxValue; | |
808 | if(localmaxima[localy-ymin][j].fRightValue > localrightvalue) | |
809 | localrightvalue = localmaxima[localy-ymin][j].fRightValue; | |
810 | if(localmaxima[localy-ymin][j].fLeftValue > localleftvalue) | |
811 | localleftvalue = localmaxima[localy-ymin][j].fLeftValue; | |
812 | localmaxima[localy-ymin][j].fMinValue = -1; | |
0bd0c1ef | 813 | found = 1; |
814 | nybins++; | |
815 | break; | |
816 | } | |
817 | else | |
818 | { | |
917e711b | 819 | if(localmaxvalue > localmaxima[localy-ymin][j].fPrevValue) |
820 | localmaxima[localy-ymin][j].fPrevValue = localmaxvalue; | |
821 | if(localmaxima[localy-ymin][j].fMaxValue > localnextvalue) | |
822 | localnextvalue = localmaxima[localy-ymin][j].fMaxValue; | |
0bd0c1ef | 823 | } |
824 | } | |
825 | } | |
917e711b | 826 | if(!found || localy == ymin)//no more local maximas to be matched, so write the final peak and break the expansion: |
0bd0c1ef | 827 | { |
917e711b | 828 | if((nybins > ysize) && ((localxend-localxstart+1) > xsize) && (localmaxvalue > localprevvalue) && (localmaxvalue > localnextvalue) && (localmaxvalue > localleftvalue) && (localmaxvalue > localrightvalue)) |
829 | // if((nybins > ysize) && ((localxend-localxstart+1) > xsize)) | |
0bd0c1ef | 830 | { |
917e711b | 831 | maxima[nmaxima].fX = ((Float_t)localxstart+(Float_t)localxend)/2.0; |
832 | maxima[nmaxima].fY = ((Float_t)ybin+(Float_t)(localy+1))/2.0; | |
833 | maxima[nmaxima].fSizeX = localxend-localxstart+1; | |
834 | maxima[nmaxima].fSizeY = nybins; | |
835 | maxima[nmaxima].fWeight = (localminvalue+localmaxvalue)/2; | |
836 | maxima[nmaxima].fStartX = localxstart; | |
837 | maxima[nmaxima].fEndX = localxend; | |
838 | maxima[nmaxima].fStartY = localy +1; | |
839 | maxima[nmaxima].fEndY = ybin; | |
0bd0c1ef | 840 | #ifdef do_mc |
917e711b | 841 | // cout<<"Peak found at: "<<((Float_t)localxstart+(Float_t)localxend)/2.0<<" "<<((Float_t)ybin+(Float_t)(localy+1))/2.0<<" "<<localminvalue<<" "<<localmaxvalue<<" "<<" with weight "<<(localminvalue+localmaxvalue)/2<<" and size "<<localxend-localxstart+1<<" by "<<nybins<<endl; |
0bd0c1ef | 842 | #endif |
de3c3890 | 843 | nmaxima++; |
0bd0c1ef | 844 | } |
845 | break; | |
846 | } | |
847 | else | |
917e711b | 848 | localy--;//Search continues... |
0bd0c1ef | 849 | } |
850 | } | |
851 | } | |
852 | ||
de3c3890 | 853 | //remove fake tracks |
de3c3890 | 854 | |
855 | for(Int_t i = 0; i < (nmaxima - 1); i++) | |
856 | { | |
1e75562a | 857 | // if(maxima[i].fWeight < 0) continue; |
de3c3890 | 858 | for(Int_t j = i + 1; j < nmaxima; j++) |
859 | { | |
1e75562a | 860 | // if(maxima[j].fWeight < 0) continue; |
861 | MergeRowPeaks(&maxima[i],&maxima[j],5.0); | |
de3c3890 | 862 | } |
863 | } | |
864 | ||
865 | //merge tracks in neighbour eta slices | |
de3c3890 | 866 | Int_t currentnpeaks = fNPeaks; |
867 | for(Int_t i = 0; i < nmaxima; i++) { | |
917e711b | 868 | if(maxima[i].fWeight < 0) continue; |
de3c3890 | 869 | Bool_t merged = kFALSE; |
870 | for(Int_t j = fN1PeaksPrevEtaSlice; j < fN2PeaksPrevEtaSlice; j++) { | |
1e75562a | 871 | // if(fWeight[j] < 0) continue; |
872 | // Merge only peaks with limited size in eta | |
873 | if((fENDETAPeaks[j]-fSTARTETAPeaks[j]) >= 2) continue; | |
874 | if((maxima[i].fStartX <= fENDXPeaks[j]+1) && (maxima[i].fEndX >= fSTARTXPeaks[j]-1) && (maxima[i].fStartY <= fENDYPeaks[j]+1) && (maxima[i].fEndY >= fSTARTYPeaks[j]-1)){ | |
875 | //merge | |
876 | merged = kTRUE; | |
877 | if(fWeight[j] > 0) { | |
917e711b | 878 | fXPeaks[fNPeaks] = (hist->GetPreciseBinCenterX(maxima[i].fX)+(fENDETAPeaks[j]-fSTARTETAPeaks[j]+1)*fXPeaks[j])/(fENDETAPeaks[j]-fSTARTETAPeaks[j]+2); |
879 | fYPeaks[fNPeaks] = (hist->GetPreciseBinCenterY(maxima[i].fY)+(fENDETAPeaks[j]-fSTARTETAPeaks[j]+1)*fYPeaks[j])/(fENDETAPeaks[j]-fSTARTETAPeaks[j]+2); | |
917e711b | 880 | fSTARTXPeaks[fNPeaks] = maxima[i].fStartX; |
881 | fSTARTYPeaks[fNPeaks] = maxima[i].fStartY; | |
882 | fENDXPeaks[fNPeaks] = maxima[i].fEndX; | |
883 | fENDYPeaks[fNPeaks] = maxima[i].fEndY; | |
1e75562a | 884 | |
885 | fWeight[fNPeaks] = abs((Int_t)maxima[i].fWeight) + abs(fWeight[j]); | |
de3c3890 | 886 | fSTARTETAPeaks[fNPeaks] = fSTARTETAPeaks[j]; |
887 | fENDETAPeaks[fNPeaks] = fCurrentEtaSlice; | |
888 | fNPeaks++; | |
de3c3890 | 889 | } |
1e75562a | 890 | fWeight[j] = -abs(fWeight[j]); |
de3c3890 | 891 | } |
892 | } | |
917e711b | 893 | fXPeaks[fNPeaks] = hist->GetPreciseBinCenterX(maxima[i].fX); |
894 | fYPeaks[fNPeaks] = hist->GetPreciseBinCenterY(maxima[i].fY); | |
de3c3890 | 895 | if(!merged) |
1e75562a | 896 | fWeight[fNPeaks] = abs((Int_t)maxima[i].fWeight); |
de3c3890 | 897 | else |
1e75562a | 898 | fWeight[fNPeaks] = -abs((Int_t)maxima[i].fWeight); |
917e711b | 899 | fSTARTXPeaks[fNPeaks] = maxima[i].fStartX; |
900 | fSTARTYPeaks[fNPeaks] = maxima[i].fStartY; | |
901 | fENDXPeaks[fNPeaks] = maxima[i].fEndX; | |
902 | fENDYPeaks[fNPeaks] = maxima[i].fEndY; | |
de3c3890 | 903 | fSTARTETAPeaks[fNPeaks] = fCurrentEtaSlice; |
904 | fENDETAPeaks[fNPeaks] = fCurrentEtaSlice; | |
905 | fNPeaks++; | |
1e75562a | 906 | } |
907 | ||
de3c3890 | 908 | fN1PeaksPrevEtaSlice = currentnpeaks; |
909 | fN2PeaksPrevEtaSlice = fNPeaks; | |
910 | ||
5929c18d | 911 | for(Int_t ybin=fNextRow[ymin]; ybin<=ymax; ybin = fNextRow[ybin+1]) |
1e75562a | 912 | delete [] localmaxima[ybin-ymin]; |
0bd0c1ef | 913 | |
917e711b | 914 | delete [] localmaxima; |
0bd0c1ef | 915 | delete [] nmaxs; |
0bd0c1ef | 916 | } |
b2a02bce | 917 | |
917e711b | 918 | void AliL3HoughMaxFinder::FindPeak1(Int_t ywindow,Int_t xbinsides) |
4fc9a6a4 | 919 | { |
b1886074 | 920 | //Testing mutliple peakfinding. |
921 | //The algorithm searches the histogram for prepreaks by looking in windows | |
917e711b | 922 | //for each bin on the xaxis. The size of these windows is controlled by ywindow. |
b1886074 | 923 | //Then the prepreaks are sorted according to their weight (sum inside window), |
924 | //and the peak positions are calculated by taking the weighted mean in both | |
917e711b | 925 | //x and y direction. The size of the peak in x-direction is controlled by xbinsides. |
b1886074 | 926 | |
4fc9a6a4 | 927 | if(!fCurrentHisto) |
928 | { | |
929 | printf("AliL3HoughMaxFinder::FindPeak1 : No input histogram\n"); | |
930 | return; | |
931 | } | |
b2a02bce | 932 | if(fCurrentHisto->GetNEntries()==0) |
933 | return; | |
934 | ||
917e711b | 935 | //Int_t ywindow=2; |
936 | //Int_t xbinsides=1; | |
b1886074 | 937 | |
afd8fed4 | 938 | //Float_t max_kappa = 0.001; |
3e87ef69 | 939 | //Float_t max_phi0 = 0.08; |
b1886074 | 940 | |
917e711b | 941 | Int_t maxsum=0; |
4fc9a6a4 | 942 | |
943 | Int_t xmin = fCurrentHisto->GetFirstXbin(); | |
944 | Int_t xmax = fCurrentHisto->GetLastXbin(); | |
945 | Int_t ymin = fCurrentHisto->GetFirstYbin(); | |
946 | Int_t ymax = fCurrentHisto->GetLastYbin(); | |
947 | Int_t nbinsx = fCurrentHisto->GetNbinsX()+1; | |
948 | ||
917e711b | 949 | AliL3AxisWindow **windowPt = new AliL3AxisWindow*[nbinsx]; |
950 | AliL3AxisWindow **anotherPt = new AliL3AxisWindow*[nbinsx]; | |
4fc9a6a4 | 951 | |
952 | for(Int_t i=0; i<nbinsx; i++) | |
953 | { | |
917e711b | 954 | windowPt[i] = new AliL3AxisWindow; |
eb86303b | 955 | #if defined(__DECCXX) |
917e711b | 956 | bzero((char *)windowPt[i],sizeof(AliL3AxisWindow)); |
eb86303b | 957 | #else |
917e711b | 958 | bzero((void*)windowPt[i],sizeof(AliL3AxisWindow)); |
eb86303b | 959 | #endif |
4fc9a6a4 | 960 | anotherPt[i] = windowPt[i]; |
961 | } | |
962 | ||
963 | for(Int_t xbin=xmin; xbin<=xmax; xbin++) | |
964 | { | |
917e711b | 965 | maxsum = 0; |
966 | for(Int_t ybin=ymin; ybin<=ymax-ywindow; ybin++) | |
4fc9a6a4 | 967 | { |
917e711b | 968 | Int_t suminwindow=0; |
969 | for(Int_t b=ybin; b<ybin+ywindow; b++) | |
4fc9a6a4 | 970 | { |
971 | //inside window | |
972 | Int_t bin = fCurrentHisto->GetBin(xbin,b); | |
917e711b | 973 | suminwindow += (Int_t)fCurrentHisto->GetBinContent(bin); |
4fc9a6a4 | 974 | } |
975 | ||
917e711b | 976 | if(suminwindow > maxsum) |
4fc9a6a4 | 977 | { |
917e711b | 978 | maxsum = suminwindow; |
979 | windowPt[xbin]->fYmin = ybin; | |
980 | windowPt[xbin]->fYmax = ybin + ywindow; | |
981 | windowPt[xbin]->fWeight = suminwindow; | |
982 | windowPt[xbin]->fXbin = xbin; | |
4fc9a6a4 | 983 | } |
984 | } | |
985 | } | |
986 | ||
987 | //Sort the windows according to the weight | |
988 | SortPeaks(windowPt,0,nbinsx); | |
989 | ||
4fc9a6a4 | 990 | Float_t top,butt; |
991 | for(Int_t i=0; i<nbinsx; i++) | |
992 | { | |
993 | top=butt=0; | |
917e711b | 994 | Int_t xbin = windowPt[i]->fXbin; |
4fc9a6a4 | 995 | |
996 | if(xbin<xmin || xbin > xmax-1) continue; | |
997 | ||
b1886074 | 998 | //Check if this is really a local maxima |
917e711b | 999 | if(anotherPt[xbin-1]->fWeight > anotherPt[xbin]->fWeight || |
1000 | anotherPt[xbin+1]->fWeight > anotherPt[xbin]->fWeight) | |
4fc9a6a4 | 1001 | continue; |
1002 | ||
917e711b | 1003 | for(Int_t j=windowPt[i]->fYmin; j<windowPt[i]->fYmax; j++) |
4fc9a6a4 | 1004 | { |
1005 | //Calculate the mean in y direction: | |
917e711b | 1006 | Int_t bin = fCurrentHisto->GetBin(windowPt[i]->fXbin,j); |
4fc9a6a4 | 1007 | top += (fCurrentHisto->GetBinCenterY(j))*(fCurrentHisto->GetBinContent(bin)); |
1008 | butt += fCurrentHisto->GetBinContent(bin); | |
1009 | } | |
afd8fed4 | 1010 | |
b2a02bce | 1011 | if(butt < fThreshold) |
1012 | continue; | |
afd8fed4 | 1013 | |
917e711b | 1014 | fXPeaks[fNPeaks] = fCurrentHisto->GetBinCenterX(windowPt[i]->fXbin); |
afd8fed4 | 1015 | fYPeaks[fNPeaks] = top/butt; |
1016 | fWeight[fNPeaks] = (Int_t)butt; | |
4ab9f8f0 | 1017 | //cout<<"mean in y "<<ypeaks[n]<<" on x "<<windowPt[i]->xbin<<" content "<<butt<<endl; |
afd8fed4 | 1018 | fNPeaks++; |
3e87ef69 | 1019 | if(fNPeaks==fNMax) |
1020 | { | |
1021 | cerr<<"AliL3HoughMaxFinder::FindPeak1 : Peak array out of range!!!"<<endl; | |
1022 | break; | |
1023 | } | |
4fc9a6a4 | 1024 | } |
1025 | ||
1026 | ||
1027 | //Improve the peaks by including the region around in x. | |
1028 | Float_t ytop,ybutt; | |
b1886074 | 1029 | Int_t prev; |
1030 | Int_t w; | |
afd8fed4 | 1031 | for(Int_t i=0; i<fNPeaks; i++) |
4fc9a6a4 | 1032 | { |
afd8fed4 | 1033 | Int_t xbin = fCurrentHisto->FindXbin(fXPeaks[i]); |
917e711b | 1034 | if(xbin - xbinsides < xmin || xbin + xbinsides > xmax) continue; |
4fc9a6a4 | 1035 | top=butt=0; |
1036 | ytop=0,ybutt=0; | |
b1886074 | 1037 | w=0; |
917e711b | 1038 | prev = xbin - xbinsides+1; |
1039 | for(Int_t j=xbin-xbinsides; j<=xbin+xbinsides; j++) | |
4fc9a6a4 | 1040 | { |
4ab9f8f0 | 1041 | /* |
b1886074 | 1042 | //Check if the windows are overlapping: |
1043 | if(anotherPt[j]->ymin > anotherPt[prev]->ymax) {prev=j; continue;} | |
1044 | if(anotherPt[j]->ymax < anotherPt[prev]->ymin) {prev=j; continue;} | |
1045 | prev = j; | |
4ab9f8f0 | 1046 | */ |
1047 | ||
917e711b | 1048 | top += fCurrentHisto->GetBinCenterX(j)*anotherPt[j]->fWeight; |
1049 | butt += anotherPt[j]->fWeight; | |
b1886074 | 1050 | |
917e711b | 1051 | for(Int_t k=anotherPt[j]->fYmin; k<anotherPt[j]->fYmax; k++) |
4fc9a6a4 | 1052 | { |
1053 | Int_t bin = fCurrentHisto->GetBin(j,k); | |
1054 | ytop += (fCurrentHisto->GetBinCenterY(k))*(fCurrentHisto->GetBinContent(bin)); | |
1055 | ybutt += fCurrentHisto->GetBinContent(bin); | |
95a00d93 | 1056 | w+=(Int_t)fCurrentHisto->GetBinContent(bin); |
4fc9a6a4 | 1057 | } |
1058 | } | |
b1886074 | 1059 | |
afd8fed4 | 1060 | fXPeaks[i] = top/butt; |
1061 | fYPeaks[i] = ytop/ybutt; | |
1062 | fWeight[i] = w; | |
3e87ef69 | 1063 | //cout<<"Setting weight "<<w<<" kappa "<<fXPeaks[i]<<" phi0 "<<fYPeaks[i]<<endl; |
b5a207b4 | 1064 | |
afd8fed4 | 1065 | /* |
b1886074 | 1066 | //Check if this peak is overlapping with a previous: |
3e87ef69 | 1067 | for(Int_t p=0; p<i; p++) |
b1886074 | 1068 | { |
3e87ef69 | 1069 | //cout<<fabs(fXPeaks[p] - fXPeaks[i])<<" "<<fabs(fYPeaks[p] - fYPeaks[i])<<endl; |
1070 | if(fabs(fXPeaks[p] - fXPeaks[i]) < max_kappa && | |
afd8fed4 | 1071 | fabs(fYPeaks[p] - fYPeaks[i]) < max_phi0) |
b1886074 | 1072 | { |
afd8fed4 | 1073 | fWeight[i]=0; |
3e87ef69 | 1074 | //break; |
b1886074 | 1075 | } |
1076 | } | |
afd8fed4 | 1077 | */ |
4fc9a6a4 | 1078 | } |
1079 | ||
1080 | for(Int_t i=0; i<nbinsx; i++) | |
1081 | delete windowPt[i]; | |
1082 | delete [] windowPt; | |
1083 | delete [] anotherPt; | |
4fc9a6a4 | 1084 | } |
1085 | ||
917e711b | 1086 | void AliL3HoughMaxFinder::SortPeaks(struct AliL3AxisWindow **a,Int_t first,Int_t last) |
4fc9a6a4 | 1087 | { |
1088 | //General sorting routine | |
1089 | //Sort according to PeakCompare() | |
1090 | ||
917e711b | 1091 | static struct AliL3AxisWindow *tmp; |
4fc9a6a4 | 1092 | static int i; // "static" to save stack space |
1093 | int j; | |
1094 | ||
1095 | while (last - first > 1) { | |
1096 | i = first; | |
1097 | j = last; | |
1098 | for (;;) { | |
1099 | while (++i < last && PeakCompare(a[i], a[first]) < 0) | |
1100 | ; | |
1101 | while (--j > first && PeakCompare(a[j], a[first]) > 0) | |
1102 | ; | |
1103 | if (i >= j) | |
1104 | break; | |
1105 | ||
1106 | tmp = a[i]; | |
1107 | a[i] = a[j]; | |
1108 | a[j] = tmp; | |
1109 | } | |
1110 | if (j == first) { | |
1111 | ++first; | |
1112 | continue; | |
1113 | } | |
1114 | tmp = a[first]; | |
1115 | a[first] = a[j]; | |
1116 | a[j] = tmp; | |
1117 | if (j - first < last - (j + 1)) { | |
1118 | SortPeaks(a, first, j); | |
1119 | first = j + 1; // QSort(j + 1, last); | |
1120 | } else { | |
1121 | SortPeaks(a, j + 1, last); | |
1122 | last = j; // QSort(first, j); | |
1123 | } | |
1124 | } | |
1125 | ||
1126 | } | |
1127 | ||
917e711b | 1128 | Int_t AliL3HoughMaxFinder::PeakCompare(struct AliL3AxisWindow *a,struct AliL3AxisWindow *b) const |
4fc9a6a4 | 1129 | { |
917e711b | 1130 | // Peak comparison based on peaks weight |
1131 | if(a->fWeight < b->fWeight) return 1; | |
1132 | if(a->fWeight > b->fWeight) return -1; | |
4fc9a6a4 | 1133 | return 0; |
4fc9a6a4 | 1134 | } |
52a2a604 | 1135 | |
7f66cda9 | 1136 | void AliL3HoughMaxFinder::FindPeak(Int_t t1,Double_t t2,Int_t t3) |
e4c21048 | 1137 | { |
1138 | //Attempt of a more sophisticated peak finder. | |
4fc9a6a4 | 1139 | //Finds the best peak in the histogram, and returns the corresponding |
1140 | //track object. | |
1141 | ||
4cafa5fc | 1142 | if(!fCurrentHisto) |
1143 | { | |
1144 | printf("AliL3HoughMaxFinder::FindPeak : No histogram!!\n"); | |
b1886074 | 1145 | return; |
4cafa5fc | 1146 | } |
1147 | AliL3Histogram *hist = fCurrentHisto; | |
b2a02bce | 1148 | if(hist->GetNEntries()==0) |
1149 | return; | |
81516384 | 1150 | |
4cafa5fc | 1151 | Int_t xmin = hist->GetFirstXbin(); |
1152 | Int_t xmax = hist->GetLastXbin(); | |
1153 | Int_t ymin = hist->GetFirstYbin(); | |
1154 | Int_t ymax = hist->GetLastYbin(); | |
1155 | Int_t nbinsx = hist->GetNbinsX()+1; | |
e4c21048 | 1156 | |
1157 | Int_t *m = new Int_t[nbinsx]; | |
917e711b | 1158 | Int_t *mlow = new Int_t[nbinsx]; |
1159 | Int_t *mup = new Int_t[nbinsx]; | |
81516384 | 1160 | |
1161 | ||
1162 | recompute: //this is a goto. | |
1163 | ||
e4c21048 | 1164 | for(Int_t i=0; i<nbinsx; i++) |
1165 | { | |
1166 | m[i]=0; | |
917e711b | 1167 | mlow[i]=0; |
1168 | mup[i]=0; | |
e4c21048 | 1169 | } |
e4c21048 | 1170 | |
917e711b | 1171 | Int_t maxx=0,sum=0,maxxbin=0,bin=0; |
e4c21048 | 1172 | |
1173 | for(Int_t xbin=xmin; xbin<=xmax; xbin++) | |
1174 | { | |
4cafa5fc | 1175 | for(Int_t ybin=ymin; ybin <= ymax - t1; ybin++) |
e4c21048 | 1176 | { |
1177 | sum = 0; | |
4cafa5fc | 1178 | for(Int_t y=ybin; y <= ybin+t1; y++) |
e4c21048 | 1179 | { |
4cafa5fc | 1180 | if(y>ymax) break; |
e4c21048 | 1181 | //Inside window |
1182 | bin = hist->GetBin(xbin,y); | |
1183 | sum += (Int_t)hist->GetBinContent(bin); | |
1184 | ||
1185 | } | |
1186 | if(sum > m[xbin]) //Max value locally in this xbin | |
1187 | { | |
1188 | m[xbin]=sum; | |
917e711b | 1189 | mlow[xbin]=ybin; |
1190 | mup[xbin]=ybin + t1; | |
e4c21048 | 1191 | } |
1192 | ||
1193 | } | |
1194 | ||
917e711b | 1195 | if(m[xbin] > maxx) //Max value globally in x-direction |
e4c21048 | 1196 | { |
917e711b | 1197 | maxxbin = xbin; |
1198 | maxx = m[xbin];//sum; | |
e4c21048 | 1199 | } |
1200 | } | |
917e711b | 1201 | //printf("maxxbin %d maxx %d mlow %d mup %d\n",maxxbin,maxx,mlow[maxxbin],mup[maxxbin]); |
1202 | //printf("ylow %f yup %f\n",hist->GetBinCenterY(mlow[maxxbin]),hist->GetBinCenterY(mup[maxxbin])); | |
e4c21048 | 1203 | |
1204 | //Determine a width in the x-direction | |
917e711b | 1205 | Int_t xlow=0,xup=0; |
4cafa5fc | 1206 | |
917e711b | 1207 | for(Int_t xbin=maxxbin-1; xbin >= xmin; xbin--) |
e4c21048 | 1208 | { |
917e711b | 1209 | if(m[xbin] < maxx*t2) |
e4c21048 | 1210 | { |
917e711b | 1211 | xlow = xbin+1; |
e4c21048 | 1212 | break; |
1213 | } | |
1214 | } | |
917e711b | 1215 | for(Int_t xbin = maxxbin+1; xbin <=xmax; xbin++) |
e4c21048 | 1216 | { |
917e711b | 1217 | if(m[xbin] < maxx*t2) |
e4c21048 | 1218 | { |
917e711b | 1219 | xup = xbin-1; |
e4c21048 | 1220 | break; |
1221 | } | |
1222 | } | |
1223 | ||
917e711b | 1224 | Double_t top=0,butt=0,value,xpeak; |
1225 | if(xup - xlow + 1 > t3) | |
81516384 | 1226 | { |
1227 | t1 -= 1; | |
917e711b | 1228 | printf("\nxrange out if limit xup %d xlow %d t1 %d\n\n",xlow,xup,t1); |
81516384 | 1229 | if(t1 > 1) |
1230 | goto recompute; | |
1231 | else | |
1232 | { | |
917e711b | 1233 | xpeak = hist->GetBinCenterX(maxxbin); |
81516384 | 1234 | goto moveon; |
1235 | } | |
1236 | } | |
1237 | ||
917e711b | 1238 | //printf("xlow %f xup %f\n",hist->GetBinCenterX(xlow),hist->GetBinCenterX(xup)); |
1239 | //printf("Spread in x %d\n",xup-xlow +1); | |
e4c21048 | 1240 | |
1241 | //Now, calculate the center of mass in x-direction | |
917e711b | 1242 | for(Int_t xbin=xlow; xbin <= xup; xbin++) |
e4c21048 | 1243 | { |
4cafa5fc | 1244 | value = hist->GetBinCenterX(xbin); |
e4c21048 | 1245 | top += value*m[xbin]; |
1246 | butt += m[xbin]; | |
1247 | } | |
917e711b | 1248 | xpeak = top/butt; |
e4c21048 | 1249 | |
81516384 | 1250 | moveon: |
e4c21048 | 1251 | |
1252 | //Find the peak in y direction: | |
917e711b | 1253 | Int_t xl = hist->FindXbin(xpeak); |
1254 | if(hist->GetBinCenterX(xl) > xpeak) | |
1255 | xl--; | |
81516384 | 1256 | |
917e711b | 1257 | Int_t xu = xl + 1; |
81516384 | 1258 | |
917e711b | 1259 | if(hist->GetBinCenterX(xl) > xpeak || hist->GetBinCenterX(xu) <= xpeak) |
1260 | printf("\nAliL3HoughMaxFinder::FindPeak : Wrong xrange %f %f %f\n\n",hist->GetBinCenterX(xl),xpeak,hist->GetBinCenterX(xu)); | |
e4c21048 | 1261 | |
917e711b | 1262 | //printf("\nxlow %f xup %f\n",hist->GetBinCenterX(xl),hist->GetBinCenterX(xu)); |
e4c21048 | 1263 | |
1264 | value=top=butt=0; | |
1265 | ||
917e711b | 1266 | //printf("ylow %f yup %f\n",hist->GetBinCenterY(mlow[xl]),hist->GetBinCenterY(mup[xl])); |
1267 | //printf("ylow %f yup %f\n",hist->GetBinCenterY(mlow[xu]),hist->GetBinCenterY(mup[xu])); | |
e4c21048 | 1268 | |
917e711b | 1269 | for(Int_t ybin=mlow[xl]; ybin <= mup[xl]; ybin++) |
e4c21048 | 1270 | { |
4cafa5fc | 1271 | value = hist->GetBinCenterY(ybin); |
917e711b | 1272 | bin = hist->GetBin(xl,ybin); |
e4c21048 | 1273 | top += value*hist->GetBinContent(bin); |
1274 | butt += hist->GetBinContent(bin); | |
1275 | } | |
917e711b | 1276 | Double_t ypeaklow = top/butt; |
e4c21048 | 1277 | |
917e711b | 1278 | //printf("ypeaklow %f\n",ypeaklow); |
e4c21048 | 1279 | |
1280 | value=top=butt=0; | |
917e711b | 1281 | for(Int_t ybin=mlow[xu]; ybin <= mup[xu]; ybin++) |
e4c21048 | 1282 | { |
4cafa5fc | 1283 | value = hist->GetBinCenterY(ybin); |
917e711b | 1284 | bin = hist->GetBin(xu,ybin); |
e4c21048 | 1285 | top += value*hist->GetBinContent(bin); |
1286 | butt += hist->GetBinContent(bin); | |
1287 | } | |
917e711b | 1288 | Double_t ypeakup = top/butt; |
e4c21048 | 1289 | |
917e711b | 1290 | //printf("ypeakup %f\n",ypeakup); |
e4c21048 | 1291 | |
917e711b | 1292 | Double_t xvalueup = hist->GetBinCenterX(xu); |
1293 | Double_t xvaluelow = hist->GetBinCenterX(xl); | |
e4c21048 | 1294 | |
917e711b | 1295 | Double_t ypeak = (ypeaklow*(xvalueup - xpeak) + ypeakup*(xpeak - xvaluelow))/(xvalueup - xvaluelow); |
e4c21048 | 1296 | |
1297 | ||
81516384 | 1298 | //Find the weight: |
917e711b | 1299 | //bin = hist->FindBin(xpeak,ypeak); |
95a00d93 | 1300 | //Int_t weight = (Int_t)hist->GetBinContent(bin); |
81516384 | 1301 | |
b1886074 | 1302 | //AliL3HoughTrack *track = new AliL3HoughTrack(); |
917e711b | 1303 | //track->SetTrackParameters(xpeak,ypeak,weight); |
1304 | fXPeaks[fNPeaks]=xpeak; | |
1305 | fYPeaks[fNPeaks]=ypeak; | |
7f66cda9 | 1306 | fWeight[fNPeaks]=(Int_t)hist->GetBinContent(bin); |
1307 | fNPeaks++; | |
52a2a604 | 1308 | |
1309 | delete [] m; | |
917e711b | 1310 | delete [] mlow; |
1311 | delete [] mup; | |
52a2a604 | 1312 | |
b1886074 | 1313 | //return track; |
e4c21048 | 1314 | } |
81516384 | 1315 | |
917e711b | 1316 | Float_t AliL3HoughMaxFinder::GetXPeakSize(Int_t i) const |
de3c3890 | 1317 | { |
917e711b | 1318 | // Get X size of a peak |
de3c3890 | 1319 | if(i<0 || i>fNMax) |
1320 | { | |
1321 | STDCERR<<"AliL3HoughMaxFinder::GetXPeakSize : Invalid index "<<i<<STDENDL; | |
1322 | return 0; | |
1323 | } | |
917e711b | 1324 | Float_t binwidth = fCurrentHisto->GetBinWidthX(); |
1325 | return binwidth*(fENDXPeaks[i]-fSTARTXPeaks[i]+1); | |
de3c3890 | 1326 | } |
1327 | ||
917e711b | 1328 | Float_t AliL3HoughMaxFinder::GetYPeakSize(Int_t i) const |
de3c3890 | 1329 | { |
917e711b | 1330 | // Get Y size of a peak |
de3c3890 | 1331 | if(i<0 || i>fNMax) |
1332 | { | |
1333 | STDCERR<<"AliL3HoughMaxFinder::GetYPeak : Invalid index "<<i<<STDENDL; | |
1334 | return 0; | |
1335 | } | |
917e711b | 1336 | Float_t binwidth = fCurrentHisto->GetBinWidthY(); |
1337 | return binwidth*(fENDYPeaks[i]-fSTARTYPeaks[i]+1); | |
de3c3890 | 1338 | } |
1e75562a | 1339 | |
1340 | Bool_t AliL3HoughMaxFinder::MergeRowPeaks(AliL3Pre2DPeak *maxima1, AliL3Pre2DPeak *maxima2, Float_t distance) | |
1341 | { | |
1342 | // Check the distance between tracks corresponding to given Hough space peaks and if the | |
1343 | // distance is smaller than some threshold value marks the smaller peak as fake | |
1344 | AliL3Histogram *hist = fCurrentHisto; | |
1345 | Int_t nxbins = hist->GetNbinsX()+2; | |
1346 | ||
1347 | Int_t xtrack1=0,xtrack2=0,ytrack1=0,ytrack2=0; | |
1348 | Int_t deltax = 9999; | |
1349 | for(Int_t ix1 = maxima1->fStartX; ix1 <= maxima1->fEndX; ix1++) { | |
1350 | for(Int_t ix2 = maxima2->fStartX; ix2 <= maxima2->fEndX; ix2++) { | |
1351 | if(abs(ix1 - ix2) < deltax) { | |
1352 | deltax = abs(ix1 - ix2); | |
1353 | xtrack1 = ix1; | |
1354 | xtrack2 = ix2; | |
1355 | } | |
1356 | } | |
1357 | } | |
1358 | Int_t deltay = 9999; | |
1359 | for(Int_t iy1 = maxima1->fStartY; iy1 <= maxima1->fEndY; iy1++) { | |
1360 | for(Int_t iy2 = maxima2->fStartY; iy2 <= maxima2->fEndY; iy2++) { | |
1361 | if(abs(iy1 - iy2) < deltay) { | |
1362 | deltay = abs(iy1 - iy2); | |
1363 | ytrack1 = iy1; | |
1364 | ytrack2 = iy2; | |
1365 | } | |
1366 | } | |
1367 | } | |
1368 | Int_t firstrow1 = fTrackFirstRow[xtrack1 + nxbins*ytrack1]; | |
1369 | Int_t lastrow1 = fTrackLastRow[xtrack1 + nxbins*ytrack1]; | |
1370 | Int_t firstrow2 = fTrackFirstRow[xtrack1 + nxbins*ytrack1]; | |
1371 | Int_t lastrow2 = fTrackLastRow[xtrack1 + nxbins*ytrack1]; | |
1372 | Int_t firstrow,lastrow; | |
1373 | if(firstrow1 < firstrow2) | |
1374 | firstrow = firstrow2; | |
1375 | else | |
1376 | firstrow = firstrow1; | |
1377 | ||
1378 | if(lastrow1 > lastrow2) | |
1379 | lastrow = lastrow2; | |
1380 | else | |
1381 | lastrow = lastrow1; | |
1382 | ||
1383 | AliL3HoughTrack track1; | |
1384 | Float_t x1 = hist->GetPreciseBinCenterX(xtrack1); | |
1385 | Float_t y1 = hist->GetPreciseBinCenterY(ytrack1); | |
1386 | Float_t psi1 = atan((x1-y1)/(AliL3HoughTransformerRow::GetBeta1()-AliL3HoughTransformerRow::GetBeta2())); | |
1387 | Float_t kappa1 = 2.0*(x1*cos(psi1)-AliL3HoughTransformerRow::GetBeta1()*sin(psi1)); | |
1388 | track1.SetTrackParameters(kappa1,psi1,1); | |
1389 | Float_t firsthit1[3]; | |
1390 | if(!track1.GetCrossingPoint(firstrow,firsthit1)) return kFALSE; | |
1391 | Float_t lasthit1[3]; | |
1392 | if(!track1.GetCrossingPoint(lastrow,lasthit1)) return kFALSE; | |
1393 | ||
1394 | AliL3HoughTrack track2; | |
1395 | Float_t x2 = hist->GetPreciseBinCenterX(xtrack2); | |
1396 | Float_t y2 = hist->GetPreciseBinCenterY(ytrack2); | |
1397 | Float_t psi2 = atan((x2-y2)/(AliL3HoughTransformerRow::GetBeta1()-AliL3HoughTransformerRow::GetBeta2())); | |
1398 | Float_t kappa2 = 2.0*(x2*cos(psi2)-AliL3HoughTransformerRow::GetBeta1()*sin(psi2)); | |
1399 | track2.SetTrackParameters(kappa2,psi2,1); | |
1400 | Float_t firsthit2[3]; | |
1401 | if(!track2.GetCrossingPoint(firstrow,firsthit2)) return kFALSE; | |
1402 | Float_t lasthit2[3]; | |
1403 | if(!track2.GetCrossingPoint(lastrow,lasthit2)) return kFALSE; | |
1404 | ||
1405 | Float_t padpitchlow = AliL3Transform::GetPadPitchWidth(AliL3Transform::GetPatch(firstrow)); | |
1406 | Float_t padpitchup = AliL3Transform::GetPadPitchWidth(AliL3Transform::GetPatch(lastrow)); | |
1407 | // check the distance between tracks at the edges | |
1408 | // cout<<"Check "<<firsthit1[1]<<" "<<firsthit2[1]<<" "<<padpitchlow<<" "<<lasthit1[1]<<" "<<lasthit2[1]<<" "<<padpitchup<<" "<<xtrack1<<" "<<ytrack1<<" "<<xtrack2<<" "<<ytrack2<<endl; | |
1409 | if((fabs(firsthit1[1]-firsthit2[1])/padpitchlow + fabs(lasthit1[1]-lasthit2[1])/padpitchup) < distance) { | |
1410 | if(maxima1->fSizeX*maxima1->fSizeY > maxima2->fSizeX*maxima2->fSizeY) | |
1411 | maxima2->fWeight = -fabs(maxima2->fWeight); | |
1412 | if(maxima1->fSizeX*maxima1->fSizeY < maxima2->fSizeX*maxima2->fSizeY) | |
1413 | maxima1->fWeight = -fabs(maxima1->fWeight); | |
1414 | if(maxima1->fSizeX*maxima1->fSizeY == maxima2->fSizeX*maxima2->fSizeY) { | |
1415 | if(maxima1->fStartX > maxima2->fStartX) | |
1416 | maxima1->fStartX = maxima2->fStartX; | |
1417 | if(maxima1->fStartY > maxima2->fStartY) | |
1418 | maxima1->fStartY = maxima2->fStartY; | |
1419 | if(maxima1->fEndX < maxima2->fEndX) | |
1420 | maxima1->fEndX = maxima2->fEndX; | |
1421 | if(maxima1->fEndY < maxima2->fEndY) | |
1422 | maxima1->fEndY = maxima2->fEndY; | |
1423 | maxima1->fX = ((Float_t)maxima1->fStartX + (Float_t)maxima1->fEndX)/2.0; | |
1424 | maxima1->fY = ((Float_t)maxima1->fStartY + (Float_t)maxima1->fEndY)/2.0; | |
1425 | maxima1->fSizeX = (maxima1->fEndX - maxima1->fStartX + 1); | |
1426 | maxima1->fSizeY = (maxima1->fEndY - maxima1->fStartY + 1); | |
1427 | maxima2->fWeight = -fabs(maxima2->fWeight); | |
1428 | } | |
1429 | return kTRUE; | |
1430 | } | |
1431 | return kFALSE; | |
1432 | } |