diff --git a/roofit/roofit/inc/RooJeffreysPrior.h b/roofit/roofit/inc/RooJeffreysPrior.h
index 8326079a1b199ded312b4fd02484ac7061998dbe..7c5177db01dcb7521b817a88ffaba2191ec634c3 100644
--- a/roofit/roofit/inc/RooJeffreysPrior.h
+++ b/roofit/roofit/inc/RooJeffreysPrior.h
@@ -17,7 +17,7 @@ class RooArgList ;
class RooJeffreysPrior : public RooAbsPdf {
public:
- RooJeffreysPrior() ;
+ RooJeffreysPrior() { };
RooJeffreysPrior(const char *name, const char *title, RooAbsPdf& nominal, const RooArgList& paramSet, const RooArgList& obsSet) ;
virtual ~RooJeffreysPrior() ;
@@ -27,22 +27,31 @@ public:
const RooArgList& lowList() const { return _obsSet ; }
const RooArgList& paramList() const { return _paramSet ; }
- Int_t getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName=0) const ;
- Double_t analyticalIntegral(Int_t code, const char* rangeName=0) const ;
-
protected:
- RooRealProxy _nominal; // The nominal value
- //RooAbsPdf* _nominal; // The nominal value
- RooArgList _ownedList ; // List of owned components
- RooListProxy _obsSet ; // Low-side variation
- RooListProxy _paramSet ; // interpolation parameters
- mutable TIterator* _paramIter ; //! Iterator over paramSet
- mutable TIterator* _obsIter ; //! Iterator over lowSet
+ RooPdfProxy _nominal; // Proxy to the PDF for this prior.
+ RooListProxy _obsSet ; // Observables of the PDF.
+ RooListProxy _paramSet ; // Parameters of the PDF.
Double_t evaluate() const;
- ClassDef(RooJeffreysPrior,1) // Sum of RooAbsReal objects
+private:
+ struct CacheElem : public RooAbsCacheElement {
+ public:
+ virtual ~CacheElem() = default;
+ // Payload
+ std::unique_ptr<RooAbsPdf> _pdf;
+ std::unique_ptr<RooArgSet> _pdfVariables;
+
+ virtual RooArgList containedArgs(Action) override {
+ RooArgList list(*_pdf);
+ list.add(*_pdfVariables, true);
+ return list;
+ }
+ };
+ mutable RooObjCacheManager _cacheMgr; //!
+
+ ClassDef(RooJeffreysPrior,2) // Sum of RooAbsReal objects
};
#endif
diff --git a/roofit/roofit/src/RooJeffreysPrior.cxx b/roofit/roofit/src/RooJeffreysPrior.cxx
index 7095346ec3dcdfc4a832874d672c91b00b454999..b604ac11a8f4c3cf3ada3e3fe1b4c2ad942123ee 100644
--- a/roofit/roofit/src/RooJeffreysPrior.cxx
+++ b/roofit/roofit/src/RooJeffreysPrior.cxx
@@ -1,17 +1,17 @@
/** \class RooJeffreysPrior
\ingroup Roofit
-RooJeffreysPrior
-**/
-
+Implementation of Jeffrey's prior. This class estimates the fisher information matrix by generating
+a binned Asimov dataset from the supplied PDFs, fitting it, retrieving the covariance matrix and inverting
+it. It returns the square root of the determinant of this matrix.
+Numerical integration is used to normalise. Since each integration step requires fits to be run,
+evaluating complicated PDFs may take long.
-#include "RooFit.h"
-
-#include "Riostream.h"
-#include "Riostream.h"
-#include <math.h>
+Check the tutorial rs302_JeffreysPriorDemo.C for a demonstration with a simple PDF.
+**/
#include "RooJeffreysPrior.h"
+
#include "RooAbsReal.h"
#include "RooAbsPdf.h"
#include "RooErrorHandler.h"
@@ -23,8 +23,7 @@ RooJeffreysPrior
#include "RooFitResult.h"
#include "RooNumIntConfig.h"
#include "RooRealVar.h"
-
-#include "TError.h"
+#include "RooHelpers.h"
using namespace std;
@@ -33,54 +32,44 @@ ClassImp(RooJeffreysPrior);
using namespace RooFit;
////////////////////////////////////////////////////////////////////////////////
-///_obsSet("!obsSet","obs-side variation",this),
+/// Construct a new JeffreysPrior.
+/// \param[in] name Name of this object.
+/// \param[in] title Title (for plotting)
+/// \param[in] nominal The PDF to base Jeffrey's prior on.
+/// \param[in] paramSet Parameters of the PDF.
+/// \param[in] obsSet Observables of the PDF.
RooJeffreysPrior::RooJeffreysPrior(const char* name, const char* title,
RooAbsPdf& nominal,
const RooArgList& paramSet,
const RooArgList& obsSet) :
RooAbsPdf(name, title),
- _nominal("nominal","nominal",this,nominal,kFALSE,kFALSE),
- _obsSet("!obsSet","obs-side variation",this,kFALSE,kFALSE),
- _paramSet("!paramSet","high-side variation",this)
+ _nominal("nominal","nominal",this, nominal, false, false),
+ _obsSet("!obsSet","Observables",this, false, false),
+ _paramSet("!paramSet","Parameters",this),
+ _cacheMgr(this, 1, true, false)
{
- _obsIter = _obsSet.createIterator() ;
- _paramIter = _paramSet.createIterator() ;
-
-
- TIterator* inputIter1 = obsSet.createIterator() ;
- RooAbsArg* comp ;
- while((comp = (RooAbsArg*)inputIter1->Next())) {
+ for (const auto comp : obsSet) {
if (!dynamic_cast<RooAbsReal*>(comp)) {
coutE(InputArguments) << "RooJeffreysPrior::ctor(" << GetName() << ") ERROR: component " << comp->GetName()
- << " in first list is not of type RooAbsReal" << endl ;
+ << " in observable list is not of type RooAbsReal" << endl ;
RooErrorHandler::softAbort() ;
}
_obsSet.add(*comp) ;
- // if (takeOwnership) {
- // _ownedList.addOwned(*comp) ;
- // }
}
- delete inputIter1 ;
- TIterator* inputIter3 = paramSet.createIterator() ;
- while((comp = (RooAbsArg*)inputIter3->Next())) {
+ for (const auto comp : paramSet) {
if (!dynamic_cast<RooAbsReal*>(comp)) {
coutE(InputArguments) << "RooJeffreysPrior::ctor(" << GetName() << ") ERROR: component " << comp->GetName()
- << " in first list is not of type RooAbsReal" << endl ;
+ << " in parameter list is not of type RooAbsReal" << endl ;
RooErrorHandler::softAbort() ;
}
_paramSet.add(*comp) ;
- // if (takeOwnership) {
- // _ownedList.addOwned(*comp) ;
- // }
}
- delete inputIter3 ;
// use a different integrator by default.
if(paramSet.getSize()==1)
this->specialIntegratorConfig(kTRUE)->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D") ;
-
}
////////////////////////////////////////////////////////////////////////////////
@@ -90,21 +79,9 @@ RooJeffreysPrior::RooJeffreysPrior(const RooJeffreysPrior& other, const char* na
RooAbsPdf(other, name),
_nominal("!nominal",this,other._nominal),
_obsSet("!obsSet",this,other._obsSet),
- _paramSet("!paramSet",this,other._paramSet)
-{
- _obsIter = _obsSet.createIterator() ;
- _paramIter = _paramSet.createIterator() ;
-
- // Member _ownedList is intentionally not copy-constructed -- ownership is not transferred
-}
-
-////////////////////////////////////////////////////////////////////////////////
-/// Default constructor
-
-RooJeffreysPrior::RooJeffreysPrior()
+ _paramSet("!paramSet",this,other._paramSet),
+ _cacheMgr(this, 1, true, false)
{
- _obsIter = NULL;
- _paramIter = NULL;
}
@@ -113,8 +90,7 @@ RooJeffreysPrior::RooJeffreysPrior()
RooJeffreysPrior::~RooJeffreysPrior()
{
- if (_obsIter) delete _obsIter ;
- if (_paramIter) delete _paramIter ;
+
}
////////////////////////////////////////////////////////////////////////////////
@@ -122,65 +98,41 @@ RooJeffreysPrior::~RooJeffreysPrior()
Double_t RooJeffreysPrior::evaluate() const
{
- RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();
- RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);
- // create Asimov dataset
- // _paramSet.Print("v");
- // return sqrt(_paramSet.getRealValue("mu"));
- *(_nominal.arg().getVariables()) = _paramSet;
- /*
- cout << "_______________"<<endl;
- _paramSet.Print("v");
- _nominal->getVariables()->Print("v");
- cout << "_______________"<<endl;
- */
- RooDataHist* data = ((RooAbsPdf&)(_nominal.arg())).generateBinned(_obsSet,ExpectedData());
- // data->Print("v");
- //RooFitResult* res = _nominal->fitTo(*data, Save(),PrintLevel(-1),Minos(kFALSE),SumW2Error(kTRUE));
- RooFitResult* res = ((RooAbsPdf&)(_nominal.arg())).fitTo(*data, Save(),PrintLevel(-1),Minos(kFALSE),SumW2Error(kFALSE));
- TMatrixDSym cov = res->covarianceMatrix();
- cov.Invert();
- double ret = sqrt(cov.Determinant());
-
- /*
- // for 1 parameter can avoid making TMatrix etc.
- // but number of params may be > 1 with others held constant
- if(_paramSet.getSize()==1){
- RooRealVar* var = (RooRealVar*) _paramSet.first();
- // also, the _paramSet proxy one does not pick up a different value
- cout << "eval at "<< ret << " " << 1/(var->getError()) << endl;
- // need to get the actual variable instance out of the pdf like below
- var = (RooRealVar*) _nominal->getVariables()->find(var->GetName());
- cout << "eval at "<< ret << " " << 1/(var->getError()) << endl;
- }
- */
+ RooHelpers::LocalChangeMsgLevel msgLvlRAII(RooFit::WARNING);
+
+
+ CacheElem* cacheElm = (CacheElem*) _cacheMgr.getObj(nullptr);
+ if (!cacheElm) {
+ //Internally, we have to enlarge the range of fit parameters to make
+ //fits converge even if we are close to the limit of a parameter. Therefore, we clone the pdf and its
+ //observables here. If something happens to the external PDF, the cache is wiped,
+ //and we start to clone again.
+ auto& pdf = _nominal.arg();
+ RooAbsPdf* clonePdf = static_cast<RooAbsPdf*>(pdf.cloneTree());
+ auto vars = clonePdf->getParameters(_obsSet);
+ for (auto varTmp : *vars) {
+ auto& var = static_cast<RooRealVar&>(*varTmp);
+ auto range = var.getRange();
+ double span = range.second - range.first;
+ var.setRange(range.first - 0.1*span, range.second + 0.1 * span);
+ }
- // res->Print();
- delete data;
- delete res;
- RooMsgService::instance().setGlobalKillBelow(msglevel);
+ cacheElm = new CacheElem;
+ cacheElm->_pdf.reset(clonePdf);
+ cacheElm->_pdfVariables.reset(vars);
- // cout << "eval at "<< ret << endl;
- // _paramSet.Print("v");
- return ret;
+ _cacheMgr.setObj(nullptr, cacheElm);
+ }
-}
+ auto& cachedPdf = *cacheElm->_pdf;
+ auto& pdfVars = *cacheElm->_pdfVariables;
+ pdfVars = _paramSet;
-////////////////////////////////////////////////////////////////////////////////
-/// if (matchArgs(allVars,analVars,x)) return 1 ;
-/// if (matchArgs(allVars,analVars,mean)) return 2 ;
-/// return 1;
+ std::unique_ptr<RooDataHist> data( cachedPdf.generateBinned(_obsSet,ExpectedData()) );
+ std::unique_ptr<RooFitResult> res( cachedPdf.fitTo(*data, Save(),PrintLevel(-1),Minos(false),SumW2Error(false)) );
+ TMatrixDSym cov = res->covarianceMatrix();
+ cov.Invert();
-Int_t RooJeffreysPrior::getAnalyticalIntegral(RooArgSet& /*allVars*/, RooArgSet& /*analVars*/, const char* /*rangeName*/) const
-{
- return 0 ;
+ return sqrt(cov.Determinant());
}
-////////////////////////////////////////////////////////////////////////////////
-
-Double_t RooJeffreysPrior::analyticalIntegral(Int_t code, const char* /*rangeName*/) const
-{
- R__ASSERT(code==1 );
- //cout << "evaluating analytic integral" << endl;
- return 1.;
-}
diff --git a/tutorials/roostats/JeffreysPriorDemo.C b/tutorials/roostats/rs302_JeffreysPriorDemo.C
similarity index 86%
rename from tutorials/roostats/JeffreysPriorDemo.C
rename to tutorials/roostats/rs302_JeffreysPriorDemo.C
index 0904ad8d5a1e027134c540992ea9484e42c01a49..c2f720a5f92e901a4a7dc711d794c5c0dcc79638 100644
--- a/tutorials/roostats/JeffreysPriorDemo.C
+++ b/tutorials/roostats/rs302_JeffreysPriorDemo.C
@@ -36,14 +36,16 @@
#include "RooWorkspace.h"
#include "RooDataHist.h"
#include "RooGenericPdf.h"
-#include "TCanvas.h"
#include "RooPlot.h"
#include "RooFitResult.h"
-#include "TMatrixDSym.h"
#include "RooRealVar.h"
#include "RooAbsPdf.h"
#include "RooNumIntConfig.h"
+
#include "TH1F.h"
+#include "TCanvas.h"
+#include "TLegend.h"
+#include "TMatrixDSym.h"
using namespace RooFit;
@@ -71,22 +73,24 @@ void JeffreysPriorDemo()
RooJeffreysPrior pi("jeffreys", "jeffreys", *w.pdf("p"), *w.set("poi"), *w.set("obs"));
- RooGenericPdf *test = new RooGenericPdf("test", "test", "1./sqrt(mu)", *w.set("poi"));
+ RooGenericPdf *test = new RooGenericPdf("Expected", "Expected = 1/#sqrt{#mu}", "1./sqrt(mu)",
+ *w.set("poi"));
TCanvas *c1 = new TCanvas;
- // The method to compute the Jeffreys prior becomes unstable at the boundaries.
- // Therefore, we don't plot it all the way.
- RooPlot *plot = w.var("mu")->frame(Range(2, 199));
+ RooPlot *plot = w.var("mu")->frame();
pi.plotOn(plot);
- test->plotOn(plot, LineColor(kRed));
+ test->plotOn(plot, LineColor(kRed), LineStyle(kDashDotted));
plot->Draw();
+
+ auto legend = plot->BuildLegend();
+ legend->DrawClone();
}
//_________________________________________________
void TestJeffreysGaussMean()
{
RooWorkspace w("w");
- w.factory("Gaussian::g(x[0,-20,20],mu[0,-5,5],sigma[1,0,10])");
+ w.factory("Gaussian::g(x[0,-20,20],mu[0,-5.,5],sigma[1,0,10])");
w.factory("n[10,.1,200]");
w.factory("ExtendPdf::p(g,n)");
w.var("sigma")->setConstant();
@@ -113,13 +117,16 @@ void TestJeffreysGaussMean()
pi.getParameters(*temp)->Print();
// return;
- RooGenericPdf *test = new RooGenericPdf("test", "test", "1", *w.set("poi"));
+ RooGenericPdf *test = new RooGenericPdf("constant", "Expected = constant", "1", *w.set("poi"));
TCanvas *c1 = new TCanvas;
RooPlot *plot = w.var("mu")->frame();
pi.plotOn(plot);
- test->plotOn(plot, LineColor(kRed), LineStyle(kDotted));
+ test->plotOn(plot, LineColor(kRed), LineStyle(kDashDotted));
plot->Draw();
+
+ auto legend = plot->BuildLegend();
+ legend->DrawClone();
}
//_________________________________________________
@@ -159,13 +166,16 @@ void TestJeffreysGaussSigma()
const RooArgSet *temp = w.set("poi");
pi.getParameters(*temp)->Print();
- RooGenericPdf *test = new RooGenericPdf("test", "test", "sqrt(2.)/sigma", *w.set("poi"));
+ RooGenericPdf *test = new RooGenericPdf("test", "Expected = #sqrt{2}/#sigma", "sqrt(2.)/sigma", *w.set("poi"));
TCanvas *c1 = new TCanvas;
RooPlot *plot = w.var("sigma")->frame();
pi.plotOn(plot);
- test->plotOn(plot, LineColor(kRed), LineStyle(kDotted));
+ test->plotOn(plot, LineColor(kRed), LineStyle(kDashDotted));
plot->Draw();
+
+ auto legend = plot->BuildLegend();
+ legend->DrawClone();
}
//_________________________________________________
@@ -177,7 +187,7 @@ void TestJeffreysGaussMeanAndSigma()
// if the Gaussian is too wide range(x) ~ sigma then PDF gets renormalized
// and the PDF falls off too fast at high sigma
RooWorkspace w("w");
- w.factory("Gaussian::g(x[0,-20,20],mu[0,-5,5],sigma[1,1,5])");
+ w.factory("Gaussian::g(x[0,-20,20],mu[0,-5,5],sigma[1,1.,5.])");
w.factory("n[100,.1,2000]");
w.factory("ExtendPdf::p(g,n)");
@@ -206,6 +216,6 @@ void TestJeffreysGaussMeanAndSigma()
// return;
TCanvas *c1 = new TCanvas;
- TH1 *Jeff2d = pi.createHistogram("2dJeffreys", *w.var("mu"), Binning(10), YVar(*w.var("sigma"), Binning(10)));
+ TH1 *Jeff2d = pi.createHistogram("2dJeffreys", *w.var("mu"), Binning(10, -5., 5), YVar(*w.var("sigma"), Binning(10, 1., 5.)));
Jeff2d->Draw("surf");
}