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