]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HLT/hough/AliL3HoughMaxFinder.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Made Constantins last checkin compile ;-)
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughMaxFinder.cxx
CommitLineData
b5a207b4 1//$Id$
2
b1886074 3// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4//*-- Copyright &copy ASV
4de874d1 5
e06900d5 6#include "AliL3StandardIncludes.h"
7
2061190a 8#ifndef no_root
208b54c5 9#include <TNtuple.h>
10#include <TFile.h>
2061190a 11#endif
208b54c5 12
e06900d5 13#include "AliL3Logging.h"
14#include "AliL3HoughMaxFinder.h"
4cafa5fc 15#include "AliL3Histogram.h"
4de874d1 16#include "AliL3TrackArray.h"
17#include "AliL3HoughTrack.h"
e06900d5 18
19#if GCCVERSION == 3
20using namespace std;
21#endif
4de874d1 22
b1886074 23//_____________________________________________________________
24// AliL3HoughMaxFinder
25//
26// Maximum finder
27
4de874d1 28ClassImp(AliL3HoughMaxFinder)
29
30
31AliL3HoughMaxFinder::AliL3HoughMaxFinder()
32{
33 //Default constructor
34 fThreshold = 0;
4de874d1 35 fHistoType=0;
3fe49b5b 36 fXPeaks=0;
37 fYPeaks=0;
38 fNPeaks=0;
39 fNMax=0;
2061190a 40#ifndef no_root
208b54c5 41 fNtuppel = 0;
2061190a 42#endif
4de874d1 43}
44
45
3fe49b5b 46AliL3HoughMaxFinder::AliL3HoughMaxFinder(Char_t *histotype,Int_t nmax,AliL3Histogram *hist)
4de874d1 47{
48 //Constructor
49
50 //fTracks = new AliL3TrackArray("AliL3HoughTrack");
51 if(strcmp(histotype,"KappaPhi")==0) fHistoType='c';
52 if(strcmp(histotype,"DPsi")==0) fHistoType='l';
53
4cafa5fc 54 if(hist)
55 fCurrentHisto = hist;
3fe49b5b 56
57 fNMax=nmax;
58 fXPeaks = new Float_t[fNMax];
59 fYPeaks = new Float_t[fNMax];
60 fWeight = new Int_t[fNMax];
2061190a 61#ifndef no_root
208b54c5 62 fNtuppel = 0;
2061190a 63#endif
208b54c5 64 fThreshold=0;
4de874d1 65}
66
67
68AliL3HoughMaxFinder::~AliL3HoughMaxFinder()
69{
70 //Destructor
3fe49b5b 71 if(fXPeaks)
72 delete [] fXPeaks;
73 if(fYPeaks)
74 delete [] fYPeaks;
75 if(fWeight)
76 delete [] fWeight;
2061190a 77#ifndef no_root
208b54c5 78 if(fNtuppel)
79 delete fNtuppel;
2061190a 80#endif
3fe49b5b 81}
4de874d1 82
3fe49b5b 83void AliL3HoughMaxFinder::Reset()
84{
85 for(Int_t i=0; i<fNMax; i++)
86 {
87 fXPeaks[i]=0;
88 fYPeaks[i]=0;
89 fWeight[i]=0;
90 }
91 fNPeaks=0;
4de874d1 92}
93
208b54c5 94void AliL3HoughMaxFinder::CreateNtuppel()
95{
2061190a 96#ifndef no_root
208b54c5 97 //content#; neighbouring bins of the peak.
2061190a 98 fNtuppel = new TNtuple("ntuppel","Peak charateristics","kappa:phi0:weigth:content3:content5:content1:content7");
99#endif
208b54c5 100}
101
102void AliL3HoughMaxFinder::WriteNtuppel(Char_t *filename)
103{
2061190a 104#ifndef no_root
208b54c5 105 TFile *file = TFile::Open(filename,"RECREATE");
106 if(!file)
107 {
108 cerr<<"AliL3HoughMaxFinder::WriteNtuppel : Error opening file "<<filename<<endl;
109 return;
110 }
111 fNtuppel->Write();
112 file->Close();
2061190a 113#endif
208b54c5 114}
115
3fe49b5b 116void AliL3HoughMaxFinder::FindAbsMaxima()
52a2a604 117{
208b54c5 118
4cafa5fc 119 if(!fCurrentHisto)
120 {
208b54c5 121 cerr<<"AliL3HoughMaxFinder::FindAbsMaxima : No histogram"<<endl;
4fc9a6a4 122 return;
4cafa5fc 123 }
124 AliL3Histogram *hist = fCurrentHisto;
4fc9a6a4 125
4cafa5fc 126 Int_t xmin = hist->GetFirstXbin();
127 Int_t xmax = hist->GetLastXbin();
128 Int_t ymin = hist->GetFirstYbin();
129 Int_t ymax = hist->GetLastYbin();
130 Int_t bin;
95a00d93 131 Double_t value,max_value=0;
3fe49b5b 132
133 Int_t max_xbin=0,max_ybin=0;
4cafa5fc 134 for(Int_t xbin=xmin; xbin<=xmax; xbin++)
135 {
136 for(Int_t ybin=ymin; ybin<=ymax; ybin++)
137 {
138 bin = hist->GetBin(xbin,ybin);
139 value = hist->GetBinContent(bin);
140 if(value>max_value)
141 {
142 max_value = value;
4fc9a6a4 143 max_xbin = xbin;
144 max_ybin = ybin;
4cafa5fc 145 }
146 }
147 }
4fc9a6a4 148
3fe49b5b 149 if(fNPeaks > fNMax)
150 {
151 cerr<<"AliL3HoughMaxFinder::FindAbsMaxima : Array out of range : "<<fNPeaks<<endl;
152 return;
153 }
154
208b54c5 155
3fe49b5b 156 Double_t max_x = hist->GetBinCenterX(max_xbin);
157 Double_t max_y = hist->GetBinCenterY(max_ybin);
158 fXPeaks[fNPeaks] = max_x;
159 fYPeaks[fNPeaks] = max_y;
160 fWeight[fNPeaks] = (Int_t)max_value;
161 fNPeaks++;
2061190a 162#ifndef no_root
208b54c5 163 if(fNtuppel)
164 {
165 Int_t bin3 = hist->GetBin(max_xbin-1,max_ybin);
166 Int_t bin5 = hist->GetBin(max_xbin+1,max_ybin);
167 Int_t bin1 = hist->GetBin(max_xbin,max_ybin-1);
168 Int_t bin7 = hist->GetBin(max_xbin,max_ybin+1);
169
170 fNtuppel->Fill(max_x,max_y,max_value,hist->GetBinContent(bin3),hist->GetBinContent(bin5),hist->GetBinContent(bin1),hist->GetBinContent(bin7));
171 }
2061190a 172#endif
4cafa5fc 173}
174
3fe49b5b 175void AliL3HoughMaxFinder::FindBigMaxima()
4cafa5fc 176{
177
3fe49b5b 178 AliL3Histogram *hist = fCurrentHisto;
4cafa5fc 179 Int_t xmin = hist->GetFirstXbin();
180 Int_t xmax = hist->GetLastXbin();
181 Int_t ymin = hist->GetFirstYbin();
182 Int_t ymax = hist->GetLastYbin();
52a2a604 183 Int_t bin[25],bin_index;
95a00d93 184 Double_t value[25];
52a2a604 185
52a2a604 186 for(Int_t xbin=xmin+2; xbin<xmax-3; xbin++)
187 {
188 for(Int_t ybin=ymin+2; ybin<ymax-3; ybin++)
189 {
190 bin_index=0;
4cafa5fc 191 for(Int_t xb=xbin-2; xb<=xbin+2; xb++)
52a2a604 192 {
4cafa5fc 193 for(Int_t yb=ybin-2; yb<=ybin+2; yb++)
52a2a604 194 {
195 bin[bin_index]=hist->GetBin(xb,yb);
196 value[bin_index]=hist->GetBinContent(bin[bin_index]);
197 bin_index++;
198 }
199
200 }
201 if(value[12]==0) continue;
202 Int_t b=0;
203 while(1)
204 {
205 if(value[b] > value[12] || b==bin_index) break;
206 b++;
3fe49b5b 207 //printf("b %d\n",b);
52a2a604 208 }
209 if(b == bin_index)
210 {
211 //Found maxima
3fe49b5b 212 if(fNPeaks > fNMax)
213 {
214 cerr<<"AliL3HoughMaxFinder::FindBigMaxima : Array out of range "<<fNPeaks<<endl;
215 return;
216 }
217
4cafa5fc 218 Double_t max_x = hist->GetBinCenterX(xbin);
219 Double_t max_y = hist->GetBinCenterY(ybin);
3fe49b5b 220 fXPeaks[fNPeaks] = max_x;
221 fYPeaks[fNPeaks] = max_y;
222 fNPeaks++;
52a2a604 223 }
224 }
225 }
52a2a604 226}
227
208b54c5 228void AliL3HoughMaxFinder::FindMaxima(Double_t grad_x,Double_t grad_y)
4de874d1 229{
230 //Locate all the maxima in input histogram.
231 //Maxima is defined as bins with more entries than the
232 //immediately neighbouring bins.
233
4ab9f8f0 234 Int_t xmin = fCurrentHisto->GetFirstXbin();
235 Int_t xmax = fCurrentHisto->GetLastXbin();
236 Int_t ymin = fCurrentHisto->GetFirstYbin();
237 Int_t ymax = fCurrentHisto->GetLastYbin();
238 Int_t bin[9];
95a00d93 239 Double_t value[9];
4de874d1 240
4de874d1 241 for(Int_t xbin=xmin+1; xbin<xmax-1; xbin++)
242 {
243 for(Int_t ybin=ymin+1; ybin<ymax-1; ybin++)
244 {
4ab9f8f0 245 bin[0] = fCurrentHisto->GetBin(xbin-1,ybin-1);
246 bin[1] = fCurrentHisto->GetBin(xbin,ybin-1);
247 bin[2] = fCurrentHisto->GetBin(xbin+1,ybin-1);
248 bin[3] = fCurrentHisto->GetBin(xbin-1,ybin);
249 bin[4] = fCurrentHisto->GetBin(xbin,ybin);
250 bin[5] = fCurrentHisto->GetBin(xbin+1,ybin);
251 bin[6] = fCurrentHisto->GetBin(xbin-1,ybin+1);
252 bin[7] = fCurrentHisto->GetBin(xbin,ybin+1);
253 bin[8] = fCurrentHisto->GetBin(xbin+1,ybin+1);
254 value[0] = fCurrentHisto->GetBinContent(bin[0]);
255 value[1] = fCurrentHisto->GetBinContent(bin[1]);
256 value[2] = fCurrentHisto->GetBinContent(bin[2]);
257 value[3] = fCurrentHisto->GetBinContent(bin[3]);
258 value[4] = fCurrentHisto->GetBinContent(bin[4]);
259 value[5] = fCurrentHisto->GetBinContent(bin[5]);
260 value[6] = fCurrentHisto->GetBinContent(bin[6]);
261 value[7] = fCurrentHisto->GetBinContent(bin[7]);
262 value[8] = fCurrentHisto->GetBinContent(bin[8]);
263
4de874d1 264
4de874d1 265
266 if(value[4]>value[0] && value[4]>value[1] && value[4]>value[2]
267 && value[4]>value[3] && value[4]>value[5] && value[4]>value[6]
268 && value[4]>value[7] && value[4]>value[8])
269 {
270 //Found a local maxima
4ab9f8f0 271 Float_t max_x = fCurrentHisto->GetBinCenterX(xbin);
272 Float_t max_y = fCurrentHisto->GetBinCenterY(ybin);
4de874d1 273
4ab9f8f0 274 if((Int_t)value[4] <= fThreshold) continue;//central bin below threshold
0309a5ee 275 if(fNPeaks >= fNMax)
4de874d1 276 {
3060f008 277 cout<<"AliL3HoughMaxFinder::FindMaxima : Array out of range "<<fNPeaks<<endl;
4ab9f8f0 278 return;
4de874d1 279 }
3060f008 280
0309a5ee 281 /*
3fe49b5b 282 //Check the gradient:
3060f008 283 if(value[3]/value[4] > 1./grad_x || value[5]/value[4] < 1./grad_x ||
284 value[1]/value[4] < 1./grad_y || value[7]/value[4] < 1./grad_y)
3fe49b5b 285 continue;
0309a5ee 286 */
3060f008 287
288
3fe49b5b 289 fXPeaks[fNPeaks] = max_x;
290 fYPeaks[fNPeaks] = max_y;
291 fWeight[fNPeaks] = (Int_t)value[4];
292 fNPeaks++;
293
294 /*
4ab9f8f0 295 //Check if the peak is overlapping with a previous:
296 Bool_t bigger = kFALSE;
297 for(Int_t p=0; p<entries; p++)
3fe49b5b 298 {
299 if(fabs(max_x - xpeaks[p]) < kappa_overlap && fabs(max_y - ypeaks[p]) < phi_overlap)
300 {
301 bigger = kTRUE;
4ab9f8f0 302 if(value[4] > weight[p]) //this peak is bigger.
303 {
304 xpeaks[p] = max_x;
305 ypeaks[p] = max_y;
306 weight[p] = (Int_t)value[4];
307 }
308 else
309 continue; //previous peak is bigger.
310 }
311 }
312 if(!bigger) //there were no overlapping peaks.
313 {
3fe49b5b 314 xpeaks[entries] = max_x;
4ab9f8f0 315 ypeaks[entries] = max_y;
316 weight[entries] = (Int_t)value[4];
317 entries++;
318 }
3fe49b5b 319 */
4de874d1 320 }
321 else
322 continue; //not a maxima
323 }
324 }
4ab9f8f0 325
4de874d1 326}
e4c21048 327
4cafa5fc 328AliL3TrackArray *AliL3HoughMaxFinder::LookForPeaks(AliL3Histogram *hist,Int_t nbins)
52a2a604 329{
330
331 AliL3TrackArray *tracks = new AliL3TrackArray("AliL3HoughTrack");
332 AliL3HoughTrack *track;
4cafa5fc 333 Int_t xmin = hist->GetFirstXbin();
334 Int_t xmax = hist->GetLastXbin();
335 Int_t ymin = hist->GetFirstYbin();
336 Int_t ymax = hist->GetLastYbin();
52a2a604 337
338 Int_t weight_loc;
52a2a604 339 for(Int_t xbin=xmin+nbins; xbin <= xmax - nbins; xbin++)
340 {
341 for(Int_t ybin=ymin+nbins; ybin <= ymax - nbins; ybin++)
342 {
343 weight_loc=0;
344 for(Int_t xbin_loc = xbin-nbins; xbin_loc <= xbin+nbins; xbin_loc++)
345 {
346 for(Int_t ybin_loc = ybin-nbins; ybin_loc <= ybin+nbins; ybin_loc++)
347 {
348 Int_t bin_loc = hist->GetBin(xbin_loc,ybin_loc);
349 weight_loc += (Int_t)hist->GetBinContent(bin_loc);
350 }
351 }
352
353 if(weight_loc > 0)
354 {
355 track = (AliL3HoughTrack*)tracks->NextTrack();
4cafa5fc 356 track->SetTrackParameters(hist->GetBinCenterX(xbin),hist->GetBinCenterY(ybin),weight_loc);
52a2a604 357 }
358
359 }
360 }
361 tracks->QSort();
362
363 AliL3HoughTrack *track1,*track2;
364
365 for(Int_t i=1; i<tracks->GetNTracks(); i++)
366 {
367 track1 = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
368 if(!track1) continue;
4cafa5fc 369 Int_t xbin1 = hist->FindXbin(track1->GetKappa());
370 Int_t ybin1 = hist->FindXbin(track1->GetPhi0());
52a2a604 371 for(Int_t j=0; j < i; j++)
372 {
373 track2 = (AliL3HoughTrack*)tracks->GetCheckedTrack(j);
374 if(!track2) continue;
4cafa5fc 375 Int_t xbin2 = hist->FindXbin(track2->GetKappa());
376 Int_t ybin2 = hist->FindYbin(track2->GetPhi0());
52a2a604 377 if(abs(xbin1-xbin2) < 10 && abs(ybin1-ybin2) < 10)
378 {
379 tracks->Remove(i);
380 break;
381 }
382 }
383
384 }
385 tracks->Compress();
386 return tracks;
387}
388
4cafa5fc 389AliL3HoughTrack *AliL3HoughMaxFinder::CalculatePeakInWindow(Int_t *maxbin,Int_t t0,Int_t t1,Double_t t2,Int_t t3)
52a2a604 390{
4cafa5fc 391 //Try to expand the area around the maxbin +- t0
52a2a604 392
4cafa5fc 393 if(!fCurrentHisto)
52a2a604 394 {
4cafa5fc 395 printf("AliL3HoughMaxFinder::LocatePeak : No histogram\n");
396 return 0;
52a2a604 397 }
4cafa5fc 398 AliL3Histogram *hist = fCurrentHisto;
399 Int_t xmin = hist->GetFirstXbin();
400 Int_t xmax = hist->GetLastXbin();
401 Int_t ymin = hist->GetFirstYbin();
402 Int_t ymax = hist->GetLastYbin();
403
404 Int_t xlow = maxbin[0]-t0;
405 if(xlow < xmin)
406 xlow = xmin;
407 Int_t xup = maxbin[0]+t0;
408 if(xup > xmax)
409 xup = xmax;
410 Int_t ylow = maxbin[1]-t0;
411 if(ylow < ymin)
412 ylow = ymin;
413 Int_t yup = maxbin[1]+t0;
414 if(yup > ymax)
415 yup = ymax;
416
417 Int_t nbinsx = hist->GetNbinsX()+1;
418
419 Int_t *m = new Int_t[nbinsx];
420 Int_t *m_low = new Int_t[nbinsx];
421 Int_t *m_up = new Int_t[nbinsx];
422
423 recompute:
52a2a604 424
425 Int_t max_x=0,sum=0,max_xbin=0,bin;
4cafa5fc 426
427 for(Int_t xbin=xlow; xbin<=xup; xbin++)
52a2a604 428 {
4cafa5fc 429 for(Int_t ybin=ylow; ybin <= yup; ybin++)
52a2a604 430 {
431 sum = 0;
4cafa5fc 432 for(Int_t y=ybin; y <= ybin+t1; y++)
52a2a604 433 {
4cafa5fc 434 if(y>yup) break; //reached the upper limit in y.
52a2a604 435 //Inside window
436 bin = hist->GetBin(xbin,y);
437 sum += (Int_t)hist->GetBinContent(bin);
438
439 }
440 if(sum > m[xbin]) //Max value locally in this xbin
441 {
442 m[xbin]=sum;
443 m_low[xbin]=ybin;
4cafa5fc 444 m_up[xbin]=ybin + t1;
52a2a604 445 }
446
447 }
448
449 if(m[xbin] > max_x) //Max value globally in x-direction
450 {
451 max_xbin = xbin;
452 max_x = m[xbin];//sum;
453 }
454 }
52a2a604 455 //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 456 //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[max_xbin]),hist->GetBinCenterY(m_up[max_xbin]));
52a2a604 457
458 //Determine a width in the x-direction
459 Int_t x_low=0,x_up=0;
460 for(Int_t xbin=max_xbin-1; xbin >= xmin; xbin--)
461 {
4cafa5fc 462 if(m[xbin] < max_x*t2)
52a2a604 463 {
464 x_low = xbin+1;
465 break;
466 }
467 }
468 for(Int_t xbin = max_xbin+1; xbin <=xmax; xbin++)
469 {
4cafa5fc 470 if(m[xbin] < max_x*t2)
52a2a604 471 {
472 x_up = xbin-1;
473 break;
474 }
475 }
4cafa5fc 476 printf("x_low %d x_up %d\n",x_low,x_up);
477
52a2a604 478 Double_t top=0,butt=0,value,x_peak;
4cafa5fc 479 if(x_up - x_low + 1 > t3)
480 {
481 t1 -= 1;
482 printf("\nxrange out if limit x_up %d x_low %d\n\n",x_low,x_up);
483 if(t1 > 1)
484 goto recompute;
485 else
486 {
487 x_peak = hist->GetBinCenterX(max_xbin);
488 goto moveon;
489 }
490 }
52a2a604 491
4cafa5fc 492 //printf("xlow %f xup %f\n",hist->GetBinCenterX(x_low),hist->GetBinCenterX(x_up));
52a2a604 493 //printf("Spread in x %d\n",x_up-x_low +1);
494
495 //Now, calculate the center of mass in x-direction
496 for(Int_t xbin=x_low; xbin <= x_up; xbin++)
497 {
4cafa5fc 498 value = hist->GetBinCenterX(xbin);
52a2a604 499 top += value*m[xbin];
500 butt += m[xbin];
501 }
502 x_peak = top/butt;
503
4cafa5fc 504 moveon:
505
52a2a604 506 //Find the peak in y direction:
4cafa5fc 507 Int_t x_l = hist->FindXbin(x_peak);
508 if(hist->GetBinCenterX(x_l) > x_peak)
52a2a604 509 x_l--;
510
511 Int_t x_u = x_l + 1;
512
4cafa5fc 513 if(hist->GetBinCenterX(x_l) > x_peak || hist->GetBinCenterX(x_u) <= x_peak)
514 printf("\nAliL3HoughMaxFinder::FindPeak : Wrong xrange %f %f %f\n\n",hist->GetBinCenterX(x_l),x_peak,hist->GetBinCenterX(x_u));
515
516 //printf("\nxlow %f xup %f\n",hist->GetBinCenterX(x_l),hist->GetBinCenterX(x_u));
52a2a604 517
518 value=top=butt=0;
519
4cafa5fc 520 //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_l]),hist->GetBinCenterY(m_up[x_l]));
521 //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_u]),hist->GetBinCenterY(m_up[x_u]));
52a2a604 522
523 for(Int_t ybin=m_low[x_l]; ybin <= m_up[x_l]; ybin++)
524 {
4cafa5fc 525 value = hist->GetBinCenterY(ybin);
52a2a604 526 bin = hist->GetBin(x_l,ybin);
527 top += value*hist->GetBinContent(bin);
528 butt += hist->GetBinContent(bin);
529 }
530 Double_t y_peak_low = top/butt;
531
532 //printf("y_peak_low %f\n",y_peak_low);
533
534 value=top=butt=0;
535 for(Int_t ybin=m_low[x_u]; ybin <= m_up[x_u]; ybin++)
536 {
4cafa5fc 537 value = hist->GetBinCenterY(ybin);
52a2a604 538 bin = hist->GetBin(x_u,ybin);
539 top += value*hist->GetBinContent(bin);
540 butt += hist->GetBinContent(bin);
541 }
542 Double_t y_peak_up = top/butt;
543
544 //printf("y_peak_up %f\n",y_peak_up);
545
4cafa5fc 546 Double_t x_value_up = hist->GetBinCenterX(x_u);
547 Double_t x_value_low = hist->GetBinCenterX(x_l);
52a2a604 548
549 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);
550
4cafa5fc 551
52a2a604 552 //Find the weight:
553 bin = hist->FindBin(x_peak,y_peak);
554 Int_t weight = (Int_t)hist->GetBinContent(bin);
555
4cafa5fc 556 AliL3HoughTrack *track = new AliL3HoughTrack();
52a2a604 557 track->SetTrackParameters(x_peak,y_peak,weight);
558
4cafa5fc 559 //Reset area around peak
560 for(Int_t xbin=x_low; xbin<=x_up; xbin++)
561 {
562 for(Int_t ybin=m_low[xbin]; ybin<=m_up[xbin]; ybin++)
563 {
564 bin = hist->GetBin(xbin,ybin);
565 hist->SetBinContent(bin,0);
566 }
567 }
568
52a2a604 569 delete [] m;
52a2a604 570 delete [] m_low;
4cafa5fc 571 delete [] m_up;
572
573 return track;
574
575
52a2a604 576}
577
4fc9a6a4 578AliL3HoughTrack *AliL3HoughMaxFinder::FindPeakLine(Double_t rho,Double_t theta)
579{
580 //Peak finder based on a second line transformation on kappa-phi space,
581 //to use as a baseline.
582
583 if(!fCurrentHisto)
584 {
585 printf("AliL3HoughTransformer::FindPeakLine : No input histogram\n");
586 return 0;
587 }
588
589 //get the line parameters:
590 Double_t a = -1./tan(theta);
591 Double_t b = rho/sin(theta);
592
593 printf("rho %f theta %f\n",rho,theta);
594 //now, start looking along the line.
595 Int_t xmin = fCurrentHisto->GetFirstXbin();
596 Int_t xmax = fCurrentHisto->GetLastXbin();
597
598 Int_t max_weight=0;
599 Double_t max_bin[2];
600 for(Int_t xbin=xmin; xbin<=xmax; xbin++)
601 {
602 Double_t x = fCurrentHisto->GetBinCenterX(xbin);
603 Double_t y = a*x + b;
604 Int_t bin = fCurrentHisto->FindBin(x,y);
605 //printf("x %f y %f weight %d\n",x,y,fCurrentHisto->GetBinContent(bin));
606 if(fCurrentHisto->GetBinContent(bin) > max_weight)
607 {
608 max_weight = (Int_t)fCurrentHisto->GetBinContent(bin);
609 max_bin[0] = x;
610 max_bin[1] = y;
611 }
612 }
613
614 AliL3HoughTrack *track = new AliL3HoughTrack();
615 track->SetTrackParameters(max_bin[0],max_bin[1],max_weight);
616 return track;
617}
618
afd8fed4 619void AliL3HoughMaxFinder::FindPeak1(Int_t y_window,Int_t x_bin_sides)
4fc9a6a4 620{
b1886074 621 //Testing mutliple peakfinding.
622 //The algorithm searches the histogram for prepreaks by looking in windows
623 //for each bin on the xaxis. The size of these windows is controlled by y_window.
624 //Then the prepreaks are sorted according to their weight (sum inside window),
625 //and the peak positions are calculated by taking the weighted mean in both
626 //x and y direction. The size of the peak in x-direction is controlled by x_bin_sides.
627
4fc9a6a4 628 if(!fCurrentHisto)
629 {
630 printf("AliL3HoughMaxFinder::FindPeak1 : No input histogram\n");
631 return;
632 }
b1886074 633 //Int_t y_window=2;
634 //Int_t x_bin_sides=1;
635
afd8fed4 636 //Float_t max_kappa = 0.001;
637 //Float_t max_phi0 = 0.05;
b1886074 638
4fc9a6a4 639 Int_t max_sum=0;
640
641 Int_t xmin = fCurrentHisto->GetFirstXbin();
642 Int_t xmax = fCurrentHisto->GetLastXbin();
643 Int_t ymin = fCurrentHisto->GetFirstYbin();
644 Int_t ymax = fCurrentHisto->GetLastYbin();
645 Int_t nbinsx = fCurrentHisto->GetNbinsX()+1;
646
647 AxisWindow **windowPt = new AxisWindow*[nbinsx];
648 AxisWindow **anotherPt = new AxisWindow*[nbinsx];
649
650 for(Int_t i=0; i<nbinsx; i++)
651 {
652 windowPt[i] = new AxisWindow;
653 bzero((void*)windowPt[i],sizeof(AxisWindow));
654 anotherPt[i] = windowPt[i];
655 }
656
657 for(Int_t xbin=xmin; xbin<=xmax; xbin++)
658 {
659 max_sum = 0;
660 for(Int_t ybin=ymin; ybin<=ymax-y_window; ybin++)
661 {
662 Int_t sum_in_window=0;
663 for(Int_t b=ybin; b<ybin+y_window; b++)
664 {
665 //inside window
666 Int_t bin = fCurrentHisto->GetBin(xbin,b);
667 sum_in_window += (Int_t)fCurrentHisto->GetBinContent(bin);
668 }
669
670 if(sum_in_window > max_sum)
671 {
672 max_sum = sum_in_window;
673 windowPt[xbin]->ymin = ybin;
674 windowPt[xbin]->ymax = ybin + y_window;
675 windowPt[xbin]->weight = sum_in_window;
676 windowPt[xbin]->xbin = xbin;
677 }
678 }
679 }
680
681 //Sort the windows according to the weight
682 SortPeaks(windowPt,0,nbinsx);
683
4fc9a6a4 684 Float_t top,butt;
685 for(Int_t i=0; i<nbinsx; i++)
686 {
687 top=butt=0;
688 Int_t xbin = windowPt[i]->xbin;
689
690 if(xbin<xmin || xbin > xmax-1) continue;
691
b1886074 692 //Check if this is really a local maxima
693 if(anotherPt[xbin-1]->weight > anotherPt[xbin]->weight ||
694 anotherPt[xbin+1]->weight > anotherPt[xbin]->weight)
4fc9a6a4 695 continue;
696
697 for(Int_t j=windowPt[i]->ymin; j<windowPt[i]->ymax; j++)
698 {
699 //Calculate the mean in y direction:
700 Int_t bin = fCurrentHisto->GetBin(windowPt[i]->xbin,j);
701 top += (fCurrentHisto->GetBinCenterY(j))*(fCurrentHisto->GetBinContent(bin));
702 butt += fCurrentHisto->GetBinContent(bin);
703 }
afd8fed4 704
705 if(butt < fThreshold)
706 continue;
707
708 fXPeaks[fNPeaks] = fCurrentHisto->GetBinCenterX(windowPt[i]->xbin);
709 fYPeaks[fNPeaks] = top/butt;
710 fWeight[fNPeaks] = (Int_t)butt;
4ab9f8f0 711 //cout<<"mean in y "<<ypeaks[n]<<" on x "<<windowPt[i]->xbin<<" content "<<butt<<endl;
afd8fed4 712 fNPeaks++;
713 if(fNPeaks==fNMax) break;
4fc9a6a4 714 }
715
716
717 //Improve the peaks by including the region around in x.
718 Float_t ytop,ybutt;
b1886074 719 Int_t prev;
720 Int_t w;
afd8fed4 721 for(Int_t i=0; i<fNPeaks; i++)
4fc9a6a4 722 {
afd8fed4 723 Int_t xbin = fCurrentHisto->FindXbin(fXPeaks[i]);
4fc9a6a4 724 if(xbin - x_bin_sides < xmin || xbin + x_bin_sides > xmax) continue;
725 top=butt=0;
726 ytop=0,ybutt=0;
b1886074 727 w=0;
728 prev = xbin - x_bin_sides+1;
4fc9a6a4 729 for(Int_t j=xbin-x_bin_sides; j<=xbin+x_bin_sides; j++)
730 {
4ab9f8f0 731 /*
b1886074 732 //Check if the windows are overlapping:
733 if(anotherPt[j]->ymin > anotherPt[prev]->ymax) {prev=j; continue;}
734 if(anotherPt[j]->ymax < anotherPt[prev]->ymin) {prev=j; continue;}
735 prev = j;
4ab9f8f0 736 */
737
4fc9a6a4 738 top += fCurrentHisto->GetBinCenterX(j)*anotherPt[j]->weight;
739 butt += anotherPt[j]->weight;
b1886074 740
4fc9a6a4 741 for(Int_t k=anotherPt[j]->ymin; k<anotherPt[j]->ymax; k++)
742 {
743 Int_t bin = fCurrentHisto->GetBin(j,k);
744 ytop += (fCurrentHisto->GetBinCenterY(k))*(fCurrentHisto->GetBinContent(bin));
745 ybutt += fCurrentHisto->GetBinContent(bin);
95a00d93 746 w+=(Int_t)fCurrentHisto->GetBinContent(bin);
4fc9a6a4 747 }
748 }
b1886074 749
afd8fed4 750 fXPeaks[i] = top/butt;
751 fYPeaks[i] = ytop/ybutt;
752 fWeight[i] = w;
4ab9f8f0 753 //cout<<"Setting weight "<<w<<" kappa "<<xpeaks[i]<<" phi0 "<<ypeaks[i]<<endl;
b5a207b4 754
afd8fed4 755 /*
b1886074 756 //Check if this peak is overlapping with a previous:
757 for(Int_t p=0; p<i-1; p++)
758 {
afd8fed4 759 if(fabs(fXPeaks[p] - fXPeaks[i]) < max_kappa ||
760 fabs(fYPeaks[p] - fYPeaks[i]) < max_phi0)
b1886074 761 {
afd8fed4 762 fWeight[i]=0;
b1886074 763 break;
764 }
765 }
afd8fed4 766 */
4fc9a6a4 767 }
768
769 for(Int_t i=0; i<nbinsx; i++)
770 delete windowPt[i];
771 delete [] windowPt;
772 delete [] anotherPt;
773
774}
775
776void AliL3HoughMaxFinder::SortPeaks(struct AxisWindow **a,Int_t first,Int_t last)
777{
778 //General sorting routine
779 //Sort according to PeakCompare()
780
781 static struct AxisWindow *tmp;
782 static int i; // "static" to save stack space
783 int j;
784
785 while (last - first > 1) {
786 i = first;
787 j = last;
788 for (;;) {
789 while (++i < last && PeakCompare(a[i], a[first]) < 0)
790 ;
791 while (--j > first && PeakCompare(a[j], a[first]) > 0)
792 ;
793 if (i >= j)
794 break;
795
796 tmp = a[i];
797 a[i] = a[j];
798 a[j] = tmp;
799 }
800 if (j == first) {
801 ++first;
802 continue;
803 }
804 tmp = a[first];
805 a[first] = a[j];
806 a[j] = tmp;
807 if (j - first < last - (j + 1)) {
808 SortPeaks(a, first, j);
809 first = j + 1; // QSort(j + 1, last);
810 } else {
811 SortPeaks(a, j + 1, last);
812 last = j; // QSort(first, j);
813 }
814 }
815
816}
817
818Int_t AliL3HoughMaxFinder::PeakCompare(struct AxisWindow *a,struct AxisWindow *b)
819{
820 if(a->weight < b->weight) return 1;
821 if(a->weight > b->weight) return -1;
822 return 0;
823
824}
52a2a604 825
7f66cda9 826void AliL3HoughMaxFinder::FindPeak(Int_t t1,Double_t t2,Int_t t3)
e4c21048 827{
828 //Attempt of a more sophisticated peak finder.
4fc9a6a4 829 //Finds the best peak in the histogram, and returns the corresponding
830 //track object.
831
4cafa5fc 832 if(!fCurrentHisto)
833 {
834 printf("AliL3HoughMaxFinder::FindPeak : No histogram!!\n");
b1886074 835 return;
4cafa5fc 836 }
837 AliL3Histogram *hist = fCurrentHisto;
81516384 838
4cafa5fc 839 Int_t xmin = hist->GetFirstXbin();
840 Int_t xmax = hist->GetLastXbin();
841 Int_t ymin = hist->GetFirstYbin();
842 Int_t ymax = hist->GetLastYbin();
843 Int_t nbinsx = hist->GetNbinsX()+1;
e4c21048 844
845 Int_t *m = new Int_t[nbinsx];
846 Int_t *m_low = new Int_t[nbinsx];
847 Int_t *m_up = new Int_t[nbinsx];
81516384 848
849
850 recompute: //this is a goto.
851
e4c21048 852 for(Int_t i=0; i<nbinsx; i++)
853 {
854 m[i]=0;
855 m_low[i]=0;
856 m_up[i]=0;
857 }
e4c21048 858
e06900d5 859 Int_t max_x=0,sum=0,max_xbin=0,bin=0;
e4c21048 860
861 for(Int_t xbin=xmin; xbin<=xmax; xbin++)
862 {
4cafa5fc 863 for(Int_t ybin=ymin; ybin <= ymax - t1; ybin++)
e4c21048 864 {
865 sum = 0;
4cafa5fc 866 for(Int_t y=ybin; y <= ybin+t1; y++)
e4c21048 867 {
4cafa5fc 868 if(y>ymax) break;
e4c21048 869 //Inside window
870 bin = hist->GetBin(xbin,y);
871 sum += (Int_t)hist->GetBinContent(bin);
872
873 }
874 if(sum > m[xbin]) //Max value locally in this xbin
875 {
876 m[xbin]=sum;
877 m_low[xbin]=ybin;
4cafa5fc 878 m_up[xbin]=ybin + t1;
e4c21048 879 }
880
881 }
882
883 if(m[xbin] > max_x) //Max value globally in x-direction
884 {
885 max_xbin = xbin;
886 max_x = m[xbin];//sum;
887 }
888 }
889 //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 890 //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[max_xbin]),hist->GetBinCenterY(m_up[max_xbin]));
e4c21048 891
892 //Determine a width in the x-direction
893 Int_t x_low=0,x_up=0;
4cafa5fc 894
e4c21048 895 for(Int_t xbin=max_xbin-1; xbin >= xmin; xbin--)
896 {
4cafa5fc 897 if(m[xbin] < max_x*t2)
e4c21048 898 {
81516384 899 x_low = xbin+1;
e4c21048 900 break;
901 }
902 }
903 for(Int_t xbin = max_xbin+1; xbin <=xmax; xbin++)
904 {
4cafa5fc 905 if(m[xbin] < max_x*t2)
e4c21048 906 {
81516384 907 x_up = xbin-1;
e4c21048 908 break;
909 }
910 }
911
81516384 912 Double_t top=0,butt=0,value,x_peak;
913 if(x_up - x_low + 1 > t3)
914 {
915 t1 -= 1;
4cafa5fc 916 printf("\nxrange out if limit x_up %d x_low %d t1 %d\n\n",x_low,x_up,t1);
81516384 917 if(t1 > 1)
918 goto recompute;
919 else
920 {
4cafa5fc 921 x_peak = hist->GetBinCenterX(max_xbin);
81516384 922 goto moveon;
923 }
924 }
925
4cafa5fc 926 //printf("xlow %f xup %f\n",hist->GetBinCenterX(x_low),hist->GetBinCenterX(x_up));
81516384 927 //printf("Spread in x %d\n",x_up-x_low +1);
e4c21048 928
929 //Now, calculate the center of mass in x-direction
e4c21048 930 for(Int_t xbin=x_low; xbin <= x_up; xbin++)
931 {
4cafa5fc 932 value = hist->GetBinCenterX(xbin);
e4c21048 933 top += value*m[xbin];
934 butt += m[xbin];
935 }
936 x_peak = top/butt;
937
81516384 938 moveon:
e4c21048 939
940 //Find the peak in y direction:
4cafa5fc 941 Int_t x_l = hist->FindXbin(x_peak);
942 if(hist->GetBinCenterX(x_l) > x_peak)
81516384 943 x_l--;
944
e4c21048 945 Int_t x_u = x_l + 1;
81516384 946
4cafa5fc 947 if(hist->GetBinCenterX(x_l) > x_peak || hist->GetBinCenterX(x_u) <= x_peak)
948 printf("\nAliL3HoughMaxFinder::FindPeak : Wrong xrange %f %f %f\n\n",hist->GetBinCenterX(x_l),x_peak,hist->GetBinCenterX(x_u));
e4c21048 949
4cafa5fc 950 //printf("\nxlow %f xup %f\n",hist->GetBinCenterX(x_l),hist->GetBinCenterX(x_u));
e4c21048 951
952 value=top=butt=0;
953
4cafa5fc 954 //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_l]),hist->GetBinCenterY(m_up[x_l]));
955 //printf("ylow %f yup %f\n",hist->GetBinCenterY(m_low[x_u]),hist->GetBinCenterY(m_up[x_u]));
e4c21048 956
957 for(Int_t ybin=m_low[x_l]; ybin <= m_up[x_l]; ybin++)
958 {
4cafa5fc 959 value = hist->GetBinCenterY(ybin);
e4c21048 960 bin = hist->GetBin(x_l,ybin);
961 top += value*hist->GetBinContent(bin);
962 butt += hist->GetBinContent(bin);
963 }
964 Double_t y_peak_low = top/butt;
965
81516384 966 //printf("y_peak_low %f\n",y_peak_low);
e4c21048 967
968 value=top=butt=0;
969 for(Int_t ybin=m_low[x_u]; ybin <= m_up[x_u]; ybin++)
970 {
4cafa5fc 971 value = hist->GetBinCenterY(ybin);
e4c21048 972 bin = hist->GetBin(x_u,ybin);
973 top += value*hist->GetBinContent(bin);
974 butt += hist->GetBinContent(bin);
975 }
976 Double_t y_peak_up = top/butt;
977
81516384 978 //printf("y_peak_up %f\n",y_peak_up);
e4c21048 979
4cafa5fc 980 Double_t x_value_up = hist->GetBinCenterX(x_u);
981 Double_t x_value_low = hist->GetBinCenterX(x_l);
e4c21048 982
983 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);
984
985
81516384 986 //Find the weight:
95a00d93 987 //bin = hist->FindBin(x_peak,y_peak);
988 //Int_t weight = (Int_t)hist->GetBinContent(bin);
81516384 989
b1886074 990 //AliL3HoughTrack *track = new AliL3HoughTrack();
991 //track->SetTrackParameters(x_peak,y_peak,weight);
7f66cda9 992 fXPeaks[fNPeaks]=x_peak;
993 fYPeaks[fNPeaks]=y_peak;
994 fWeight[fNPeaks]=(Int_t)hist->GetBinContent(bin);
995 fNPeaks++;
52a2a604 996
997 delete [] m;
998 delete [] m_low;
999 delete [] m_up;
1000
b1886074 1001 //return track;
e4c21048 1002
52a2a604 1003
e4c21048 1004}
81516384 1005