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