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