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 
22 Multivariate 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 
40 using namespace std;
41 
43  ;
44 
45 ////////////////////////////////////////////////////////////////////////////////
46 
47 RooMultiVarGaussian::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 
60  _det = _cov.Determinant() ;
61 
62  // Invert covariance matrix
63  _covI.Invert() ;
64 }
65 
66 
67 ////////////////////////////////////////////////////////////////////////////////
68 
69 RooMultiVarGaussian::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 {
78  _det = _cov.Determinant() ;
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 
106 RooMultiVarGaussian::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 
129 RooMultiVarGaussian::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++) {
141  _mu.add(RooFit::RooConst(0)) ;
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 
195 Int_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 
289 Double_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 
372 Int_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 
594 void 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 
617 void 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 
RooAbsArg::Clone
virtual TObject * Clone(const char *newname=0) const
Make a clone of an object using the Streamer facility.
Definition: RooAbsArg.h:85
RooMultiVarGaussian::GenData::omap
std::vector< int > omap
Definition: RooMultiVarGaussian.h:65
TVectorT::GetNrows
Int_t GetNrows() const
Definition: TVectorT.h:81
RooMultiVarGaussian.h
TMatrixTSym::ResizeTo
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1)
RooRealVar::setVal
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:226
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:91
RooAbsReal.h
RooMultiVarGaussian::generateEvent
void generateEvent(Int_t code)
Retrieve generator config from cache.
Definition: RooMultiVarGaussian.cxx:432
RooAbsRealLValue::getMax
virtual Double_t getMax(const char *name=0) const
Get maximum of currently defined range.
Definition: RooAbsRealLValue.h:83
RooFit.h
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:364
Form
char * Form(const char *fmt,...)
RooMultiVarGaussian::BitBlock::b1
Int_t b1
Definition: RooMultiVarGaussian.h:81
RooMultiVarGaussian::analyticalIntegral
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Handle full integral here.
Definition: RooMultiVarGaussian.cxx:289
RooMultiVarGaussian
Definition: RooMultiVarGaussian.h:31
coutW
#define coutW(a)
Definition: RooMsgService.h:32
RooArgList
Definition: RooArgList.h:21
RooMultiVarGaussian::evaluate
Double_t evaluate() const
Do not persist.
Definition: RooMultiVarGaussian.cxx:176
RooMultiVarGaussian::getAnalyticalIntegral
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
Definition: RooMultiVarGaussian.cxx:195
exp
double exp(double)
RooAbsReal::getVal
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
RooMultiVarGaussian::BitBlock::setBit
void setBit(Int_t ibit)
Definition: RooMultiVarGaussian.cxx:642
RooAbsCollection::remove
virtual Bool_t remove(const RooAbsArg &var, Bool_t silent=kFALSE, Bool_t matchByNameOnly=kFALSE)
Remove the specified argument from our list.
Definition: RooAbsCollection.cxx:584
RooAbsCollection::find
RooAbsArg * find(const char *name) const
Find object with given name in list.
Definition: RooAbsCollection.cxx:813
x
Double_t x[n]
Definition: legend1.C:17
RooArgSet::add
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
RooAbsReal
Definition: RooAbsReal.h:61
TMatrixTSym
Definition: TMatrixDSymfwd.h:22
RooMultiVarGaussian::decodeCode
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...
Definition: RooMultiVarGaussian.cxx:594
RooArgList::at
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:88
TVectorT::ResizeTo
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition: TVectorT.cxx:294
RooMultiVarGaussian::BitBlock::b0
Int_t b0
Definition: RooMultiVarGaussian.h:80
TMatrixT
Definition: TMatrixDfwd.h:22
b
#define b(i)
Definition: RSha256.hxx:118
RooMultiVarGaussian::GenData::UT
TMatrixD UT
Definition: RooMultiVarGaussian.h:64
bool
RooAbsRealLValue::getMin
virtual Double_t getMin(const char *name=0) const
Get miniminum of currently defined range.
Definition: RooAbsRealLValue.h:80
RooFitResult
Definition: RooFitResult.h:40
TDecompChol
Definition: TDecompChol.h:24
RooMultiVarGaussian::GenData::pmap
std::vector< int > pmap
Definition: RooMultiVarGaussian.h:66
RooMultiVarGaussian::GenData::mu1
TVectorD mu1
Definition: RooMultiVarGaussian.h:67
TMatrixT::kInverted
@ kInverted
Definition: TMatrixT.h:58
TDecompChol.h
RooMultiVarGaussian::_anaIntCache
std::map< int, AnaIntData > _anaIntCache
Definition: RooMultiVarGaussian.h:94
ROOT::Math::fabs
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Definition: UnaryOperators.h:131
RooAbsCollection::contains
Bool_t contains(const RooAbsArg &var) const
Check if collection contains an argument with the same name as var.
Definition: RooAbsCollection.h:102
RooMultiVarGaussian::_covI
TMatrixDSym _covI
Definition: RooMultiVarGaussian.h:102
RooMultiVarGaussian::RooMultiVarGaussian
RooMultiVarGaussian()
Definition: RooMultiVarGaussian.h:34
TDecompChol::GetU
const TMatrixD & GetU() const
Definition: TDecompChol.h:51
RooMultiVarGaussian::_x
RooListProxy _x
Definition: RooMultiVarGaussian.h:99
TMatrixT::kTransposed
@ kTransposed
Definition: TMatrixT.h:58
RooMultiVarGaussian::anaIntData
AnaIntData & anaIntData(Int_t code) const
Check if cache entry was previously created.
Definition: RooMultiVarGaussian.cxx:318
TMatrixT::Invert
TMatrixT< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant.
Definition: TMatrixT.cxx:1397
RooMultiVarGaussian::GenData
Definition: RooMultiVarGaussian.h:62
RooMultiVarGaussian::genData
GenData & genData(Int_t code) const
WVE – CHECK THAT GENDATA IS VALID GIVEN CURRENT VALUES OF _MU.
Definition: RooMultiVarGaussian.cxx:499
TMatrixT::ResizeTo
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1)
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
RooRandom.h
RooMultiVarGaussian::_muVec
TVectorD _muVec
Definition: RooMultiVarGaussian.h:107
RooMultiVarGaussian::getGenerator
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, Bool_t staticInitOK=kTRUE) const
Special case: generate all observables.
Definition: RooMultiVarGaussian.cxx:372
RooMultiVarGaussian::BitBlock
Definition: RooMultiVarGaussian.h:72
RooMultiVarGaussian::BitBlock::getBit
Bool_t getBit(Int_t ibit)
Definition: RooMultiVarGaussian.cxx:650
RooMultiVarGaussian::_det
Double_t _det
Definition: RooMultiVarGaussian.h:103
RooMultiVarGaussian::BitBlock::b2
Int_t b2
Definition: RooMultiVarGaussian.h:82
RooListProxy::addOwned
virtual Bool_t addOwned(RooAbsArg &var, Bool_t silent=kFALSE)
Reimplementation of standard RooArgList::addOwned()
Definition: RooListProxy.cxx:114
RooListProxy::add
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Reimplementation of standard RooArgList::add()
Definition: RooListProxy.cxx:104
sqrt
double sqrt(double)
RooMultiVarGaussian::AnaIntData
Definition: RooMultiVarGaussian.h:54
RooMultiVarGaussian::GenData::S12S22I
TMatrixD S12S22I
Definition: RooMultiVarGaussian.h:69
RooMultiVarGaussian::_aicMap
std::vector< BitBlock > _aicMap
Definition: RooMultiVarGaussian.h:97
RooRealVar.h
RooMultiVarGaussian::AnaIntData::pmap
std::vector< int > pmap
Definition: RooMultiVarGaussian.h:58
RooFitResult.h
unsigned int
RooMultiVarGaussian::GenData::mu2
TVectorD mu2
Definition: RooMultiVarGaussian.h:68
RooConstVar.h
RooMultiVarGaussian::AnaIntData::S22det
Double_t S22det
Definition: RooMultiVarGaussian.h:57
RooMultiVarGaussian::BitBlock::operator==
Bool_t operator==(const BitBlock &other)
Definition: RooMultiVarGaussian.cxx:659
RooGlobalFunc.h
RooMultiVarGaussian::AnaIntData::S22bar
TMatrixD S22bar
Definition: RooMultiVarGaussian.h:56
TVectorT
Definition: TMatrixTBase.h:78
Double_t
double Double_t
Definition: RtypesCore.h:59
TGeant4Unit::pi
static constexpr double pi
Definition: TGeant4SystemOfUnits.h:73
TMatrixTSym::Invert
TMatrixTSym< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant Notice that the LU decomposition is used instead of B...
Definition: TMatrixTSym.cxx:961
RooRandom::gaussian
static Double_t gaussian(TRandom *generator=randomGenerator())
Return a Gaussian random variable with mean 0 and variance 1.
Definition: RooRandom.cxx:110
RooMultiVarGaussian::_mu
RooListProxy _mu
Definition: RooMultiVarGaussian.h:100
RooMultiVarGaussian::blockDecompose
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.
Definition: RooMultiVarGaussian.cxx:617
name
char name[80]
Definition: TGX11.cxx:110
RooAbsRealLValue::setConstant
void setConstant(Bool_t value=kTRUE)
Definition: RooAbsRealLValue.h:109
x2
static const double x2[5]
Definition: RooGaussKronrodIntegrator1D.cxx:364
RooMultiVarGaussian::syncMuVec
void syncMuVec() const
Definition: RooMultiVarGaussian.cxx:164
RooFitResult::floatParsFinal
const RooArgList & floatParsFinal() const
Definition: RooFitResult.h:110
RooAbsArg
Definition: RooAbsArg.h:73
TMatrixTSym::Determinant
virtual Double_t Determinant() const
Definition: TMatrixTSym.cxx:935
RooMultiVarGaussian::AnaIntData::nint
Int_t nint
Definition: RooMultiVarGaussian.h:59
RooAbsPdf
Definition: RooAbsPdf.h:40
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:53
RooMultiVarGaussian::initGenerator
void initGenerator(Int_t code)
Clear the GenData cache as its content is not invariant under changes in the mu vector.
Definition: RooMultiVarGaussian.cxx:420
RooMultiVarGaussian::_genCache
std::map< int, GenData > _genCache
Definition: RooMultiVarGaussian.h:95
RooFit::Integration
@ Integration
Definition: RooGlobalFunc.h:67
RooRealVar
Definition: RooRealVar.h:35
Riostream.h
pow
double pow(double, double)
TDecompChol::Decompose
virtual Bool_t Decompose()
Matrix A is decomposed in component U so that A = U^T * U If the decomposition succeeds,...
Definition: TDecompChol.cxx:106
RooMath.h
RooMultiVarGaussian::_cov
TMatrixDSym _cov
Definition: RooMultiVarGaussian.h:101
RooAbsCollection::getSize
Int_t getSize() const
Definition: RooAbsCollection.h:171
RooArgSet
Definition: RooArgSet.h:28
cxcoutD
#define cxcoutD(a)
Definition: RooMsgService.h:81
RooMultiVarGaussian::_z
Double_t _z
Definition: RooMultiVarGaussian.h:104
int
RooMultiVarGaussian::BitBlock::b3
Int_t b3
Definition: RooMultiVarGaussian.h:83
RooFit::RooConst
RooConstVar & RooConst(Double_t val)
Definition: RooGlobalFunc.cxx:341