Logo ROOT  
Reference Guide
RooJeffreysPrior.cxx
Go to the documentation of this file.
1/** \class RooJeffreysPrior
2\ingroup Roofit
3
4Implementation of Jeffrey's prior. This class estimates the fisher information matrix by generating
5a binned Asimov dataset from the supplied PDFs, fitting it, retrieving the covariance matrix and inverting
6it. It returns the square root of the determinant of this matrix.
7Numerical integration is used to normalise. Since each integration step requires fits to be run,
8evaluating complicated PDFs may take long.
9
10Check the tutorial rs302_JeffreysPriorDemo.C for a demonstration with a simple PDF.
11**/
12
13#include "RooJeffreysPrior.h"
14
15#include "RooAbsReal.h"
16#include "RooAbsPdf.h"
17#include "RooErrorHandler.h"
18#include "RooArgSet.h"
19#include "RooMsgService.h"
20#include "RooFitResult.h"
21#include "TMatrixDSym.h"
22#include "RooDataHist.h"
23#include "RooNumIntConfig.h"
24#include "RooRealVar.h"
25#include "RooHelpers.h"
26
27using namespace std;
28
30
31using namespace RooFit;
32
33////////////////////////////////////////////////////////////////////////////////
34/// Construct a new JeffreysPrior.
35/// \param[in] name Name of this object.
36/// \param[in] title Title (for plotting)
37/// \param[in] nominal The PDF to base Jeffrey's prior on.
38/// \param[in] paramSet Parameters of the PDF.
39/// \param[in] obsSet Observables of the PDF.
40
41RooJeffreysPrior::RooJeffreysPrior(const char* name, const char* title,
42 RooAbsPdf& nominal,
43 const RooArgList& paramSet,
44 const RooArgList& obsSet) :
45 RooAbsPdf(name, title),
46 _nominal("nominal","nominal",this, nominal, false, false),
47 _obsSet("!obsSet","Observables",this, false, false),
48 _paramSet("!paramSet","Parameters",this),
49 _cacheMgr(this, 1, true, false)
50{
51 for (const auto comp : obsSet) {
52 if (!dynamic_cast<RooAbsReal*>(comp)) {
53 coutE(InputArguments) << "RooJeffreysPrior::ctor(" << GetName() << ") ERROR: component " << comp->GetName()
54 << " in observable list is not of type RooAbsReal" << endl ;
56 }
57 _obsSet.add(*comp) ;
58 }
59
60 for (const auto comp : paramSet) {
61 if (!dynamic_cast<RooAbsReal*>(comp)) {
62 coutE(InputArguments) << "RooJeffreysPrior::ctor(" << GetName() << ") ERROR: component " << comp->GetName()
63 << " in parameter list is not of type RooAbsReal" << endl ;
65 }
66 _paramSet.add(*comp) ;
67 }
68
69 // use a different integrator by default.
70 if(paramSet.getSize()==1)
71 this->specialIntegratorConfig(kTRUE)->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D") ;
72}
73
74////////////////////////////////////////////////////////////////////////////////
75/// Copy constructor
76
78 RooAbsPdf(other, name),
79 _nominal("!nominal",this,other._nominal),
80 _obsSet("!obsSet",this,other._obsSet),
81 _paramSet("!paramSet",this,other._paramSet),
82 _cacheMgr(this, 1, true, false)
83{
84
85}
86
87////////////////////////////////////////////////////////////////////////////////
88/// Destructor
89
91{
92
93}
94
95////////////////////////////////////////////////////////////////////////////////
96/// Calculate and return current value of self
97
99{
101
102
103 CacheElem* cacheElm = (CacheElem*) _cacheMgr.getObj(nullptr);
104 if (!cacheElm) {
105 //Internally, we have to enlarge the range of fit parameters to make
106 //fits converge even if we are close to the limit of a parameter. Therefore, we clone the pdf and its
107 //observables here. If something happens to the external PDF, the cache is wiped,
108 //and we start to clone again.
109 auto& pdf = _nominal.arg();
110 RooAbsPdf* clonePdf = static_cast<RooAbsPdf*>(pdf.cloneTree());
111 auto vars = clonePdf->getParameters(_obsSet);
112 for (auto varTmp : *vars) {
113 auto& var = static_cast<RooRealVar&>(*varTmp);
114 auto range = var.getRange();
115 double span = range.second - range.first;
116 var.setRange(range.first - 0.1*span, range.second + 0.1 * span);
117 }
118
119 cacheElm = new CacheElem;
120 cacheElm->_pdf.reset(clonePdf);
121 cacheElm->_pdfVariables.reset(vars);
122
123 _cacheMgr.setObj(nullptr, cacheElm);
124 }
125
126 auto& cachedPdf = *cacheElm->_pdf;
127 auto& pdfVars = *cacheElm->_pdfVariables;
128 pdfVars.assign(_paramSet);
129
130 std::unique_ptr<RooDataHist> data( cachedPdf.generateBinned(_obsSet,ExpectedData()) );
131 std::unique_ptr<RooFitResult> res( cachedPdf.fitTo(*data, Save(),PrintLevel(-1),Minos(false),SumW2Error(false)) );
132 TMatrixDSym cov = res->covarianceMatrix();
133 cov.Invert();
134
135 return sqrt(cov.Determinant());
136}
137
#define coutE(a)
Definition: RooMsgService.h:33
const Bool_t kTRUE
Definition: RtypesCore.h:100
#define ClassImp(name)
Definition: Rtypes.h:364
char name[80]
Definition: TGX11.cxx:110
RooArgSet * getParameters(const RooAbsData *data, bool stripDisconnected=true) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don't match any of...
Definition: RooAbsArg.cxx:578
Int_t getSize() const
friend class CacheElem
The cache manager.
Definition: RooAbsPdf.h:377
std::pair< double, double > getRange(const char *name=0) const
Get low and high bound of the variable.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:63
RooNumIntConfig * specialIntegratorConfig() const
Returns the specialized integrator configuration for this RooAbsReal.
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:22
T * getObj(const RooArgSet *nset, Int_t *sterileIndex=0, const TNamed *isetRangeName=0)
Int_t setObj(const RooArgSet *nset, T *obj, const TNamed *isetRangeName=0)
virtual Bool_t setLabel(const char *label, bool printError=true) override
Set value by specifying the name of the desired state.
static void softAbort()
Switches the message service to a different level while the instance is alive.
Definition: RooHelpers.h:42
Implementation of Jeffrey's prior.
virtual ~RooJeffreysPrior()
Destructor.
RooTemplateProxy< RooAbsPdf > _nominal
Double_t evaluate() const
Calculate and return current value of self.
RooObjCacheManager _cacheMgr
RooListProxy _obsSet
RooListProxy _paramSet
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE) override
Reimplementation of standard RooArgList::add()
RooCategory & method1D()
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:39
const T & arg() const
Return reference to object held in proxy.
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 const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
RooCmdArg SumW2Error(Bool_t flag)
RooCmdArg Save(Bool_t flag=kTRUE)
RooCmdArg PrintLevel(Int_t code)
RooCmdArg Minos(Bool_t flag=kTRUE)
RooCmdArg ExpectedData(Bool_t flag=kTRUE)
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition: Common.h:18
@ InputArguments
Definition: RooGlobalFunc.h:61
std::unique_ptr< RooAbsPdf > _pdf
std::unique_ptr< RooArgSet > _pdfVariables