ROOT   6.10/09 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 *
11  * *
12  * Redistribution and use in source and binary forms, *
13  * with or without modification, are permitted according to the terms *
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 {
57
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) ;
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()) ;
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 {
116
117  for (Int_t i=0 ; i<mu.GetNrows() ; 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 {
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
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) {
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 ;
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 ;
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()) {
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) ;
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
const int nx
Definition: kalman.C:16
virtual Double_t getMin(const char *name=0) const
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
std::map< int, GenData > _genCache
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition: TVectorT.cxx:292
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, Bool_t staticInitOK=kTRUE) const
Special case: generate all observables.
const RooArgList & floatParsFinal() const
Definition: RooFitResult.h:108
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 Double_t getMax(const char *name=0) const
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:86
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported...
static Double_t gaussian(TRandom *generator=randomGenerator())
Return a Gaussian random variable with mean 0 and variance 1.
Definition: RooRandom.cxx:111
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.
const double pi
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
#define cxcoutD(a)
Definition: RooMsgService.h:79
Int_t GetNcols() const
Definition: TMatrixTBase.h:125
Int_t GetNrows() const
Definition: TVectorT.h:75
Bool_t operator==(const BitBlock &other)
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
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...
virtual TObject * Clone(const char *newname=0) const
Make a clone of an object using the Streamer facility.
Definition: RooAbsArg.h:75
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Handle full integral here.
STL namespace.
#define coutW(a)
Definition: RooMsgService.h:33
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
Double_t evaluate() const
Do not persist.
std::vector< BitBlock > _aicMap
double sqrt(double)
std::map< int, AnaIntData > _anaIntCache
static const double x2[5]
Double_t x[n]
Definition: legend1.C:17
double pow(double, double)
void initGenerator(Int_t code)
Clear the GenData cache as its content is not invariant under changes in the mu vector.
const TMatrixD & GetU() const
Definition: TDecompChol.h:45
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
TMatrixT< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant.
Definition: TMatrixT.cxx:1396
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
virtual void setVal(Double_t value)
Set value of variable to &#39;value&#39;.
Definition: RooRealVar.cxx:205
Int_t getSize() const
Cholesky Decomposition class.
Definition: TDecompChol.h:24
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
virtual Double_t Determinant() const
RooAbsArg * at(Int_t idx) const
Definition: RooArgList.h:84
void setConstant(Bool_t value=kTRUE)
void generateEvent(Int_t code)
Retrieve generator config from cache.
unsigned int UInt_t
Definition: RtypesCore.h:42
char * Form(const char *fmt,...)
GenData & genData(Int_t code) const
WVE – CHECK THAT GENDATA IS VALID GIVEN CURRENT VALUES OF _MU.
const Bool_t kFALSE
Definition: RtypesCore.h:92
Multivariate Gaussian p.d.f.
#define ClassImp(name)
Definition: Rtypes.h:336
RooAbsArg * find(const char *name) const
Find object with given name in list.
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
virtual Bool_t addOwned(RooAbsArg &var, Bool_t silent=kFALSE)
virtual Bool_t Decompose()
Matrix A is decomposed in component U so that A = U^T * U If the decomposition succeeds, bit kDecomposed is set , otherwise kSingular.
AnaIntData & anaIntData(Int_t code) const
Check if cache entry was previously created.
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 Bool_t remove(const RooAbsArg &var, Bool_t silent=kFALSE, Bool_t matchByNameOnly=kFALSE)
Remove the specified argument from our list.
RooConstVar & RooConst(Double_t val)
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
Bool_t contains(const RooAbsArg &var) const
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:66
double exp(double)
const Bool_t kTRUE
Definition: RtypesCore.h:91