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