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 "RooAbsPdf.h"
16#include "RooErrorHandler.h"
17#include "RooArgSet.h"
18#include "RooMsgService.h"
19#include "RooFitResult.h"
20#include "TMatrixDSym.h"
21#include "RooDataHist.h"
22#include "RooNumIntConfig.h"
23#include "RooRealVar.h"
24#include "RooHelpers.h"
25
26using namespace std;
27
29
30using namespace RooFit;
31
32////////////////////////////////////////////////////////////////////////////////
33/// Construct a new JeffreysPrior.
34/// \param[in] name Name of this object.
35/// \param[in] title Title (for plotting)
36/// \param[in] nominal The PDF to base Jeffrey's prior on.
37/// \param[in] paramSet Parameters of the PDF.
38/// \param[in] obsSet Observables of the PDF.
39
40RooJeffreysPrior::RooJeffreysPrior(const char* name, const char* title,
41 RooAbsPdf& nominal,
42 const RooArgList& paramSet,
43 const RooArgList& obsSet) :
44 RooAbsPdf(name, title),
45 _nominal("nominal","nominal",this, nominal, false, false),
46 _obsSet("!obsSet","Observables",this, false, false),
47 _paramSet("!paramSet","Parameters",this),
48 _cacheMgr(this, 1, true, false)
49{
50 for (const auto comp : obsSet) {
51 if (!dynamic_cast<RooAbsReal*>(comp)) {
52 coutE(InputArguments) << "RooJeffreysPrior::ctor(" << GetName() << ") ERROR: component " << comp->GetName()
53 << " in observable list is not of type RooAbsReal" << endl ;
55 }
56 _obsSet.add(*comp) ;
57 }
58
59 for (const auto comp : paramSet) {
60 if (!dynamic_cast<RooAbsReal*>(comp)) {
61 coutE(InputArguments) << "RooJeffreysPrior::ctor(" << GetName() << ") ERROR: component " << comp->GetName()
62 << " in parameter list is not of type RooAbsReal" << endl ;
64 }
65 _paramSet.add(*comp) ;
66 }
67
68 // use a different integrator by default.
69 if(paramSet.getSize()==1)
70 this->specialIntegratorConfig(true)->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D") ;
71}
72
73////////////////////////////////////////////////////////////////////////////////
74/// Copy constructor
75
77 RooAbsPdf(other, name),
78 _nominal("!nominal",this,other._nominal),
79 _obsSet("!obsSet",this,other._obsSet),
80 _paramSet("!paramSet",this,other._paramSet),
81 _cacheMgr(this, 1, true, false)
82{
83
84}
85
86////////////////////////////////////////////////////////////////////////////////
87/// Destructor
88
90{
91
92}
93
94////////////////////////////////////////////////////////////////////////////////
95/// Calculate and return current value of self
96
98{
100
101
102 CacheElem* cacheElm = (CacheElem*) _cacheMgr.getObj(nullptr);
103 if (!cacheElm) {
104 //Internally, we have to enlarge the range of fit parameters to make
105 //fits converge even if we are close to the limit of a parameter. Therefore, we clone the pdf and its
106 //observables here. If something happens to the external PDF, the cache is wiped,
107 //and we start to clone again.
108 auto& pdf = _nominal.arg();
109 RooAbsPdf* clonePdf = static_cast<RooAbsPdf*>(pdf.cloneTree());
110 auto vars = clonePdf->getParameters(_obsSet);
111 for (auto varTmp : *vars) {
112 auto& var = static_cast<RooRealVar&>(*varTmp);
113 auto range = var.getRange();
114 double span = range.second - range.first;
115 var.setRange(range.first - 0.1*span, range.second + 0.1 * span);
116 }
117
118 cacheElm = new CacheElem;
119 cacheElm->_pdf.reset(clonePdf);
120 cacheElm->_pdfVariables.reset(vars);
121
122 _cacheMgr.setObj(nullptr, cacheElm);
123 }
124
125 auto& cachedPdf = *cacheElm->_pdf;
126 auto& pdfVars = *cacheElm->_pdfVariables;
127 pdfVars.assign(_paramSet);
128
129 std::unique_ptr<RooDataHist> data( cachedPdf.generateBinned(_obsSet,ExpectedData()) );
130 std::unique_ptr<RooFitResult> res( cachedPdf.fitTo(*data, Save(),PrintLevel(-1),Minos(false),SumW2Error(false)) );
131 TMatrixDSym cov = res->covarianceMatrix();
132 cov.Invert();
133
134 return sqrt(cov.Determinant());
135}
136
#define coutE(a)
Definition: RooMsgService.h:37
#define ClassImp(name)
Definition: Rtypes.h:375
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
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:541
Int_t getSize() const
Return the number of elements in the collection.
std::pair< double, double > getRange(const char *name=nullptr) 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:62
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
Int_t setObj(const RooArgSet *nset, T *obj, const TNamed *isetRangeName=nullptr)
Setter function without integration set.
T * getObj(const RooArgSet *nset, Int_t *sterileIndex=nullptr, const TNamed *isetRangeName=nullptr)
Getter function without integration set.
bool setLabel(const char *label, bool printError=true) override
Set value by specifying the name of the desired state.
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...
static void softAbort()
Soft abort function that interrupts macro execution but doesn't kill ROOT.
Switches the message service to a different level while the instance is alive.
Definition: RooHelpers.h:42
Implementation of Jeffrey's prior.
RooTemplateProxy< RooAbsPdf > _nominal
~RooJeffreysPrior() override
Destructor.
double evaluate() const override
Calculate and return current value of self.
RooObjCacheManager _cacheMgr
RooListProxy _obsSet
RooListProxy _paramSet
RooCategory & method1D()
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:40
const T & arg() const
Return reference to object held in proxy.
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:47
RooCmdArg Save(bool flag=true)
RooCmdArg SumW2Error(bool flag)
RooCmdArg Minos(bool flag=true)
RooCmdArg PrintLevel(Int_t code)
RooCmdArg ExpectedData(bool flag=true)
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:62
std::unique_ptr< RooAbsPdf > _pdf
std::unique_ptr< RooArgSet > _pdfVariables