Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooMultiVarGaussian.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitModels *
4 * @(#)root/roofit:$Id$
5 * Authors: *
6 * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7 * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8 * *
9 * Copyright (c) 2000-2005, Regents of the University of California *
10 * and Stanford University. All rights reserved. *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
14 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15 *****************************************************************************/
16
17/**
18\file RooMultiVarGaussian.cxx
19\class RooMultiVarGaussian
20\ingroup Roofitcore
21
22Multivariate Gaussian p.d.f. with correlations
23**/
24
25#include "Riostream.h"
26#include <cmath>
27
28#include "RooMultiVarGaussian.h"
29#include "RooAbsReal.h"
30#include "RooRealVar.h"
31#include "RooRandom.h"
32#include "RooMath.h"
33#include "RooGlobalFunc.h"
34#include "RooConstVar.h"
35#include "TDecompChol.h"
36#include "RooFitResult.h"
37
38using std::string, std::list, std::map, std::vector;
39
40
41////////////////////////////////////////////////////////////////////////////////
42
43RooMultiVarGaussian::RooMultiVarGaussian(const char *name, const char *title,
44 const RooArgList& xvec, const RooArgList& mu, const TMatrixDBase& cov) :
45 RooAbsPdf(name,title),
46 _x("x","Observables",this,true,false),
47 _mu("mu","Offset vector",this,true,false),
48 _cov{cov.GetNrows()},
49 _covI{cov.GetNrows()},
50 _z(4)
51{
52 if(!cov.IsSymmetric()) {
53 std::stringstream errorMsg;
54 errorMsg << "RooMultiVarGaussian::RooMultiVarGaussian(" << GetName()
55 << ") input covariance matrix is not symmetric!";
56 coutE(InputArguments) << errorMsg.str() << std::endl;
57 throw std::invalid_argument(errorMsg.str().c_str());
58 }
59
60 _cov.SetSub(0, cov);
61 _covI.SetSub(0, cov);
62
63 _x.add(xvec) ;
64
65 _mu.add(mu) ;
66
68
69 // Invert covariance matrix
70 _covI.Invert() ;
71}
72
73
74////////////////////////////////////////////////////////////////////////////////
75
76RooMultiVarGaussian::RooMultiVarGaussian(const char *name, const char *title, const RooArgList &xvec,
77 const RooFitResult &fr, bool reduceToConditional)
78 : RooAbsPdf(name, title),
79 _x("x", "Observables", this, true, false),
80 _mu("mu", "Offset vector", this, true, false),
81 _cov(reduceToConditional ? fr.conditionalCovarianceMatrix(xvec) : fr.reducedCovarianceMatrix(xvec)),
82 _covI(_cov),
83 _det(_cov.Determinant()),
84 _z(4)
85{
86
87 // Fill mu vector with constant RooRealVars
89 const RooArgList& fpf = fr.floatParsFinal() ;
90 for (std::size_t i=0 ; i<fpf.size() ; i++) {
91 if (xvec.find(fpf.at(i)->GetName())) {
92 std::unique_ptr<RooRealVar> parclone{static_cast<RooRealVar*>(fpf.at(i)->Clone(Form("%s_centralvalue",fpf.at(i)->GetName())))};
93 parclone->setConstant(true) ;
94 _mu.addOwned(std::move(parclone));
95 munames.push_back(fpf.at(i)->GetName()) ;
96 }
97 }
98
99 // Fill X vector in same order as mu vector
100 for (list<string>::iterator iter=munames.begin() ; iter!=munames.end() ; ++iter) {
101 RooRealVar* xvar = static_cast<RooRealVar*>(xvec.find(iter->c_str())) ;
102 _x.add(*xvar) ;
103 }
104
105 // Invert covariance matrix
106 _covI.Invert() ;
107
108}
109
110namespace {
111
113{
114 RooArgList out;
115 for (int i = 0; i < mu.GetNrows(); i++) {
116 out.add(RooFit::RooConst(mu(i)));
117 }
118 return out;
119}
120
121RooArgList muConstants(std::size_t n)
122{
123 out;
124 for (std::size_t i = 0; i < n; i++) {
125 RooFit::RooConst(0));
126 }
127 return out;
128}
129
130} // namespace
131
132
133////////////////////////////////////////////////////////////////////////////////
134
135RooMultiVarGaussian::RooMultiVarGaussian(const char *name, const char *title, const RooArgList &xvec,
136 const TVectorD &mu, const TMatrixDBase &cov)
137 : RooMultiVarGaussian{name, title, xvec, muFromVector(mu), cov}
138{
139}
140
141////////////////////////////////////////////////////////////////////////////////
142
143RooMultiVarGaussian::RooMultiVarGaussian(const char *name, const char *title, const RooArgList &xvec,
144 const TMatrixDBase &cov)
145 : RooMultiVarGaussian{name, title, xvec, muConstants(xvec.size()), cov}
146{
147}
148
149
150////////////////////////////////////////////////////////////////////////////////
151
153 RooAbsPdf(other,name), _aicMap(other._aicMap), _x("x",this,other._x), _mu("mu",this,other._mu),
154 _cov(other._cov), _covI(other._covI), _det(other._det), _z(other._z)
155{
156}
157
158
159
160////////////////////////////////////////////////////////////////////////////////
161
163{
165 for (std::size_t i=0 ; i<_mu.size() ; i++) {
166 _muVec[i] = static_cast<RooAbsReal*>(_mu.at(i))->getVal() ;
167 }
168}
169
170
171////////////////////////////////////////////////////////////////////////////////
172/// Represent observables as vector
173
175{
176 TVectorD x(_x.size()) ;
177 for (std::size_t i=0 ; i<_x.size() ; i++) {
178 x[i] = static_cast<RooAbsReal*>(_x.at(i))->getVal() ;
179 }
180
181 // Calculate return value
182 syncMuVec() ;
184
185 double alpha = x_min_mu * (_covI * x_min_mu) ;
186 return exp(-0.5*alpha) ;
187}
188
189////////////////////////////////////////////////////////////////////////////////
190
192{
193 RooArgSet allVars(allVarsIn) ;
194
195 // If allVars contains x_i it cannot contain mu_i
196 for (std::size_t i=0 ; i<_x.size() ; i++) {
197 if (allVars.contains(*_x.at(i))) {
198 allVars.remove(*_mu.at(i),true,true) ;
199 }
200 }
201
202
203 // Analytical integral known over all observables
204 if (allVars.size()==_x.size() && !rangeName) {
205 analVars.add(allVars) ;
206 return -1 ;
207 }
208
209 Int_t code(0) ;
210
211 Int_t nx = _x.size() ;
212 if (nx>127) {
213 // Warn that analytical integration is only provided for the first 127 observables
214 coutW(Integration) << "RooMultiVarGaussian::getAnalyticalIntegral(" << GetName() << ") WARNING: p.d.f. has " << _x.size()
215 << " observables, analytical integration is only implemented for the first 127 observables" << std::endl ;
216 nx=127 ;
217 }
218
219 // Advertise partial analytical integral over all observables for which is wide enough to
220 // use asymptotic integral calculation
221 BitBlock bits ;
222 bool anyBits(false) ;
223 syncMuVec() ;
224 for (std::size_t i=0 ; i<_x.size() ; i++) {
225
226 // Check if integration over observable #i is requested
227 if (allVars.find(_x.at(i)->GetName())) {
228 // Check if range is wider than Z sigma
229 RooRealVar* xi = static_cast<RooRealVar*>(_x.at(i)) ;
230 if (xi->getMin(rangeName)<_muVec(i)-_z*sqrt(_cov(i,i)) && xi->getMax(rangeName) > _muVec(i)+_z*sqrt(_cov(i,i))) {
231 cxcoutD(Integration) << "RooMultiVarGaussian::getAnalyticalIntegral(" << GetName()
232 << ") Advertising analytical integral over " << xi->GetName() << " as range is >" << _z << " sigma" << std::endl ;
233 bits.setBit(i) ;
234 anyBits = true ;
235 analVars.add(*allVars.find(_x.at(i)->GetName())) ;
236 } else {
237 cxcoutD(Integration) << "RooMultiVarGaussian::getAnalyticalIntegral(" << GetName() << ") Range of " << xi->GetName() << " is <"
238 << _z << " sigma, relying on numeric integral" << std::endl ;
239 }
240 }
241
242 // Check if integration over parameter #i is requested
243 if (allVars.find(_mu.at(i)->GetName())) {
244 // Check if range is wider than Z sigma
245 RooRealVar* pi = static_cast<RooRealVar*>(_mu.at(i)) ;
246 if (pi->getMin(rangeName)<_muVec(i)-_z*sqrt(_cov(i,i)) && pi->getMax(rangeName) > _muVec(i)+_z*sqrt(_cov(i,i))) {
247 cxcoutD(Integration) << "RooMultiVarGaussian::getAnalyticalIntegral(" << GetName()
248 << ") Advertising analytical integral over " << pi->GetName() << " as range is >" << _z << " sigma" << std::endl ;
249 bits.setBit(i) ;
250 anyBits = true ;
251 analVars.add(*allVars.find(_mu.at(i)->GetName())) ;
252 } else {
253 cxcoutD(Integration) << "RooMultiVarGaussian::getAnalyticalIntegral(" << GetName() << ") Range of " << pi->GetName() << " is <"
254 << _z << " sigma, relying on numeric integral" << std::endl ;
255 }
256 }
257
258
259 }
260
261 // Full numeric integration over requested observables maps always to code zero
262 if (!anyBits) {
263 return 0 ;
264 }
265
266 // Map BitBlock into return code
267 for (UInt_t i=0 ; i<_aicMap.size() ; i++) {
268 if (_aicMap[i]==bits) {
269 code = i+1 ;
270 }
271 }
272 if (code==0) {
273 _aicMap.push_back(bits) ;
274 code = _aicMap.size() ;
275 }
276
277 return code ;
278}
279
280
281
282////////////////////////////////////////////////////////////////////////////////
283/// Handle full integral here
284
285double RooMultiVarGaussian::analyticalIntegral(Int_t code, const char* /*rangeName*/) const
286{
287 if (code==-1) {
288 return pow(2*3.14159268,_x.size()/2.)*sqrt(std::abs(_det)) ;
289 }
290
291 // Handle partial integrals here
292
293 // Retrieve |S22|, S22bar from cache
294 AnaIntData& aid = anaIntData(code) ;
295
296 // Fill position vector for non-integrated observables
297 syncMuVec() ;
298 TVectorD u(aid.pmap.size()) ;
299 for (UInt_t i=0 ; i<aid.pmap.size() ; i++) {
300 u(i) = (static_cast<RooAbsReal*>(_x.at(aid.pmap[i])))->getVal() - _muVec(aid.pmap[i]) ;
301 }
302
303 // Calculate partial integral
304 double ret = pow(2*3.14159268,aid.nint/2.)/sqrt(std::abs(aid.S22det))*exp(-0.5*u*(aid.S22bar*u)) ;
305
306 return ret ;
307}
308
309
310////////////////////////////////////////////////////////////////////////////////
311/// Check if cache entry was previously created
312
314{
315 map<int,AnaIntData>::iterator iter = _anaIntCache.find(code) ;
316 if (iter != _anaIntCache.end()) {
317 return iter->second ;
318 }
319
320 // Calculate cache contents
321
322 // Decode integration code
325 decodeCode(code,map1,map2) ;
326
327 // Rearrange observables so that all non-integrated observables
328 // go first (preserving relative order) and all integrated observables
329 // go last (preserving relative order)
335
336 // Begin calculation of partial integrals
337 // ___
338 // sqrt(2pi)^(#intObs) (-0.5 * u1T S22 u1 )
339 // I = ------------------- * e
340 // sqrt(|det(S22)|)
341 // ___
342 // Where S22 is the sub-matrix of covI for the integrated observables and S22
343 // is the Schur complement of S22
344 // ___ -1
345 // S22 = S11 - S12 * S22 * S21
346 //
347 // and u1 is the vector of non-integrated observables
348
349 // Calculate Schur complement S22bar
351 S22inv.Invert() ;
352 TMatrixD S22bar = S11 - S12*S22inv*S21 ;
353
354 // Create new cache entry
356 cacheData.S22bar.ResizeTo(S22bar) ;
357 cacheData.S22bar=S22bar ;
358 cacheData.S22det= S22.Determinant() ;
359 cacheData.pmap = map1 ;
360 cacheData.nint = map2.size() ;
361
362 return cacheData ;
363}
364
365
366
367////////////////////////////////////////////////////////////////////////////////
368/// Special case: generate all observables
369
371{
372 if (directVars.size()==_x.size()) {
374 return -1 ;
375 }
376
377 Int_t nx = _x.size() ;
378 if (nx>127) {
379 // Warn that analytical integration is only provided for the first 127 observables
380 coutW(Integration) << "RooMultiVarGaussian::getGenerator(" << GetName() << ") WARNING: p.d.f. has " << _x.size()
381 << " observables, partial internal generation is only implemented for the first 127 observables" << std::endl ;
382 nx=127 ;
383 }
384
385 // Advertise partial generation over all permutations of observables
386 Int_t code(0) ;
387 BitBlock bits ;
388 for (std::size_t i=0 ; i<_x.size() ; i++) {
389 RooAbsArg* arg = directVars.find(_x.at(i)->GetName()) ;
390 if (arg) {
391 bits.setBit(i) ;
392// code |= (1<<i) ;
393 generateVars.add(*arg) ;
394 }
395 }
396
397 // Map BitBlock into return code
398 for (UInt_t i=0 ; i<_aicMap.size() ; i++) {
399 if (_aicMap[i]==bits) {
400 code = i+1 ;
401 }
402 }
403 if (code==0) {
404 _aicMap.push_back(bits) ;
405 code = _aicMap.size() ;
406 }
407
408
409 return code ;
410}
411
412
413
414////////////////////////////////////////////////////////////////////////////////
415/// Clear the GenData cache as its content is not invariant under changes in
416/// the mu vector.
417
419{
420 _genCache.clear() ;
421
422}
423
424
425
426
427////////////////////////////////////////////////////////////////////////////////
428/// Retrieve generator config from cache
429
431{
432 GenData& gd = genData(code) ;
433 TMatrixD& TU = gd.UT ;
434 Int_t nobs = TU.GetNcols() ;
435 vector<int>& omap = gd.omap ;
436
437 while(true) {
438
439 // Create unit Gaussian vector
441 for(Int_t k= 0; k <nobs; k++) {
443 }
444
445 // Apply transformation matrix
446 xgen *= TU ;
447
448 // Apply shift
449 if (code == -1) {
450
451 // Simple shift if we generate all observables
452 xgen += gd.mu1 ;
453
454 } else {
455
456 // Non-generated observable dependent shift for partial generations
457
458 // mubar = mu1 + S12 S22Inv ( x2 - mu2)
459 TVectorD mubar(gd.mu1) ;
460 TVectorD x2(gd.pmap.size()) ;
461 for (UInt_t i=0 ; i<gd.pmap.size() ; i++) {
462 x2(i) = (static_cast<RooAbsReal*>(_x.at(gd.pmap[i])))->getVal() ;
463 }
464 mubar += gd.S12S22I * (x2 - gd.mu2) ;
465
466 xgen += mubar ;
467
468 }
469
470 // Transfer values and check if values are in range
471 bool ok(true) ;
472 for (int i=0 ; i<nobs ; i++) {
473 RooRealVar* xi = static_cast<RooRealVar*>(_x.at(omap[i])) ;
474 if (xgen(i)<xi->getMin() || xgen(i)>xi->getMax()) {
475 ok = false ;
476 break ;
477 } else {
478 xi->setVal(xgen(i)) ;
479 }
480 }
481
482 // If all values are in range, accept event and return
483 // otherwise retry
484 if (ok) {
485 break ;
486 }
487 }
488
489 return;
490}
491
492
493
494////////////////////////////////////////////////////////////////////////////////
495/// WVE -- CHECK THAT GENDATA IS VALID GIVEN CURRENT VALUES OF _MU
496
498{
499 // Check if cache entry was previously created
500 map<int,GenData>::iterator iter = _genCache.find(code) ;
501 if (iter != _genCache.end()) {
502 return iter->second ;
503 }
504
505 // Create new entry
506 GenData& cacheData = _genCache[code] ;
507
508 if (code==-1) {
509
510 // Do eigen value decomposition
512 tdc.Decompose() ;
513 TMatrixD U = tdc.GetU() ;
515
516 // Fill cache data
517 cacheData.UT.ResizeTo(TU) ;
518 cacheData.UT = TU ;
519 cacheData.omap.resize(_x.size()) ;
520 for (std::size_t i=0 ; i<_x.size() ; i++) {
521 cacheData.omap[i] = i ;
522 }
523 syncMuVec() ;
524 cacheData.mu1.ResizeTo(_muVec) ;
525 cacheData.mu1 = _muVec ;
526
527 } else {
528
529 // Construct observables: map1 = generated, map2 = given
532 decodeCode(code,map2,map1) ;
533
534 // Do block decomposition of covariance matrix
540
541 // Constructed conditional matrix form
542 // -1
543 // F(X1|X2) --> CovI --> S22bar = S11 - S12 S22 S21
544 // -1
545 // --> mu --> mubar = mu1 + S12 S22 ( x2 - mu2)
546
547 // Do eigenvalue decomposition
549 TMatrixD S22bar = S11 - S12 * (S22Inv * S21) ;
550
551 // Do eigen value decomposition of S22bar
552 TDecompChol tdc(S22bar) ;
553 tdc.Decompose() ;
554 TMatrixD U = tdc.GetU() ;
556
557 // Split mu vector into mu1 and mu2
558 TVectorD mu1(map1.size());
559 TVectorD mu2(map2.size());
560 syncMuVec() ;
561 for (UInt_t i=0 ; i<map1.size() ; i++) {
562 mu1(i) = _muVec(map1[i]) ;
563 }
564 for (UInt_t i=0 ; i<map2.size() ; i++) {
565 mu2(i) = _muVec(map2[i]) ;
566 }
567
568 // Calculate rotation matrix for mu vector
570
571 // Fill cache data
572 cacheData.UT.ResizeTo(TU) ;
573 cacheData.UT = TU ;
574 cacheData.omap = map1 ;
575 cacheData.pmap = map2 ;
576 cacheData.mu1.ResizeTo(mu1) ;
577 cacheData.mu2.ResizeTo(mu2) ;
578 cacheData.mu1 = mu1 ;
579 cacheData.mu2 = mu2 ;
580 cacheData.S12S22I.ResizeTo(S12S22Inv) ;
581 cacheData.S12S22I = S12S22Inv ;
582
583 }
584
585
586 return cacheData ;
587}
588
589
590
591
592////////////////////////////////////////////////////////////////////////////////
593/// Decode analytical integration/generation code into index map of integrated/generated (map2)
594/// and non-integrated/generated observables (map1)
595
597{
598 if (code<0 || code> (Int_t)_aicMap.size()) {
599 std::cout << "RooMultiVarGaussian::decodeCode(" << GetName() << ") ERROR don't have bit pattern for code " << code << std::endl ;
600 throw string("RooMultiVarGaussian::decodeCode() ERROR don't have bit pattern for code") ;
601 }
602
603 BitBlock b = _aicMap[code-1] ;
604 map1.clear() ;
605 map2.clear() ;
606 for (std::size_t i=0 ; i<_x.size() ; i++) {
607 if (b.getBit(i)) {
608 map2.push_back(i) ;
609 } else {
610 map1.push_back(i) ;
611 }
612 }
613}
614
615
616////////////////////////////////////////////////////////////////////////////////
617/// Block decomposition of covI according to given maps of observables
618
620{
621 // Allocate and fill reordered covI matrix in 2x2 block structure
622
623 S11.ResizeTo(map1.size(),map1.size()) ;
624 S12.ResizeTo(map1.size(),map2.size()) ;
625 S21.ResizeTo(map2.size(),map1.size()) ;
626 S22.ResizeTo(map2.size(),map2.size()) ;
627
628 for (UInt_t i=0 ; i<map1.size() ; i++) {
629 for (UInt_t j=0 ; j<map1.size() ; j++)
630 S11(i,j) = input(map1[i],map1[j]) ;
631 for (UInt_t j=0 ; j<map2.size() ; j++)
632 S12(i,j) = input(map1[i],map2[j]) ;
633 }
634 for (UInt_t i=0 ; i<map2.size() ; i++) {
635 for (UInt_t j=0 ; j<map1.size() ; j++)
636 S21(i,j) = input(map2[i],map1[j]) ;
637 for (UInt_t j=0 ; j<map2.size() ; j++)
638 S22(i,j) = input(map2[i],map2[j]) ;
639 }
640
641}
642
643
645{
646 if (ibit<32) { b0 |= (1<<ibit) ; return ; }
647 if (ibit<64) { b1 |= (1<<(ibit-32)) ; return ; }
648 if (ibit<96) { b2 |= (1<<(ibit-64)) ; return ; }
649 if (ibit<128) { b3 |= (1<<(ibit-96)) ; return ; }
650}
651
653{
654 if (ibit<32) return (b0 & (1<<ibit)) ;
655 if (ibit<64) return (b1 & (1<<(ibit-32))) ;
656 if (ibit<96) return (b2 & (1<<(ibit-64))) ;
657 if (ibit<128) return (b3 & (1<<(ibit-96))) ;
658 return false ;
659}
660
662{
663 if (lhs.b0 != rhs.b0)
664 return false;
665 if (lhs.b1 != rhs.b1)
666 return false;
667 if (lhs.b2 != rhs.b2)
668 return false;
669 if (lhs.b3 != rhs.b3)
670 return false;
671 return true;
672}
#define b(i)
Definition RSha256.hxx:100
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define cxcoutD(a)
#define coutW(a)
#define coutE(a)
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void input
Option_t Option_t TPoint TPoint const char x2
char name[80]
Definition TGX11.cxx:110
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition TString.cxx:2489
const_iterator begin() const
const_iterator end() const
Common abstract base class for objects that represent a value and a "shape" in RooFit.
Definition RooAbsArg.h:77
virtual bool remove(const RooAbsArg &var, bool silent=false, bool matchByNameOnly=false)
Remove the specified argument from our list.
bool contains(const RooAbsArg &var) const
Check if collection contains an argument with the same name as var.
Storage_t::size_type size() const
RooAbsArg * find(const char *name) const
Find object with given name in list.
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
virtual double getMax(const char *name=nullptr) const
Get maximum of currently defined range.
virtual double getMin(const char *name=nullptr) const
Get minimum of currently defined range.
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
bool operator==(double value) const
Equality operator comparing to a double.
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:24
bool addOwned(RooAbsArg &var, bool silent=false) override
Overloaded RooCollection_t::addOwned() method insert object into owning set and registers object as s...
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...
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
const RooArgList & floatParsFinal() const
Return list of floating parameters after fit.
Multivariate Gaussian p.d.f.
double analyticalIntegral(Int_t code, const char *rangeName=nullptr) const override
Handle full integral here.
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=nullptr) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
void decodeCode(Int_t code, std::vector< int > &map1, std::vector< int > &map2) const
Decode analytical integration/generation code into index map of integrated/generated (map2) and non-i...
std::vector< BitBlock > _aicMap
!
static void blockDecompose(const TMatrixD &input, const std::vector< int > &map1, const std::vector< int > &map2, TMatrixDSym &S11, TMatrixD &S12, TMatrixD &S21, TMatrixDSym &S22)
Block decomposition of covI according to given maps of observables.
AnaIntData & anaIntData(Int_t code) const
Check if cache entry was previously created.
double evaluate() const override
Do not persist.
std::map< int, GenData > _genCache
!
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, bool staticInitOK=true) const override
Special case: generate all observables.
void initGenerator(Int_t code) override
Clear the GenData cache as its content is not invariant under changes in the mu vector.
GenData & genData(Int_t code) const
WVE – CHECK THAT GENDATA IS VALID GIVEN CURRENT VALUES OF _MU.
std::map< int, AnaIntData > _anaIntCache
!
void generateEvent(Int_t code) override
Retrieve generator config from cache.
static double gaussian(TRandom *generator=randomGenerator())
Return a Gaussian random variable with mean 0 and variance 1.
Variable that can be changed from the outside.
Definition RooRealVar.h:37
void setVal(double value) override
Set value of variable to 'value'.
Cholesky Decomposition class.
Definition TDecompChol.h:25
TMatrixTSym< Element > & SetSub(Int_t row_lwb, const TMatrixTBase< Element > &source)
Insert matrix source starting at [row_lwb][row_lwb], thereby overwriting the part [row_lwb....
Double_t Determinant() const override
TMatrixTSym< Element > & Invert(Double_t *det=nullptr)
Invert the matrix and calculate its determinant Notice that the LU decomposition is used instead of B...
const char * GetName() const override
Returns name of object.
Definition TNamed.h:49
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition TVectorT.cxx:294
Int_t GetNrows() const
Definition TVectorT.h:75
RooConstVar & RooConst(double val)
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16