Logo ROOT  
Reference Guide
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 "RooFit.h"
26
27#include "Riostream.h"
28#include <math.h>
29
30#include "RooMultiVarGaussian.h"
31#include "RooAbsReal.h"
32#include "RooRealVar.h"
33#include "RooRandom.h"
34#include "RooMath.h"
35#include "RooGlobalFunc.h"
36#include "RooConstVar.h"
37#include "TDecompChol.h"
38#include "RooFitResult.h"
39
40using namespace std;
41
43 ;
44
45////////////////////////////////////////////////////////////////////////////////
46
47RooMultiVarGaussian::RooMultiVarGaussian(const char *name, const char *title,
48 const RooArgList& xvec, const RooArgList& mu, const TMatrixDSym& cov) :
49 RooAbsPdf(name,title),
50 _x("x","Observables",this,kTRUE,kFALSE),
51 _mu("mu","Offset vector",this,kTRUE,kFALSE),
52 _cov(cov),
53 _covI(cov),
54 _z(4)
55{
56 _x.add(xvec) ;
57
58 _mu.add(mu) ;
59
61
62 // Invert covariance matrix
63 _covI.Invert() ;
64}
65
66
67////////////////////////////////////////////////////////////////////////////////
68
69RooMultiVarGaussian::RooMultiVarGaussian(const char *name, const char *title,
70 const RooArgList& xvec, const RooFitResult& fr, Bool_t reduceToConditional) :
71 RooAbsPdf(name,title),
72 _x("x","Observables",this,kTRUE,kFALSE),
73 _mu("mu","Offset vector",this,kTRUE,kFALSE),
74 _cov(reduceToConditional ? fr.conditionalCovarianceMatrix(xvec) : fr.reducedCovarianceMatrix(xvec)),
75 _covI(_cov),
76 _z(4)
77{
79
80 // Fill mu vector with constant RooRealVars
81 list<string> munames ;
82 const RooArgList& fpf = fr.floatParsFinal() ;
83 for (Int_t i=0 ; i<fpf.getSize() ; i++) {
84 if (xvec.find(fpf.at(i)->GetName())) {
85 RooRealVar* parclone = (RooRealVar*) fpf.at(i)->Clone(Form("%s_centralvalue",fpf.at(i)->GetName())) ;
86 parclone->setConstant(kTRUE) ;
87 _mu.addOwned(*parclone) ;
88 munames.push_back(fpf.at(i)->GetName()) ;
89 }
90 }
91
92 // Fill X vector in same order as mu vector
93 for (list<string>::iterator iter=munames.begin() ; iter!=munames.end() ; ++iter) {
94 RooRealVar* xvar = (RooRealVar*) xvec.find(iter->c_str()) ;
95 _x.add(*xvar) ;
96 }
97
98 // Invert covariance matrix
99 _covI.Invert() ;
100
101}
102
103
104////////////////////////////////////////////////////////////////////////////////
105
106RooMultiVarGaussian::RooMultiVarGaussian(const char *name, const char *title,
107 const RooArgList& xvec, const TVectorD& mu, const TMatrixDSym& cov) :
108 RooAbsPdf(name,title),
109 _x("x","Observables",this,kTRUE,kFALSE),
110 _mu("mu","Offset vector",this,kTRUE,kFALSE),
111 _cov(cov),
112 _covI(cov),
113 _z(4)
114{
115 _x.add(xvec) ;
116
117 for (Int_t i=0 ; i<mu.GetNrows() ; i++) {
118 _mu.add(RooFit::RooConst(mu(i))) ;
119 }
120
121 _det = _cov.Determinant() ;
122
123 // Invert covariance matrix
124 _covI.Invert() ;
125}
126
127////////////////////////////////////////////////////////////////////////////////
128
129RooMultiVarGaussian::RooMultiVarGaussian(const char *name, const char *title,
130 const RooArgList& xvec, const TMatrixDSym& cov) :
131 RooAbsPdf(name,title),
132 _x("x","Observables",this,kTRUE,kFALSE),
133 _mu("mu","Offset vector",this,kTRUE,kFALSE),
134 _cov(cov),
135 _covI(cov),
136 _z(4)
137{
138 _x.add(xvec) ;
139
140 for (Int_t i=0 ; i<xvec.getSize() ; i++) {
142 }
143
144 _det = _cov.Determinant() ;
145
146 // Invert covariance matrix
147 _covI.Invert() ;
148}
149
150
151
152////////////////////////////////////////////////////////////////////////////////
153
155 RooAbsPdf(other,name), _aicMap(other._aicMap), _x("x",this,other._x), _mu("mu",this,other._mu),
156 _cov(other._cov), _covI(other._covI), _det(other._det), _z(other._z)
157{
158}
159
160
161
162////////////////////////////////////////////////////////////////////////////////
163
165{
167 for (Int_t i=0 ; i<_mu.getSize() ; i++) {
168 _muVec[i] = ((RooAbsReal*)_mu.at(i))->getVal() ;
169 }
170}
171
172
173////////////////////////////////////////////////////////////////////////////////
174/// Represent observables as vector
175
177{
178 TVectorD x(_x.getSize()) ;
179 for (int i=0 ; i<_x.getSize() ; i++) {
180 x[i] = ((RooAbsReal*)_x.at(i))->getVal() ;
181 }
182
183 // Calculate return value
184 syncMuVec() ;
185 TVectorD x_min_mu = x - _muVec ;
186
187 Double_t alpha = x_min_mu * (_covI * x_min_mu) ;
188 return exp(-0.5*alpha) ;
189}
190
191
192
193////////////////////////////////////////////////////////////////////////////////
194
195Int_t RooMultiVarGaussian::getAnalyticalIntegral(RooArgSet& allVarsIn, RooArgSet& analVars, const char* rangeName) const
196{
197 RooArgSet allVars(allVarsIn) ;
198
199 // If allVars contains x_i it cannot contain mu_i
200 for (Int_t i=0 ; i<_x.getSize() ; i++) {
201 if (allVars.contains(*_x.at(i))) {
202 allVars.remove(*_mu.at(i),kTRUE,kTRUE) ;
203 }
204 }
205
206
207 // Analytical integral known over all observables
208 if (allVars.getSize()==_x.getSize() && !rangeName) {
209 analVars.add(allVars) ;
210 return -1 ;
211 }
212
213 Int_t code(0) ;
214
215 Int_t nx = _x.getSize() ;
216 if (nx>127) {
217 // Warn that analytical integration is only provided for the first 127 observables
218 coutW(Integration) << "RooMultiVarGaussian::getAnalyticalIntegral(" << GetName() << ") WARNING: p.d.f. has " << _x.getSize()
219 << " observables, analytical integration is only implemented for the first 127 observables" << endl ;
220 nx=127 ;
221 }
222
223 // Advertise partial analytical integral over all observables for which is wide enough to
224 // use asymptotic integral calculation
225 BitBlock bits ;
226 Bool_t anyBits(kFALSE) ;
227 syncMuVec() ;
228 for (int i=0 ; i<_x.getSize() ; i++) {
229
230 // Check if integration over observable #i is requested
231 if (allVars.find(_x.at(i)->GetName())) {
232 // Check if range is wider than Z sigma
233 RooRealVar* xi = (RooRealVar*)_x.at(i) ;
234 if (xi->getMin(rangeName)<_muVec(i)-_z*sqrt(_cov(i,i)) && xi->getMax(rangeName) > _muVec(i)+_z*sqrt(_cov(i,i))) {
235 cxcoutD(Integration) << "RooMultiVarGaussian::getAnalyticalIntegral(" << GetName()
236 << ") Advertising analytical integral over " << xi->GetName() << " as range is >" << _z << " sigma" << endl ;
237 bits.setBit(i) ;
238 anyBits = kTRUE ;
239 analVars.add(*allVars.find(_x.at(i)->GetName())) ;
240 } else {
241 cxcoutD(Integration) << "RooMultiVarGaussian::getAnalyticalIntegral(" << GetName() << ") Range of " << xi->GetName() << " is <"
242 << _z << " sigma, relying on numeric integral" << endl ;
243 }
244 }
245
246 // Check if integration over parameter #i is requested
247 if (allVars.find(_mu.at(i)->GetName())) {
248 // Check if range is wider than Z sigma
249 RooRealVar* pi = (RooRealVar*)_mu.at(i) ;
250 if (pi->getMin(rangeName)<_muVec(i)-_z*sqrt(_cov(i,i)) && pi->getMax(rangeName) > _muVec(i)+_z*sqrt(_cov(i,i))) {
251 cxcoutD(Integration) << "RooMultiVarGaussian::getAnalyticalIntegral(" << GetName()
252 << ") Advertising analytical integral over " << pi->GetName() << " as range is >" << _z << " sigma" << endl ;
253 bits.setBit(i) ;
254 anyBits = kTRUE ;
255 analVars.add(*allVars.find(_mu.at(i)->GetName())) ;
256 } else {
257 cxcoutD(Integration) << "RooMultiVarGaussian::getAnalyticalIntegral(" << GetName() << ") Range of " << pi->GetName() << " is <"
258 << _z << " sigma, relying on numeric integral" << endl ;
259 }
260 }
261
262
263 }
264
265 // Full numeric integration over requested observables maps always to code zero
266 if (!anyBits) {
267 return 0 ;
268 }
269
270 // Map BitBlock into return code
271 for (UInt_t i=0 ; i<_aicMap.size() ; i++) {
272 if (_aicMap[i]==bits) {
273 code = i+1 ;
274 }
275 }
276 if (code==0) {
277 _aicMap.push_back(bits) ;
278 code = _aicMap.size() ;
279 }
280
281 return code ;
282}
283
284
285
286////////////////////////////////////////////////////////////////////////////////
287/// Handle full integral here
288
289Double_t RooMultiVarGaussian::analyticalIntegral(Int_t code, const char* /*rangeName*/) const
290{
291 if (code==-1) {
292 return pow(2*3.14159268,_x.getSize()/2.)*sqrt(fabs(_det)) ;
293 }
294
295 // Handle partial integrals here
296
297 // Retrieve |S22|, S22bar from cache
298 AnaIntData& aid = anaIntData(code) ;
299
300 // Fill position vector for non-integrated observables
301 syncMuVec() ;
302 TVectorD u(aid.pmap.size()) ;
303 for (UInt_t i=0 ; i<aid.pmap.size() ; i++) {
304 u(i) = ((RooAbsReal*)_x.at(aid.pmap[i]))->getVal() - _muVec(aid.pmap[i]) ;
305 }
306
307 // Calculate partial integral
308 Double_t ret = pow(2*3.14159268,aid.nint/2.)/sqrt(fabs(aid.S22det))*exp(-0.5*u*(aid.S22bar*u)) ;
309
310 return ret ;
311}
312
313
314
315////////////////////////////////////////////////////////////////////////////////
316/// Check if cache entry was previously created
317
319{
320 map<int,AnaIntData>::iterator iter = _anaIntCache.find(code) ;
321 if (iter != _anaIntCache.end()) {
322 return iter->second ;
323 }
324
325 // Calculate cache contents
326
327 // Decode integration code
328 vector<int> map1,map2 ;
329 decodeCode(code,map1,map2) ;
330
331 // Rearrage observables so that all non-integrated observables
332 // go first (preserving relative order) and all integrated observables
333 // go last (preserving relative order)
334 TMatrixDSym S11, S22 ;
335 TMatrixD S12, S21 ;
336 blockDecompose(_covI,map1,map2,S11,S12,S21,S22) ;
337
338 // Begin calculation of partial integrals
339 // ___
340 // sqrt(2pi)^(#intObs) (-0.5 * u1T S22 u1 )
341 // I = ------------------- * e
342 // sqrt(|det(S22)|)
343 // ___
344 // Where S22 is the sub-matrix of covI for the integrated observables and S22
345 // is the Schur complement of S22
346 // ___ -1
347 // S22 = S11 - S12 * S22 * S21
348 //
349 // and u1 is the vector of non-integrated observables
350
351 // Calculate Schur complement S22bar
352 TMatrixD S22inv(S22) ;
353 S22inv.Invert() ;
354 TMatrixD S22bar = S11 - S12*S22inv*S21 ;
355
356 // Create new cache entry
357 AnaIntData& cacheData = _anaIntCache[code] ;
358 cacheData.S22bar.ResizeTo(S22bar) ;
359 cacheData.S22bar=S22bar ;
360 cacheData.S22det= S22.Determinant() ;
361 cacheData.pmap = map1 ;
362 cacheData.nint = map2.size() ;
363
364 return cacheData ;
365}
366
367
368
369////////////////////////////////////////////////////////////////////////////////
370/// Special case: generate all observables
371
372Int_t RooMultiVarGaussian::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
373{
374 if (directVars.getSize()==_x.getSize()) {
375 generateVars.add(directVars) ;
376 return -1 ;
377 }
378
379 Int_t nx = _x.getSize() ;
380 if (nx>127) {
381 // Warn that analytical integration is only provided for the first 127 observables
382 coutW(Integration) << "RooMultiVarGaussian::getGenerator(" << GetName() << ") WARNING: p.d.f. has " << _x.getSize()
383 << " observables, partial internal generation is only implemented for the first 127 observables" << endl ;
384 nx=127 ;
385 }
386
387 // Advertise partial generation over all permutations of observables
388 Int_t code(0) ;
389 BitBlock bits ;
390 for (int i=0 ; i<_x.getSize() ; i++) {
391 RooAbsArg* arg = directVars.find(_x.at(i)->GetName()) ;
392 if (arg) {
393 bits.setBit(i) ;
394// code |= (1<<i) ;
395 generateVars.add(*arg) ;
396 }
397 }
398
399 // Map BitBlock into return code
400 for (UInt_t i=0 ; i<_aicMap.size() ; i++) {
401 if (_aicMap[i]==bits) {
402 code = i+1 ;
403 }
404 }
405 if (code==0) {
406 _aicMap.push_back(bits) ;
407 code = _aicMap.size() ;
408 }
409
410
411 return code ;
412}
413
414
415
416////////////////////////////////////////////////////////////////////////////////
417/// Clear the GenData cache as its content is not invariant under changes in
418/// the mu vector.
419
421{
422 _genCache.clear() ;
423
424}
425
426
427
428
429////////////////////////////////////////////////////////////////////////////////
430/// Retrieve generator config from cache
431
433{
434 GenData& gd = genData(code) ;
435 TMatrixD& TU = gd.UT ;
436 Int_t nobs = TU.GetNcols() ;
437 vector<int>& omap = gd.omap ;
438
439 while(1) {
440
441 // Create unit Gaussian vector
442 TVectorD xgen(nobs);
443 for(Int_t k= 0; k <nobs; k++) {
444 xgen(k)= RooRandom::gaussian();
445 }
446
447 // Apply transformation matrix
448 xgen *= TU ;
449
450 // Apply shift
451 if (code == -1) {
452
453 // Simple shift if we generate all observables
454 xgen += gd.mu1 ;
455
456 } else {
457
458 // Non-generated observable dependent shift for partial generations
459
460 // mubar = mu1 + S12 S22Inv ( x2 - mu2)
461 TVectorD mubar(gd.mu1) ;
462 TVectorD x2(gd.pmap.size()) ;
463 for (UInt_t i=0 ; i<gd.pmap.size() ; i++) {
464 x2(i) = ((RooAbsReal*)_x.at(gd.pmap[i]))->getVal() ;
465 }
466 mubar += gd.S12S22I * (x2 - gd.mu2) ;
467
468 xgen += mubar ;
469
470 }
471
472 // Transfer values and check if values are in range
473 Bool_t ok(kTRUE) ;
474 for (int i=0 ; i<nobs ; i++) {
475 RooRealVar* xi = (RooRealVar*)_x.at(omap[i]) ;
476 if (xgen(i)<xi->getMin() || xgen(i)>xi->getMax()) {
477 ok = kFALSE ;
478 break ;
479 } else {
480 xi->setVal(xgen(i)) ;
481 }
482 }
483
484 // If all values are in range, accept event and return
485 // otherwise retry
486 if (ok) {
487 break ;
488 }
489 }
490
491 return;
492}
493
494
495
496////////////////////////////////////////////////////////////////////////////////
497/// WVE -- CHECK THAT GENDATA IS VALID GIVEN CURRENT VALUES OF _MU
498
500{
501 // Check if cache entry was previously created
502 map<int,GenData>::iterator iter = _genCache.find(code) ;
503 if (iter != _genCache.end()) {
504 return iter->second ;
505 }
506
507 // Create new entry
508 GenData& cacheData = _genCache[code] ;
509
510 if (code==-1) {
511
512 // Do eigen value decomposition
513 TDecompChol tdc(_cov) ;
514 tdc.Decompose() ;
515 TMatrixD U = tdc.GetU() ;
517
518 // Fill cache data
519 cacheData.UT.ResizeTo(TU) ;
520 cacheData.UT = TU ;
521 cacheData.omap.resize(_x.getSize()) ;
522 for (int i=0 ; i<_x.getSize() ; i++) {
523 cacheData.omap[i] = i ;
524 }
525 syncMuVec() ;
526 cacheData.mu1.ResizeTo(_muVec) ;
527 cacheData.mu1 = _muVec ;
528
529 } else {
530
531 // Construct observables: map1 = generated, map2 = given
532 vector<int> map1, map2 ;
533 decodeCode(code,map2,map1) ;
534
535 // Do block decomposition of covariance matrix
536 TMatrixDSym S11, S22 ;
537 TMatrixD S12, S21 ;
538 blockDecompose(_cov,map1,map2,S11,S12,S21,S22) ;
539
540 // Constructed conditional matrix form
541 // -1
542 // F(X1|X2) --> CovI --> S22bar = S11 - S12 S22 S21
543 // -1
544 // --> mu --> mubar = mu1 + S12 S22 ( x2 - mu2)
545
546 // Do eigenvalue decomposition
547 TMatrixD S22Inv(TMatrixD::kInverted,S22) ;
548 TMatrixD S22bar = S11 - S12 * (S22Inv * S21) ;
549
550 // Do eigen value decomposition of S22bar
551 TDecompChol tdc(S22bar) ;
552 tdc.Decompose() ;
553 TMatrixD U = tdc.GetU() ;
555
556 // Split mu vector into mu1 and mu2
557 TVectorD mu1(map1.size()),mu2(map2.size()) ;
558 syncMuVec() ;
559 for (UInt_t i=0 ; i<map1.size() ; i++) {
560 mu1(i) = _muVec(map1[i]) ;
561 }
562 for (UInt_t i=0 ; i<map2.size() ; i++) {
563 mu2(i) = _muVec(map2[i]) ;
564 }
565
566 // Calculate rotation matrix for mu vector
567 TMatrixD S12S22Inv = S12 * S22Inv ;
568
569 // Fill cache data
570 cacheData.UT.ResizeTo(TU) ;
571 cacheData.UT = TU ;
572 cacheData.omap = map1 ;
573 cacheData.pmap = map2 ;
574 cacheData.mu1.ResizeTo(mu1) ;
575 cacheData.mu2.ResizeTo(mu2) ;
576 cacheData.mu1 = mu1 ;
577 cacheData.mu2 = mu2 ;
578 cacheData.S12S22I.ResizeTo(S12S22Inv) ;
579 cacheData.S12S22I = S12S22Inv ;
580
581 }
582
583
584 return cacheData ;
585}
586
587
588
589
590////////////////////////////////////////////////////////////////////////////////
591/// Decode analytical integration/generation code into index map of integrated/generated (map2)
592/// and non-integrated/generated observables (map1)
593
594void RooMultiVarGaussian::decodeCode(Int_t code, vector<int>& map1, vector<int>& map2) const
595{
596 if (code<0 || code> (Int_t)_aicMap.size()) {
597 cout << "RooMultiVarGaussian::decodeCode(" << GetName() << ") ERROR don't have bit pattern for code " << code << endl ;
598 throw string("RooMultiVarGaussian::decodeCode() ERROR don't have bit pattern for code") ;
599 }
600
601 BitBlock b = _aicMap[code-1] ;
602 map1.clear() ;
603 map2.clear() ;
604 for (int i=0 ; i<_x.getSize() ; i++) {
605 if (b.getBit(i)) {
606 map2.push_back(i) ;
607 } else {
608 map1.push_back(i) ;
609 }
610 }
611}
612
613
614////////////////////////////////////////////////////////////////////////////////
615/// Block decomposition of covI according to given maps of observables
616
617void RooMultiVarGaussian::blockDecompose(const TMatrixD& input, const vector<int>& map1, const vector<int>& map2, TMatrixDSym& S11, TMatrixD& S12, TMatrixD& S21, TMatrixDSym& S22)
618{
619 // Allocate and fill reordered covI matrix in 2x2 block structure
620
621 S11.ResizeTo(map1.size(),map1.size()) ;
622 S12.ResizeTo(map1.size(),map2.size()) ;
623 S21.ResizeTo(map2.size(),map1.size()) ;
624 S22.ResizeTo(map2.size(),map2.size()) ;
625
626 for (UInt_t i=0 ; i<map1.size() ; i++) {
627 for (UInt_t j=0 ; j<map1.size() ; j++)
628 S11(i,j) = input(map1[i],map1[j]) ;
629 for (UInt_t j=0 ; j<map2.size() ; j++)
630 S12(i,j) = input(map1[i],map2[j]) ;
631 }
632 for (UInt_t i=0 ; i<map2.size() ; i++) {
633 for (UInt_t j=0 ; j<map1.size() ; j++)
634 S21(i,j) = input(map2[i],map1[j]) ;
635 for (UInt_t j=0 ; j<map2.size() ; j++)
636 S22(i,j) = input(map2[i],map2[j]) ;
637 }
638
639}
640
641
643{
644 if (ibit<32) { b0 |= (1<<ibit) ; return ; }
645 if (ibit<64) { b1 |= (1<<(ibit-32)) ; return ; }
646 if (ibit<96) { b2 |= (1<<(ibit-64)) ; return ; }
647 if (ibit<128) { b3 |= (1<<(ibit-96)) ; return ; }
648}
649
651{
652 if (ibit<32) return (b0 & (1<<ibit)) ;
653 if (ibit<64) return (b1 & (1<<(ibit-32))) ;
654 if (ibit<96) return (b2 & (1<<(ibit-64))) ;
655 if (ibit<128) return (b3 & (1<<(ibit-96))) ;
656 return kFALSE ;
657}
658
660{
661 if (b0 != other.b0) return kFALSE ;
662 if (b1 != other.b1) return kFALSE ;
663 if (b2 != other.b2) return kFALSE ;
664 if (b3 != other.b3) return kFALSE ;
665 return kTRUE ;
666}
667
668
669
670
671
#define b(i)
Definition: RSha256.hxx:100
static const double x2[5]
#define cxcoutD(a)
Definition: RooMsgService.h:82
#define coutW(a)
Definition: RooMsgService.h:33
int Int_t
Definition: RtypesCore.h:41
unsigned int UInt_t
Definition: RtypesCore.h:42
const Bool_t kFALSE
Definition: RtypesCore.h:88
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
const Bool_t kTRUE
Definition: RtypesCore.h:87
#define ClassImp(name)
Definition: Rtypes.h:365
char name[80]
Definition: TGX11.cxx:109
double pow(double, double)
double sqrt(double)
double exp(double)
char * Form(const char *fmt,...)
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:71
virtual TObject * Clone(const char *newname=0) const
Make a clone of an object using the Streamer facility.
Definition: RooAbsArg.h:83
Int_t getSize() const
Bool_t contains(const RooAbsArg &var) const
virtual Bool_t remove(const RooAbsArg &var, Bool_t silent=kFALSE, Bool_t matchByNameOnly=kFALSE)
Remove the specified argument from our list.
RooAbsArg * find(const char *name) const
Find object with given name in list.
virtual Double_t getMax(const char *name=0) const
Get maximum of currently defined range.
void setConstant(Bool_t value=kTRUE)
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:59
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:87
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:74
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
virtual Bool_t add(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling add() for each element in the source coll...
Definition: RooArgSet.h:88
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Definition: RooFitResult.h:40
const RooArgList & floatParsFinal() const
Definition: RooFitResult.h:110
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Reimplementation of standard RooArgList::add()
virtual Bool_t addOwned(RooAbsArg &var, Bool_t silent=kFALSE)
Reimplementation of standard RooArgList::addOwned()
Bool_t operator==(const BitBlock &other)
Multivariate Gaussian p.d.f.
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
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
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.
std::map< int, GenData > _genCache
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Handle full integral here.
void generateEvent(Int_t code)
Retrieve generator config from cache.
GenData & genData(Int_t code) const
WVE – CHECK THAT GENDATA IS VALID GIVEN CURRENT VALUES OF _MU.
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, Bool_t staticInitOK=kTRUE) const
Special case: generate all observables.
std::map< int, AnaIntData > _anaIntCache
void initGenerator(Int_t code)
Clear the GenData cache as its content is not invariant under changes in the mu vector.
Double_t evaluate() const
Do not persist.
static Double_t gaussian(TRandom *generator=randomGenerator())
Return a Gaussian random variable with mean 0 and variance 1.
Definition: RooRandom.cxx:111
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:278
Cholesky Decomposition class.
Definition: TDecompChol.h:25
const TMatrixD & GetU() const
Definition: TDecompChol.h:45
virtual Bool_t Decompose()
Matrix A is decomposed in component U so that A = U^T * U If the decomposition succeeds,...
Int_t GetNcols() const
Definition: TMatrixTBase.h:127
virtual Double_t Determinant() const
TMatrixTSym< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant Notice that the LU decomposition is used instead of B...
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1)
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1)
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
Definition: TMatrixT.cxx:1210
TMatrixT< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant.
Definition: TMatrixT.cxx:1396
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition: TVectorT.cxx:292
Int_t GetNrows() const
Definition: TVectorT.h:75
Double_t x[n]
Definition: legend1.C:17
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
@ Integration
Definition: RooGlobalFunc.h:67
RooConstVar & RooConst(Double_t val)
static constexpr double pi