Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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_t rho, Double_t nSigma, Bool_t rotate, Bool_t sortInput)
73 : RooAbsPdf(name, title), _varList("varList", "List of variables", this),
74 _rhoList("rhoList", "List of rho parameters", this), _data(&data), _options(options), _widthFactor(rho),
75 _nSigma(nSigma), _rotate(rotate), _sortInput(sortInput), _nAdpt(1)
76{
77 // Constructor
78
79 for (const auto var : varList) {
80 if (!dynamic_cast<const RooAbsReal*>(var)) {
81 coutE(InputArguments) << "RooNDKeysPdf::ctor(" << GetName() << ") ERROR: variable " << var->GetName()
82 << " is not of type RooAbsReal" << endl ;
83 R__ASSERT(0) ;
84 }
85 _varList.add(*var) ;
86 _varName.push_back(var->GetName());
87 }
88
89 createPdf();
90}
91
92////////////////////////////////////////////////////////////////////////////////
93/// Constructor
94
95RooNDKeysPdf::RooNDKeysPdf(const char *name, const char *title, const RooArgList &varList, const TH1 &hist,
96 TString options, Double_t rho, Double_t nSigma, Bool_t rotate, Bool_t sortInput)
97 : RooAbsPdf(name, title), _varList("varList", "List of variables", this),
98 _rhoList("rhoList", "List of rho parameters", this), _ownedData(createDatasetFromHist(varList, hist)), _data(_ownedData.get()),
99 _options(options), _widthFactor(rho), _nSigma(nSigma), _rotate(rotate),
100 _sortInput(sortInput), _nAdpt(1)
101{
102 for (const auto var : varList) {
103 if (!dynamic_cast<const RooAbsReal *>(var)) {
104 coutE(InputArguments) << "RooNDKeysPdf::ctor(" << GetName() << ") ERROR: variable " << var->GetName()
105 << " is not of type RooAbsReal" << endl;
106 assert(0);
107 }
108 _varList.add(*var);
109 _varName.push_back(var->GetName());
110 }
111
112 createPdf();
113}
114
115////////////////////////////////////////////////////////////////////////////////
116/// Constructor
117
118RooNDKeysPdf::RooNDKeysPdf(const char *name, const char *title, const RooArgList &varList, const RooDataSet &data,
119 const TVectorD &rho, TString options, Double_t nSigma, Bool_t rotate, Bool_t sortInput)
120 : RooAbsPdf(name, title), _varList("varList", "List of variables", this),
121 _rhoList("rhoList", "List of rho parameters", this), _data(&data), _options(options), _widthFactor(-1.0),
122 _nSigma(nSigma), _rotate(rotate), _sortInput(sortInput), _nAdpt(1)
123{
124 for (const auto var : varList) {
125 if (!dynamic_cast<const RooAbsReal*>(var)) {
126 coutE(InputArguments) << "RooNDKeysPdf::ctor(" << GetName() << ") ERROR: variable " << var->GetName()
127 << " is not of type RooAbsReal" << endl;
128 R__ASSERT(0);
129 }
130 _varList.add(*var) ;
131 _varName.push_back(var->GetName());
132 }
133
134 // copy rho widths
135 if( _varList.getSize() != rho.GetNrows() ) {
136 coutE(InputArguments)
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_t nSigma, Bool_t rotate, Bool_t sortInput)
159 : RooAbsPdf(name, title), _varList("varList", "List of variables", this),
160 _rhoList("rhoList", "List of rho parameters", this), _data(&data), _options(options), _widthFactor(-1.0),
161 _nSigma(nSigma), _rotate(rotate), _sortInput(sortInput), _nAdpt(1)
162{
163 for (const auto var : varList) {
164 if (!dynamic_cast<RooAbsReal *>(var)) {
165 coutE(InputArguments) << "RooNDKeysPdf::ctor(" << GetName() << ") ERROR: variable " << var->GetName()
166 << " is not of type RooAbsReal" << endl;
167 assert(0);
168 }
169 _varList.add(*var);
170 _varName.push_back(var->GetName());
171 }
172
173 _rho.resize(rhoList.getSize(), 1.0);
174
175 for (unsigned int i=0; i < rhoList.size(); ++i) {
176 const auto rho = rhoList.at(i);
177 if (!dynamic_cast<const RooAbsReal *>(rho)) {
178 coutE(InputArguments) << "RooNDKeysPdf::ctor(" << GetName() << ") ERROR: parameter " << rho->GetName()
179 << " is not of type RooRealVar" << endl;
180 assert(0);
181 }
182 _rhoList.add(*rho);
183 _rho[i] = (dynamic_cast<RooAbsReal *>(rho))->getVal();
184 }
185
186 // copy rho widths
187 if ((_varList.getSize() != _rhoList.getSize())) {
188 coutE(InputArguments) << "ERROR: RooNDKeysPdf::RooNDKeysPdf() : The size of rhoList is different from varList."
189 << "Unable to create the PDF." << endl;
190 assert(_varList.getSize() == _rhoList.getSize());
191 }
192
193 // keep track of changes in rho parameters
194 _tracker = new RooChangeTracker("tracker", "track rho parameters", _rhoList, true); // check for value updates
195 (void)_tracker->hasChanged(true); // first evaluation always true for new parameters (?)
196
197 createPdf();
198}
199
200////////////////////////////////////////////////////////////////////////////////
201/// Constructor
202
203RooNDKeysPdf::RooNDKeysPdf(const char *name, const char *title, const RooArgList &varList, const TH1 &hist,
204 const RooArgList &rhoList, TString options, Double_t nSigma, Bool_t rotate, Bool_t sortInput)
205 : RooAbsPdf(name, title), _varList("varList", "List of variables", this),
206 _rhoList("rhoList", "List of rho parameters", this), _ownedData(createDatasetFromHist(varList, hist)), _data(_ownedData.get()),
207 _options(options), _widthFactor(-1), _nSigma(nSigma), _rotate(rotate), _sortInput(sortInput),
208 _nAdpt(1)
209{
210 for (const auto var : varList) {
211 if (!dynamic_cast<RooAbsReal *>(var)) {
212 coutE(InputArguments) << "RooNDKeysPdf::ctor(" << GetName() << ") ERROR: variable " << var->GetName()
213 << " is not of type RooAbsReal" << endl;
214 assert(0);
215 }
216 _varList.add(*var);
217 _varName.push_back(var->GetName());
218 }
219
220 // copy rho widths
221 _rho.resize(rhoList.getSize(), 1.0);
222
223 for (unsigned int i=0; i < rhoList.size(); ++i) {
224 const auto rho = rhoList.at(i);
225 if (!dynamic_cast<RooAbsReal *>(rho)) {
226 coutE(InputArguments) << "RooNDKeysPdf::ctor(" << GetName() << ") ERROR: parameter " << rho->GetName()
227 << " is not of type RooRealVar" << endl;
228 assert(0);
229 }
230 _rhoList.add(*rho);
231 _rho[i] = (dynamic_cast<RooAbsReal *>(rho))->getVal();
232 }
233
234 if ((_varList.getSize() != _rhoList.getSize())) {
235 coutE(InputArguments) << "ERROR: RooNDKeysPdf::RooNDKeysPdf() : The size of rhoList is different from varList."
236 << "Unable to create the PDF." << endl;
237 assert(_varList.getSize() == _rhoList.getSize());
238 }
239
240 // keep track of changes in rho parameters
241 _tracker = new RooChangeTracker("tracker", "track rho parameters", _rhoList, true); // check for value updates
242 (void)_tracker->hasChanged(true); // first evaluation always true for new parameters (?)
243
244 createPdf();
245}
246
247////////////////////////////////////////////////////////////////////////////////
248/// Constructor
249
250RooNDKeysPdf::RooNDKeysPdf(const char *name, const char *title, RooAbsReal &x, const RooDataSet &data, Mirror mirror,
251 Double_t rho, Double_t nSigma, Bool_t rotate, Bool_t sortInput)
252 : RooAbsPdf(name, title), _varList("varList", "List of variables", this),
253 _rhoList("rhoList", "List of rho parameters", this), _data(&data), _options("a"), _widthFactor(rho),
254 _nSigma(nSigma), _rotate(rotate), _sortInput(sortInput), _nAdpt(1)
255{
256 _varList.add(x);
257 _varName.push_back(x.GetName());
258
259 if (mirror != NoMirror) {
260 if (mirror != MirrorBoth)
261 coutW(InputArguments) << "RooNDKeysPdf::RooNDKeysPdf() : Warning : asymmetric mirror(s) no longer supported."
262 << endl;
263 _options = "m";
264 }
265
266 createPdf();
267}
268
269////////////////////////////////////////////////////////////////////////////////
270/// Backward compatibility constructor for Roo2DKeysPdf. If you are a new user,
271/// please use the first constructor form.
272
273RooNDKeysPdf::RooNDKeysPdf(const char *name, const char *title, RooAbsReal &x, RooAbsReal &y, const RooDataSet &data,
274 TString options, Double_t rho, Double_t nSigma, Bool_t rotate, Bool_t sortInput)
275 : RooAbsPdf(name, title), _varList("varList", "List of variables", this),
276 _rhoList("rhoList", "List of rho parameters", this), _data(&data), _options(options), _widthFactor(rho),
277 _nSigma(nSigma), _rotate(rotate), _sortInput(sortInput), _nAdpt(1)
278{
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 != NULL ? new RooChangeTracker(*other._tracker) : NULL);
297 // if (_tracker!=NULL) { _tracker->hasChanged(true); }
298
299 _fixedShape = other._fixedShape;
300 _mirror = other._mirror;
301 _debug = other._debug;
302 _verbose = other._verbose;
303 _nDim = other._nDim;
304 _nEvents = other._nEvents;
305 _nEventsM = other._nEventsM;
306 _nEventsW = other._nEventsW;
307 _d = other._d;
308 _n = other._n;
309 _dataPts = other._dataPts;
310 _dataPtsR = other._dataPtsR;
311 _weights0 = other._weights0;
312 _weights1 = other._weights1;
314 _sortTVIdcs = other._sortTVIdcs;
315 _varName = other._varName;
316 _rho = other._rho;
317 _x = other._x;
318 _x0 = other._x0;
319 _x1 = other._x1;
320 _x2 = other._x2;
321 _xDatLo = other._xDatLo;
322 _xDatHi = other._xDatHi;
323 _xDatLo3s = other._xDatLo3s;
324 _xDatHi3s = other._xDatHi3s;
325 _mean = other._mean;
326 _sigma = other._sigma;
327
328 // BoxInfo
329 _netFluxZ = other._netFluxZ;
330 _nEventsBW = other._nEventsBW;
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
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{
461
463 _nEventsBW = 0;
464 _nEventsBMSW = 0;
465
466 if(_nDim==0) {
467 coutE(InputArguments) << "ERROR: RooNDKeysPdf::initialize() : The observable list is empty. "
468 << "Unable to begin generating the PDF." << endl;
469 R__ASSERT (_nDim!=0);
470 }
471
472 if(_nEvents==0) {
473 coutE(InputArguments) << "ERROR: RooNDKeysPdf::initialize() : The input data set is empty. "
474 << "Unable to begin generating the PDF." << endl;
475 R__ASSERT (_nEvents!=0);
476 }
477
478 _d = static_cast<Double_t>(_nDim);
479
480 vector<Double_t> dummy(_nDim,0.);
481 _dataPts.resize(_nEvents,dummy);
482 _weights0.resize(_nEvents,dummy);
483
484 //rdh _rho.resize(_nDim,_widthFactor);
485
486 if (_widthFactor>0) { _rho.resize(_nDim,_widthFactor); }
487 // else: _rho has been provided as external input
488
489 _x.resize(_nDim,0.);
490 _x0.resize(_nDim,0.);
491 _x1.resize(_nDim,0.);
492 _x2.resize(_nDim,0.);
493
494 _mean.resize(_nDim,0.);
495 _sigma.resize(_nDim,0.);
496
497 _xDatLo.resize(_nDim,0.);
498 _xDatHi.resize(_nDim,0.);
499 _xDatLo3s.resize(_nDim,0.);
500 _xDatHi3s.resize(_nDim,0.);
501
502 boxInfoInit(&_fullBoxInfo,0,0xFFFF);
503
504 _minWeight=0;
505 _maxWeight=0;
506 _wMap.clear();
507
508 _covMat = 0;
509 _corrMat= 0;
510 _rotMat = 0;
511 _sigmaR = 0;
512 _dx = new TVectorD(_nDim); _dx->Zero();
513 _dataPtsR.resize(_nEvents,*_dx);
514
515 for(unsigned int j=0; j < _varList.size(); ++j) {
516 auto& var = static_cast<const RooRealVar&>(_varList[j]);
517 _xDatLo[j] = var.getMin();
518 _xDatHi[j] = var.getMax();
519 }
520}
521
522////////////////////////////////////////////////////////////////////////////////
523/// copy the dataset and calculate some useful variables
524
526{
527 // first some initialization
528 _nEventsW=0.;
529
530 TMatrixD mat(_nDim,_nDim);
531 if (!_covMat) _covMat = new TMatrixDSym(_nDim);
533 if (!_rotMat) _rotMat = new TMatrixD(_nDim,_nDim);
534 if (!_sigmaR) _sigmaR = new TVectorD(_nDim);
535
536 mat.Zero();
537 _covMat->Zero();
538 _corrMat->Zero();
539 _rotMat->Zero();
540 _sigmaR->Zero();
541
542 const RooArgSet* values= _data->get();
543 vector<RooRealVar*> dVars(_nDim);
544 for (Int_t j=0; j<_nDim; j++) {
545 dVars[j] = (RooRealVar*)values->find(_varName[j].c_str());
546 _x0[j]=_x1[j]=_x2[j]=0.;
547 }
548
549 _idx.clear();
550 for (Int_t i=0; i<_nEvents; i++) {
551
552 _data->get(i); // fills dVars
553 _idx.push_back(i);
554 vector<Double_t>& point = _dataPts[i];
555 TVectorD& pointV = _dataPtsR[i];
556
557 Double_t myweight = _data->weight(); // default is one?
558 if ( TMath::Abs(myweight)>_maxWeight ) { _maxWeight = TMath::Abs(myweight); }
559 _nEventsW += myweight;
560
561 for (Int_t j=0; j<_nDim; j++) {
562 for (Int_t k=0; k<_nDim; k++) {
563 mat(j,k) += dVars[j]->getVal() * dVars[k]->getVal() * myweight;
564 }
565 // only need to do once
566 if (firstCall)
567 point[j] = pointV[j] = dVars[j]->getVal();
568
569 _x0[j] += 1. * myweight;
570 _x1[j] += point[j] * myweight ;
571 _x2[j] += point[j] * point[j] * myweight ;
572 if (_x2[j]!=_x2[j]) exit(3);
573
574 // only need to do once
575 if (firstCall) {
576 if (point[j]<_xDatLo[j]) { _xDatLo[j]=point[j]; }
577 if (point[j]>_xDatHi[j]) { _xDatHi[j]=point[j]; }
578 }
579 }
580 }
581
582 _n = TMath::Power(4./(_nEventsW*(_d+2.)), 1./(_d+4.)) ;
583 // = (4/[n(dim(R) + 2)])^1/(dim(R)+4); dim(R) = 2
584 _minWeight = (0.5 - TMath::Erf(_nSigma/sqrt(2.))/2.) * _maxWeight;
585
586 for (Int_t j=0; j<_nDim; j++) {
587 _mean[j] = _x1[j]/_x0[j];
588 _sigma[j] = sqrt(_x2[j]/_x0[j]-_mean[j]*_mean[j]);
589 }
590
591 for (Int_t j=0; j<_nDim; j++) {
592 for (Int_t k=0; k<_nDim; k++) {
593 (*_covMat)(j,k) = mat(j,k)/_x0[j] - _mean[j]*_mean[k];
594 }
595 }
596
597 for (Int_t j=0; j<_nDim; j++) {
598 for (Int_t k=0; k<_nDim; k++)
599 (*_corrMat)(j,k) = (*_covMat)(j,k)/(_sigma[j]*_sigma[k]) ;
600 }
601
602 // use raw sigmas (without rho) for sigmaAvgR
603 TMatrixDSymEigen evCalculator(*_covMat);
604 TVectorD sigmaRraw = evCalculator.GetEigenValues();
605 for (Int_t j=0; j<_nDim; j++) { sigmaRraw[j] = sqrt(sigmaRraw[j]); }
606
607 _sigmaAvgR=1.;
608 for (Int_t j=0; j<_nDim; j++) { _sigmaAvgR *= sigmaRraw[j]; }
610
611 // find decorrelation matrix and eigenvalues (R)
612 if (_nDim > 1 && _rotate) {
613 // new: rotation matrix now independent of rho evaluation
614 *_rotMat = evCalculator.GetEigenVectors();
615 *_rotMat = _rotMat->T(); // transpose
616 } else {
617 TMatrixD haar(_nDim, _nDim);
618 TMatrixD unit(TMatrixD::kUnit, haar);
619 *_rotMat = unit;
620 }
621
622 // update sigmas (rho dependent)
623 updateRho();
624
625 //// rho no longer used after this.
626 //// Now set rho = 1 because sigmaR now contains rho
627 // for (Int_t j=0; j<_nDim; j++) { _rho[j] = 1.; } // reset: important!
628
629 if (_verbose) {
630 //_covMat->Print();
631 _rotMat->Print();
632 _corrMat->Print();
633 _sigmaR->Print();
634 }
635
636 if (_nDim > 1 && _rotate) {
637 // apply rotation
638 for (Int_t i = 0; i < _nEvents; i++) {
639 TVectorD &pointR = _dataPtsR[i];
640 pointR *= *_rotMat;
641 }
642 }
643
644 coutI(Contents) << "RooNDKeysPdf::loadDataSet(" << this << ")"
645 << "\n Number of events in dataset: " << _nEvents
646 << "\n Weighted number of events in dataset: " << _nEventsW << endl;
647}
648
649////////////////////////////////////////////////////////////////////////////////
650/// determine mirror dataset.
651/// mirror points are added around the physical boundaries of the dataset
652/// Two steps:
653/// 1. For each entry, determine if it should be mirrored (the mirror configuration).
654/// 2. For each mirror configuration, make the mirror points.
655
657{
658 for (Int_t j=0; j<_nDim; j++) {
659 _xDatLo3s[j] = _xDatLo[j] + _nSigma * (_n * _sigma[j]);
660 _xDatHi3s[j] = _xDatHi[j] - _nSigma * (_n * _sigma[j]);
661
662 // cout<<"xDatLo3s["<<j<<"]="<<_xDatLo3s[j]<<endl;
663 // cout<<"xDatHi3s["<<j<<"]="<<_xDatHi3s[j]<<endl;
664 }
665
666 vector<Double_t> dummy(_nDim,0.);
667
668 // 1.
669 for (Int_t i=0; i<_nEvents; i++) {
670 vector<Double_t>& x = _dataPts[i];
671
672 Int_t size = 1;
673 vector<vector<Double_t> > mpoints(size,dummy);
674 vector<vector<Int_t> > mjdcs(size);
675
676 // create all mirror configurations for event i
677 for (Int_t j=0; j<_nDim; j++) {
678
679 vector<Int_t>& mjdxK = mjdcs[0];
680 vector<Double_t>& mpointK = mpoints[0];
681
682 // single mirror *at physical boundaries*
683 if ((x[j]>_xDatLo[j] && x[j]<_xDatLo3s[j]) && x[j]<(_xDatLo[j]+_xDatHi[j])/2.) {
684 mpointK[j] = 2.*_xDatLo[j]-x[j];
685 mjdxK.push_back(j);
686 } else if ((x[j]<_xDatHi[j] && x[j]>_xDatHi3s[j]) && x[j]>(_xDatLo[j]+_xDatHi[j])/2.) {
687 mpointK[j] = 2.*_xDatHi[j]-x[j];
688 mjdxK.push_back(j);
689 }
690 }
691
692 vector<Int_t>& mjdx0 = mjdcs[0];
693 // no mirror point(s) for this event
694 if (size==1 && mjdx0.size()==0) continue;
695
696 // 2.
697 // generate all mirror points for event i
698 vector<Int_t>& mjdx = mjdcs[0];
699 vector<Double_t>& mpoint = mpoints[0];
700
701 // number of mirror points for this mirror configuration
702 Int_t eMir = 1 << mjdx.size();
703 vector<vector<Double_t> > epoints(eMir,x);
704
705 for (Int_t m=0; m<Int_t(mjdx.size()); m++) {
706 Int_t size1 = 1 << m;
707 Int_t size2 = 1 << (m+1);
708 // copy all previous mirror points
709 for (Int_t l=size1; l<size2; ++l) {
710 epoints[l] = epoints[l-size1];
711 // fill high mirror points
712 vector<Double_t>& epoint = epoints[l];
713 epoint[mjdx[Int_t(mjdx.size()-1)-m]] = mpoint[mjdx[Int_t(mjdx.size()-1)-m]];
714 }
715 }
716
717 // remove duplicate mirror points
718 // note that: first epoint == x
719 epoints.erase(epoints.begin());
720
721 // add mirror points of event i to total dataset
722 TVectorD pointR(_nDim);
723
724 for (Int_t m=0; m<Int_t(epoints.size()); m++) {
725 _idx.push_back(i);
726 _dataPts.push_back(epoints[m]);
727 //_weights0.push_back(_weights0[i]);
728 for (Int_t j=0; j<_nDim; j++) { pointR[j] = (epoints[m])[j]; }
729 if (_nDim > 1 && _rotate) {
730 pointR *= *_rotMat;
731 }
732 _dataPtsR.push_back(pointR);
733 }
734
735 epoints.clear();
736 mpoints.clear();
737 mjdcs.clear();
738 } // end of event loop
739
740 _nEventsM = Int_t(_dataPts.size());
741}
742
743////////////////////////////////////////////////////////////////////////////////
744
746{
747 _wMap.clear();
748
749 for (Int_t i=0; i<_nEventsM; i++) {
750 _data->get(_idx[i]);
751 Double_t myweight = _data->weight();
752 //if ( TMath::Abs(myweight)>_minWeight ) {
753 _wMap[i] = myweight;
754 //}
755 }
756
757 coutI(Contents) << "RooNDKeysPdf::loadWeightSet(" << this << ") : Number of weighted events : " << _wMap.size() << endl;
758}
759
760////////////////////////////////////////////////////////////////////////////////
761/// determine points in +/- nSigma shell around the box determined by the variable
762/// ranges. These points are needed in the normalization, to determine probability
763/// leakage in and out of the box.
764
766{
767 for (Int_t j=0; j<_nDim; j++) {
768 if (bi->xVarLo[j]==_xDatLo[j] && bi->xVarHi[j]==_xDatHi[j]) {
769 bi->netFluxZ = bi->netFluxZ && kTRUE;
770 } else { bi->netFluxZ = kFALSE; }
771
772 bi->xVarLoM3s[j] = bi->xVarLo[j] - _nSigma * (_n * _sigma[j]);
773 bi->xVarLoP3s[j] = bi->xVarLo[j] + _nSigma * (_n * _sigma[j]);
774 bi->xVarHiM3s[j] = bi->xVarHi[j] - _nSigma * (_n * _sigma[j]);
775 bi->xVarHiP3s[j] = bi->xVarHi[j] + _nSigma * (_n * _sigma[j]);
776
777 //cout<<"bi->xVarLoM3s["<<j<<"]="<<bi->xVarLoM3s[j]<<endl;
778 //cout<<"bi->xVarLoP3s["<<j<<"]="<<bi->xVarLoP3s[j]<<endl;
779 //cout<<"bi->xVarHiM3s["<<j<<"]="<<bi->xVarHiM3s[j]<<endl;
780 //cout<<"bi->xVarHiM3s["<<j<<"]="<<bi->xVarHiM3s[j]<<endl;
781 }
782
783 //for (Int_t i=0; i<_nEventsM; i++) {
784 for (const auto& wMapItr : _wMap) {
785 Int_t i = wMapItr.first;
786
787 const vector<Double_t>& x = _dataPts[i];
788 Bool_t inVarRange(kTRUE);
789 Bool_t inVarRangePlusShell(kTRUE);
790
791 for (Int_t j=0; j<_nDim; j++) {
792
793 if (x[j]>bi->xVarLo[j] && x[j]<bi->xVarHi[j]) {
794 inVarRange = inVarRange && kTRUE;
795 } else { inVarRange = inVarRange && kFALSE; }
796
797 if (x[j]>bi->xVarLoM3s[j] && x[j]<bi->xVarHiP3s[j]) {
798 inVarRangePlusShell = inVarRangePlusShell && kTRUE;
799 } else { inVarRangePlusShell = inVarRangePlusShell && kFALSE; }
800 }
801
802 // event in range?
803 if (inVarRange) {
804 bi->bIdcs.push_back(i);
805 }
806
807 // event in shell?
808 if (inVarRangePlusShell) {
809 bi->bpsIdcs[i] = kTRUE;
810 Bool_t inShell(kFALSE);
811 for (Int_t j=0; j<_nDim; j++) {
812 if ((x[j]>bi->xVarLoM3s[j] && x[j]<bi->xVarLoP3s[j]) && x[j]<(bi->xVarLo[j]+bi->xVarHi[j])/2.) {
813 inShell = kTRUE;
814 } else if ((x[j]>bi->xVarHiM3s[j] && x[j]<bi->xVarHiP3s[j]) && x[j]>(bi->xVarLo[j]+bi->xVarHi[j])/2.) {
815 inShell = kTRUE;
816 }
817 }
818 if (inShell) bi->sIdcs.push_back(i); // needed for normalization
819 else {
820 bi->bmsIdcs.push_back(i); // idem
821 }
822 }
823 }
824
825 coutI(Contents) << "RooNDKeysPdf::calculateShell() : "
826 << "\n Events in shell " << bi->sIdcs.size()
827 << "\n Events in box " << bi->bIdcs.size()
828 << "\n Events in box and shell " << bi->bpsIdcs.size()
829 << endl;
830}
831
832////////////////////////////////////////////////////////////////////////////////
833///bi->nEventsBMSW=0.;
834///bi->nEventsBW=0.;
835
837{
838 // box minus shell
839 for (Int_t i=0; i<Int_t(bi->bmsIdcs.size()); i++)
840 bi->nEventsBMSW += _wMap.at(bi->bmsIdcs[i]);
841
842 // box
843 for (Int_t i=0; i<Int_t(bi->bIdcs.size()); i++)
844 bi->nEventsBW += _wMap.at(bi->bIdcs[i]);
845
846 cxcoutD(Eval) << "RooNDKeysPdf::calculatePreNorm() : "
847 << "\n nEventsBMSW " << bi->nEventsBMSW
848 << "\n nEventsBW " << bi->nEventsBW
849 << endl;
850}
851
852////////////////////////////////////////////////////////////////////////////////
853/// sort entries, as needed for loopRange()
854
856{
857 // will loop over all events by default
858 if (!_sortInput) {
859 _ibNoSort.clear();
860 for (unsigned int i = 0; i < _dataPtsR.size(); ++i) {
861 _ibNoSort[i] = kTRUE;
862 }
863 return;
864 }
865
866 itVec itrVecR;
867 vector<TVectorD>::iterator dpRItr = _dataPtsR.begin();
868 for (Int_t i = 0; dpRItr != _dataPtsR.end(); ++dpRItr, ++i) {
869 if (bi) {
870 if (bi->bpsIdcs.find(i) != bi->bpsIdcs.end())
871 // if (_wMap.find(i)!=_wMap.end())
872 itrVecR.push_back(itPair(i, dpRItr));
873 } else
874 itrVecR.push_back(itPair(i, dpRItr));
875 }
876
877 _sortTVIdcs.resize(_nDim);
878 for (Int_t j=0; j<_nDim; j++) {
879 _sortTVIdcs[j].clear();
880 sort(itrVecR.begin(), itrVecR.end(), [=](const itPair& a, const itPair& b) {
881 return (*a.second)[j] < (*b.second)[j];
882 });
883 _sortTVIdcs[j] = itrVecR;
884 }
885
886 for (Int_t j=0; j<_nDim; j++) {
887 cxcoutD(Eval) << "RooNDKeysPdf::sortDataIndices() : Number of sorted events : " << _sortTVIdcs[j].size() << endl;
888 }
889}
890
891////////////////////////////////////////////////////////////////////////////////
892
894{
895 cxcoutD(Eval) << "RooNDKeysPdf::calculateBandWidth()" << endl;
896
897 const bool adaptive = _options.Contains("a");
898 if (_weights != &_weights1 || _weights != &_weights0) {
899 _weights = adaptive ? &_weights1 : &_weights0;
900 }
901
902 // non-adaptive bandwidth
903 // (default, and needed to calculate adaptive bandwidth)
904
905 if(!adaptive) {
906 cxcoutD(Eval) << "RooNDKeysPdf::calculateBandWidth() Using static bandwidth." << endl;
907 }
908
909 // fixed width approximation
910 for (Int_t i=0; i<_nEvents; i++) {
911 vector<Double_t>& weight = _weights0[i];
912 for (Int_t j = 0; j < _nDim; j++) {
913 weight[j] = _n * (*_sigmaR)[j];
914 // cout<<"j: "<<j<<", _n: "<<_n<<", sigmaR="<<(*_sigmaR)[j]<<", weight="<<weight[j]<<endl;
915 }
916 }
917
918 // adaptive width
919 if (adaptive) {
920 cxcoutD(Eval) << "RooNDKeysPdf::calculateBandWidth() Using adaptive bandwidth." << endl;
921
922 double sqrt12 = sqrt(12.);
923 double sqrtSigmaAvgR = sqrt(_sigmaAvgR);
924
925 vector<Double_t> dummy(_nDim, 0.);
926 _weights1.resize(_nEvents, dummy);
927
928 std::vector<std::vector<Double_t>> *weights_prev(0);
929 std::vector<std::vector<Double_t>> *weights_new(0);
930
931 // cout << "Number of adaptive iterations: " << _nAdpt << endl;
932
933 for (Int_t k = 1; k <= _nAdpt; ++k) {
934
935 // cout << " Cycle: " << k << endl;
936
937 // if multiple adaptive iterations, need to swap weight sets
938 if (k % 2) {
939 weights_prev = &_weights0;
940 weights_new = &_weights1;
941 } else {
942 weights_prev = &_weights1;
943 weights_new = &_weights0;
944 }
945
946 for (Int_t i = 0; i < _nEvents; ++i) {
947 vector<Double_t> &x = _dataPts[i];
948 Double_t f = TMath::Power(gauss(x, *weights_prev) / _nEventsW, -1. / (2. * _d));
949
950 vector<Double_t> &weight = (*weights_new)[i];
951 for (Int_t j = 0; j < _nDim; j++) {
952 Double_t norm = (_n * (*_sigmaR)[j]) / sqrtSigmaAvgR;
953 weight[j] = norm * f / sqrt12; // note additional factor of sqrt(12) compared with HEP-EX/0011057
954 }
955 }
956 }
957 // this is the latest updated weights set
958 _weights = weights_new;
959 }
960}
961
962////////////////////////////////////////////////////////////////////////////////
963/// loop over all closest point to x, as determined by loopRange()
964
965Double_t RooNDKeysPdf::gauss(vector<Double_t>& x, vector<vector<Double_t> >& weights) const
966{
967 if(_nEvents==0) return 0.;
968
969 const double sqrt2pi = std::sqrt(TMath::TwoPi());
970
971 Double_t z=0.;
972 map<Int_t,Bool_t> ibMap;
973
974 // determine input loop range for event x
975 if (_sortInput) {
976 if (_sortTVIdcs.empty())
977 const_cast<RooNDKeysPdf*>(this)->sortDataIndices();
978
979 loopRange(x, ibMap);
980 }
981
982 for (const auto& ibMapItr : _sortInput ? ibMap : _ibNoSort) {
983 Int_t i = ibMapItr.first;
984
985 Double_t g(1.);
986
987 if (i >= (Int_t)_idx.size()) {
988 continue;
989 } //---> 1.myline
990
991 const vector<Double_t> &point = _dataPts[i];
992 const vector<Double_t> &weight = weights[_idx[i]];
993
994 for (Int_t j = 0; j < _nDim; j++) {
995 (*_dx)[j] = x[j] - point[j];
996 }
997
998 if (_nDim > 1 && _rotate) {
999 *_dx *= *_rotMat; // rotate to decorrelated frame!
1000 }
1001
1002 for (Int_t j = 0; j < _nDim; j++) {
1003 Double_t r = (*_dx)[j]; // x[j] - point[j];
1004 Double_t c = 1. / (2. * weight[j] * weight[j]);
1005
1006 // cout << "j = " << j << " x[j] = " << point[j] << " w = " << weight[j] << endl;
1007
1008 g *= exp(-c * r * r);
1009 g *= 1. / (sqrt2pi * weight[j]);
1010 }
1011 z += (g * _wMap.at(_idx[i]));
1012 }
1013 return z;
1014}
1015
1016////////////////////////////////////////////////////////////////////////////////
1017/// determine closest points to x, to loop over in evaluate()
1018
1019void RooNDKeysPdf::loopRange(vector<Double_t>& x, map<Int_t,Bool_t>& ibMap) const
1020{
1021 ibMap.clear();
1022 TVectorD xRm(_nDim);
1023 TVectorD xRp(_nDim);
1024
1025 for (Int_t j = 0; j < _nDim; j++) {
1026 xRm[j] = xRp[j] = x[j];
1027 }
1028
1029 if (_nDim > 1 && _rotate) {
1030 xRm *= *_rotMat;
1031 xRp *= *_rotMat;
1032 }
1033 for (Int_t j = 0; j < _nDim; j++) {
1034 xRm[j] -= _nSigma * (_n * (*_sigmaR)[j]);
1035 xRp[j] += _nSigma * (_n * (*_sigmaR)[j]);
1036 // cout<<"xRm["<<j<<"]="<<xRm[j]<<endl;
1037 // cout<<"xRp["<<j<<"]="<<xRp[j]<<endl;
1038 }
1039
1040 vector<TVectorD> xvecRm(1,xRm);
1041 vector<TVectorD> xvecRp(1,xRp);
1042
1043 map<Int_t,Bool_t> ibMapRT;
1044
1045 for (Int_t j=0; j<_nDim; j++) {
1046 ibMap.clear();
1047 auto comp = [=](const itPair& a, const itPair& b) {
1048 return (*a.second)[j] < (*b.second)[j];
1049 };
1050
1051 auto lo = lower_bound(_sortTVIdcs[j].begin(), _sortTVIdcs[j].end(), itPair(0, xvecRm.begin()), comp);
1052 auto hi = upper_bound(_sortTVIdcs[j].begin(), _sortTVIdcs[j].end(), itPair(0, xvecRp.begin()), comp);
1053 auto it = lo;
1054
1055 if (j==0) {
1056 if (_nDim==1) { for (it=lo; it!=hi; ++it) ibMap[(*it).first] = kTRUE; }
1057 else { for (it=lo; it!=hi; ++it) ibMapRT[(*it).first] = kTRUE; }
1058 continue;
1059 }
1060
1061 for (it=lo; it!=hi; ++it)
1062 if (ibMapRT.find((*it).first)!=ibMapRT.end()) { ibMap[(*it).first] = kTRUE; }
1063
1064 ibMapRT.clear();
1065 if (j!=_nDim-1) { ibMapRT = ibMap; }
1066 }
1067}
1068
1069////////////////////////////////////////////////////////////////////////////////
1070
1071void RooNDKeysPdf::boxInfoInit(BoxInfo* bi, const char* rangeName, Int_t /*code*/) const
1072{
1073 vector<Bool_t> doInt(_nDim,kTRUE);
1074
1075 bi->filled = kFALSE;
1076
1077 bi->xVarLo.resize(_nDim,0.);
1078 bi->xVarHi.resize(_nDim,0.);
1079 bi->xVarLoM3s.resize(_nDim,0.);
1080 bi->xVarLoP3s.resize(_nDim,0.);
1081 bi->xVarHiM3s.resize(_nDim,0.);
1082 bi->xVarHiP3s.resize(_nDim,0.);
1083
1084 bi->netFluxZ = kTRUE;
1085 bi->bpsIdcs.clear();
1086 bi->bIdcs.clear();
1087 bi->sIdcs.clear();
1088 bi->bmsIdcs.clear();
1089
1090 bi->nEventsBMSW=0.;
1091 bi->nEventsBW=0.;
1092
1093 for (unsigned int j=0; j < _varList.size(); ++j) {
1094 auto var = static_cast<const RooAbsRealLValue*>(_varList.at(j));
1095 if (doInt[j]) {
1096 bi->xVarLo[j] = var->getMin(rangeName);
1097 bi->xVarHi[j] = var->getMax(rangeName);
1098 } else {
1099 bi->xVarLo[j] = var->getVal() ;
1100 bi->xVarHi[j] = var->getVal() ;
1101 }
1102 }
1103}
1104
1105////////////////////////////////////////////////////////////////////////////////
1106
1108{
1109 if ( (_tracker && _tracker->hasChanged(kTRUE)) || (_weights != &_weights0 && _weights != &_weights1) ) {
1110 updateRho(); // update internal rho parameters
1111 // redetermine static and/or adaptive bandwidth
1112 const_cast<RooNDKeysPdf*>(this)->calculateBandWidth();
1113 }
1114
1115 const RooArgSet *nset = _varList.nset();
1116 for (unsigned int j=0; j < _varList.size(); ++j) {
1117 auto var = static_cast<const RooAbsReal*>(_varList.at(j));
1118 _x[j] = var->getVal(nset);
1119 }
1120
1121 Double_t val = gauss(_x,*_weights);
1122 //cout<<"returning "<<val<<endl;
1123
1124 if (val>=1E-20)
1125 return val ;
1126 else
1127 return (1E-20) ;
1128}
1129
1130////////////////////////////////////////////////////////////////////////////////
1131
1132Int_t RooNDKeysPdf::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const
1133{
1134 if (rangeName) return 0 ;
1135
1136 Int_t code=0;
1137 if (matchArgs(allVars,analVars,RooArgSet(_varList))) { code=1; }
1138
1139 return code;
1140
1141}
1142
1143////////////////////////////////////////////////////////////////////////////////
1144
1145Double_t RooNDKeysPdf::analyticalIntegral(Int_t code, const char* rangeName) const
1146{
1148
1149 cxcoutD(Eval) << "Calling RooNDKeysPdf::analyticalIntegral(" << GetName() << ") with code " << code
1150 << " and rangeName " << (rangeName?rangeName:"<none>") << endl;
1151
1152 // determine which observables need to be integrated over ...
1153 Int_t nComb = 1 << (_nDim);
1154 R__ASSERT(code>=1 && code<nComb) ;
1155
1156 vector<Bool_t> doInt(_nDim,kTRUE);
1157
1158 // get BoxInfo
1159 BoxInfo* bi(0);
1160
1161 if (rangeName) {
1162 string rangeNameStr(rangeName) ;
1163 bi = _rangeBoxInfo[make_pair(rangeNameStr,code)] ;
1164 if (!bi) {
1165 bi = new BoxInfo ;
1166 _rangeBoxInfo[make_pair(rangeNameStr,code)] = bi ;
1167 boxInfoInit(bi,rangeName,code);
1168 }
1169 } else bi= &_fullBoxInfo ;
1170
1171 // have boundaries changed?
1172 Bool_t newBounds(kFALSE);
1173 for (unsigned int j=0; j < _varList.size(); ++j) {
1174 auto var = static_cast<const RooAbsRealLValue*>(_varList.at(j));
1175 if ((var->getMin(rangeName)-bi->xVarLo[j]!=0) ||
1176 (var->getMax(rangeName)-bi->xVarHi[j]!=0)) {
1177 newBounds = kTRUE;
1178 }
1179 }
1180
1181 // reset
1182 if (newBounds) {
1183 cxcoutD(Eval) << "RooNDKeysPdf::analyticalIntegral() : Found new boundaries ... " << (rangeName?rangeName:"<none>") << endl;
1184 boxInfoInit(bi,rangeName,code);
1185 }
1186
1187 // recalculates netFluxZero and nEventsIR
1188 if (!bi->filled || newBounds) {
1189 // Fill box info with contents
1190 calculateShell(bi);
1191 calculatePreNorm(bi);
1192 bi->filled = kTRUE;
1193 const_cast<RooNDKeysPdf*>(this)->sortDataIndices(bi);
1194 }
1195
1196 // first guess
1197 Double_t norm=bi->nEventsBW;
1198
1199 if (_mirror && bi->netFluxZ) {
1200 // KEYS expression is self-normalized
1201 cxcoutD(Eval) << "RooNDKeysPdf::analyticalIntegral() : Using mirrored normalization : " << bi->nEventsBW << endl;
1202 return bi->nEventsBW;
1203 }
1204 // calculate leakage in and out of variable range box
1205 else
1206 {
1207 norm = bi->nEventsBMSW;
1208 if (norm<0.) norm=0.;
1209
1210 for (Int_t i=0; i<Int_t(bi->sIdcs.size()); ++i) {
1211 Double_t prob=1.;
1212 const vector<Double_t>& x = _dataPts[bi->sIdcs[i]];
1213 const vector<Double_t>& weight = (*_weights)[_idx[bi->sIdcs[i]]];
1214
1215 vector<Double_t> chi(_nDim,100.);
1216
1217 for (Int_t j=0; j<_nDim; j++) {
1218 if(!doInt[j]) continue;
1219
1220 if ((x[j]>bi->xVarLoM3s[j] && x[j]<bi->xVarLoP3s[j]) && x[j]<(bi->xVarLo[j]+bi->xVarHi[j])/2.)
1221 chi[j] = (x[j]-bi->xVarLo[j])/weight[j];
1222 else if ((x[j]>bi->xVarHiM3s[j] && x[j]<bi->xVarHiP3s[j]) && x[j]>(bi->xVarLo[j]+bi->xVarHi[j])/2.)
1223 chi[j] = (bi->xVarHi[j]-x[j])/weight[j];
1224
1225 if (chi[j]>0) // inVarRange
1226 prob *= (0.5 + TMath::Erf(fabs(chi[j])/sqrt(2.))/2.);
1227 else // outside Var range
1228 prob *= (0.5 - TMath::Erf(fabs(chi[j])/sqrt(2.))/2.);
1229 }
1230
1231 norm += prob * _wMap.at(_idx[bi->sIdcs[i]]);
1232 }
1233
1234 cxcoutD(Eval) << "RooNDKeysPdf::analyticalIntegral() : Final normalization : " << norm << " " << bi->nEventsBW << endl;
1235 return norm;
1236 }
1237}
1238
1239
1240////////////////////////////////////////////////////////////////////////////////
1241
1243{
1244 std::vector<RooRealVar *> varVec;
1245 RooArgSet varsAndWeightSet;
1246
1247 for (const auto var : varList) {
1248 if (!dynamic_cast<RooRealVar *>(var)) {
1249 coutE(InputArguments) << "RooNDKeysPdf::createDatasetFromHist(" << GetName() << ") WARNING: variable "
1250 << var->GetName() << " is not of type RooRealVar. Skip." << endl;
1251 continue;
1252 }
1253 varsAndWeightSet.add(*var); // used for dataset creation
1254 varVec.push_back(static_cast<RooRealVar *>(var)); // used for setting the variables.
1255 }
1256
1257 /// Add weight
1258 RooRealVar weight("weight", "event weight", 0);
1259 varsAndWeightSet.add(weight);
1260
1261 /// determine histogram dimensionality
1262 unsigned int histndim(0);
1263 std::string classname = hist.ClassName();
1264 if (classname.find("TH1") == 0) {
1265 histndim = 1;
1266 } else if (classname.find("TH2") == 0) {
1267 histndim = 2;
1268 } else if (classname.find("TH3") == 0) {
1269 histndim = 3;
1270 }
1271 assert(histndim == varVec.size());
1272
1273 if (histndim > 3 || histndim <= 0) {
1274 coutE(InputArguments) << "RooNDKeysPdf::createDatasetFromHist(" << GetName()
1275 << ") ERROR: input histogram dimension not between [1-3]: " << histndim << endl;
1276 assert(0);
1277 }
1278
1279 /// dataset creation
1280 RooDataSet *dataFromHist = new RooDataSet("datasetFromHist", "datasetFromHist", varsAndWeightSet, weight.GetName());
1281
1282 /// dataset filling
1283 for (int i = 1; i <= hist.GetXaxis()->GetNbins(); ++i) {
1284 // 1 or more dimension
1285
1286 Double_t xval = hist.GetXaxis()->GetBinCenter(i);
1287 varVec[0]->setVal(xval);
1288
1289 if (varVec.size() == 1) {
1290 Double_t fval = hist.GetBinContent(i);
1291 weight.setVal(fval);
1292 dataFromHist->add(varsAndWeightSet, fval);
1293 } else { // 2 or more dimensions
1294
1295 for (int j = 1; j <= hist.GetYaxis()->GetNbins(); ++j) {
1296 Double_t yval = hist.GetYaxis()->GetBinCenter(j);
1297 varVec[1]->setVal(yval);
1298
1299 if (varVec.size() == 2) {
1300 Double_t fval = hist.GetBinContent(i, j);
1301 weight.setVal(fval);
1302 dataFromHist->add(varsAndWeightSet, fval);
1303 } else { // 3 dimensions
1304
1305 for (int k = 1; k <= hist.GetZaxis()->GetNbins(); ++k) {
1306 Double_t zval = hist.GetZaxis()->GetBinCenter(k);
1307 varVec[2]->setVal(zval);
1308
1309 Double_t fval = hist.GetBinContent(i, j, k);
1310 weight.setVal(fval);
1311 dataFromHist->add(varsAndWeightSet, fval);
1312 }
1313 }
1314 }
1315 }
1316 }
1317
1318 return dataFromHist;
1319}
1320
1321
1322////////////////////////////////////////////////////////////////////////////////
1323/// Return evaluated weights
1325{
1327
1328 TMatrixD mref(_nEvents, _nDim + 1);
1329
1330 cxcoutD(Eval) << "RooNDKeysPdf::getWeights() Return evaluated weights." << endl;
1331
1332 for (Int_t i = 0; i < _nEvents; ++i) {
1333 const vector<Double_t>& x = _dataPts[i];
1334 for (Int_t j = 0; j < _nDim; j++) {
1335 mref(i, j) = x[j];
1336 }
1337
1338 const vector<Double_t>& weight = (*_weights)[i];
1339 mref(i, _nDim) = weight[k];
1340 }
1341
1342 return mref;
1343}
1344
1345 ////////////////////////////////////////////////////////////////////////////////
1347{
1348 for (unsigned int j = 0; j < _rhoList.size(); ++j) {
1349 auto rho = static_cast<const RooAbsReal*>(_rhoList.at(j));
1350 _rho[j] = rho->getVal();
1351 }
1352
1353 if (_nDim > 1 && _rotate) {
1354 TMatrixDSym covMatRho(_nDim); // covariance matrix times rho parameters
1355 for (Int_t j = 0; j < _nDim; j++) {
1356 for (Int_t k = 0; k < _nDim; k++) {
1357 covMatRho(j, k) = (*_covMat)(j, k) * _rho[j] * _rho[k];
1358 }
1359 }
1360 // find decorrelation matrix and eigenvalues (R)
1361 TMatrixDSymEigen evCalculatorRho(covMatRho);
1362 *_sigmaR = evCalculatorRho.GetEigenValues();
1363 for (Int_t j = 0; j < _nDim; j++) {
1364 (*_sigmaR)[j] = sqrt((*_sigmaR)[j]);
1365 }
1366 } else {
1367 for (Int_t j = 0; j < _nDim; j++) {
1368 (*_sigmaR)[j] = (_sigma[j] * _rho[j]);
1369 } // * rho
1370 }
1371}
typedef void(GLAPIENTRYP _GLUfuncptr)(void)
ROOT::R::TRInterface & r
Definition Object.C:4
#define b(i)
Definition RSha256.hxx:100
#define f(i)
Definition RSha256.hxx:104
#define c(i)
Definition RSha256.hxx:101
#define g(i)
Definition RSha256.hxx:105
#define a(i)
Definition RSha256.hxx:99
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define coutI(a)
#define cxcoutD(a)
#define coutW(a)
#define coutE(a)
std::vector< itPair > itVec
std::pair< Int_t, VecTVecDouble::iterator > itPair
int Int_t
Definition RtypesCore.h:45
const Bool_t kFALSE
Definition RtypesCore.h:101
const Bool_t kTRUE
Definition RtypesCore.h:100
#define ClassImp(name)
Definition Rtypes.h:364
#define R__ASSERT(e)
Definition TError.h:118
char name[80]
Definition TGX11.cxx:110
#define hi
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:687
friend class RooArgSet
Definition RooAbsArg.h:642
Int_t getSize() const
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
Storage_t::size_type size() const
const char * GetName() const
Returns name of object.
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.
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...
virtual Double_t getMin(const char *name=0) const
Get miniminum of currently defined range.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition RooAbsReal.h:64
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:94
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:35
RooChangeTracker is a meta object that tracks value changes in a given set of RooAbsArgs by registeri...
Bool_t hasChanged(Bool_t clearState)
Returns true if state has changed since last call with clearState=kTRUE.
RooDataSet is a container class to hold unbinned data.
Definition RooDataSet.h:36
virtual const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
virtual void add(const RooArgSet &row, Double_t weight=1.0, Double_t weightError=0) override
Add a data point, with its coordinates specified in the 'data' argset, to the data set.
virtual Double_t weight() const override
Return event weight of current event.
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE) override
Reimplementation of standard RooArgList::add()
Generic N-dimensional implementation of a kernel estimation p.d.f.
BoxInfo _fullBoxInfo
RooListProxy _rhoList
Double_t _nSigma
TVectorD * _dx
Double_t _minWeight
std::vector< Int_t > _sIdcs
void initialize()
initialization
std::vector< Double_t > _xVarLo
void calculatePreNorm(BoxInfo *bi) const
bi->nEventsBMSW=0.; bi->nEventsBW=0.;
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
TVectorD * _sigmaR
std::vector< Double_t > _xDatLo3s
Double_t _nEventsBW
std::vector< Double_t > _x0
virtual ~RooNDKeysPdf()
std::vector< Double_t > _xVarLoP3s
Bool_t _fixedShape
Bool_t _sortInput
std::vector< Double_t > _xDatHi3s
const RooDataSet * _data
std::vector< Int_t > _bIdcs
std::vector< TVectorD > _dataPtsR
void loopRange(std::vector< Double_t > &x, std::map< Int_t, Bool_t > &ibMap) const
determine closest points to x, to loop over in evaluate()
std::vector< Double_t > _xVarHi
void createPdf(Bool_t firstCall=kTRUE)
evaluation order of constructor.
Double_t _nEventsBMSW
void updateRho() const
Double_t _nEventsW
void calculateShell(BoxInfo *bi) const
determine points in +/- nSigma shell around the box determined by the variable ranges.
void calculateBandWidth()
std::vector< Double_t > _x1
Double_t _maxWeight
std::vector< Double_t > _mean
RooListProxy _varList
TMatrixDSym * _corrMat
std::vector< Int_t > _bmsIdcs
std::vector< Double_t > _rho
std::vector< Double_t > _xDatHi
std::vector< Double_t > _x
std::vector< itVec > _sortTVIdcs
Weights to be used. Points either to _weights0 or _weights1.
void boxInfoInit(BoxInfo *bi, const char *rangeName, Int_t code) const
std::map< std::pair< std::string, int >, BoxInfo * > _rangeBoxInfo
std::vector< Double_t > _xVarLoM3s
Double_t _sigmaAvgR
Double_t _widthFactor
std::vector< Int_t > _idx
std::vector< Double_t > _xDatLo
void sortDataIndices(BoxInfo *bi=0)
sort entries, as needed for loopRange()
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
void loadDataSet(Bool_t firstCall)
copy the dataset and calculate some useful variables
TString _options
do not persist
std::vector< std::vector< Double_t > > * _weights
TMatrixD getWeights(const int &k) const
Return evaluated weights.
std::vector< Double_t > _xVarHiM3s
Double_t evaluate() const
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
std::map< Int_t, Bool_t > _ibNoSort
std::vector< std::vector< Double_t > > _weights0
void mirrorDataSet()
determine mirror dataset.
void setOptions()
set the configuration
std::vector< std::vector< Double_t > > _dataPts
RooDataSet * createDatasetFromHist(const RooArgList &varList, const TH1 &hist) const
RooNDKeysPdf()=default
std::vector< Double_t > _sigma
std::vector< std::string > _varName
Double_t gauss(std::vector< Double_t > &x, std::vector< std::vector< Double_t > > &weights) const
loop over all closest point to x, as determined by loopRange()
void checkInitWeights() const
std::map< Int_t, Double_t > _wMap
std::vector< std::vector< Double_t > > _weights1
TMatrixD * _rotMat
std::vector< Double_t > _x2
TMatrixDSym * _covMat
std::vector< Double_t > _xVarHiP3s
std::map< Int_t, Bool_t > _bpsIdcs
RooChangeTracker * _tracker
RooRealVar represents a variable that can be changed from the outside.
Definition RooRealVar.h:39
virtual void setVal(Double_t value)
Set value of variable to 'value'.
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition TAxis.cxx:478
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:58
TAxis * GetZaxis()
Definition TH1.h:322
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition TH1.h:320
TAxis * GetYaxis()
Definition TH1.h:321
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:4994
TMatrixDSymEigen.
const TVectorD & GetEigenValues() const
const TMatrixD & GetEigenVectors() const
void Print(Option_t *name="") const
Print the matrix as a table of elements.
virtual TMatrixTBase< Element > & Zero()
Set matrix elements to zero.
TMatrixT< Element > & T()
Definition TMatrixT.h:150
virtual const char * GetName() const
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:200
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
Print the vector as a list of elements.
void box(Int_t pat, Double_t x1, Double_t y1, Double_t x2, Double_t y2)
Definition fillpatterns.C:1
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
Double_t Erf(Double_t x)
Computation of the error function erf(x).
Definition TMath.cxx:184
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition TMath.h:685
Short_t Abs(Short_t d)
Definition TMathBase.h:120
constexpr Double_t TwoPi()
Definition TMath.h:44
std::vector< Double_t > xVarHiP3s
std::vector< Double_t > xVarHi
std::map< Int_t, Bool_t > bpsIdcs
std::vector< Int_t > bIdcs
std::vector< Double_t > xVarLoM3s
std::vector< Double_t > xVarHiM3s
std::vector< Double_t > xVarLoP3s
std::vector< Double_t > xVarLo
std::vector< Int_t > sIdcs
std::vector< Int_t > bmsIdcs
auto * m
Definition textangle.C:8
auto * l
Definition textangle.C:4