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