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