Logo ROOT  
Reference Guide
RooNDKeysPdf.cxx
Go to the documentation of this file.
1 /*****************************************************************************
2  * Project: RooFit *
3  * Package: RooFitModels *
4  * File: $Id: RooNDKeysPdf.cxx 31258 2009-11-17 22:41:06Z wouter $
5  * Authors: *
6  * Max Baak, CERN, mbaak@cern.ch *
7  * *
8  * Redistribution and use in source and binary forms, *
9  * with or without modification, are permitted according to the terms *
10  * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
11  *****************************************************************************/
12 
13 /** \class RooNDKeysPdf
14  \ingroup Roofit
15 
16 Generic N-dimensional implementation of a kernel estimation p.d.f. This p.d.f. models the distribution
17 of an arbitrary input dataset as a superposition of Gaussian kernels, one for each data point,
18 each contributing 1/N to the total integral of the p.d.f.
19 If the 'adaptive mode' is enabled, the width of the Gaussian is adaptively calculated from the
20 local density of events, i.e. narrow for regions with high event density to preserve details and
21 wide for regions with log event density to promote smoothness. The details of the general algorithm
22 are described in the following paper:
23 Cranmer KS, Kernel Estimation in High-Energy Physics.
24  Computer Physics Communications 136:198-207,2001 - e-Print Archive: hep ex/0011057
25 For multi-dimensional datasets, the kernels are modeled by multidimensional Gaussians. The kernels are
26 constructed such that they reflect the correlation coefficients between the observables
27 in the input dataset.
28 **/
29 
30 #include <iostream>
31 #include <algorithm>
32 #include <string>
33 
34 #include "TMath.h"
35 #include "TMatrixDSymEigen.h"
36 #include "RooNDKeysPdf.h"
37 #include "RooAbsReal.h"
38 #include "RooRealVar.h"
39 #include "RooRandom.h"
40 #include "RooHist.h"
41 #include "RooMsgService.h"
42 #include "RooChangeTracker.h"
43 
44 #include "TError.h"
45 
46 using namespace std;
47 
49 
50 ////////////////////////////////////////////////////////////////////////////////
51 /// Construct N-dimensional kernel estimation p.d.f. in observables 'varList'
52 /// from dataset 'data'. Options can be
53 ///
54 /// - 'a' = Use adaptive kernels (width varies with local event density)
55 /// - 'm' = Mirror data points over observable boundaries. Improves modeling
56 /// behavior at edges for distributions that are not close to zero
57 /// at edge
58 /// - 'd' = Debug flag
59 /// - 'v' = Verbose flag
60 ///
61 /// The parameter rho (default = 1) provides an overall scale factor that can
62 /// be applied to the bandwith calculated for each kernel. The nSigma parameter
63 /// determines the size of the box that is used to search for contributing kernels
64 /// around a given point in observable space. The nSigma parameters is used
65 /// in case of non-adaptive bandwidths and for the 1st non-adaptive pass for
66 /// the calculation of adaptive keys p.d.f.s.
67 ///
68 /// The optional weight arguments allows to specify an observable or function
69 /// expression in observables that specifies the weight of each event.
70 
71 RooNDKeysPdf::RooNDKeysPdf(const char *name, const char *title, const RooArgList &varList, const RooDataSet &data,
72  TString options, Double_t rho, Double_t nSigma, Bool_t rotate, Bool_t sortInput)
73  : RooAbsPdf(name, title), _varList("varList", "List of variables", this),
74  _rhoList("rhoList", "List of rho parameters", this), _data(&data), _options(options), _widthFactor(rho),
75  _nSigma(nSigma), _rotate(rotate), _sortInput(sortInput), _nAdpt(1)
76 {
77  // Constructor
78 
79  for (const auto var : varList) {
80  if (!dynamic_cast<const RooAbsReal*>(var)) {
81  coutE(InputArguments) << "RooNDKeysPdf::ctor(" << GetName() << ") ERROR: variable " << var->GetName()
82  << " is not of type RooAbsReal" << endl ;
83  R__ASSERT(0) ;
84  }
85  _varList.add(*var) ;
86  _varName.push_back(var->GetName());
87  }
88 
89  createPdf();
90 }
91 
92 ////////////////////////////////////////////////////////////////////////////////
93 /// Constructor
94 
95 RooNDKeysPdf::RooNDKeysPdf(const char *name, const char *title, const RooArgList &varList, const TH1 &hist,
96  TString options, Double_t rho, Double_t nSigma, Bool_t rotate, Bool_t sortInput)
97  : RooAbsPdf(name, title), _varList("varList", "List of variables", this),
98  _rhoList("rhoList", "List of rho parameters", this), _ownedData(createDatasetFromHist(varList, hist)), _data(_ownedData.get()),
99  _options(options), _widthFactor(rho), _nSigma(nSigma), _rotate(rotate),
100  _sortInput(sortInput), _nAdpt(1)
101 {
102  for (const auto var : varList) {
103  if (!dynamic_cast<const RooAbsReal *>(var)) {
104  coutE(InputArguments) << "RooNDKeysPdf::ctor(" << GetName() << ") ERROR: variable " << var->GetName()
105  << " is not of type RooAbsReal" << endl;
106  assert(0);
107  }
108  _varList.add(*var);
109  _varName.push_back(var->GetName());
110  }
111 
112  createPdf();
113 }
114 
115 ////////////////////////////////////////////////////////////////////////////////
116 /// Constructor
117 
118 RooNDKeysPdf::RooNDKeysPdf(const char *name, const char *title, const RooArgList &varList, const RooDataSet &data,
119  const TVectorD &rho, TString options, Double_t nSigma, Bool_t rotate, Bool_t sortInput)
120  : RooAbsPdf(name, title), _varList("varList", "List of variables", this),
121  _rhoList("rhoList", "List of rho parameters", this), _data(&data), _options(options), _widthFactor(-1.0),
122  _nSigma(nSigma), _rotate(rotate), _sortInput(sortInput), _nAdpt(1)
123 {
124  for (const auto var : varList) {
125  if (!dynamic_cast<const RooAbsReal*>(var)) {
126  coutE(InputArguments) << "RooNDKeysPdf::ctor(" << GetName() << ") ERROR: variable " << var->GetName()
127  << " is not of type RooAbsReal" << endl;
128  R__ASSERT(0);
129  }
130  _varList.add(*var) ;
131  _varName.push_back(var->GetName());
132  }
133 
134  // copy rho widths
135  if( _varList.getSize() != rho.GetNrows() ) {
137  << "ERROR: RooNDKeysPdf::RooNDKeysPdf() : The vector-size of rho is different from that of varList."
138  << "Unable to create the PDF." << endl;
139  R__ASSERT(_varList.getSize() == rho.GetNrows());
140  }
141 
142  // negative width factor will serve as a switch in initialize()
143  // negative value means that a vector has been provided as input,
144  // and that _rho has already been set ...
145  _rho.resize( rho.GetNrows() );
146  for (Int_t j = 0; j < rho.GetNrows(); j++) {
147  _rho[j] = rho[j]; /*cout<<"RooNDKeysPdf ctor, _rho["<<j<<"]="<<_rho[j]<<endl;*/
148  }
149 
150  createPdf(); // calls initialize ...
151 }
152 
153 ////////////////////////////////////////////////////////////////////////////////
154 /// Backward compatibility constructor for (1-dim) RooKeysPdf. If you are a new user,
155 /// please use the first constructor form.
156 
157 RooNDKeysPdf::RooNDKeysPdf(const char *name, const char *title, const RooArgList &varList, const RooDataSet &data,
158  const RooArgList &rhoList, TString options, Double_t nSigma, Bool_t rotate, Bool_t sortInput)
159  : RooAbsPdf(name, title), _varList("varList", "List of variables", this),
160  _rhoList("rhoList", "List of rho parameters", this), _data(&data), _options(options), _widthFactor(-1.0),
161  _nSigma(nSigma), _rotate(rotate), _sortInput(sortInput), _nAdpt(1)
162 {
163  for (const auto var : varList) {
164  if (!dynamic_cast<RooAbsReal *>(var)) {
165  coutE(InputArguments) << "RooNDKeysPdf::ctor(" << GetName() << ") ERROR: variable " << var->GetName()
166  << " is not of type RooAbsReal" << endl;
167  assert(0);
168  }
169  _varList.add(*var);
170  _varName.push_back(var->GetName());
171  }
172 
173  _rho.resize(rhoList.getSize(), 1.0);
174 
175  for (unsigned int i=0; i < rhoList.size(); ++i) {
176  const auto rho = rhoList.at(i);
177  if (!dynamic_cast<const RooAbsReal *>(rho)) {
178  coutE(InputArguments) << "RooNDKeysPdf::ctor(" << GetName() << ") ERROR: parameter " << rho->GetName()
179  << " is not of type RooRealVar" << endl;
180  assert(0);
181  }
182  _rhoList.add(*rho);
183  _rho[i] = (dynamic_cast<RooAbsReal *>(rho))->getVal();
184  }
185 
186  // copy rho widths
187  if ((_varList.getSize() != _rhoList.getSize())) {
188  coutE(InputArguments) << "ERROR: RooNDKeysPdf::RooNDKeysPdf() : The size of rhoList is different from varList."
189  << "Unable to create the PDF." << endl;
190  assert(_varList.getSize() == _rhoList.getSize());
191  }
192 
193  // keep track of changes in rho parameters
194  _tracker = new RooChangeTracker("tracker", "track rho parameters", _rhoList, true); // check for value updates
195  (void)_tracker->hasChanged(true); // first evaluation always true for new parameters (?)
196 
197  createPdf();
198 }
199 
200 ////////////////////////////////////////////////////////////////////////////////
201 /// Constructor
202 
203 RooNDKeysPdf::RooNDKeysPdf(const char *name, const char *title, const RooArgList &varList, const TH1 &hist,
204  const RooArgList &rhoList, TString options, Double_t nSigma, Bool_t rotate, Bool_t sortInput)
205  : RooAbsPdf(name, title), _varList("varList", "List of variables", this),
206  _rhoList("rhoList", "List of rho parameters", this), _ownedData(createDatasetFromHist(varList, hist)), _data(_ownedData.get()),
207  _options(options), _widthFactor(-1), _nSigma(nSigma), _rotate(rotate), _sortInput(sortInput),
208  _nAdpt(1)
209 {
210  for (const auto var : varList) {
211  if (!dynamic_cast<RooAbsReal *>(var)) {
212  coutE(InputArguments) << "RooNDKeysPdf::ctor(" << GetName() << ") ERROR: variable " << var->GetName()
213  << " is not of type RooAbsReal" << endl;
214  assert(0);
215  }
216  _varList.add(*var);
217  _varName.push_back(var->GetName());
218  }
219 
220  // copy rho widths
221  _rho.resize(rhoList.getSize(), 1.0);
222 
223  for (unsigned int i=0; i < rhoList.size(); ++i) {
224  const auto rho = rhoList.at(i);
225  if (!dynamic_cast<RooAbsReal *>(rho)) {
226  coutE(InputArguments) << "RooNDKeysPdf::ctor(" << GetName() << ") ERROR: parameter " << rho->GetName()
227  << " is not of type RooRealVar" << endl;
228  assert(0);
229  }
230  _rhoList.add(*rho);
231  _rho[i] = (dynamic_cast<RooAbsReal *>(rho))->getVal();
232  }
233 
234  if ((_varList.getSize() != _rhoList.getSize())) {
235  coutE(InputArguments) << "ERROR: RooNDKeysPdf::RooNDKeysPdf() : The size of rhoList is different from varList."
236  << "Unable to create the PDF." << endl;
237  assert(_varList.getSize() == _rhoList.getSize());
238  }
239 
240  // keep track of changes in rho parameters
241  _tracker = new RooChangeTracker("tracker", "track rho parameters", _rhoList, true); // check for value updates
242  (void)_tracker->hasChanged(true); // first evaluation always true for new parameters (?)
243 
244  createPdf();
245 }
246 
247 ////////////////////////////////////////////////////////////////////////////////
248 /// Constructor
249 
250 RooNDKeysPdf::RooNDKeysPdf(const char *name, const char *title, RooAbsReal &x, const RooDataSet &data, Mirror mirror,
251  Double_t rho, Double_t nSigma, Bool_t rotate, Bool_t sortInput)
252  : RooAbsPdf(name, title), _varList("varList", "List of variables", this),
253  _rhoList("rhoList", "List of rho parameters", this), _data(&data), _options("a"), _widthFactor(rho),
254  _nSigma(nSigma), _rotate(rotate), _sortInput(sortInput), _nAdpt(1)
255 {
256  _varList.add(x);
257  _varName.push_back(x.GetName());
258 
259  if (mirror != NoMirror) {
260  if (mirror != MirrorBoth)
261  coutW(InputArguments) << "RooNDKeysPdf::RooNDKeysPdf() : Warning : asymmetric mirror(s) no longer supported."
262  << endl;
263  _options = "m";
264  }
265 
266  createPdf();
267 }
268 
269 ////////////////////////////////////////////////////////////////////////////////
270 /// Backward compatibility constructor for Roo2DKeysPdf. If you are a new user,
271 /// please use the first constructor form.
272 
273 RooNDKeysPdf::RooNDKeysPdf(const char *name, const char *title, RooAbsReal &x, RooAbsReal &y, const RooDataSet &data,
274  TString options, Double_t rho, Double_t nSigma, Bool_t rotate, Bool_t sortInput)
275  : RooAbsPdf(name, title), _varList("varList", "List of variables", this),
276  _rhoList("rhoList", "List of rho parameters", this), _data(&data), _options(options), _widthFactor(rho),
277  _nSigma(nSigma), _rotate(rotate), _sortInput(sortInput), _nAdpt(1)
278 {
279  _varList.add(RooArgSet(x, y));
280  _varName.push_back(x.GetName());
281  _varName.push_back(y.GetName());
282 
283  createPdf();
284 }
285 
286 ////////////////////////////////////////////////////////////////////////////////
287 /// Constructor
288 
289 RooNDKeysPdf::RooNDKeysPdf(const RooNDKeysPdf &other, const char *name)
290  : RooAbsPdf(other, name), _varList("varList", this, other._varList), _rhoList("rhoList", this, other._rhoList),
291  _ownedData(other._ownedData ? new RooDataSet(*other._ownedData) : nullptr),
292  _data(_ownedData ? _ownedData.get() : other._data), _options(other._options), _widthFactor(other._widthFactor),
293  _nSigma(other._nSigma), _rotate(other._rotate), _sortInput(other._sortInput),
294  _nAdpt(other._nAdpt)
295 {
296  _tracker = (other._tracker != NULL ? new RooChangeTracker(*other._tracker) : NULL);
297  // if (_tracker!=NULL) { _tracker->hasChanged(true); }
298 
299  _fixedShape = other._fixedShape;
300  _mirror = other._mirror;
301  _debug = other._debug;
302  _verbose = other._verbose;
303  _nDim = other._nDim;
304  _nEvents = other._nEvents;
305  _nEventsM = other._nEventsM;
306  _nEventsW = other._nEventsW;
307  _d = other._d;
308  _n = other._n;
309  _dataPts = other._dataPts;
310  _dataPtsR = other._dataPtsR;
311  _weights0 = other._weights0;
312  _weights1 = other._weights1;
314  _sortTVIdcs = other._sortTVIdcs;
315  _varName = other._varName;
316  _rho = other._rho;
317  _x = other._x;
318  _x0 = other._x0;
319  _x1 = other._x1;
320  _x2 = other._x2;
321  _xDatLo = other._xDatLo;
322  _xDatHi = other._xDatHi;
323  _xDatLo3s = other._xDatLo3s;
324  _xDatHi3s = other._xDatHi3s;
325  _mean = other._mean;
326  _sigma = other._sigma;
327 
328  // BoxInfo
329  _netFluxZ = other._netFluxZ;
330  _nEventsBW = other._nEventsBW;
331  _nEventsBMSW = other._nEventsBMSW;
332  _xVarLo = other._xVarLo;
333  _xVarHi = other._xVarHi;
334  _xVarLoM3s = other._xVarLoM3s;
335  _xVarLoP3s = other._xVarLoP3s;
336  _xVarHiM3s = other._xVarHiM3s;
337  _xVarHiP3s = other._xVarHiP3s;
338  _bpsIdcs = other._bpsIdcs;
339  _ibNoSort = other._ibNoSort;
340  _sIdcs = other._sIdcs;
341  _bIdcs = other._bIdcs;
342  _bmsIdcs = other._bmsIdcs;
343 
345  _fullBoxInfo = other._fullBoxInfo;
346 
347  _idx = other._idx;
348  _minWeight = other._minWeight;
349  _maxWeight = other._maxWeight;
350  _wMap = other._wMap;
351 
352  _covMat = new TMatrixDSym(*other._covMat);
353  _corrMat = new TMatrixDSym(*other._corrMat);
354  _rotMat = new TMatrixD(*other._rotMat);
355  _sigmaR = new TVectorD(*other._sigmaR);
356  _dx = new TVectorD(*other._dx);
357  _sigmaAvgR = other._sigmaAvgR;
358 }
359 
360 ////////////////////////////////////////////////////////////////////////////////
361 
363 {
364  if (_covMat) delete _covMat;
365  if (_corrMat) delete _corrMat;
366  if (_rotMat) delete _rotMat;
367  if (_sigmaR) delete _sigmaR;
368  if (_dx) delete _dx;
369  if (_tracker)
370  delete _tracker;
371 
372  // delete all the boxinfos map
373  while ( !_rangeBoxInfo.empty() ) {
374  map<pair<string,int>,BoxInfo*>::iterator iter = _rangeBoxInfo.begin();
375  BoxInfo* box= (*iter).second;
376  _rangeBoxInfo.erase(iter);
377  delete box;
378  }
379 }
380 
381 ////////////////////////////////////////////////////////////////////////////////
382 /// evaluation order of constructor.
383 
385 {
386  if (firstCall) {
387  // set options
388  setOptions();
389  // initialization
390  initialize();
391  }
392 
393 
394  // copy dataset, calculate sigma_i's, determine min and max event weight
395  loadDataSet(firstCall);
396 
397  // mirror dataset around dataset boundaries -- does not depend on event weights
398  if (_mirror) mirrorDataSet();
399 
400  // store indices and weights of events with high enough weights
401  loadWeightSet();
402 
403  // store indices of events in variable boundaries and box shell.
404 //calculateShell(&_fullBoxInfo);
405  // calculate normalization needed in analyticalIntegral()
406 //calculatePreNorm(&_fullBoxInfo);
407 
408  // lookup table for determining which events contribute to a certain coordinate
409  const_cast<RooNDKeysPdf*>(this)->sortDataIndices();
410 
411  // determine static and/or adaptive bandwidth
413 }
414 
415 ////////////////////////////////////////////////////////////////////////////////
416 /// set the configuration
417 
419 {
420  _options.ToLower();
421 
422  if( _options.Contains("a") ) { _weights = &_weights1; }
423  else { _weights = &_weights0; }
424  if( _options.Contains("m") ) _mirror = true;
425  else _mirror = false;
426  if( _options.Contains("d") ) _debug = true;
427  else _debug = false;
428  if( _options.Contains("v") ) { _debug = true; _verbose = true; }
429  else { _debug = false; _verbose = false; }
430 
431  cxcoutD(InputArguments) << "RooNDKeysPdf::setOptions() options = " << _options
432  << "\n\tbandWidthType = " << _options.Contains("a")
433  << "\n\tmirror = " << _mirror
434  << "\n\tdebug = " << _debug
435  << "\n\tverbose = " << _verbose
436  << endl;
437 
438  if (_nSigma<2.0) {
439  coutW(InputArguments) << "RooNDKeysPdf::setOptions() : Warning : nSigma = " << _nSigma << " < 2.0. "
440  << "Calculated normalization could be too large."
441  << endl;
442  }
443 
444  // number of adaptive width iterations. Default is 1.
445  if (_options.Contains("a")) {
446  if (!sscanf(_options.Data(), "%d%*s", &_nAdpt)) {
447  _nAdpt = 1;
448  }
449  }
450 }
451 
452 ////////////////////////////////////////////////////////////////////////////////
453 /// initialization
454 
456 {
457  _nDim = _varList.getSize();
461 
462  _netFluxZ = kFALSE;
463  _nEventsBW = 0;
464  _nEventsBMSW = 0;
465 
466  if(_nDim==0) {
467  coutE(InputArguments) << "ERROR: RooNDKeysPdf::initialize() : The observable list is empty. "
468  << "Unable to begin generating the PDF." << endl;
469  R__ASSERT (_nDim!=0);
470  }
471 
472  if(_nEvents==0) {
473  coutE(InputArguments) << "ERROR: RooNDKeysPdf::initialize() : The input data set is empty. "
474  << "Unable to begin generating the PDF." << endl;
475  R__ASSERT (_nEvents!=0);
476  }
477 
478  _d = static_cast<Double_t>(_nDim);
479 
480  vector<Double_t> dummy(_nDim,0.);
481  _dataPts.resize(_nEvents,dummy);
482  _weights0.resize(_nEvents,dummy);
483 
484  //rdh _rho.resize(_nDim,_widthFactor);
485 
486  if (_widthFactor>0) { _rho.resize(_nDim,_widthFactor); }
487  // else: _rho has been provided as external input
488 
489  _x.resize(_nDim,0.);
490  _x0.resize(_nDim,0.);
491  _x1.resize(_nDim,0.);
492  _x2.resize(_nDim,0.);
493 
494  _mean.resize(_nDim,0.);
495  _sigma.resize(_nDim,0.);
496 
497  _xDatLo.resize(_nDim,0.);
498  _xDatHi.resize(_nDim,0.);
499  _xDatLo3s.resize(_nDim,0.);
500  _xDatHi3s.resize(_nDim,0.);
501 
502  boxInfoInit(&_fullBoxInfo,0,0xFFFF);
503 
504  _minWeight=0;
505  _maxWeight=0;
506  _wMap.clear();
507 
508  _covMat = 0;
509  _corrMat= 0;
510  _rotMat = 0;
511  _sigmaR = 0;
512  _dx = new TVectorD(_nDim); _dx->Zero();
513  _dataPtsR.resize(_nEvents,*_dx);
514 
515  for(unsigned int j=0; j < _varList.size(); ++j) {
516  auto& var = static_cast<const RooRealVar&>(_varList[j]);
517  _xDatLo[j] = var.getMin();
518  _xDatHi[j] = var.getMax();
519  }
520 }
521 
522 ////////////////////////////////////////////////////////////////////////////////
523 /// copy the dataset and calculate some useful variables
524 
526 {
527  // first some initialization
528  _nEventsW=0.;
529 
530  TMatrixD mat(_nDim,_nDim);
531  if (!_covMat) _covMat = new TMatrixDSym(_nDim);
532  if (!_corrMat) _corrMat= new TMatrixDSym(_nDim);
533  if (!_rotMat) _rotMat = new TMatrixD(_nDim,_nDim);
534  if (!_sigmaR) _sigmaR = new TVectorD(_nDim);
535 
536  mat.Zero();
537  _covMat->Zero();
538  _corrMat->Zero();
539  _rotMat->Zero();
540  _sigmaR->Zero();
541 
542  const RooArgSet* values= _data->get();
543  vector<RooRealVar*> dVars(_nDim);
544  for (Int_t j=0; j<_nDim; j++) {
545  dVars[j] = (RooRealVar*)values->find(_varName[j].c_str());
546  _x0[j]=_x1[j]=_x2[j]=0.;
547  }
548 
549  _idx.clear();
550  for (Int_t i=0; i<_nEvents; i++) {
551 
552  _data->get(i); // fills dVars
553  _idx.push_back(i);
554  vector<Double_t>& point = _dataPts[i];
555  TVectorD& pointV = _dataPtsR[i];
556 
557  Double_t myweight = _data->weight(); // default is one?
558  if ( TMath::Abs(myweight)>_maxWeight ) { _maxWeight = TMath::Abs(myweight); }
559  _nEventsW += myweight;
560 
561  for (Int_t j=0; j<_nDim; j++) {
562  for (Int_t k=0; k<_nDim; k++) {
563  mat(j,k) += dVars[j]->getVal() * dVars[k]->getVal() * myweight;
564  }
565  // only need to do once
566  if (firstCall)
567  point[j] = pointV[j] = dVars[j]->getVal();
568 
569  _x0[j] += 1. * myweight;
570  _x1[j] += point[j] * myweight ;
571  _x2[j] += point[j] * point[j] * myweight ;
572  if (_x2[j]!=_x2[j]) exit(3);
573 
574  // only need to do once
575  if (firstCall) {
576  if (point[j]<_xDatLo[j]) { _xDatLo[j]=point[j]; }
577  if (point[j]>_xDatHi[j]) { _xDatHi[j]=point[j]; }
578  }
579  }
580  }
581 
582  _n = TMath::Power(4./(_nEventsW*(_d+2.)), 1./(_d+4.)) ;
583  // = (4/[n(dim(R) + 2)])^1/(dim(R)+4); dim(R) = 2
584  _minWeight = (0.5 - TMath::Erf(_nSigma/sqrt(2.))/2.) * _maxWeight;
585 
586  for (Int_t j=0; j<_nDim; j++) {
587  _mean[j] = _x1[j]/_x0[j];
588  _sigma[j] = sqrt(_x2[j]/_x0[j]-_mean[j]*_mean[j]);
589  }
590 
591  for (Int_t j=0; j<_nDim; j++) {
592  for (Int_t k=0; k<_nDim; k++) {
593  (*_covMat)(j,k) = mat(j,k)/_x0[j] - _mean[j]*_mean[k];
594  }
595  }
596 
597  for (Int_t j=0; j<_nDim; j++) {
598  for (Int_t k=0; k<_nDim; k++)
599  (*_corrMat)(j,k) = (*_covMat)(j,k)/(_sigma[j]*_sigma[k]) ;
600  }
601 
602  // use raw sigmas (without rho) for sigmaAvgR
603  TMatrixDSymEigen evCalculator(*_covMat);
604  TVectorD sigmaRraw = evCalculator.GetEigenValues();
605  for (Int_t j=0; j<_nDim; j++) { sigmaRraw[j] = sqrt(sigmaRraw[j]); }
606 
607  _sigmaAvgR=1.;
608  for (Int_t j=0; j<_nDim; j++) { _sigmaAvgR *= sigmaRraw[j]; }
610 
611  // find decorrelation matrix and eigenvalues (R)
612  if (_nDim > 1 && _rotate) {
613  // new: rotation matrix now independent of rho evaluation
614  *_rotMat = evCalculator.GetEigenVectors();
615  *_rotMat = _rotMat->T(); // transpose
616  } else {
617  TMatrixD haar(_nDim, _nDim);
618  TMatrixD unit(TMatrixD::kUnit, haar);
619  *_rotMat = unit;
620  }
621 
622  // update sigmas (rho dependent)
623  updateRho();
624 
625  //// rho no longer used after this.
626  //// Now set rho = 1 because sigmaR now contains rho
627  // for (Int_t j=0; j<_nDim; j++) { _rho[j] = 1.; } // reset: important!
628 
629  if (_verbose) {
630  //_covMat->Print();
631  _rotMat->Print();
632  _corrMat->Print();
633  _sigmaR->Print();
634  }
635 
636  if (_nDim > 1 && _rotate) {
637  // apply rotation
638  for (Int_t i = 0; i < _nEvents; i++) {
639  TVectorD &pointR = _dataPtsR[i];
640  pointR *= *_rotMat;
641  }
642  }
643 
644  coutI(Contents) << "RooNDKeysPdf::loadDataSet(" << this << ")"
645  << "\n Number of events in dataset: " << _nEvents
646  << "\n Weighted number of events in dataset: " << _nEventsW << endl;
647 }
648 
649 ////////////////////////////////////////////////////////////////////////////////
650 /// determine mirror dataset.
651 /// mirror points are added around the physical boundaries of the dataset
652 /// Two steps:
653 /// 1. For each entry, determine if it should be mirrored (the mirror configuration).
654 /// 2. For each mirror configuration, make the mirror points.
655 
657 {
658  for (Int_t j=0; j<_nDim; j++) {
659  _xDatLo3s[j] = _xDatLo[j] + _nSigma * (_n * _sigma[j]);
660  _xDatHi3s[j] = _xDatHi[j] - _nSigma * (_n * _sigma[j]);
661 
662  // cout<<"xDatLo3s["<<j<<"]="<<_xDatLo3s[j]<<endl;
663  // cout<<"xDatHi3s["<<j<<"]="<<_xDatHi3s[j]<<endl;
664  }
665 
666  vector<Double_t> dummy(_nDim,0.);
667 
668  // 1.
669  for (Int_t i=0; i<_nEvents; i++) {
670  vector<Double_t>& x = _dataPts[i];
671 
672  Int_t size = 1;
673  vector<vector<Double_t> > mpoints(size,dummy);
674  vector<vector<Int_t> > mjdcs(size);
675 
676  // create all mirror configurations for event i
677  for (Int_t j=0; j<_nDim; j++) {
678 
679  vector<Int_t>& mjdxK = mjdcs[0];
680  vector<Double_t>& mpointK = mpoints[0];
681 
682  // single mirror *at physical boundaries*
683  if ((x[j]>_xDatLo[j] && x[j]<_xDatLo3s[j]) && x[j]<(_xDatLo[j]+_xDatHi[j])/2.) {
684  mpointK[j] = 2.*_xDatLo[j]-x[j];
685  mjdxK.push_back(j);
686  } else if ((x[j]<_xDatHi[j] && x[j]>_xDatHi3s[j]) && x[j]>(_xDatLo[j]+_xDatHi[j])/2.) {
687  mpointK[j] = 2.*_xDatHi[j]-x[j];
688  mjdxK.push_back(j);
689  }
690  }
691 
692  vector<Int_t>& mjdx0 = mjdcs[0];
693  // no mirror point(s) for this event
694  if (size==1 && mjdx0.size()==0) continue;
695 
696  // 2.
697  // generate all mirror points for event i
698  vector<Int_t>& mjdx = mjdcs[0];
699  vector<Double_t>& mpoint = mpoints[0];
700 
701  // number of mirror points for this mirror configuration
702  Int_t eMir = 1 << mjdx.size();
703  vector<vector<Double_t> > epoints(eMir,x);
704 
705  for (Int_t m=0; m<Int_t(mjdx.size()); m++) {
706  Int_t size1 = 1 << m;
707  Int_t size2 = 1 << (m+1);
708  // copy all previous mirror points
709  for (Int_t l=size1; l<size2; ++l) {
710  epoints[l] = epoints[l-size1];
711  // fill high mirror points
712  vector<Double_t>& epoint = epoints[l];
713  epoint[mjdx[Int_t(mjdx.size()-1)-m]] = mpoint[mjdx[Int_t(mjdx.size()-1)-m]];
714  }
715  }
716 
717  // remove duplicate mirror points
718  // note that: first epoint == x
719  epoints.erase(epoints.begin());
720 
721  // add mirror points of event i to total dataset
722  TVectorD pointR(_nDim);
723 
724  for (Int_t m=0; m<Int_t(epoints.size()); m++) {
725  _idx.push_back(i);
726  _dataPts.push_back(epoints[m]);
727  //_weights0.push_back(_weights0[i]);
728  for (Int_t j=0; j<_nDim; j++) { pointR[j] = (epoints[m])[j]; }
729  if (_nDim > 1 && _rotate) {
730  pointR *= *_rotMat;
731  }
732  _dataPtsR.push_back(pointR);
733  }
734 
735  epoints.clear();
736  mpoints.clear();
737  mjdcs.clear();
738  } // end of event loop
739 
740  _nEventsM = Int_t(_dataPts.size());
741 }
742 
743 ////////////////////////////////////////////////////////////////////////////////
744 
746 {
747  _wMap.clear();
748 
749  for (Int_t i=0; i<_nEventsM; i++) {
750  _data->get(_idx[i]);
751  Double_t myweight = _data->weight();
752  //if ( TMath::Abs(myweight)>_minWeight ) {
753  _wMap[i] = myweight;
754  //}
755  }
756 
757  coutI(Contents) << "RooNDKeysPdf::loadWeightSet(" << this << ") : Number of weighted events : " << _wMap.size() << endl;
758 }
759 
760 ////////////////////////////////////////////////////////////////////////////////
761 /// determine points in +/- nSigma shell around the box determined by the variable
762 /// ranges. These points are needed in the normalization, to determine probability
763 /// leakage in and out of the box.
764 
766 {
767  for (Int_t j=0; j<_nDim; j++) {
768  if (bi->xVarLo[j]==_xDatLo[j] && bi->xVarHi[j]==_xDatHi[j]) {
769  bi->netFluxZ = bi->netFluxZ && kTRUE;
770  } else { bi->netFluxZ = kFALSE; }
771 
772  bi->xVarLoM3s[j] = bi->xVarLo[j] - _nSigma * (_n * _sigma[j]);
773  bi->xVarLoP3s[j] = bi->xVarLo[j] + _nSigma * (_n * _sigma[j]);
774  bi->xVarHiM3s[j] = bi->xVarHi[j] - _nSigma * (_n * _sigma[j]);
775  bi->xVarHiP3s[j] = bi->xVarHi[j] + _nSigma * (_n * _sigma[j]);
776 
777  //cout<<"bi->xVarLoM3s["<<j<<"]="<<bi->xVarLoM3s[j]<<endl;
778  //cout<<"bi->xVarLoP3s["<<j<<"]="<<bi->xVarLoP3s[j]<<endl;
779  //cout<<"bi->xVarHiM3s["<<j<<"]="<<bi->xVarHiM3s[j]<<endl;
780  //cout<<"bi->xVarHiM3s["<<j<<"]="<<bi->xVarHiM3s[j]<<endl;
781  }
782 
783  //for (Int_t i=0; i<_nEventsM; i++) {
784  for (const auto& wMapItr : _wMap) {
785  Int_t i = wMapItr.first;
786 
787  const vector<Double_t>& x = _dataPts[i];
788  Bool_t inVarRange(kTRUE);
789  Bool_t inVarRangePlusShell(kTRUE);
790 
791  for (Int_t j=0; j<_nDim; j++) {
792 
793  if (x[j]>bi->xVarLo[j] && x[j]<bi->xVarHi[j]) {
794  inVarRange = inVarRange && kTRUE;
795  } else { inVarRange = inVarRange && kFALSE; }
796 
797  if (x[j]>bi->xVarLoM3s[j] && x[j]<bi->xVarHiP3s[j]) {
798  inVarRangePlusShell = inVarRangePlusShell && kTRUE;
799  } else { inVarRangePlusShell = inVarRangePlusShell && kFALSE; }
800  }
801 
802  // event in range?
803  if (inVarRange) {
804  bi->bIdcs.push_back(i);
805  }
806 
807  // event in shell?
808  if (inVarRangePlusShell) {
809  bi->bpsIdcs[i] = kTRUE;
810  Bool_t inShell(kFALSE);
811  for (Int_t j=0; j<_nDim; j++) {
812  if ((x[j]>bi->xVarLoM3s[j] && x[j]<bi->xVarLoP3s[j]) && x[j]<(bi->xVarLo[j]+bi->xVarHi[j])/2.) {
813  inShell = kTRUE;
814  } else if ((x[j]>bi->xVarHiM3s[j] && x[j]<bi->xVarHiP3s[j]) && x[j]>(bi->xVarLo[j]+bi->xVarHi[j])/2.) {
815  inShell = kTRUE;
816  }
817  }
818  if (inShell) bi->sIdcs.push_back(i); // needed for normalization
819  else {
820  bi->bmsIdcs.push_back(i); // idem
821  }
822  }
823  }
824 
825  coutI(Contents) << "RooNDKeysPdf::calculateShell() : "
826  << "\n Events in shell " << bi->sIdcs.size()
827  << "\n Events in box " << bi->bIdcs.size()
828  << "\n Events in box and shell " << bi->bpsIdcs.size()
829  << endl;
830 }
831 
832 ////////////////////////////////////////////////////////////////////////////////
833 ///bi->nEventsBMSW=0.;
834 ///bi->nEventsBW=0.;
835 
837 {
838  // box minus shell
839  for (Int_t i=0; i<Int_t(bi->bmsIdcs.size()); i++)
840  bi->nEventsBMSW += _wMap.at(bi->bmsIdcs[i]);
841 
842  // box
843  for (Int_t i=0; i<Int_t(bi->bIdcs.size()); i++)
844  bi->nEventsBW += _wMap.at(bi->bIdcs[i]);
845 
846  cxcoutD(Eval) << "RooNDKeysPdf::calculatePreNorm() : "
847  << "\n nEventsBMSW " << bi->nEventsBMSW
848  << "\n nEventsBW " << bi->nEventsBW
849  << endl;
850 }
851 
852 ////////////////////////////////////////////////////////////////////////////////
853 /// sort entries, as needed for loopRange()
854 
856 {
857  // will loop over all events by default
858  if (!_sortInput) {
859  _ibNoSort.clear();
860  for (unsigned int i = 0; i < _dataPtsR.size(); ++i) {
861  _ibNoSort[i] = kTRUE;
862  }
863  return;
864  }
865 
866  itVec itrVecR;
867  vector<TVectorD>::iterator dpRItr = _dataPtsR.begin();
868  for (Int_t i = 0; dpRItr != _dataPtsR.end(); ++dpRItr, ++i) {
869  if (bi) {
870  if (bi->bpsIdcs.find(i) != bi->bpsIdcs.end())
871  // if (_wMap.find(i)!=_wMap.end())
872  itrVecR.push_back(itPair(i, dpRItr));
873  } else
874  itrVecR.push_back(itPair(i, dpRItr));
875  }
876 
877  _sortTVIdcs.resize(_nDim);
878  for (Int_t j=0; j<_nDim; j++) {
879  _sortTVIdcs[j].clear();
880  sort(itrVecR.begin(), itrVecR.end(), [=](const itPair& a, const itPair& b) {
881  return (*a.second)[j] < (*b.second)[j];
882  });
883  _sortTVIdcs[j] = itrVecR;
884  }
885 
886  for (Int_t j=0; j<_nDim; j++) {
887  cxcoutD(Eval) << "RooNDKeysPdf::sortDataIndices() : Number of sorted events : " << _sortTVIdcs[j].size() << endl;
888  }
889 }
890 
891 ////////////////////////////////////////////////////////////////////////////////
892 
894 {
895  cxcoutD(Eval) << "RooNDKeysPdf::calculateBandWidth()" << endl;
896 
897  const bool adaptive = _options.Contains("a");
898  if (_weights != &_weights1 || _weights != &_weights0) {
899  _weights = adaptive ? &_weights1 : &_weights0;
900  }
901 
902  // non-adaptive bandwidth
903  // (default, and needed to calculate adaptive bandwidth)
904 
905  if(!adaptive) {
906  cxcoutD(Eval) << "RooNDKeysPdf::calculateBandWidth() Using static bandwidth." << endl;
907  }
908 
909  // fixed width approximation
910  for (Int_t i=0; i<_nEvents; i++) {
911  vector<Double_t>& weight = _weights0[i];
912  for (Int_t j = 0; j < _nDim; j++) {
913  weight[j] = _n * (*_sigmaR)[j];
914  // cout<<"j: "<<j<<", _n: "<<_n<<", sigmaR="<<(*_sigmaR)[j]<<", weight="<<weight[j]<<endl;
915  }
916  }
917 
918  // adaptive width
919  if (adaptive) {
920  cxcoutD(Eval) << "RooNDKeysPdf::calculateBandWidth() Using adaptive bandwidth." << endl;
921 
922  double sqrt12 = sqrt(12.);
923  double sqrtSigmaAvgR = sqrt(_sigmaAvgR);
924 
925  vector<Double_t> dummy(_nDim, 0.);
926  _weights1.resize(_nEvents, dummy);
927 
928  std::vector<std::vector<Double_t>> *weights_prev(0);
929  std::vector<std::vector<Double_t>> *weights_new(0);
930 
931  // cout << "Number of adaptive iterations: " << _nAdpt << endl;
932 
933  for (Int_t k = 1; k <= _nAdpt; ++k) {
934 
935  // cout << " Cycle: " << k << endl;
936 
937  // if multiple adaptive iterations, need to swap weight sets
938  if (k % 2) {
939  weights_prev = &_weights0;
940  weights_new = &_weights1;
941  } else {
942  weights_prev = &_weights1;
943  weights_new = &_weights0;
944  }
945 
946  for (Int_t i = 0; i < _nEvents; ++i) {
947  vector<Double_t> &x = _dataPts[i];
948  Double_t f = TMath::Power(gauss(x, *weights_prev) / _nEventsW, -1. / (2. * _d));
949 
950  vector<Double_t> &weight = (*weights_new)[i];
951  for (Int_t j = 0; j < _nDim; j++) {
952  Double_t norm = (_n * (*_sigmaR)[j]) / sqrtSigmaAvgR;
953  weight[j] = norm * f / sqrt12; // note additional factor of sqrt(12) compared with HEP-EX/0011057
954  }
955  }
956  }
957  // this is the latest updated weights set
958  _weights = weights_new;
959  }
960 }
961 
962 ////////////////////////////////////////////////////////////////////////////////
963 /// loop over all closest point to x, as determined by loopRange()
964 
965 Double_t RooNDKeysPdf::gauss(vector<Double_t>& x, vector<vector<Double_t> >& weights) const
966 {
967  if(_nEvents==0) return 0.;
968 
969  const double sqrt2pi = std::sqrt(TMath::TwoPi());
970 
971  Double_t z=0.;
972  map<Int_t,Bool_t> ibMap;
973 
974  // determine input loop range for event x
975  if (_sortInput) {
976  if (_sortTVIdcs.empty())
977  const_cast<RooNDKeysPdf*>(this)->sortDataIndices();
978 
979  loopRange(x, ibMap);
980  }
981 
982  for (const auto& ibMapItr : _sortInput ? ibMap : _ibNoSort) {
983  Int_t i = ibMapItr.first;
984 
985  Double_t g(1.);
986 
987  if (i >= (Int_t)_idx.size()) {
988  continue;
989  } //---> 1.myline
990 
991  const vector<Double_t> &point = _dataPts[i];
992  const vector<Double_t> &weight = weights[_idx[i]];
993 
994  for (Int_t j = 0; j < _nDim; j++) {
995  (*_dx)[j] = x[j] - point[j];
996  }
997 
998  if (_nDim > 1 && _rotate) {
999  *_dx *= *_rotMat; // rotate to decorrelated frame!
1000  }
1001 
1002  for (Int_t j = 0; j < _nDim; j++) {
1003  Double_t r = (*_dx)[j]; // x[j] - point[j];
1004  Double_t c = 1. / (2. * weight[j] * weight[j]);
1005 
1006  // cout << "j = " << j << " x[j] = " << point[j] << " w = " << weight[j] << endl;
1007 
1008  g *= exp(-c * r * r);
1009  g *= 1. / (sqrt2pi * weight[j]);
1010  }
1011  z += (g * _wMap.at(_idx[i]));
1012  }
1013  return z;
1014 }
1015 
1016 ////////////////////////////////////////////////////////////////////////////////
1017 /// determine closest points to x, to loop over in evaluate()
1018 
1019 void RooNDKeysPdf::loopRange(vector<Double_t>& x, map<Int_t,Bool_t>& ibMap) const
1020 {
1021  ibMap.clear();
1022  TVectorD xRm(_nDim);
1023  TVectorD xRp(_nDim);
1024 
1025  for (Int_t j = 0; j < _nDim; j++) {
1026  xRm[j] = xRp[j] = x[j];
1027  }
1028 
1029  if (_nDim > 1 && _rotate) {
1030  xRm *= *_rotMat;
1031  xRp *= *_rotMat;
1032  }
1033  for (Int_t j = 0; j < _nDim; j++) {
1034  xRm[j] -= _nSigma * (_n * (*_sigmaR)[j]);
1035  xRp[j] += _nSigma * (_n * (*_sigmaR)[j]);
1036  // cout<<"xRm["<<j<<"]="<<xRm[j]<<endl;
1037  // cout<<"xRp["<<j<<"]="<<xRp[j]<<endl;
1038  }
1039 
1040  vector<TVectorD> xvecRm(1,xRm);
1041  vector<TVectorD> xvecRp(1,xRp);
1042 
1043  map<Int_t,Bool_t> ibMapRT;
1044 
1045  for (Int_t j=0; j<_nDim; j++) {
1046  ibMap.clear();
1047  auto comp = [=](const itPair& a, const itPair& b) {
1048  return (*a.second)[j] < (*b.second)[j];
1049  };
1050 
1051  auto lo = lower_bound(_sortTVIdcs[j].begin(), _sortTVIdcs[j].end(), itPair(0, xvecRm.begin()), comp);
1052  auto hi = upper_bound(_sortTVIdcs[j].begin(), _sortTVIdcs[j].end(), itPair(0, xvecRp.begin()), comp);
1053  auto it = lo;
1054 
1055  if (j==0) {
1056  if (_nDim==1) { for (it=lo; it!=hi; ++it) ibMap[(*it).first] = kTRUE; }
1057  else { for (it=lo; it!=hi; ++it) ibMapRT[(*it).first] = kTRUE; }
1058  continue;
1059  }
1060 
1061  for (it=lo; it!=hi; ++it)
1062  if (ibMapRT.find((*it).first)!=ibMapRT.end()) { ibMap[(*it).first] = kTRUE; }
1063 
1064  ibMapRT.clear();
1065  if (j!=_nDim-1) { ibMapRT = ibMap; }
1066  }
1067 }
1068 
1069 ////////////////////////////////////////////////////////////////////////////////
1070 
1071 void RooNDKeysPdf::boxInfoInit(BoxInfo* bi, const char* rangeName, Int_t /*code*/) const
1072 {
1073  vector<Bool_t> doInt(_nDim,kTRUE);
1074 
1075  bi->filled = kFALSE;
1076 
1077  bi->xVarLo.resize(_nDim,0.);
1078  bi->xVarHi.resize(_nDim,0.);
1079  bi->xVarLoM3s.resize(_nDim,0.);
1080  bi->xVarLoP3s.resize(_nDim,0.);
1081  bi->xVarHiM3s.resize(_nDim,0.);
1082  bi->xVarHiP3s.resize(_nDim,0.);
1083 
1084  bi->netFluxZ = kTRUE;
1085  bi->bpsIdcs.clear();
1086  bi->bIdcs.clear();
1087  bi->sIdcs.clear();
1088  bi->bmsIdcs.clear();
1089 
1090  bi->nEventsBMSW=0.;
1091  bi->nEventsBW=0.;
1092 
1093  for (unsigned int j=0; j < _varList.size(); ++j) {
1094  auto var = static_cast<const RooAbsRealLValue*>(_varList.at(j));
1095  if (doInt[j]) {
1096  bi->xVarLo[j] = var->getMin(rangeName);
1097  bi->xVarHi[j] = var->getMax(rangeName);
1098  } else {
1099  bi->xVarLo[j] = var->getVal() ;
1100  bi->xVarHi[j] = var->getVal() ;
1101  }
1102  }
1103 }
1104 
1105 ////////////////////////////////////////////////////////////////////////////////
1106 
1108 {
1109  if ( (_tracker && _tracker->hasChanged(kTRUE)) || (_weights != &_weights0 && _weights != &_weights1) ) {
1110  updateRho(); // update internal rho parameters
1111  // redetermine static and/or adaptive bandwidth
1112  const_cast<RooNDKeysPdf*>(this)->calculateBandWidth();
1113  }
1114 
1115  const RooArgSet *nset = _varList.nset();
1116  for (unsigned int j=0; j < _varList.size(); ++j) {
1117  auto var = static_cast<const RooAbsReal*>(_varList.at(j));
1118  _x[j] = var->getVal(nset);
1119  }
1120 
1121  Double_t val = gauss(_x,*_weights);
1122  //cout<<"returning "<<val<<endl;
1123 
1124  if (val>=1E-20)
1125  return val ;
1126  else
1127  return (1E-20) ;
1128 }
1129 
1130 ////////////////////////////////////////////////////////////////////////////////
1131 
1132 Int_t RooNDKeysPdf::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const
1133 {
1134  if (rangeName) return 0 ;
1135 
1136  Int_t code=0;
1137  if (matchArgs(allVars,analVars,RooArgSet(_varList))) { code=1; }
1138 
1139  return code;
1140 
1141 }
1142 
1143 ////////////////////////////////////////////////////////////////////////////////
1144 
1145 Double_t RooNDKeysPdf::analyticalIntegral(Int_t code, const char* rangeName) const
1146 {
1147  checkInitWeights();
1148 
1149  cxcoutD(Eval) << "Calling RooNDKeysPdf::analyticalIntegral(" << GetName() << ") with code " << code
1150  << " and rangeName " << (rangeName?rangeName:"<none>") << endl;
1151 
1152  // determine which observables need to be integrated over ...
1153  Int_t nComb = 1 << (_nDim);
1154  R__ASSERT(code>=1 && code<nComb) ;
1155 
1156  vector<Bool_t> doInt(_nDim,kTRUE);
1157 
1158  // get BoxInfo
1159  BoxInfo* bi(0);
1160 
1161  if (rangeName) {
1162  string rangeNameStr(rangeName) ;
1163  bi = _rangeBoxInfo[make_pair(rangeNameStr,code)] ;
1164  if (!bi) {
1165  bi = new BoxInfo ;
1166  _rangeBoxInfo[make_pair(rangeNameStr,code)] = bi ;
1167  boxInfoInit(bi,rangeName,code);
1168  }
1169  } else bi= &_fullBoxInfo ;
1170 
1171  // have boundaries changed?
1172  Bool_t newBounds(kFALSE);
1173  for (unsigned int j=0; j < _varList.size(); ++j) {
1174  auto var = static_cast<const RooAbsRealLValue*>(_varList.at(j));
1175  if ((var->getMin(rangeName)-bi->xVarLo[j]!=0) ||
1176  (var->getMax(rangeName)-bi->xVarHi[j]!=0)) {
1177  newBounds = kTRUE;
1178  }
1179  }
1180 
1181  // reset
1182  if (newBounds) {
1183  cxcoutD(Eval) << "RooNDKeysPdf::analyticalIntegral() : Found new boundaries ... " << (rangeName?rangeName:"<none>") << endl;
1184  boxInfoInit(bi,rangeName,code);
1185  }
1186 
1187  // recalculates netFluxZero and nEventsIR
1188  if (!bi->filled || newBounds) {
1189  // Fill box info with contents
1190  calculateShell(bi);
1191  calculatePreNorm(bi);
1192  bi->filled = kTRUE;
1193  const_cast<RooNDKeysPdf*>(this)->sortDataIndices(bi);
1194  }
1195 
1196  // first guess
1197  Double_t norm=bi->nEventsBW;
1198 
1199  if (_mirror && bi->netFluxZ) {
1200  // KEYS expression is self-normalized
1201  cxcoutD(Eval) << "RooNDKeysPdf::analyticalIntegral() : Using mirrored normalization : " << bi->nEventsBW << endl;
1202  return bi->nEventsBW;
1203  }
1204  // calculate leakage in and out of variable range box
1205  else
1206  {
1207  norm = bi->nEventsBMSW;
1208  if (norm<0.) norm=0.;
1209 
1210  for (Int_t i=0; i<Int_t(bi->sIdcs.size()); ++i) {
1211  Double_t prob=1.;
1212  const vector<Double_t>& x = _dataPts[bi->sIdcs[i]];
1213  const vector<Double_t>& weight = (*_weights)[_idx[bi->sIdcs[i]]];
1214 
1215  vector<Double_t> chi(_nDim,100.);
1216 
1217  for (Int_t j=0; j<_nDim; j++) {
1218  if(!doInt[j]) continue;
1219 
1220  if ((x[j]>bi->xVarLoM3s[j] && x[j]<bi->xVarLoP3s[j]) && x[j]<(bi->xVarLo[j]+bi->xVarHi[j])/2.)
1221  chi[j] = (x[j]-bi->xVarLo[j])/weight[j];
1222  else if ((x[j]>bi->xVarHiM3s[j] && x[j]<bi->xVarHiP3s[j]) && x[j]>(bi->xVarLo[j]+bi->xVarHi[j])/2.)
1223  chi[j] = (bi->xVarHi[j]-x[j])/weight[j];
1224 
1225  if (chi[j]>0) // inVarRange
1226  prob *= (0.5 + TMath::Erf(fabs(chi[j])/sqrt(2.))/2.);
1227  else // outside Var range
1228  prob *= (0.5 - TMath::Erf(fabs(chi[j])/sqrt(2.))/2.);
1229  }
1230 
1231  norm += prob * _wMap.at(_idx[bi->sIdcs[i]]);
1232  }
1233 
1234  cxcoutD(Eval) << "RooNDKeysPdf::analyticalIntegral() : Final normalization : " << norm << " " << bi->nEventsBW << endl;
1235  return norm;
1236  }
1237 }
1238 
1239 
1240 ////////////////////////////////////////////////////////////////////////////////
1241 
1243 {
1244  std::vector<RooRealVar *> varVec;
1245  RooArgSet varsAndWeightSet;
1246 
1247  for (const auto var : varList) {
1248  if (!dynamic_cast<RooRealVar *>(var)) {
1249  coutE(InputArguments) << "RooNDKeysPdf::createDatasetFromHist(" << GetName() << ") WARNING: variable "
1250  << var->GetName() << " is not of type RooRealVar. Skip." << endl;
1251  continue;
1252  }
1253  varsAndWeightSet.add(*var); // used for dataset creation
1254  varVec.push_back(static_cast<RooRealVar *>(var)); // used for setting the variables.
1255  }
1256 
1257  /// Add weight
1258  RooRealVar weight("weight", "event weight", 0);
1259  varsAndWeightSet.add(weight);
1260 
1261  /// determine histogram dimensionality
1262  unsigned int histndim(0);
1263  std::string classname = hist.ClassName();
1264  if (classname.find("TH1") == 0) {
1265  histndim = 1;
1266  } else if (classname.find("TH2") == 0) {
1267  histndim = 2;
1268  } else if (classname.find("TH3") == 0) {
1269  histndim = 3;
1270  }
1271  assert(histndim == varVec.size());
1272 
1273  if (histndim > 3 || histndim <= 0) {
1274  coutE(InputArguments) << "RooNDKeysPdf::createDatasetFromHist(" << GetName()
1275  << ") ERROR: input histogram dimension not between [1-3]: " << histndim << endl;
1276  assert(0);
1277  }
1278 
1279  /// dataset creation
1280  RooDataSet *dataFromHist = new RooDataSet("datasetFromHist", "datasetFromHist", varsAndWeightSet, weight.GetName());
1281 
1282  /// dataset filling
1283  for (int i = 1; i <= hist.GetXaxis()->GetNbins(); ++i) {
1284  // 1 or more dimension
1285 
1286  Double_t xval = hist.GetXaxis()->GetBinCenter(i);
1287  varVec[0]->setVal(xval);
1288 
1289  if (varVec.size() == 1) {
1290  Double_t fval = hist.GetBinContent(i);
1291  weight.setVal(fval);
1292  dataFromHist->add(varsAndWeightSet, fval);
1293  } else { // 2 or more dimensions
1294 
1295  for (int j = 1; j <= hist.GetYaxis()->GetNbins(); ++j) {
1296  Double_t yval = hist.GetYaxis()->GetBinCenter(j);
1297  varVec[1]->setVal(yval);
1298 
1299  if (varVec.size() == 2) {
1300  Double_t fval = hist.GetBinContent(i, j);
1301  weight.setVal(fval);
1302  dataFromHist->add(varsAndWeightSet, fval);
1303  } else { // 3 dimensions
1304 
1305  for (int k = 1; k <= hist.GetZaxis()->GetNbins(); ++k) {
1306  Double_t zval = hist.GetZaxis()->GetBinCenter(k);
1307  varVec[2]->setVal(zval);
1308 
1309  Double_t fval = hist.GetBinContent(i, j, k);
1310  weight.setVal(fval);
1311  dataFromHist->add(varsAndWeightSet, fval);
1312  }
1313  }
1314  }
1315  }
1316  }
1317 
1318  return dataFromHist;
1319 }
1320 
1321 
1322 ////////////////////////////////////////////////////////////////////////////////
1323 /// Return evaluated weights
1325 {
1326  checkInitWeights();
1327 
1328  TMatrixD mref(_nEvents, _nDim + 1);
1329 
1330  cxcoutD(Eval) << "RooNDKeysPdf::getWeights() Return evaluated weights." << endl;
1331 
1332  for (Int_t i = 0; i < _nEvents; ++i) {
1333  const vector<Double_t>& x = _dataPts[i];
1334  for (Int_t j = 0; j < _nDim; j++) {
1335  mref(i, j) = x[j];
1336  }
1337 
1338  const vector<Double_t>& weight = (*_weights)[i];
1339  mref(i, _nDim) = weight[k];
1340  }
1341 
1342  return mref;
1343 }
1344 
1345  ////////////////////////////////////////////////////////////////////////////////
1347 {
1348  for (unsigned int j = 0; j < _rhoList.size(); ++j) {
1349  auto rho = static_cast<const RooAbsReal*>(_rhoList.at(j));
1350  _rho[j] = rho->getVal();
1351  }
1352 
1353  if (_nDim > 1 && _rotate) {
1354  TMatrixDSym covMatRho(_nDim); // covariance matrix times rho parameters
1355  for (Int_t j = 0; j < _nDim; j++) {
1356  for (Int_t k = 0; k < _nDim; k++) {
1357  covMatRho(j, k) = (*_covMat)(j, k) * _rho[j] * _rho[k];
1358  }
1359  }
1360  // find decorrelation matrix and eigenvalues (R)
1361  TMatrixDSymEigen evCalculatorRho(covMatRho);
1362  *_sigmaR = evCalculatorRho.GetEigenValues();
1363  for (Int_t j = 0; j < _nDim; j++) {
1364  (*_sigmaR)[j] = sqrt((*_sigmaR)[j]);
1365  }
1366  } else {
1367  for (Int_t j = 0; j < _nDim; j++) {
1368  (*_sigmaR)[j] = (_sigma[j] * _rho[j]);
1369  } // * rho
1370  }
1371 }
RooNDKeysPdf::createPdf
void createPdf(Bool_t firstCall=kTRUE)
evaluation order of constructor.
Definition: RooNDKeysPdf.cxx:384
c
#define c(i)
Definition: RSha256.hxx:119
TVectorT::GetNrows
Int_t GetNrows() const
Definition: TVectorT.h:81
l
auto * l
Definition: textangle.C:4
m
auto * m
Definition: textangle.C:8
RooNDKeysPdf::BoxInfo::filled
Bool_t filled
Definition: RooNDKeysPdf.h:101
RooRealVar::setVal
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:216
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:91
RooAbsReal.h
RooMsgService.h
f
#define f(i)
Definition: RSha256.hxx:122
RooNDKeysPdf::_d
Double_t _d
Definition: RooNDKeysPdf.h:156
RooNDKeysPdf::calculateBandWidth
void calculateBandWidth()
Definition: RooNDKeysPdf.cxx:893
RooNDKeysPdf::~RooNDKeysPdf
virtual ~RooNDKeysPdf()
Definition: RooNDKeysPdf.cxx:362
RooNDKeysPdf::_bIdcs
std::vector< Int_t > _bIdcs
Definition: RooNDKeysPdf.h:185
RooFit::InputArguments
@ InputArguments
Definition: RooGlobalFunc.h:68
TString::Data
const char * Data() const
Definition: TString.h:369
RooNDKeysPdf::_rotMat
TMatrixD * _rotMat
Definition: RooNDKeysPdf.h:199
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:364
RooNDKeysPdf::initialize
void initialize()
initialization
Definition: RooNDKeysPdf.cxx:455
RooNDKeysPdf::BoxInfo::xVarHiM3s
std::vector< Double_t > xVarHiM3s
Definition: RooNDKeysPdf.h:106
r
ROOT::R::TRInterface & r
Definition: Object.C:4
RooNDKeysPdf::BoxInfo
Definition: RooNDKeysPdf.h:100
RooNDKeysPdf::_nEventsW
Double_t _nEventsW
Definition: RooNDKeysPdf.h:155
RooNDKeysPdf::Mirror
Mirror
Definition: RooNDKeysPdf.h:52
RooNDKeysPdf::_xDatLo3s
std::vector< Double_t > _xDatLo3s
Definition: RooNDKeysPdf.h:175
coutE
#define coutE(a)
Definition: RooMsgService.h:33
coutW
#define coutW(a)
Definition: RooMsgService.h:32
RooArgList
Definition: RooArgList.h:21
RooNDKeysPdf::_weights1
std::vector< std::vector< Double_t > > _weights1
Definition: RooNDKeysPdf.h:163
exp
double exp(double)
RooNDKeysPdf::_xVarLo
std::vector< Double_t > _xVarLo
Definition: RooNDKeysPdf.h:180
Int_t
int Int_t
Definition: RtypesCore.h:45
TAxis::GetBinCenter
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:478
box
void box(Int_t pat, Double_t x1, Double_t y1, Double_t x2, Double_t y2)
Definition: fillpatterns.C:1
RooNDKeysPdf::BoxInfo::xVarLoM3s
std::vector< Double_t > xVarLoM3s
Definition: RooNDKeysPdf.h:106
RooNDKeysPdf::BoxInfo::netFluxZ
Bool_t netFluxZ
Definition: RooNDKeysPdf.h:102
RooAbsCollection::find
RooAbsArg * find(const char *name) const
Find object with given name in list.
Definition: RooAbsCollection.cxx:813
TString::Contains
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:624
x
Double_t x[n]
Definition: legend1.C:17
RooNDKeysPdf::loopRange
void loopRange(std::vector< Double_t > &x, std::map< Int_t, Bool_t > &ibMap) const
determine closest points to x, to loop over in evaluate()
Definition: RooNDKeysPdf.cxx:1019
coutI
#define coutI(a)
Definition: RooMsgService.h:30
TMath::Erf
Double_t Erf(Double_t x)
Computation of the error function erf(x).
Definition: TMath.cxx:184
RooArgSet::add
virtual Bool_t add(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling add() for each element in the source coll...
Definition: RooArgSet.h:88
RooNDKeysPdf::_netFluxZ
Bool_t _netFluxZ
Definition: RooNDKeysPdf.h:177
RooAbsReal
Definition: RooAbsReal.h:61
TMatrixTSym
Definition: TMatrixDSymfwd.h:22
RooArgList::at
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:88
TMath::Abs
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
RooNDKeysPdf::getAnalyticalIntegral
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
Definition: RooNDKeysPdf.cxx:1132
RooNDKeysPdf::boxInfoInit
void boxInfoInit(BoxInfo *bi, const char *rangeName, Int_t code) const
Definition: RooNDKeysPdf.cxx:1071
TVectorT::Zero
TVectorT< Element > & Zero()
Set vector elements to zero.
Definition: TVectorT.cxx:453
TString
Definition: TString.h:136
TMatrixT
Definition: TMatrixDfwd.h:22
RooNDKeysPdf::_rhoList
RooListProxy _rhoList
Definition: RooNDKeysPdf.h:116
RooNDKeysPdf::getWeights
TMatrixD getWeights(const int &k) const
Return evaluated weights.
Definition: RooNDKeysPdf.cxx:1324
RooNDKeysPdf::_sigmaAvgR
Double_t _sigmaAvgR
Definition: RooNDKeysPdf.h:202
RooNDKeysPdf::_nSigma
Double_t _nSigma
Definition: RooNDKeysPdf.h:145
TMatrixT::T
TMatrixT< Element > & T()
Definition: TMatrixT.h:150
RooNDKeysPdf::_nEventsM
Int_t _nEventsM
Definition: RooNDKeysPdf.h:154
TH1::GetZaxis
TAxis * GetZaxis()
Definition: TH1.h:319
b
#define b(i)
Definition: RSha256.hxx:118
RooNDKeysPdf::_covMat
TMatrixDSym * _covMat
Definition: RooNDKeysPdf.h:197
RooNDKeysPdf::_nEvents
Int_t _nEvents
Definition: RooNDKeysPdf.h:153
RooNDKeysPdf::mirrorDataSet
void mirrorDataSet()
determine mirror dataset.
Definition: RooNDKeysPdf.cxx:656
bool
RooNDKeysPdf::_xDatLo
std::vector< Double_t > _xDatLo
Definition: RooNDKeysPdf.h:174
RooDataSet::get
virtual const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
Definition: RooDataSet.cxx:1038
RooNDKeysPdf::BoxInfo::xVarHiP3s
std::vector< Double_t > xVarHiP3s
Definition: RooNDKeysPdf.h:106
RooChangeTracker
Definition: RooChangeTracker.h:26
RooNDKeysPdf::BoxInfo::sIdcs
std::vector< Int_t > sIdcs
Definition: RooNDKeysPdf.h:108
itPair
std::pair< Int_t, VecTVecDouble::iterator > itPair
Definition: RooNDKeysPdf.h:42
RooNDKeysPdf::_dx
TVectorD * _dx
Definition: RooNDKeysPdf.h:201
RooNDKeysPdf::_maxWeight
Double_t _maxWeight
Definition: RooNDKeysPdf.h:194
RooNDKeysPdf::_xVarHiP3s
std::vector< Double_t > _xVarHiP3s
Definition: RooNDKeysPdf.h:181
RooNDKeysPdf::_nEventsBMSW
Double_t _nEventsBMSW
Definition: RooNDKeysPdf.h:179
hi
float type_of_call hi(const int &, const int &)
RooNDKeysPdf::_x1
std::vector< Double_t > _x1
Definition: RooNDKeysPdf.h:172
RooNDKeysPdf::_nDim
Int_t _nDim
Definition: RooNDKeysPdf.h:152
RooNDKeysPdf::_widthFactor
Double_t _widthFactor
Definition: RooNDKeysPdf.h:144
RooDataSet::add
virtual void add(const RooArgSet &row, Double_t weight=1.0, Double_t weightError=0) override
Add a data point, with its coordinates specified in the 'data' argset, to the data set.
Definition: RooDataSet.cxx:1164
RooNDKeysPdf::BoxInfo::xVarLo
std::vector< Double_t > xVarLo
Definition: RooNDKeysPdf.h:105
TMatrixDSymEigen.h
RooNDKeysPdf::_nEventsBW
Double_t _nEventsBW
Definition: RooNDKeysPdf.h:178
ROOT::Math::fabs
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Definition: UnaryOperators.h:131
RooNDKeysPdf::_verbose
Bool_t _verbose
Definition: RooNDKeysPdf.h:150
RooNDKeysPdf::_fullBoxInfo
BoxInfo _fullBoxInfo
Definition: RooNDKeysPdf.h:190
RooNDKeysPdf::calculatePreNorm
void calculatePreNorm(BoxInfo *bi) const
bi->nEventsBMSW=0.
Definition: RooNDKeysPdf.cxx:836
RooNDKeysPdf::_xVarLoM3s
std::vector< Double_t > _xVarLoM3s
Definition: RooNDKeysPdf.h:181
itVec
std::vector< itPair > itVec
Definition: RooNDKeysPdf.h:43
TH1::GetBinContent
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4906
RooNDKeysPdf::_debug
Bool_t _debug
Definition: RooNDKeysPdf.h:149
RooNDKeysPdf::_ibNoSort
std::map< Int_t, Bool_t > _ibNoSort
Definition: RooNDKeysPdf.h:183
RooNDKeysPdf::BoxInfo::xVarHi
std::vector< Double_t > xVarHi
Definition: RooNDKeysPdf.h:105
TH1::GetYaxis
TAxis * GetYaxis()
Definition: TH1.h:318
RooNDKeysPdf::BoxInfo::bpsIdcs
std::map< Int_t, Bool_t > bpsIdcs
Definition: RooNDKeysPdf.h:107
RooNDKeysPdf::_rangeBoxInfo
std::map< std::pair< std::string, int >, BoxInfo * > _rangeBoxInfo
Definition: RooNDKeysPdf.h:189
RooNDKeysPdf::_bpsIdcs
std::map< Int_t, Bool_t > _bpsIdcs
Definition: RooNDKeysPdf.h:182
RooNDKeysPdf::_xDatHi
std::vector< Double_t > _xDatHi
Definition: RooNDKeysPdf.h:174
RooNDKeysPdf::analyticalIntegral
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Definition: RooNDKeysPdf.cxx:1145
TMatrixDSymEigen::GetEigenVectors
const TMatrixD & GetEigenVectors() const
Definition: TMatrixDSymEigen.h:59
RooNDKeysPdf::loadWeightSet
void loadWeightSet()
Definition: RooNDKeysPdf.cxx:745
a
auto * a
Definition: textangle.C:12
RooNDKeysPdf::RooNDKeysPdf
RooNDKeysPdf()=default
RooNDKeysPdf::BoxInfo::nEventsBMSW
Double_t nEventsBMSW
Definition: RooNDKeysPdf.h:104
RooNDKeysPdf::loadDataSet
void loadDataSet(Bool_t firstCall)
copy the dataset and calculate some useful variables
Definition: RooNDKeysPdf.cxx:525
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
RooRandom.h
RooNDKeysPdf::_weights
std::vector< std::vector< Double_t > > * _weights
Definition: RooNDKeysPdf.h:164
RooNDKeysPdf::_xDatHi3s
std::vector< Double_t > _xDatHi3s
Definition: RooNDKeysPdf.h:175
TMath::Power
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:735
RooNDKeysPdf::_tracker
RooChangeTracker * _tracker
Definition: RooNDKeysPdf.h:208
RooNDKeysPdf::BoxInfo::xVarLoP3s
std::vector< Double_t > xVarLoP3s
Definition: RooNDKeysPdf.h:106
RooNDKeysPdf::_sigmaR
TVectorD * _sigmaR
Definition: RooNDKeysPdf.h:200
RooNDKeysPdf::_wMap
std::map< Int_t, Double_t > _wMap
Definition: RooNDKeysPdf.h:195
RooNDKeysPdf::_xVarHiM3s
std::vector< Double_t > _xVarHiM3s
Definition: RooNDKeysPdf.h:181
RooAbsArg::RooArgSet
friend class RooArgSet
Definition: RooAbsArg.h:579
RooNDKeysPdf::_xVarHi
std::vector< Double_t > _xVarHi
Definition: RooNDKeysPdf.h:180
RooNDKeysPdf::_idx
std::vector< Int_t > _idx
Definition: RooNDKeysPdf.h:192
RooAbsCollection::size
Storage_t::size_type size() const
Definition: RooAbsCollection.h:154
RooListProxy::add
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Reimplementation of standard RooArgList::add()
Definition: RooListProxy.cxx:104
RooAbsData::numEntries
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Definition: RooAbsData.cxx:307
RooNDKeysPdf::_weights0
std::vector< std::vector< Double_t > > _weights0
Definition: RooNDKeysPdf.h:162
y
Double_t y[n]
Definition: legend1.C:17
sqrt
double sqrt(double)
RooAbsProxy::nset
const RooArgSet * nset() const
Definition: RooAbsProxy.h:59
RooNDKeysPdf::MirrorBoth
@ MirrorBoth
Definition: RooNDKeysPdf.h:52
RooNDKeysPdf::_nAdpt
Int_t _nAdpt
Definition: RooNDKeysPdf.h:206
RooNDKeysPdf::_x
std::vector< Double_t > _x
Definition: RooNDKeysPdf.h:171
RooChangeTracker::hasChanged
Bool_t hasChanged(Bool_t clearState)
Returns true if state has changed since last call with clearState=kTRUE.
Definition: RooChangeTracker.cxx:120
RooRealVar.h
RooNDKeysPdf::_fixedShape
Bool_t _fixedShape
Definition: RooNDKeysPdf.h:147
RooNDKeysPdf::createDatasetFromHist
RooDataSet * createDatasetFromHist(const RooArgList &varList, const TH1 &hist) const
Definition: RooNDKeysPdf.cxx:1242
RooNDKeysPdf::_dataPtsR
std::vector< TVectorD > _dataPtsR
Definition: RooNDKeysPdf.h:161
void
typedef void((*Func_t)())
RooHist.h
TMath::TwoPi
constexpr Double_t TwoPi()
Definition: TMath.h:50
RooNDKeysPdf::_mirror
Bool_t _mirror
Definition: RooNDKeysPdf.h:148
RooNDKeysPdf::_sIdcs
std::vector< Int_t > _sIdcs
Definition: RooNDKeysPdf.h:184
RooNDKeysPdf::NoMirror
@ NoMirror
Definition: RooNDKeysPdf.h:52
RooNDKeysPdf::checkInitWeights
void checkInitWeights() const
Definition: RooNDKeysPdf.h:135
TVectorT
Definition: TMatrixTBase.h:78
RooNDKeysPdf::_xVarLoP3s
std::vector< Double_t > _xVarLoP3s
Definition: RooNDKeysPdf.h:181
Double_t
double Double_t
Definition: RtypesCore.h:59
RooAbsArg::RooDataSet
friend class RooDataSet
Definition: RooAbsArg.h:627
RooNDKeysPdf::_sigma
std::vector< Double_t > _sigma
Definition: RooNDKeysPdf.h:173
RooNDKeysPdf::_n
Double_t _n
Definition: RooNDKeysPdf.h:157
R__ASSERT
#define R__ASSERT(e)
Definition: TError.h:120
RooChangeTracker.h
RooNDKeysPdf::gauss
Double_t gauss(std::vector< Double_t > &x, std::vector< std::vector< Double_t >> &weights) const
loop over all closest point to x, as determined by loopRange()
Definition: RooNDKeysPdf.cxx:965
RooNDKeysPdf::_bmsIdcs
std::vector< Int_t > _bmsIdcs
Definition: RooNDKeysPdf.h:186
TVectorD
TVectorT< Double_t > TVectorD
Definition: TVectorDfwd.h:22
TMatrixD
TMatrixT< Double_t > TMatrixD
Definition: TMatrixDfwd.h:22
dummy
static RooMathCoreReg dummy
Definition: RooMathCoreReg.cxx:27
RooNDKeysPdf::calculateShell
void calculateShell(BoxInfo *bi) const
determine points in +/- nSigma shell around the box determined by the variable ranges.
Definition: RooNDKeysPdf.cxx:765
RooNDKeysPdf::_data
const RooDataSet * _data
Definition: RooNDKeysPdf.h:142
RooNDKeysPdf::evaluate
Double_t evaluate() const
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
Definition: RooNDKeysPdf.cxx:1107
RooNDKeysPdf::_sortInput
Bool_t _sortInput
Definition: RooNDKeysPdf.h:205
RooNDKeysPdf::BoxInfo::nEventsBW
Double_t nEventsBW
Definition: RooNDKeysPdf.h:103
RooNDKeysPdf::_varList
RooListProxy _varList
Definition: RooNDKeysPdf.h:115
RooNDKeysPdf::_x0
std::vector< Double_t > _x0
Definition: RooNDKeysPdf.h:172
RooNDKeysPdf::_minWeight
Double_t _minWeight
Definition: RooNDKeysPdf.h:193
TH1
Definition: TH1.h:57
RooNDKeysPdf::_sortTVIdcs
std::vector< itVec > _sortTVIdcs
Weights to be used. Points either to _weights0 or _weights1.
Definition: RooNDKeysPdf.h:166
TMatrixDSymEigen::GetEigenValues
const TVectorD & GetEigenValues() const
Definition: TMatrixDSymEigen.h:60
RooNDKeysPdf.h
name
char name[80]
Definition: TGX11.cxx:110
TVectorT::Print
void Print(Option_t *option="") const
Print the vector as a list of elements.
Definition: TVectorT.cxx:1363
RooNDKeysPdf::_x2
std::vector< Double_t > _x2
Definition: RooNDKeysPdf.h:172
RooNDKeysPdf::BoxInfo::bmsIdcs
std::vector< Int_t > bmsIdcs
Definition: RooNDKeysPdf.h:110
RooDataSet
Definition: RooDataSet.h:33
RooNDKeysPdf::_rotate
Bool_t _rotate
Definition: RooNDKeysPdf.h:204
RooNDKeysPdf::_varName
std::vector< std::string > _varName
Definition: RooNDKeysPdf.h:168
RooNDKeysPdf
Definition: RooNDKeysPdf.h:48
RooNDKeysPdf::updateRho
void updateRho() const
Definition: RooNDKeysPdf.cxx:1346
RooAbsReal::matchArgs
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
Definition: RooAbsReal.cxx:3400
RooNDKeysPdf::sortDataIndices
void sortDataIndices(BoxInfo *bi=0)
sort entries, as needed for loopRange()
Definition: RooNDKeysPdf.cxx:855
RooNDKeysPdf::_options
TString _options
do not persist
Definition: RooNDKeysPdf.h:143
RooAbsPdf
Definition: RooAbsPdf.h:40
RooNDKeysPdf::_mean
std::vector< Double_t > _mean
Definition: RooNDKeysPdf.h:173
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:53
TMatrixT::kUnit
@ kUnit
Definition: TMatrixT.h:58
TMatrixDSym
TMatrixTSym< Double_t > TMatrixDSym
Definition: TMatrixDSymfwd.h:22
TString::ToLower
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1125
RooRealVar
Definition: RooRealVar.h:36
TH1::GetXaxis
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:317
RooFit::Eval
@ Eval
Definition: RooGlobalFunc.h:68
RooNDKeysPdf::BoxInfo::bIdcs
std::vector< Int_t > bIdcs
Definition: RooNDKeysPdf.h:109
RooNDKeysPdf::_corrMat
TMatrixDSym * _corrMat
Definition: RooNDKeysPdf.h:198
RooFit::Contents
@ Contents
Definition: RooGlobalFunc.h:69
RooNDKeysPdf::_rho
std::vector< Double_t > _rho
Definition: RooNDKeysPdf.h:169
RooAbsRealLValue
Definition: RooAbsRealLValue.h:31
RooNDKeysPdf::setOptions
void setOptions()
set the configuration
Definition: RooNDKeysPdf.cxx:418
TObject::ClassName
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:130
TMath::E
constexpr Double_t E()
Base of natural log:
Definition: TMath.h:102
TMatrixDSymEigen
Definition: TMatrixDSymEigen.h:27
RooAbsCollection::getSize
Int_t getSize() const
Definition: RooAbsCollection.h:171
RooNDKeysPdf::_dataPts
std::vector< std::vector< Double_t > > _dataPts
Definition: RooNDKeysPdf.h:160
TMath.h
RooArgSet
Definition: RooArgSet.h:28
RooDataSet::weight
virtual Double_t weight() const override
Return event weight of current event.
Definition: RooDataSet.cxx:994
cxcoutD
#define cxcoutD(a)
Definition: RooMsgService.h:81
int
TError.h
g
#define g(i)
Definition: RSha256.hxx:123