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
16Generic N-dimensional implementation of a kernel estimation p.d.f. This p.d.f. models the distribution
17of an arbitrary input dataset as a superposition of Gaussian kernels, one for each data point,
18each contributing 1/N to the total integral of the p.d.f.
19If the 'adaptive mode' is enabled, the width of the Gaussian is adaptively calculated from the
20local density of events, i.e. narrow for regions with high event density to preserve details and
21wide for regions with log event density to promote smoothness. The details of the general algorithm
22are described in the following paper:
23Cranmer KS, Kernel Estimation in High-Energy Physics.
24 Computer Physics Communications 136:198-207,2001 - e-Print Archive: hep ex/0011057
25For multi-dimensional datasets, the kernels are modeled by multidimensional Gaussians. The kernels are
26constructed such that they reflect the correlation coefficients between the observables
27in 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
46using 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
71RooNDKeysPdf::RooNDKeysPdf(const char *name, const char *title, const RooArgList &varList, const RooDataSet &data,
72 TString options, double rho, double nSigma, bool rotate, bool 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
95RooNDKeysPdf::RooNDKeysPdf(const char *name, const char *title, const RooArgList &varList, const TH1 &hist,
96 TString options, double rho, double nSigma, bool rotate, bool 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
118RooNDKeysPdf::RooNDKeysPdf(const char *name, const char *title, const RooArgList &varList, const RooDataSet &data,
119 const TVectorD &rho, TString options, double nSigma, bool rotate, bool 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;
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
157RooNDKeysPdf::RooNDKeysPdf(const char *name, const char *title, const RooArgList &varList, const RooDataSet &data,
158 const RooArgList &rhoList, TString options, double nSigma, bool rotate, bool 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
203RooNDKeysPdf::RooNDKeysPdf(const char *name, const char *title, const RooArgList &varList, const TH1 &hist,
204 const RooArgList &rhoList, TString options, double nSigma, bool rotate, bool 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
250RooNDKeysPdf::RooNDKeysPdf(const char *name, const char *title, RooAbsReal &x, const RooDataSet &data, Mirror mirror,
251 double rho, double nSigma, bool rotate, bool 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
273RooNDKeysPdf::RooNDKeysPdf(const char *name, const char *title, RooAbsReal &x, RooAbsReal &y, const RooDataSet &data,
274 TString options, double rho, double nSigma, bool rotate, bool 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{
280 _varName.push_back(x.GetName());
281 _varName.push_back(y.GetName());
282
283 createPdf();
284}
285
286////////////////////////////////////////////////////////////////////////////////
287/// Constructor
288
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 != nullptr ? new RooChangeTracker(*other._tracker) : nullptr);
297 // if (_tracker!=nullptr) { _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;
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
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
384void RooNDKeysPdf::createPdf(bool firstCall)
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
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{
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{
460 _fixedShape= false;
461
462 _netFluxZ = false;
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>(_nDim);
479
480 vector<double> 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
525void RooNDKeysPdf::loadDataSet(bool firstCall)
526{
527 // first some initialization
528 _nEventsW=0.;
529
530 TMatrixD mat(_nDim,_nDim);
531 if (!_covMat) _covMat = 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>& point = _dataPts[i];
555 TVectorD& pointV = _dataPtsR[i];
556
557 double 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> dummy(_nDim,0.);
667
668 // 1.
669 for (Int_t i=0; i<_nEvents; i++) {
670 vector<double>& x = _dataPts[i];
671
672 Int_t size = 1;
673 vector<vector<double> > 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>& 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.empty()) continue;
695
696 // 2.
697 // generate all mirror points for event i
698 vector<Int_t>& mjdx = mjdcs[0];
699 vector<double>& mpoint = mpoints[0];
700
701 // number of mirror points for this mirror configuration
702 Int_t eMir = 1 << mjdx.size();
703 vector<vector<double> > 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>& 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 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 && true;
770 } else { bi->netFluxZ = false; }
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>& x = _dataPts[i];
788 bool inVarRange(true);
789 bool inVarRangePlusShell(true);
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 && true;
795 } else { inVarRange = inVarRange && false; }
796
797 if (x[j]>bi->xVarLoM3s[j] && x[j]<bi->xVarHiP3s[j]) {
798 inVarRangePlusShell = inVarRangePlusShell && true;
799 } else { inVarRangePlusShell = inVarRangePlusShell && false; }
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] = true;
810 bool inShell(false);
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 = true;
814 } else if ((x[j]>bi->xVarHiM3s[j] && x[j]<bi->xVarHiP3s[j]) && x[j]>(bi->xVarLo[j]+bi->xVarHi[j])/2.) {
815 inShell = true;
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] = true;
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>& 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> dummy(_nDim, 0.);
926 _weights1.resize(_nEvents, dummy);
927
928 std::vector<std::vector<double>> *weights_prev(0);
929 std::vector<std::vector<double>> *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> &x = _dataPts[i];
948 double f = TMath::Power(gauss(x, *weights_prev) / _nEventsW, -1. / (2. * _d));
949
950 vector<double> &weight = (*weights_new)[i];
951 for (Int_t j = 0; j < _nDim; j++) {
952 double 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
965double RooNDKeysPdf::gauss(vector<double>& x, vector<vector<double> >& weights) const
966{
967 if(_nEvents==0) return 0.;
968
969 const double sqrt2pi = std::sqrt(TMath::TwoPi());
970
971 double z=0.;
972 std::vector<int> indices;
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, indices);
980 } else {
981 indices.reserve(_ibNoSort.size());
982 for (const auto& ibMapItr : _ibNoSort) {
983 indices.push_back(ibMapItr.first);
984 }
985 }
986
987 for (const auto& i : indices) {
988 double g(1.);
989
990 if (i >= (Int_t)_idx.size()) {
991 continue;
992 } //---> 1.myline
993
994 const vector<double> &point = _dataPts[i];
995 const vector<double> &weight = weights[_idx[i]];
996
997 for (Int_t j = 0; j < _nDim; j++) {
998 (*_dx)[j] = x[j] - point[j];
999 }
1000
1001 if (_nDim > 1 && _rotate) {
1002 *_dx *= *_rotMat; // rotate to decorrelated frame!
1003 }
1004
1005 for (Int_t j = 0; j < _nDim; j++) {
1006 double r = (*_dx)[j]; // x[j] - point[j];
1007 double c = 1. / (2. * weight[j] * weight[j]);
1008
1009 // cout << "j = " << j << " x[j] = " << point[j] << " w = " << weight[j] << endl;
1010
1011 g *= exp(-c * r * r);
1012 g *= 1. / (sqrt2pi * weight[j]);
1013 }
1014 z += (g * _wMap.at(_idx[i]));
1015 }
1016 return z;
1017}
1018
1019////////////////////////////////////////////////////////////////////////////////
1020/// determine closest points to x, to loop over in evaluate()
1021
1022void RooNDKeysPdf::loopRange(vector<double>& x, std::vector<Int_t>& ibMap) const
1023{
1024 ibMap.clear();
1025 TVectorD xRm(_nDim);
1026 TVectorD xRp(_nDim);
1027
1028 for (Int_t j = 0; j < _nDim; j++) {
1029 xRm[j] = xRp[j] = x[j];
1030 }
1031
1032 if (_nDim > 1 && _rotate) {
1033 xRm *= *_rotMat;
1034 xRp *= *_rotMat;
1035 }
1036 for (Int_t j = 0; j < _nDim; j++) {
1037 xRm[j] -= _nSigma * (_n * (*_sigmaR)[j]);
1038 xRp[j] += _nSigma * (_n * (*_sigmaR)[j]);
1039 }
1040
1041 std::vector<TVectorD> xvecRm(1,xRm);
1042 std::vector<TVectorD> xvecRp(1,xRp);
1043
1044 // To remember the indices that were selected in the previous dimensions
1045 std::vector<Int_t> ibMapRT;
1046
1047 for (Int_t j=0; j<_nDim; j++) {
1048 ibMap.clear();
1049 auto comp = [=](const itPair& a, const itPair& b) {
1050 return (*a.second)[j] < (*b.second)[j];
1051 };
1052
1053 // Figure out which elements are to be selected in the current dimension
1054 auto lo = std::lower_bound(_sortTVIdcs[j].begin(), _sortTVIdcs[j].end(), itPair(0, xvecRm.begin()), comp);
1055 auto hi = std::upper_bound(_sortTVIdcs[j].begin(), _sortTVIdcs[j].end(), itPair(0, xvecRp.begin()), comp);
1056 auto it = lo;
1057
1058 if (j==0) {
1059 // In the first dimension, we can skip straight to filling the remembered
1060 // indices for the previous dimension, or in the 1-d case, just fill the
1061 // output map.
1062 auto& m = _nDim==1 ? ibMap : ibMapRT;
1063 m.reserve(std::distance(lo, hi));
1064 for (it=lo; it!=hi; ++it) {
1065 m.push_back(it->first);
1066 }
1067 // Sort by indices so we can use binary search
1068 std::sort(m.begin(), m.end());
1069 continue;
1070 }
1071
1072 // Iterate over elements in current dimension and add them to the list if
1073 // they were also selected in the previous dimension.
1074 for (it=lo; it!=hi; ++it) {
1075 // Look up in ibMapRT with binary search
1076 auto found = std::lower_bound(ibMapRT.begin(), ibMapRT.end(), it->first);
1077 if (found != ibMapRT.end() && !(it->first < *found)) {
1078 ibMap.push_back(it->first);
1079 }
1080 }
1081 // Sort by indices so we can use binary search
1082 std::sort(ibMap.begin(), ibMap.end());
1083
1084 ibMapRT.clear();
1085 if (j!=_nDim-1) { ibMapRT = std::move(ibMap); }
1086 }
1087}
1088
1089////////////////////////////////////////////////////////////////////////////////
1090
1091void RooNDKeysPdf::boxInfoInit(BoxInfo* bi, const char* rangeName, Int_t /*code*/) const
1092{
1093 vector<bool> doInt(_nDim,true);
1094
1095 bi->filled = false;
1096
1097 bi->xVarLo.resize(_nDim,0.);
1098 bi->xVarHi.resize(_nDim,0.);
1099 bi->xVarLoM3s.resize(_nDim,0.);
1100 bi->xVarLoP3s.resize(_nDim,0.);
1101 bi->xVarHiM3s.resize(_nDim,0.);
1102 bi->xVarHiP3s.resize(_nDim,0.);
1103
1104 bi->netFluxZ = true;
1105 bi->bpsIdcs.clear();
1106 bi->bIdcs.clear();
1107 bi->sIdcs.clear();
1108 bi->bmsIdcs.clear();
1109
1110 bi->nEventsBMSW=0.;
1111 bi->nEventsBW=0.;
1112
1113 for (unsigned int j=0; j < _varList.size(); ++j) {
1114 auto var = static_cast<const RooAbsRealLValue*>(_varList.at(j));
1115 if (doInt[j]) {
1116 bi->xVarLo[j] = var->getMin(rangeName);
1117 bi->xVarHi[j] = var->getMax(rangeName);
1118 } else {
1119 bi->xVarLo[j] = var->getVal() ;
1120 bi->xVarHi[j] = var->getVal() ;
1121 }
1122 }
1123}
1124
1125////////////////////////////////////////////////////////////////////////////////
1126
1128{
1129 if ( (_tracker && _tracker->hasChanged(true)) || (_weights != &_weights0 && _weights != &_weights1) ) {
1130 updateRho(); // update internal rho parameters
1131 // redetermine static and/or adaptive bandwidth
1132 const_cast<RooNDKeysPdf*>(this)->calculateBandWidth();
1133 }
1134
1135 const RooArgSet *nset = _varList.nset();
1136 for (unsigned int j=0; j < _varList.size(); ++j) {
1137 auto var = static_cast<const RooAbsReal*>(_varList.at(j));
1138 _x[j] = var->getVal(nset);
1139 }
1140
1141 double val = gauss(_x,*_weights);
1142 //cout<<"returning "<<val<<endl;
1143
1144 if (val>=1E-20)
1145 return val ;
1146 else
1147 return (1E-20) ;
1148}
1149
1150////////////////////////////////////////////////////////////////////////////////
1151
1152Int_t RooNDKeysPdf::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const
1153{
1154 if (rangeName) return 0 ;
1155
1156 Int_t code=0;
1157 if (matchArgs(allVars,analVars,RooArgSet(_varList))) { code=1; }
1158
1159 return code;
1160
1161}
1162
1163////////////////////////////////////////////////////////////////////////////////
1164
1165double RooNDKeysPdf::analyticalIntegral(Int_t code, const char* rangeName) const
1166{
1168
1169 cxcoutD(Eval) << "Calling RooNDKeysPdf::analyticalIntegral(" << GetName() << ") with code " << code
1170 << " and rangeName " << (rangeName?rangeName:"<none>") << endl;
1171
1172 // determine which observables need to be integrated over ...
1173 Int_t nComb = 1 << (_nDim);
1174 R__ASSERT(code>=1 && code<nComb) ;
1175
1176 vector<bool> doInt(_nDim,true);
1177
1178 // get BoxInfo
1179 BoxInfo* bi(0);
1180
1181 if (rangeName) {
1182 string rangeNameStr(rangeName) ;
1183 bi = _rangeBoxInfo[make_pair(rangeNameStr,code)] ;
1184 if (!bi) {
1185 bi = new BoxInfo ;
1186 _rangeBoxInfo[make_pair(rangeNameStr,code)] = bi ;
1187 boxInfoInit(bi,rangeName,code);
1188 }
1189 } else bi= &_fullBoxInfo ;
1190
1191 // have boundaries changed?
1192 bool newBounds(false);
1193 for (unsigned int j=0; j < _varList.size(); ++j) {
1194 auto var = static_cast<const RooAbsRealLValue*>(_varList.at(j));
1195 if ((var->getMin(rangeName)-bi->xVarLo[j]!=0) ||
1196 (var->getMax(rangeName)-bi->xVarHi[j]!=0)) {
1197 newBounds = true;
1198 }
1199 }
1200
1201 // reset
1202 if (newBounds) {
1203 cxcoutD(Eval) << "RooNDKeysPdf::analyticalIntegral() : Found new boundaries ... " << (rangeName?rangeName:"<none>") << endl;
1204 boxInfoInit(bi,rangeName,code);
1205 }
1206
1207 // recalculates netFluxZero and nEventsIR
1208 if (!bi->filled || newBounds) {
1209 // Fill box info with contents
1210 calculateShell(bi);
1211 calculatePreNorm(bi);
1212 bi->filled = true;
1213 const_cast<RooNDKeysPdf*>(this)->sortDataIndices(bi);
1214 }
1215
1216 // first guess
1217 double norm=bi->nEventsBW;
1218
1219 if (_mirror && bi->netFluxZ) {
1220 // KEYS expression is self-normalized
1221 cxcoutD(Eval) << "RooNDKeysPdf::analyticalIntegral() : Using mirrored normalization : " << bi->nEventsBW << endl;
1222 return bi->nEventsBW;
1223 }
1224 // calculate leakage in and out of variable range box
1225 else
1226 {
1227 norm = bi->nEventsBMSW;
1228 if (norm<0.) norm=0.;
1229
1230 for (Int_t i=0; i<Int_t(bi->sIdcs.size()); ++i) {
1231 double prob=1.;
1232 const vector<double>& x = _dataPts[bi->sIdcs[i]];
1233 const vector<double>& weight = (*_weights)[_idx[bi->sIdcs[i]]];
1234
1235 vector<double> chi(_nDim,100.);
1236
1237 for (Int_t j=0; j<_nDim; j++) {
1238 if(!doInt[j]) continue;
1239
1240 if ((x[j]>bi->xVarLoM3s[j] && x[j]<bi->xVarLoP3s[j]) && x[j]<(bi->xVarLo[j]+bi->xVarHi[j])/2.)
1241 chi[j] = (x[j]-bi->xVarLo[j])/weight[j];
1242 else if ((x[j]>bi->xVarHiM3s[j] && x[j]<bi->xVarHiP3s[j]) && x[j]>(bi->xVarLo[j]+bi->xVarHi[j])/2.)
1243 chi[j] = (bi->xVarHi[j]-x[j])/weight[j];
1244
1245 if (chi[j]>0) // inVarRange
1246 prob *= (0.5 + TMath::Erf(fabs(chi[j])/sqrt(2.))/2.);
1247 else // outside Var range
1248 prob *= (0.5 - TMath::Erf(fabs(chi[j])/sqrt(2.))/2.);
1249 }
1250
1251 norm += prob * _wMap.at(_idx[bi->sIdcs[i]]);
1252 }
1253
1254 cxcoutD(Eval) << "RooNDKeysPdf::analyticalIntegral() : Final normalization : " << norm << " " << bi->nEventsBW << endl;
1255 return norm;
1256 }
1257}
1258
1259
1260////////////////////////////////////////////////////////////////////////////////
1261
1263{
1264 std::vector<RooRealVar *> varVec;
1265 RooArgSet varsAndWeightSet;
1266
1267 for (const auto var : varList) {
1268 if (!dynamic_cast<RooRealVar *>(var)) {
1269 coutE(InputArguments) << "RooNDKeysPdf::createDatasetFromHist(" << GetName() << ") WARNING: variable "
1270 << var->GetName() << " is not of type RooRealVar. Skip." << endl;
1271 continue;
1272 }
1273 varsAndWeightSet.add(*var); // used for dataset creation
1274 varVec.push_back(static_cast<RooRealVar *>(var)); // used for setting the variables.
1275 }
1276
1277 /// Add weight
1278 RooRealVar weight("weight", "event weight", 0);
1279 varsAndWeightSet.add(weight);
1280
1281 /// determine histogram dimensionality
1282 unsigned int histndim(0);
1283 std::string classname = hist.ClassName();
1284 if (classname.find("TH1") == 0) {
1285 histndim = 1;
1286 } else if (classname.find("TH2") == 0) {
1287 histndim = 2;
1288 } else if (classname.find("TH3") == 0) {
1289 histndim = 3;
1290 }
1291 assert(histndim == varVec.size());
1292
1293 if (histndim > 3 || histndim <= 0) {
1294 coutE(InputArguments) << "RooNDKeysPdf::createDatasetFromHist(" << GetName()
1295 << ") ERROR: input histogram dimension not between [1-3]: " << histndim << endl;
1296 assert(0);
1297 }
1298
1299 /// dataset creation
1300 RooDataSet *dataFromHist = new RooDataSet("datasetFromHist", "datasetFromHist", varsAndWeightSet, weight.GetName());
1301
1302 /// dataset filling
1303 for (int i = 1; i <= hist.GetXaxis()->GetNbins(); ++i) {
1304 // 1 or more dimension
1305
1306 double xval = hist.GetXaxis()->GetBinCenter(i);
1307 varVec[0]->setVal(xval);
1308
1309 if (varVec.size() == 1) {
1310 double fval = hist.GetBinContent(i);
1311 weight.setVal(fval);
1312 dataFromHist->add(varsAndWeightSet, fval);
1313 } else { // 2 or more dimensions
1314
1315 for (int j = 1; j <= hist.GetYaxis()->GetNbins(); ++j) {
1316 double yval = hist.GetYaxis()->GetBinCenter(j);
1317 varVec[1]->setVal(yval);
1318
1319 if (varVec.size() == 2) {
1320 double fval = hist.GetBinContent(i, j);
1321 weight.setVal(fval);
1322 dataFromHist->add(varsAndWeightSet, fval);
1323 } else { // 3 dimensions
1324
1325 for (int k = 1; k <= hist.GetZaxis()->GetNbins(); ++k) {
1326 double zval = hist.GetZaxis()->GetBinCenter(k);
1327 varVec[2]->setVal(zval);
1328
1329 double fval = hist.GetBinContent(i, j, k);
1330 weight.setVal(fval);
1331 dataFromHist->add(varsAndWeightSet, fval);
1332 }
1333 }
1334 }
1335 }
1336 }
1337
1338 return dataFromHist;
1339}
1340
1341
1342////////////////////////////////////////////////////////////////////////////////
1343/// Return evaluated weights
1345{
1347
1348 TMatrixD mref(_nEvents, _nDim + 1);
1349
1350 cxcoutD(Eval) << "RooNDKeysPdf::getWeights() Return evaluated weights." << endl;
1351
1352 for (Int_t i = 0; i < _nEvents; ++i) {
1353 const vector<double>& x = _dataPts[i];
1354 for (Int_t j = 0; j < _nDim; j++) {
1355 mref(i, j) = x[j];
1356 }
1357
1358 const vector<double>& weight = (*_weights)[i];
1359 mref(i, _nDim) = weight[k];
1360 }
1361
1362 return mref;
1363}
1364
1365 ////////////////////////////////////////////////////////////////////////////////
1367{
1368 for (unsigned int j = 0; j < _rhoList.size(); ++j) {
1369 auto rho = static_cast<const RooAbsReal*>(_rhoList.at(j));
1370 _rho[j] = rho->getVal();
1371 }
1372
1373 if (_nDim > 1 && _rotate) {
1374 TMatrixDSym covMatRho(_nDim); // covariance matrix times rho parameters
1375 for (Int_t j = 0; j < _nDim; j++) {
1376 for (Int_t k = 0; k < _nDim; k++) {
1377 covMatRho(j, k) = (*_covMat)(j, k) * _rho[j] * _rho[k];
1378 }
1379 }
1380 // find decorrelation matrix and eigenvalues (R)
1381 TMatrixDSymEigen evCalculatorRho(covMatRho);
1382 *_sigmaR = evCalculatorRho.GetEigenValues();
1383 for (Int_t j = 0; j < _nDim; j++) {
1384 (*_sigmaR)[j] = sqrt((*_sigmaR)[j]);
1385 }
1386 } else {
1387 for (Int_t j = 0; j < _nDim; j++) {
1388 (*_sigmaR)[j] = (_sigma[j] * _rho[j]);
1389 } // * rho
1390 }
1391}
#define f(i)
Definition: RSha256.hxx:104
#define c(i)
Definition: RSha256.hxx:101
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define coutI(a)
Definition: RooMsgService.h:34
#define cxcoutD(a)
Definition: RooMsgService.h:85
#define coutW(a)
Definition: RooMsgService.h:36
#define coutE(a)
Definition: RooMsgService.h:37
std::vector< itPair > itVec
Definition: RooNDKeysPdf.h:43
std::pair< Int_t, VecTVecDouble::iterator > itPair
Definition: RooNDKeysPdf.h:42
int Int_t
Definition: RtypesCore.h:45
#define ClassImp(name)
Definition: Rtypes.h:375
#define R__ASSERT(e)
Definition: TError.h:118
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t b
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t g
char name[80]
Definition: TGX11.cxx:110
#define hi
Definition: THbookFile.cxx:128
TMatrixTSym< Double_t > TMatrixDSym
TMatrixT< Double_t > TMatrixD
Definition: TMatrixDfwd.h:23
TVectorT< Double_t > TVectorD
Definition: TVectorDfwd.h:23
friend class RooDataSet
Definition: RooAbsArg.h:690
friend class RooArgSet
Definition: RooAbsArg.h:645
Int_t getSize() const
Return the number of elements in the collection.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Storage_t::size_type size() const
RooAbsArg * find(const char *name) const
Find object with given name in list.
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Definition: RooAbsData.cxx:379
const RooArgSet * nset() const
Definition: RooAbsProxy.h:45
RooAbsRealLValue is the common abstract base class for objects that represent a real value that may a...
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:64
bool matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:22
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:110
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:57
RooChangeTracker is a meta object that tracks value changes in a given set of RooAbsArgs by registeri...
bool hasChanged(bool clearState)
Returns true if state has changed since last call with clearState=true.
bool add(const RooAbsArg &var, bool valueServer, bool shapeServer, bool silent)
Overloaded RooCollection_t::add() method insert object into set and registers object as server to own...
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:55
const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
void add(const RooArgSet &row, double weight=1.0, double weightError=0) override
Add one ore more rows of data.
double weight() const override
Return event weight of current event.
Definition: RooDataSet.cxx:931
Generic N-dimensional implementation of a kernel estimation p.d.f.
Definition: RooNDKeysPdf.h:48
BoxInfo _fullBoxInfo
Definition: RooNDKeysPdf.h:190
RooListProxy _rhoList
Definition: RooNDKeysPdf.h:116
std::vector< double > _xVarLoM3s
Definition: RooNDKeysPdf.h:181
TVectorD * _dx
Definition: RooNDKeysPdf.h:201
void loadWeightSet()
std::map< Int_t, bool > _ibNoSort
Definition: RooNDKeysPdf.h:183
std::vector< Int_t > _sIdcs
Definition: RooNDKeysPdf.h:184
std::vector< std::vector< double > > _weights0
Definition: RooNDKeysPdf.h:162
void initialize()
initialization
~RooNDKeysPdf() override
void calculatePreNorm(BoxInfo *bi) const
bi->nEventsBMSW=0.; bi->nEventsBW=0.;
std::vector< double > _xDatHi
Definition: RooNDKeysPdf.h:174
void loadDataSet(bool firstCall)
copy the dataset and calculate some useful variables
std::vector< std::vector< double > > * _weights
Definition: RooNDKeysPdf.h:164
TVectorD * _sigmaR
Definition: RooNDKeysPdf.h:200
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
void loopRange(std::vector< double > &x, std::vector< Int_t > &indices) const
determine closest points to x, to loop over in evaluate()
std::vector< double > _xDatLo
Definition: RooNDKeysPdf.h:174
std::vector< double > _x
Definition: RooNDKeysPdf.h:171
std::map< Int_t, double > _wMap
Definition: RooNDKeysPdf.h:195
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
double _maxWeight
Definition: RooNDKeysPdf.h:194
double _nEventsBMSW
Definition: RooNDKeysPdf.h:179
std::vector< double > _xVarHi
Definition: RooNDKeysPdf.h:180
const RooDataSet * _data
Definition: RooNDKeysPdf.h:142
std::vector< Int_t > _bIdcs
Definition: RooNDKeysPdf.h:185
std::vector< TVectorD > _dataPtsR
Definition: RooNDKeysPdf.h:161
void updateRho() const
std::vector< double > _mean
Definition: RooNDKeysPdf.h:173
std::vector< double > _xVarLo
Definition: RooNDKeysPdf.h:180
void calculateShell(BoxInfo *bi) const
determine points in +/- nSigma shell around the box determined by the variable ranges.
void calculateBandWidth()
std::vector< double > _xDatLo3s
Definition: RooNDKeysPdf.h:175
std::vector< double > _x1
Definition: RooNDKeysPdf.h:172
double _widthFactor
Definition: RooNDKeysPdf.h:144
std::vector< std::vector< double > > _dataPts
Definition: RooNDKeysPdf.h:160
double _sigmaAvgR
Definition: RooNDKeysPdf.h:202
std::vector< double > _xVarHiM3s
Definition: RooNDKeysPdf.h:181
std::vector< double > _x0
Definition: RooNDKeysPdf.h:172
std::vector< double > _x2
Definition: RooNDKeysPdf.h:172
double gauss(std::vector< double > &x, std::vector< std::vector< double > > &weights) const
loop over all closest point to x, as determined by loopRange()
std::vector< double > _rho
Definition: RooNDKeysPdf.h:169
RooListProxy _varList
Definition: RooNDKeysPdf.h:115
TMatrixDSym * _corrMat
Definition: RooNDKeysPdf.h:198
std::vector< Int_t > _bmsIdcs
Definition: RooNDKeysPdf.h:186
std::vector< itVec > _sortTVIdcs
Weights to be used. Points either to _weights0 or _weights1.
Definition: RooNDKeysPdf.h:166
void boxInfoInit(BoxInfo *bi, const char *rangeName, Int_t code) const
std::map< std::pair< std::string, int >, BoxInfo * > _rangeBoxInfo
Definition: RooNDKeysPdf.h:189
std::vector< Int_t > _idx
Definition: RooNDKeysPdf.h:192
void sortDataIndices(BoxInfo *bi=0)
sort entries, as needed for loopRange()
std::vector< double > _xVarLoP3s
Definition: RooNDKeysPdf.h:181
std::vector< std::vector< double > > _weights1
Definition: RooNDKeysPdf.h:163
std::vector< double > _xDatHi3s
Definition: RooNDKeysPdf.h:175
TString _options
do not persist
Definition: RooNDKeysPdf.h:143
std::vector< double > _xVarHiP3s
Definition: RooNDKeysPdf.h:181
TMatrixD getWeights(const int &k) const
Return evaluated weights.
std::vector< double > _sigma
Definition: RooNDKeysPdf.h:173
void mirrorDataSet()
determine mirror dataset.
void setOptions()
set the configuration
RooDataSet * createDatasetFromHist(const RooArgList &varList, const TH1 &hist) const
double _nSigma
Definition: RooNDKeysPdf.h:145
RooNDKeysPdf()=default
std::vector< std::string > _varName
Definition: RooNDKeysPdf.h:168
void checkInitWeights() const
Definition: RooNDKeysPdf.h:135
double _nEventsBW
Definition: RooNDKeysPdf.h:178
TMatrixD * _rotMat
Definition: RooNDKeysPdf.h:199
void createPdf(bool firstCall=true)
evaluation order of constructor.
std::map< Int_t, bool > _bpsIdcs
Definition: RooNDKeysPdf.h:182
TMatrixDSym * _covMat
Definition: RooNDKeysPdf.h:197
double analyticalIntegral(Int_t code, const char *rangeName=0) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
double _nEventsW
Definition: RooNDKeysPdf.h:155
double _minWeight
Definition: RooNDKeysPdf.h:193
RooChangeTracker * _tracker
Definition: RooNDKeysPdf.h:208
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:40
void setVal(double value) override
Set value of variable to 'value'.
Definition: RooRealVar.cxx:281
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:479
TH1 is the base class of all histogram classes in ROOT.
Definition: TH1.h:58
TAxis * GetZaxis()
Definition: TH1.h:321
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:319
TAxis * GetYaxis()
Definition: TH1.h:320
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:5035
TMatrixDSymEigen.
const TVectorD & GetEigenValues() const
const TMatrixD & GetEigenVectors() const
virtual TMatrixTBase< Element > & Zero()
Set matrix elements to zero.
void Print(Option_t *name="") const override
Print the matrix as a table of elements.
TMatrixT< Element > & T()
Definition: TMatrixT.h:150
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:130
Basic string class.
Definition: TString.h:136
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1150
const char * Data() const
Definition: TString.h:369
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:624
TVectorT< Element > & Zero()
Set vector elements to zero.
Definition: TVectorT.cxx:453
Int_t GetNrows() const
Definition: TVectorT.h:75
void Print(Option_t *option="") const override
Print the vector as a list of elements.
Definition: TVectorT.cxx:1363
void box(Int_t pat, Double_t x1, Double_t y1, Double_t x2, Double_t y2)
Definition: fillpatterns.C:1
RVec< PromoteType< T > > exp(const RVec< T > &v)
Definition: RVec.hxx:1744
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
void(off) SmallVectorTemplateBase< T
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
@ InputArguments
Definition: RooGlobalFunc.h:64
Double_t Erf(Double_t x)
Computation of the error function erf(x).
Definition: TMath.cxx:190
constexpr Double_t E()
Base of natural log: .
Definition: TMath.h:93
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition: TMath.h:718
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition: TMathBase.h:123
constexpr Double_t TwoPi()
Definition: TMath.h:44
std::vector< double > xVarHiM3s
Definition: RooNDKeysPdf.h:106
std::vector< Int_t > bIdcs
Definition: RooNDKeysPdf.h:109
std::vector< double > xVarHiP3s
Definition: RooNDKeysPdf.h:106
std::vector< double > xVarLo
Definition: RooNDKeysPdf.h:105
std::vector< double > xVarHi
Definition: RooNDKeysPdf.h:105
std::map< Int_t, bool > bpsIdcs
Definition: RooNDKeysPdf.h:107
std::vector< Int_t > sIdcs
Definition: RooNDKeysPdf.h:108
std::vector< double > xVarLoM3s
Definition: RooNDKeysPdf.h:106
std::vector< double > xVarLoP3s
Definition: RooNDKeysPdf.h:106
std::vector< Int_t > bmsIdcs
Definition: RooNDKeysPdf.h:110
auto * m
Definition: textangle.C:8
auto * l
Definition: textangle.C:4
auto * a
Definition: textangle.C:12