meerkat is hosted by Hepforge, IPPP Durham
Meerkat  v1r3
Multidimensional kernel density estimation package
BinnedKernelDensity.cpp
Go to the documentation of this file.
1 #include <stdio.h>
2 #include <vector>
3 #include <stdlib.h>
4 
5 #include "Compression.h"
6 #include "TMath.h"
7 #include "TFile.h"
8 #include "TTree.h"
9 #include "TRandom3.h"
10 
11 #include "AbsPhaseSpace.hh"
12 #include "AbsDensity.hh"
13 #include "BinnedKernelDensity.hh"
14 
15 #include "Logger.hh"
16 
17 BinnedKernelDensity::BinnedKernelDensity(const char* pdfName,
18  AbsPhaseSpace* thePhaseSpace,
19  TTree* tree,
20  std::vector<TString> &vars,
21  std::vector<UInt_t> &binning,
22  std::vector<Double_t> &width,
23  AbsDensity* d,
24  UInt_t toyEvents,
25  UInt_t maxEvents,
26  UInt_t skipEvents
27  ) : AbsDensity(pdfName) {
28  init(thePhaseSpace, tree, vars, binning, width, d, toyEvents, maxEvents, skipEvents);
29 }
30 
31 BinnedKernelDensity::BinnedKernelDensity(const char* pdfName,
32  AbsPhaseSpace* thePhaseSpace,
33  TTree* tree,
34  const char *vars1,
35  UInt_t bins1,
36  Double_t width1,
37  AbsDensity* d,
38  UInt_t toyEvents,
39  UInt_t maxEvents,
40  UInt_t skipEvents
41  ) : AbsDensity(pdfName) {
42 
43  std::vector<TString> vars;
44  std::vector<UInt_t> binning;
45  std::vector<Double_t> width;
46 
47  vars.push_back(TString(vars1));
48  binning.push_back(bins1);
49  width.push_back(width1);
50 
51  init(thePhaseSpace, tree, vars, binning, width, d, toyEvents, maxEvents, skipEvents);
52 }
53 
54 BinnedKernelDensity::BinnedKernelDensity(const char* pdfName,
55  AbsPhaseSpace* thePhaseSpace,
56  TTree* tree,
57  const char *vars1,
58  const char *weight,
59  UInt_t bins1,
60  Double_t width1,
61  AbsDensity* d,
62  UInt_t toyEvents,
63  UInt_t maxEvents,
64  UInt_t skipEvents
65  ) : AbsDensity(pdfName) {
66 
67  std::vector<TString> vars;
68  std::vector<UInt_t> binning;
69  std::vector<Double_t> width;
70 
71  vars.push_back(TString(vars1));
72  vars.push_back(TString(weight));
73  binning.push_back(bins1);
74  width.push_back(width1);
75 
76  init(thePhaseSpace, tree, vars, binning, width, d, toyEvents, maxEvents, skipEvents);
77 }
78 
79 BinnedKernelDensity::BinnedKernelDensity(const char* pdfname,
80  AbsPhaseSpace* thephaseSpace,
81  TTree* tree,
82  const char *vars1, const char *vars2,
83  UInt_t bins1, UInt_t bins2,
84  Double_t width1, Double_t width2,
85  AbsDensity* d,
86  UInt_t toyEvents,
87  UInt_t maxEvents,
88  UInt_t skipEvents
89  ) : AbsDensity(pdfname) {
90 
91  std::vector<TString> vars;
92  std::vector<UInt_t> binning;
93  std::vector<Double_t> width;
94 
95  vars.push_back(TString(vars1));
96  vars.push_back(TString(vars2));
97  binning.push_back(bins1);
98  binning.push_back(bins2);
99  width.push_back(width1);
100  width.push_back(width2);
101 
102  init(thephaseSpace, tree, vars, binning, width, d, toyEvents, maxEvents, skipEvents);
103 }
104 
105 BinnedKernelDensity::BinnedKernelDensity(const char* pdfname,
106  AbsPhaseSpace* thephaseSpace,
107  TTree* tree,
108  const char *vars1, const char *vars2,
109  const char *weight,
110  UInt_t bins1, UInt_t bins2,
111  Double_t width1, Double_t width2,
112  AbsDensity* d,
113  UInt_t toyEvents,
114  UInt_t maxEvents,
115  UInt_t skipEvents
116  ) : AbsDensity(pdfname) {
117 
118  std::vector<TString> vars;
119  std::vector<UInt_t> binning;
120  std::vector<Double_t> width;
121 
122  vars.push_back(TString(vars1));
123  vars.push_back(TString(vars2));
124  vars.push_back(TString(weight));
125  binning.push_back(bins1);
126  binning.push_back(bins2);
127  width.push_back(width1);
128  width.push_back(width2);
129 
130  init(thephaseSpace, tree, vars, binning, width, d, toyEvents, maxEvents, skipEvents);
131 }
132 
133 BinnedKernelDensity::BinnedKernelDensity(const char* pdfname,
134  AbsPhaseSpace* thephaseSpace,
135  TTree* tree,
136  const char *vars1, const char *vars2, const char *vars3,
137  UInt_t bins1, UInt_t bins2, UInt_t bins3,
138  Double_t width1, Double_t width2, Double_t width3,
139  AbsDensity* d,
140  UInt_t toyEvents,
141  UInt_t maxEvents,
142  UInt_t skipEvents
143  ) : AbsDensity(pdfname) {
144 
145  std::vector<TString> vars;
146  std::vector<UInt_t> binning;
147  std::vector<Double_t> width;
148 
149  vars.push_back(TString(vars1));
150  vars.push_back(TString(vars2));
151  vars.push_back(TString(vars3));
152  binning.push_back(bins1);
153  binning.push_back(bins2);
154  binning.push_back(bins3);
155  width.push_back(width1);
156  width.push_back(width2);
157  width.push_back(width3);
158 
159  init(thephaseSpace, tree, vars, binning, width, d, toyEvents, maxEvents, skipEvents);
160 }
161 
162 BinnedKernelDensity::BinnedKernelDensity(const char* pdfname,
163  AbsPhaseSpace* thephaseSpace,
164  TTree* tree,
165  const char *vars1, const char *vars2, const char *vars3,
166  const char *weight,
167  UInt_t bins1, UInt_t bins2, UInt_t bins3,
168  Double_t width1, Double_t width2, Double_t width3,
169  AbsDensity* d,
170  UInt_t toyEvents,
171  UInt_t maxEvents,
172  UInt_t skipEvents
173  ) : AbsDensity(pdfname) {
174 
175  std::vector<TString> vars;
176  std::vector<UInt_t> binning;
177  std::vector<Double_t> width;
178 
179  vars.push_back(TString(vars1));
180  vars.push_back(TString(vars2));
181  vars.push_back(TString(vars3));
182  vars.push_back(TString(weight));
183  binning.push_back(bins1);
184  binning.push_back(bins2);
185  binning.push_back(bins3);
186  width.push_back(width1);
187  width.push_back(width2);
188  width.push_back(width3);
189 
190  init(thephaseSpace, tree, vars, binning, width, d, toyEvents, maxEvents, skipEvents);
191 }
192 
193 BinnedKernelDensity::BinnedKernelDensity(const char* pdfname,
194  AbsPhaseSpace* thephaseSpace,
195  TTree* tree,
196  const char *vars1, const char *vars2, const char *vars3, const char *vars4,
197  UInt_t bins1, UInt_t bins2, UInt_t bins3, UInt_t bins4,
198  Double_t width1, Double_t width2, Double_t width3, Double_t width4,
199  AbsDensity* d,
200  UInt_t toyEvents,
201  UInt_t maxEvents,
202  UInt_t skipEvents
203  ) : AbsDensity(pdfname) {
204 
205  std::vector<TString> vars;
206  std::vector<UInt_t> binning;
207  std::vector<Double_t> width;
208 
209  vars.push_back(TString(vars1));
210  vars.push_back(TString(vars2));
211  vars.push_back(TString(vars3));
212  vars.push_back(TString(vars4));
213  binning.push_back(bins1);
214  binning.push_back(bins2);
215  binning.push_back(bins3);
216  binning.push_back(bins4);
217  width.push_back(width1);
218  width.push_back(width2);
219  width.push_back(width3);
220  width.push_back(width4);
221 
222  init(thephaseSpace, tree, vars, binning, width, d, toyEvents, maxEvents, skipEvents);
223 }
224 
225 BinnedKernelDensity::BinnedKernelDensity(const char* pdfname,
226  AbsPhaseSpace* thephaseSpace,
227  TTree* tree,
228  const char *vars1, const char *vars2, const char *vars3, const char *vars4,
229  const char *weight,
230  UInt_t bins1, UInt_t bins2, UInt_t bins3, UInt_t bins4,
231  Double_t width1, Double_t width2, Double_t width3, Double_t width4,
232  AbsDensity* d,
233  UInt_t toyEvents,
234  UInt_t maxEvents,
235  UInt_t skipEvents
236  ) : AbsDensity(pdfname) {
237 
238  std::vector<TString> vars;
239  std::vector<UInt_t> binning;
240  std::vector<Double_t> width;
241 
242  vars.push_back(TString(vars1));
243  vars.push_back(TString(vars2));
244  vars.push_back(TString(vars3));
245  vars.push_back(TString(vars4));
246  vars.push_back(TString(weight));
247  binning.push_back(bins1);
248  binning.push_back(bins2);
249  binning.push_back(bins3);
250  binning.push_back(bins4);
251  width.push_back(width1);
252  width.push_back(width2);
253  width.push_back(width3);
254  width.push_back(width4);
255 
256  init(thephaseSpace, tree, vars, binning, width, d, toyEvents, maxEvents, skipEvents);
257 }
258 
259 BinnedKernelDensity::BinnedKernelDensity(const char* pdfname,
260  AbsPhaseSpace* thephaseSpace,
261  TTree* tree,
262  const char *vars1, const char *vars2, const char *vars3, const char *vars4, const char *vars5,
263  UInt_t bins1, UInt_t bins2, UInt_t bins3, UInt_t bins4, UInt_t bins5,
264  Double_t width1, Double_t width2, Double_t width3, Double_t width4, Double_t width5,
265  AbsDensity* d,
266  UInt_t toyEvents,
267  UInt_t maxEvents,
268  UInt_t skipEvents
269  ) : AbsDensity(pdfname) {
270 
271  std::vector<TString> vars;
272  std::vector<UInt_t> binning;
273  std::vector<Double_t> width;
274 
275  vars.push_back(TString(vars1));
276  vars.push_back(TString(vars2));
277  vars.push_back(TString(vars3));
278  vars.push_back(TString(vars4));
279  vars.push_back(TString(vars5));
280  binning.push_back(bins1);
281  binning.push_back(bins2);
282  binning.push_back(bins3);
283  binning.push_back(bins4);
284  binning.push_back(bins5);
285  width.push_back(width1);
286  width.push_back(width2);
287  width.push_back(width3);
288  width.push_back(width4);
289  width.push_back(width5);
290 
291  init(thephaseSpace, tree, vars, binning, width, d, toyEvents, maxEvents, skipEvents);
292 }
293 
294 BinnedKernelDensity::BinnedKernelDensity(const char* pdfname,
295  AbsPhaseSpace* thephaseSpace,
296  TTree* tree,
297  const char *vars1, const char *vars2, const char *vars3, const char *vars4, const char *vars5,
298  const char *weight,
299  UInt_t bins1, UInt_t bins2, UInt_t bins3, UInt_t bins4, UInt_t bins5,
300  Double_t width1, Double_t width2, Double_t width3, Double_t width4, Double_t width5,
301  AbsDensity* d,
302  UInt_t toyEvents,
303  UInt_t maxEvents,
304  UInt_t skipEvents
305  ) : AbsDensity(pdfname) {
306 
307  std::vector<TString> vars;
308  std::vector<UInt_t> binning;
309  std::vector<Double_t> width;
310 
311  vars.push_back(TString(vars1));
312  vars.push_back(TString(vars2));
313  vars.push_back(TString(vars3));
314  vars.push_back(TString(vars4));
315  vars.push_back(TString(vars5));
316  vars.push_back(TString(weight));
317  binning.push_back(bins1);
318  binning.push_back(bins2);
319  binning.push_back(bins3);
320  binning.push_back(bins4);
321  binning.push_back(bins5);
322  width.push_back(width1);
323  width.push_back(width2);
324  width.push_back(width3);
325  width.push_back(width4);
326  width.push_back(width5);
327 
328  init(thephaseSpace, tree, vars, binning, width, d, toyEvents, maxEvents, skipEvents);
329 }
330 
331 BinnedKernelDensity::~BinnedKernelDensity() {
332 
333 }
334 
336 void BinnedKernelDensity::init(AbsPhaseSpace* thephaseSpace,
337  TTree* tree,
338  std::vector<TString> &vars,
339  std::vector<UInt_t> &binning,
340  std::vector<Double_t> &width,
341  AbsDensity* d,
342  UInt_t toyEvents,
343  UInt_t maxEvents,
344  UInt_t skipEvents
345  ) {
346 
347  m_phaseSpace = thephaseSpace;
348  m_binning = binning;
349  m_width = width;
350  m_approxDensity = d;
351  m_dim = m_phaseSpace->dimensionality();
352 
353  m_fractionalMode = false;
354 
355  Logger::print(0,"%20.20s INFO: Creating binned kernel density over %dD phase space\n", m_name, m_dim );
356 
357  if (m_binning.size() != m_dim) {
358  Logger::print(0,"%20.20s ERROR: Dimensionality of phase space (%d) does not match binning vector size (%d)\n",
359  m_name, m_dim, (UInt_t)m_binning.size());
360  abort();
361  }
362 
363  if (m_width.size() != m_dim) {
364  Logger::print(2,"%20.20s ERROR: Dimensionality of phase space (%d) does not match kernel width vector size (%d)\n",
365  m_name, m_dim, (UInt_t)m_width.size());
366  abort();
367  }
368 
369  UInt_t size = 1;
370  std::vector<UInt_t>::iterator i;
371  for (i=m_binning.begin(); i!=m_binning.end(); i++) {
372  size *= (*i);
373  }
374 
375  Logger::print(0,"%20.20s INFO: Map size=%d\n", m_name, size);
376  if (size > 20000000) {
377  Logger::print(2,"%20.20s ERROR: Map size (%d) too large!\n", m_name, size);
378  abort();
379  }
380 
381  m_map.resize(size);
382  m_approxMap.resize(size);
383 
384  fillMapFromTree(tree, vars, maxEvents, skipEvents);
385  fillMapFromDensity(m_approxDensity, toyEvents);
386 
387  normalise();
388 
389 }
390 
392 UInt_t BinnedKernelDensity::iterToIndex( std::vector<UInt_t> &iter ) {
393  UInt_t index = 0;
394  Int_t n;
395  for (n=m_dim-1; n>=0; n--) {
396  if (n==(Int_t)m_dim-1) {
397  index = iter[n];
398  } else {
399  index = index*m_binning[n] + iter[n];
400  }
401  }
402  return index;
403 }
404 
405 void BinnedKernelDensity::addToMap(std::vector<Double_t> &map, std::vector<Double_t> &point, Double_t weight) {
406 
407  // Fill the map
408  std::vector<UInt_t> initBin(m_dim);
409  std::vector<UInt_t> finalBin(m_dim);
410  std::vector<Double_t> lowLimit(m_dim);
411  std::vector<Double_t> coeff(m_dim);
412 
413  UInt_t size = map.size();
414 
415  // Calculate the initial and final N-dim bins
416  Int_t n;
417  for (n=0; n<(Int_t)m_dim; n++) {
418  Double_t x1 = point[n] - m_width[n];
419  Double_t x2 = point[n] + m_width[n];
420 
421  Double_t low = m_phaseSpace->lowerLimit(n);
422  Double_t up = m_phaseSpace->upperLimit(n);
423  Int_t i1 = (int)TMath::Ceil( (x1-low)/(up-low)*(m_binning[n]-1) );
424  if (i1 < 0) i1 = 0;
425  if (i1 >= (Int_t)m_binning[n]) i1 = m_binning[n] - 1;
426  Int_t i2 = (int)TMath::Floor( (x2-low)/(up-low)*(m_binning[n]-1) );
427  if (i2 < 0) i2 = 0;
428  if (i2 >= (Int_t)m_binning[n]) i2 = m_binning[n] - 1;
429 
430  if (i1 > i2) {
431  Logger::print(2,"%20.20s ERROR: i1(%d)>i2(%d) for n=%d, bin size is larger than kernel width!\n", m_name, i1, i2, n);
432  abort();
433  }
434 
435  initBin[n] = i1;
436  finalBin[n] = i2;
437 
438  lowLimit[n] = phaseSpace()->lowerLimit(n);
439  coeff[n] = ( phaseSpace()->upperLimit(n) - lowLimit[n] ) / ((Double_t)m_binning[n]-1);
440  lowLimit[n] -= point[n];
441 
442  coeff[n] /= m_width[n];
443  lowLimit[n] /= m_width[n];
444  }
445 
446  // Loop through the map and fill bins
447  std::vector<UInt_t> iter(m_dim);
448  for (n=0; n<(Int_t)m_dim; n++) iter[n] = initBin[n];
449 
450  do {
451 
452  UInt_t index = iterToIndex( iter );
453 
454  if (index >= size) {
455  Logger::print(2,"%20.20s ERROR: index (%d) is larger than array size (%d)\n", m_name, index, size);
456  abort();
457  } else {
458  Double_t sqsum = 0;
459  for (n=0; n<(Int_t)m_dim; n++) {
460 // Double_t low = m_phaseSpace->lowerLimit(n);
461 // Double_t up = m_phaseSpace->upperLimit(n);
462 // Double_t binx = low + (Double_t)iter[n]/((Double_t)m_binning[n]-1)*(up-low);
463 // Double_t dx = binx/m_width[n];
464  Double_t dx = lowLimit[n] + (Double_t)iter[n]*coeff[n];
465  if (fabs(dx) < 1.) sqsum += dx*dx;
466  }
467  if (sqsum < 1.) map[index] += weight*(1.-sqsum);
468 
469  }
470 
471  // Increase iterator
472  Bool_t run = 0;
473  for (n=0; n<(Int_t)m_dim; n++) {
474  if (iter[n] < finalBin[n]) {
475  iter[n]++;
476  run = 1;
477  break;
478  } else {
479  iter[n] = initBin[n];
480  }
481  }
482  if (!run) break;
483 
484  } while(1);
485 
486 }
487 
488 void BinnedKernelDensity::fillMapFromTree( TTree* tree, std::vector<TString> &vars, UInt_t maxEvents, UInt_t skipEvents) {
489 
490  if (vars.size() != m_dim && vars.size() != m_dim + 1) {
491  Logger::print(2,"%20.20s ERROR: Number of TTree variables (%d) in tree \"%s\" does not correspond to phase space dimensionality (%d)\n",
492  m_name, (UInt_t)vars.size(), tree->GetName(), m_phaseSpace->dimensionality() );
493  abort();
494  }
495 
496  if (vars.size() == m_dim + 1) {
497  Logger::print(0,"%20.20s INFO: Using variable \"%s\" as weight\n",
498  m_name, vars[m_dim].Data() );
499  }
500 
501  UInt_t nvars = vars.size();
502 
503  tree->ResetBranchAddresses();
504 
505  Long64_t nentries = tree->GetEntries();
506  if (maxEvents > 0 && skipEvents + maxEvents < nentries) nentries = skipEvents + maxEvents;
507 
508  Logger::print(0,"%20.20s INFO: Will read %lld events (skipping first %d)\n",
509  m_name, nentries-skipEvents, skipEvents );
510 
511  Long64_t i;
512 
513  std::vector<Float_t> varArray(nvars);
514 
515  Int_t n;
516  for (n=0; n < (Int_t)nvars; n++) {
517  Logger::print(0,"%20.20s INFO: Will read branch \"%s\" from tree \"%s\"\n", m_name, vars[n].Data(), tree->GetName());
518  Int_t status = tree->SetBranchAddress(vars[n], &( varArray[n] ));
519  if (status < 0) {
520  if (n+1 == (Int_t)nvars) {
521  Logger::print(1,"%20.20s WARNING: Weight branch \"%s\" not found, reverting to unweighted mode!\n", m_name, vars[n].Data() );
522  nvars--;
523  break;
524  }
525  Logger::print(2,"%20.20s ERROR: Error setting branch, status=%d\n", m_name, status);
526  abort();
527  }
528  }
529 
530  std::vector<Double_t> point(m_dim);
531 
532  Long64_t nout = 0;
533 
534  Logger::setTimer();
535 
536  for(i=skipEvents; i<nentries; i++) {
537  tree->GetEntry(i);
538  for (n=0; n<(Int_t)m_dim; n++) {
539  point[n] = varArray[n];
540  }
541  Double_t weight = 1.;
542  if (m_dim + 1 == nvars)
543  weight = varArray[m_dim];
544 
545  if (!m_phaseSpace->withinLimits( point )) {
546  nout ++;
547 // Logger::print(0,"%20.20s WARNING: Ntuple point (", m_name);
548 // for (n=0; n<(Int_t)nvars; n++) {
549 // Logger::print(0,"%f ", point[n]);
550 // }
551 // Logger::print(0,") outside phase space\n");
552  } else {
553  addToMap(m_map, point, weight);
554  }
555 
556  if (i % 100 == 0 && Logger::timer(2)) {
557  Logger::print(0,"%20.20s INFO: Read %lld/%lld events (%f%%), %lld out\n", m_name, i-skipEvents, nentries-skipEvents,
558  100.*float(i-skipEvents)/float(nentries-skipEvents), nout);
559  }
560  }
561 
562  Logger::print(0,"%20.20s INFO: %lld events read in from \"%s\", %lld out\n", m_name, nentries-nout, tree->GetName(), nout );
563 
564 }
565 
566 void BinnedKernelDensity::fillMapFromDensity(AbsDensity* theDensity, UInt_t toyEvents) {
567 
568  std::vector<Double_t> x(m_dim);
569 // UInt_t size = m_approxMap.size();
570 
571  if (theDensity == 0) {
572  Logger::print(0,"%20.20s INFO: Will use uniform density for approximation\n", m_name);
573  } else {
574  Logger::print(0,"%20.20s INFO: Will use density \"%s\" for approximation\n", m_name, theDensity->name());
575  }
576 
577  if (toyEvents == 0) {
578  // Fill map in nodes of the binning
579  Logger::print(0,"%20.20s INFO: Convolution of approx. density using rectangular grid\n", m_name);
580 
581  // Create iterator vector
582  std::vector<UInt_t> iter(m_dim);
583  Int_t j;
584  for (j=0; j<(Int_t)m_dim; j++) {
585  iter[j] = 0;
586  }
587 
588  // Loop through the nodes
589  Int_t index = 0;
590 
591  Logger::setTimer();
592  do {
593  index++;
594  for (j=m_dim-1; j>=0; j--) {
595  Double_t low = m_phaseSpace->lowerLimit(j);
596  Double_t up = m_phaseSpace->upperLimit(j);
597  x[j] = low + (Double_t)iter[j]/((Double_t)m_binning[j]-1)*(up-low);
598  }
599 
600  Double_t e = 1.;
601  if (theDensity) e = theDensity->density(x);
602 
603  if ((index % 100) == 0 && Logger::timer(2))
604  Logger::print(0,"%20.20s INFO: Index %d, density=%f\n", m_name, index, e);
605 
606  if (m_phaseSpace->withinLimits(x)) addToMap(m_approxMap, x, e);
607 
608  // Increase iterator
609  Bool_t run = 0;
610  for (j=0; j<(Int_t)m_dim; j++) {
611  if (iter[j] < m_binning[j]-1) {
612  iter[j]++;
613  run = 1;
614  break;
615  } else {
616  iter[j] = 0;
617  }
618  }
619  if (!run) break;
620 
621  } while(1);
622 
623  } else {
624 
625  // Fill map from random points
626 
627  Logger::print(0,"%20.20s INFO: Convolution of approx. density using MC with %d events\n", m_name, toyEvents);
628  Int_t i;
629  std::vector<Double_t> lower(m_dim);
630  std::vector<Double_t> coeff(m_dim);
631  for (i=0; i<(Int_t)m_dim; i++) {
632  lower[i] = m_phaseSpace->lowerLimit(i);
633  coeff[i] = m_phaseSpace->upperLimit(i) - lower[i];
634  }
635 
636  Logger::setTimer();
637  for (i=0; i<(Int_t)toyEvents; i++) {
638 
639  Bool_t success = 0;
640  UInt_t t;
641  for (t = 0; t < m_maxTries; t++) {
642 
643  // Generate random point
644  UInt_t var;
645  for (var = 0; var < m_dim; var++) {
646  x[var] = lower[var] + m_rnd.Rndm()*coeff[var];
647  }
648 
649  Bool_t inPhsp = m_phaseSpace->withinLimits(x);
650  if (inPhsp) {
651  success = 1;
652  break;
653  }
654  }
655  if (!success) {
656  Logger::print(1,"%20.20s WARNING: failed to generate a point within phase space after %d tries\n", m_name, m_maxTries);
657  continue;
658  }
659 
660  Double_t e = 1;
661  if (theDensity) e = theDensity->density(x);
662 
663  if ((i % 100) == 0 && Logger::timer(1))
664  Logger::print(0,"%20.20s INFO: toy event %d/%d (%f%%), density=%f\n", m_name, i, toyEvents, 100*float(i)/float(toyEvents), e);
665 
666  addToMap(m_approxMap, x, e);
667 
668  }
669  }
670 }
671 
673 void BinnedKernelDensity::writeToFile(const char* filename) {
674  if (TString(filename).EndsWith(".root") ) {
675  writeToRootFile(filename);
676  } else {
677  writeToTextFile(filename);
678  }
679 }
680 
682 void BinnedKernelDensity::writeToTextFile(const char* filename) {
683 
684  Logger::print(0,"%20.20s INFO: Writing binned density to text file \"%s\"\n", m_name, filename );
685 
686  FILE* file = fopen(filename, "w+");
687  fprintf(file, "%d\n", m_dim);
688 
689  Int_t j;
690  for (j=0; j<(Int_t)m_dim; j++) {
691  fprintf(file, "%d\n", m_binning[j]);
692  }
693 
694  // Zero iterator vector
695  std::vector<Double_t> x(m_dim);
696  std::vector<UInt_t> iter(m_dim);
697  for (j=0; j<(Int_t)m_dim; j++) {
698  iter[j] = 0;
699  }
700 
701  UInt_t size = m_map.size();
702 
703  // Loop through the bins
704  do {
705 
706  UInt_t index = 0;
707  for (j=m_dim-1; j>=0; j--) {
708  Double_t low = m_phaseSpace->lowerLimit(j);
709  Double_t up = m_phaseSpace->upperLimit(j);
710  x[j] = low + (Double_t)iter[j]/((Double_t)m_binning[j]-1)*(up-low);
711  if (j==(Int_t)m_dim-1) {
712  index = iter[j];
713  } else {
714  index = index*m_binning[j] + iter[j];
715  }
716 // if (iter[0] == 0)
717 // Logger::print(0," Dim%d, bin%d, x=%f\n", j, iter[j], x[j]);
718  }
719 
720  if (index >= size) {
721  Logger::print(2,"%20.20s ERROR: index (%d) is larger than array size (%d)\n", m_name, index, size);
722  abort();
723  } else {
724  fprintf(file, "%f %d\n", density(x), m_phaseSpace->withinLimits(x) );
725  }
726 
727  Bool_t run = 0;
728  for (j=0; j<(Int_t)m_dim; j++) {
729  if (iter[j] < m_binning[j]-1) {
730  iter[j]++;
731  run = 1;
732  break;
733  } else {
734  iter[j] = 0;
735  }
736  }
737  if (!run) break;
738 
739  } while(1);
740  fclose(file);
741 }
742 
744 void BinnedKernelDensity::writeToRootFile(const char* filename) {
745 
746  Logger::print(0,"%20.20s INFO: Writing binned density to ROOT file \"%s\"\n", m_name, filename );
747 
748  TDirectory* curr_dir = gDirectory;
749 
750  TFile* file = TFile::Open(filename, "RECREATE");
751  file->SetCompressionSettings(ROOT::CompressionSettings(ROOT::kLZMA, 9)); // Maximum compression
752  TTree dimTree("DimTree", "DimTree");
753 
754  Int_t bins;
755  dimTree.Branch("bins",&bins,"bins/I");
756 
757  Int_t j;
758  for (j=0; j<(Int_t)m_dim; j++) {
759  bins = m_binning[j];
760  dimTree.Fill();
761  }
762  dimTree.Write();
763 
764  // Zero iterator vector
765  std::vector<Double_t> x(m_dim);
766  std::vector<UInt_t> iter(m_dim);
767  for (j=0; j<(Int_t)m_dim; j++) {
768  iter[j] = 0;
769  }
770 
771  UInt_t size = m_map.size();
772 
773  TTree mapTree("MapTree", "MapTree");
774 
775  Char_t inphsp;
776  Float_t dens;
777  mapTree.Branch("dens", &dens,"dens/F");
778  mapTree.Branch("inphsp",&inphsp,"inphsp/B");
779 
780  // Loop through the bins
781  do {
782 
783  UInt_t index = 0;
784  for (j=m_dim-1; j>=0; j--) {
785  Double_t low = m_phaseSpace->lowerLimit(j);
786  Double_t up = m_phaseSpace->upperLimit(j);
787  x[j] = low + (Double_t)iter[j]/((Double_t)m_binning[j]-1)*(up-low);
788  if (j==(Int_t)m_dim-1) {
789  index = iter[j];
790  } else {
791  index = index*m_binning[j] + iter[j];
792  }
793 // if (iter[0] == 0)
794 // Logger::print(0," Dim%d, bin%d, x=%f\n", j, iter[j], x[j]);
795  }
796 
797  if (index >= size) {
798  Logger::print(2,"%20.20s ERROR: index (%d) is larger than array size (%d)\n", m_name, index, size);
799  abort();
800  } else {
801  dens = density(x);
802  inphsp = m_phaseSpace->withinLimits(x);
803  mapTree.Fill();
804  }
805 
806  Bool_t run = 0;
807  for (j=0; j<(Int_t)m_dim; j++) {
808  if (iter[j] < m_binning[j]-1) {
809  iter[j]++;
810  run = 1;
811  break;
812  } else {
813  iter[j] = 0;
814  }
815  }
816  if (!run) break;
817 
818  } while(1);
819 
820  mapTree.Write();
821  file->Close();
822 
823  gDirectory = curr_dir;
824 }
825 
827 void BinnedKernelDensity::normalise(void) {
828 
829  Logger::print(0,"%20.20s INFO: Normalising density\n", m_name);
830 
831  Double_t sum = 0;
832  UInt_t num = 0;
833 
834  // Zero iterator vector
835  std::vector<Double_t> x(m_dim);
836  std::vector<UInt_t> iter(m_dim);
837  UInt_t size = m_map.size();
838 
839  Int_t j;
840 
841  for (j=0; j<(Int_t)m_dim; j++) {
842  iter[j] = 0;
843  }
844 
845  // Loop through the bins
846  do {
847 
848  UInt_t index = 0;
849  for (j=m_dim-1; j>=0; j--) {
850  Double_t low = m_phaseSpace->lowerLimit(j);
851  Double_t up = m_phaseSpace->upperLimit(j);
852  x[j] = low + (Double_t)iter[j]/((Double_t)m_binning[j]-1)*(up-low);
853  if (j==(Int_t)m_dim-1) {
854  index = iter[j];
855  } else {
856  index = index*m_binning[j] + iter[j];
857  }
858 // if (iter[0] == 0)
859 // Logger::print(0," Dim%d, bin%d, x=%f\n", j, iter[j], x[j]);
860  }
861 
862  if (index >= size) {
863  Logger::print(2,"%20.20s ERROR: index (%d) is larger than array size (%d)\n", m_name, index, size);
864  abort();
865  } else {
866  if (m_phaseSpace->withinLimits(x)) {
867  sum += density(x);
868  num ++;
869  }
870  }
871 
872  Bool_t run = 0;
873  for (j=0; j<(Int_t)m_dim; j++) {
874  if (iter[j] < m_binning[j]-1) {
875  iter[j]++;
876  run = 1;
877  break;
878  } else {
879  iter[j] = 0;
880  }
881  }
882  if (!run) break;
883  } while(1);
884 
885  sum /= (Double_t)num;
886 
887  Logger::print(0,"%20.20s INFO: Average PDF value before normalisation is %f\n", m_name, sum);
888 
889  // Loop through the map and scale its entries
890  for (j=0; j<(Int_t)size; j++) m_map[j] /= sum;
891 
892 }
893 
894 
895 Double_t BinnedKernelDensity::mapDensity(std::vector<Double_t> &map, std::vector<Double_t> &x) {
896 
897  Int_t j;
898  std::vector<UInt_t> ivect(m_dim);
899 
900  for (j=0; j<(Int_t)m_dim; j++) {
901  Double_t low = m_phaseSpace->lowerLimit(j);
902  Double_t up = m_phaseSpace->upperLimit(j);
903  Double_t xj = x[j];
904  if (xj < low || xj > up) {
905  return 0.;
906  }
907  Int_t ij = (Int_t)floor((xj-low)/(up-low)*(m_binning[j]-1));
908 
909  if (ij == (Int_t)m_binning[j]-1) ij--;
910 
911  ivect[j] = ij;
912  }
913 
914  Double_t e = 0.;
915  Double_t wsum = 0.;
916 
917  std::vector<UInt_t> iter(m_dim);
918 
919  for (j=0; j<(Int_t)m_dim; j++) {
920  iter[j] = 0;
921  }
922 
923  // Loop through the vertices of the N-dim cube
924  do {
925 
926  // Calculate weight
927  Double_t weight = 1;
928  for (j=0; j<(Int_t)m_dim; j++) {
929  Double_t low = m_phaseSpace->lowerLimit(j);
930  Double_t up = m_phaseSpace->upperLimit(j);
931 
932  Double_t xj1 = low + ((Double_t)ivect[j] )/((Double_t)m_binning[j]-1.)*(up-low);
933  Double_t xj2 = low + ((Double_t)ivect[j]+1.)/((Double_t)m_binning[j]-1.)*(up-low);
934 
935 // if (x[j] < xj1 || x[j] > xj2) {
936 // Logger::print(0,"%20.20s WARNING: dim %d: x=%f, x1=%f, x2=%f\n", m_name, j, x[j], xj1, xj2);
937 // }
938 
939  Double_t fweight;
940  if (iter[j] == 0) {
941  fweight = 1. - (x[j]-xj1)/(xj2-xj1);
942  } else {
943  fweight = (x[j]-xj1)/(xj2-xj1);
944  }
945  if (fweight < 0.) {
946  Logger::print(1,"%20.20s WARNING: dim %d: x=%f, weight=%f\n", m_name, j, x[j], fweight);
947  fweight = 0.;
948  }
949  if (fweight > 1.) {
950  Logger::print(1,"%20.20s WARNING: dim %d: x=%f, weight=%f\n", m_name, j, x[j], fweight);
951  fweight = 1.;
952  }
953 
954  weight *= fweight;
955 
956 // Logger::print(0,"DEBUG: Weight fraction: dim%d, i=%d, x=%f, x1=%f, x2=%f, fweight=%f\n", j, iter[j], x[j], xj1, xj2, fweight);
957 
958  }
959 
960  UInt_t index = 0;
961  for (j=m_dim-1; j>=0; j--) {
962  Int_t ij = ivect[j] + iter[j];
963  if (j==(Int_t)m_dim-1) {
964  index = ij;
965  } else {
966  index = index*m_binning[j] + ij;
967  }
968  }
969 
970  e += weight*map[index];
971  wsum += weight;
972 
973 // Logger::print(0,"DEBUG: Weight=%f, index=%d, density=%f\n", weight, index, m_map[index]);
974 
975  // Increment iterator
976  Bool_t run = 0;
977  for (j=0; j<(Int_t)m_dim; j++) {
978  if (iter[j] == 0) {
979  iter[j]++;
980  run = 1;
981  break;
982  } else {
983  iter[j] = 0;
984  }
985  }
986  if (!run) break;
987 
988  } while(1);
989 
990 // Logger::print(0,"density=%f, wsum=%f\n", e, wsum);
991 
992  return e;
993 
994 }
995 
996 Double_t BinnedKernelDensity::density(std::vector<Double_t> &x) {
997  Double_t a = mapDensity(m_approxMap, x);
998  if (a>0.) {
999  if (m_approxDensity && !m_fractionalMode) {
1000  return mapDensity(m_map, x)/a*m_approxDensity->density(x);
1001  } else {
1002  return mapDensity(m_map, x)/a;
1003  }
1004  } else {
1005  return 0.;
1006  }
1007 }
AbsDensity::name
const char * name(void)
Return the name of the PDF.
Definition: AbsDensity.hh:106
BinnedKernelDensity::writeToRootFile
void writeToRootFile(const char *fileName)
Write the binned density into a ROOT file.
Definition: BinnedKernelDensity.cpp:744
BinnedKernelDensity::m_phaseSpace
AbsPhaseSpace * m_phaseSpace
Reference to phase space.
Definition: BinnedKernelDensity.hh:466
BinnedKernelDensity::m_width
std::vector< Double_t > m_width
Vector of kernel widths.
Definition: BinnedKernelDensity.hh:472
AbsPhaseSpace::upperLimit
virtual Double_t upperLimit(UInt_t var)=0
Return upper allowed limit of the variable.
AbsDensity
Abstract class which defines probability density interface.
Definition: AbsDensity.hh:16
BinnedKernelDensity::m_map
std::vector< Double_t > m_map
Bin map of estimated PDF.
Definition: BinnedKernelDensity.hh:457
BinnedKernelDensity::BinnedKernelDensity
BinnedKernelDensity(const char *pdfName, AbsPhaseSpace *thePhaseSpace, TTree *tree, const char *vars1, UInt_t bins1, Double_t width1, AbsDensity *approx=0, UInt_t toyEvents=0, UInt_t maxEvents=0, UInt_t skipEvents=0)
Constructor for 1-dimensional kernel PDF with binned interpolation from the sample of points in an NT...
Definition: BinnedKernelDensity.cpp:31
Logger::print
void print(int level, const char *format,...)
Definition: Logger.cpp:27
AbsDensity::m_name
char m_name[256]
PDF name.
Definition: AbsDensity.hh:118
BinnedKernelDensity::normalise
void normalise(void)
Normalise the PDF such that the average PDF value over the allowed phase space equals to 1...
Definition: BinnedKernelDensity.cpp:827
BinnedKernelDensity::m_fractionalMode
Bool_t m_fractionalMode
Fractional mode flag.
Definition: BinnedKernelDensity.hh:478
BinnedKernelDensity::m_binning
std::vector< UInt_t > m_binning
Vector of bin numbers in each variable.
Definition: BinnedKernelDensity.hh:463
AbsDensity::m_rnd
TRandom3 m_rnd
Random number generator.
Definition: AbsDensity.hh:127
BinnedKernelDensity::iterToIndex
UInt_t iterToIndex(std::vector< UInt_t > &iter)
Convert an N-dimensional iterator vector into a linear bin index in the bin map.
Definition: BinnedKernelDensity.cpp:392
BinnedKernelDensity::density
Double_t density(std::vector< Double_t > &x)
Calculate PDF density at the point. Multilinear binned interpolation is used which is calculated in t...
Definition: BinnedKernelDensity.cpp:996
AbsPhaseSpace
Abstract class which defines phase space interface.
Definition: AbsPhaseSpace.hh:10
AbsPhaseSpace::dimensionality
virtual UInt_t dimensionality()=0
Get dimensionality of phase space.
AbsDensity::m_maxTries
UInt_t m_maxTries
Maximum number of tries for accept-reject method.
Definition: AbsDensity.hh:124
Logger::setTimer
void setTimer(void)
Reset time counter.
Definition: Logger.cpp:9
AbsPhaseSpace::lowerLimit
virtual Double_t lowerLimit(UInt_t var)=0
Return lower allowed limit of the variable.
BinnedKernelDensity::fillMapFromTree
void fillMapFromTree(TTree *tree, std::vector< TString > &vars, UInt_t maxEvents=0, UInt_t skipEvents=0)
Definition: BinnedKernelDensity.cpp:488
Logger::timer
int timer(int secs)
Return 1 if more than a certain number of seconds have passed since the last call of this function or...
Definition: Logger.cpp:13
BinnedKernelDensity::addToMap
void addToMap(std::vector< Double_t > &map, std::vector< Double_t > &point, Double_t weight=1.)
Add a kernel density of a single data point to the binned map.
Definition: BinnedKernelDensity.cpp:405
BinnedKernelDensity::~BinnedKernelDensity
virtual ~BinnedKernelDensity()
Destructor.
Definition: BinnedKernelDensity.cpp:331
BinnedKernelDensity::init
void init(AbsPhaseSpace *thePhaseSpace, TTree *tree, std::vector< TString > &vars, std::vector< UInt_t > &binning, std::vector< Double_t > &width, AbsDensity *approx=0, UInt_t toyEvents=0, UInt_t maxEvents=0, UInt_t skipEvents=0)
Common initialise method used by all constructors.
Definition: BinnedKernelDensity.cpp:336
BinnedKernelDensity::mapDensity
Double_t mapDensity(std::vector< Double_t > &map, std::vector< Double_t > &x)
Calculate the raw density using the binned map (estimated or approximation) at a given point...
Definition: BinnedKernelDensity.cpp:895
AbsPhaseSpace::withinLimits
virtual Bool_t withinLimits(std::vector< Double_t > &x)=0
Check if the point is within phase space limits.
AbsDensity::density
virtual Double_t density(std::vector< Double_t > &x)=0
Calculate PDF value at the given point.
BinnedKernelDensity::writeToTextFile
void writeToTextFile(const char *fileName)
Write the binned PDF to a text file.
Definition: BinnedKernelDensity.cpp:682
BinnedKernelDensity::writeToFile
void writeToFile(const char *fileName)
Write density map to a file, ROOT or text depending on extension.
Definition: BinnedKernelDensity.cpp:673
BinnedKernelDensity::fillMapFromDensity
void fillMapFromDensity(AbsDensity *density, UInt_t toyEvents=0)
Definition: BinnedKernelDensity.cpp:566
BinnedKernelDensity::phaseSpace
AbsPhaseSpace * phaseSpace()
Return the phase space.
Definition: BinnedKernelDensity.hh:394
BinnedKernelDensity::m_dim
UInt_t m_dim
Cached dimensionality of the phase space.
Definition: BinnedKernelDensity.hh:475
BinnedKernelDensity::m_approxMap
std::vector< Double_t > m_approxMap
Bin map of approximation PDF convolved with the kernel.
Definition: BinnedKernelDensity.hh:460
BinnedKernelDensity::m_approxDensity
AbsDensity * m_approxDensity
Reference to approximation PDF.
Definition: BinnedKernelDensity.hh:469
Generated by   doxygen 1.8.9.1