Skip to content

Commit

Permalink
[RF] Add RooWrapperPdf.
Browse files Browse the repository at this point in the history
RooFit contains a number of functions that cannot be used as PDFs since
they don't have automatic normalisation. When wrapped into the wrapper
PDF, functions can be used in the same way as PDFs.
  • Loading branch information
hageboeck committed Sep 13, 2019
1 parent 0cfc1b6 commit 2712117
Show file tree
Hide file tree
Showing 6 changed files with 267 additions and 1 deletion.
4 changes: 3 additions & 1 deletion roofit/roofitcore/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -232,6 +232,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore
RooSpan.h
BatchData.h
RooVDTHeaders.h
RooWrapperPdf.h
SOURCES
src/BidirMMapPipe.cxx
src/BidirMMapPipe.h
Expand Down Expand Up @@ -450,7 +451,8 @@ ROOT_STANDARD_LIBRARY_PACKAGE(RooFitCore
src/RooXYChi2Var.cxx
src/RooHelpers.cxx
src/BatchData.cxx
DICTIONARY_OPTIONS
src/RooWrapperPdf.cxx
DICTIONARY_OPTIONS
"-writeEmptyRootPCM"
DEPENDENCIES
Core
Expand Down
1 change: 1 addition & 0 deletions roofit/roofitcore/inc/LinkDef1.h
Original file line number Diff line number Diff line change
Expand Up @@ -90,6 +90,7 @@
#pragma link C++ class RooEffProd+ ;
#pragma link C++ class RooExtendPdf+ ;
#pragma link off class RooErrorHandler+ ;
#pragma link C++ class RooWrapperPdf+;
#endif


109 changes: 109 additions & 0 deletions roofit/roofitcore/inc/RooWrapperPdf.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,109 @@
// Author: Stephan Hageboeck, CERN
/*****************************************************************************
* Project: RooFit *
* Package: RooFitCore *
* Authors: *
* WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
* DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
* *
* Copyright (c) 2000-2018, Regents of the University of California *
* and Stanford University. All rights reserved. *
* *
* Redistribution and use in source and binary forms, *
* with or without modification, are permitted according to the terms *
* listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
*****************************************************************************/
#ifndef ROO_WRAPPER_PDF
#define ROO_WRAPPER_PDF

#include "RooAbsReal.h"
#include "RooRealProxy.h"
#include "RooAbsPdf.h"

class RooWrapperPdf final : public RooAbsPdf {
public:

RooWrapperPdf() { };
/// Construct a new RooWrapperPdf.
/// \param[in] name A name to identify this object.
/// \param[in] title Title (for e.g. plotting)
/// \param[in] inputFunction Any RooAbsReal that should be converted into a PDF. Although it's possible
/// to pass a PDF, it only makes sense for non-PDF functions.
RooWrapperPdf(const char *name, const char *title, RooAbsReal& inputFunction) :
RooAbsPdf(name, title),
_func("inputFunction", "Function to be converted into a PDF", this, inputFunction) { }
virtual ~RooWrapperPdf() {};

RooWrapperPdf(const RooWrapperPdf& other, const char* name = 0) :
RooAbsPdf(other, name),
_func("inputFunction", this, other._func) { }

virtual TObject* clone(const char* newname) const override {
return new RooWrapperPdf(*this, newname);
}

// Analytical Integration handling
Bool_t forceAnalyticalInt(const RooAbsArg& dep) const override {
return _func.arg().forceAnalyticalInt(dep);
}
Int_t getAnalyticalIntegralWN(RooArgSet& allVars, RooArgSet& analVars, const RooArgSet* normSet,
const char* rangeName=0) const override {
return _func.arg().getAnalyticalIntegralWN(allVars, analVars, normSet, rangeName);
}
Int_t getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& numVars,
const char* rangeName=0) const override {
return _func.arg().getAnalyticalIntegral(allVars, numVars, rangeName);
}
double analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* rangeName) const override {
return _func.arg().analyticalIntegralWN(code, normSet, rangeName);
}
double analyticalIntegral(Int_t code, const char* rangeName=0) const override {
return _func.arg().analyticalIntegral(code, rangeName);
}


// Internal toy generation. Since our _func is not a PDF (if it is, it doesn't make sense to use this wrapper),
// we cannot do anything.
/// Get specialised generator. Since the underlying function is not a PDF, this will always return zero.
// Int_t getGenerator(const RooArgSet& /*directVars*/, RooArgSet& /*generateVars*/,
// bool /*staticInitOK = true*/) const override { return 0; }
// void initGenerator(Int_t /*code*/) override { }
// void generateEvent(Int_t /*code*/) override { }
// Bool_t isDirectGenSafe(const RooAbsArg& /*arg*/) const override { return false; }


// Hints for optimized brute-force sampling
Int_t getMaxVal(const RooArgSet& vars) const override {
return _func.arg().getMaxVal(vars);
}
Double_t maxVal(Int_t code) const override {
return _func.arg().maxVal(code);
}
Int_t minTrialSamples(const RooArgSet& arGenObs) const override {
return _func.arg().minTrialSamples(arGenObs);
}

// Plotting and binning hints
Bool_t isBinnedDistribution(const RooArgSet& obs) const override {
return _func.arg().isBinnedDistribution(obs);
}
std::list<Double_t>* binBoundaries(RooAbsRealLValue& obs, Double_t xlo, Double_t xhi) const override {
return _func.arg().binBoundaries(obs, xlo, xhi);
}
std::list<Double_t>* plotSamplingHint(RooAbsRealLValue& obs, Double_t xlo, Double_t xhi) const override {
return _func.arg().plotSamplingHint(obs, xlo, xhi);
}



private:
RooRealProxy _func;

double evaluate() const override {
return _func;
}

ClassDefOverride(RooWrapperPdf,1)
};

#endif
34 changes: 34 additions & 0 deletions roofit/roofitcore/src/RooWrapperPdf.cxx
Original file line number Diff line number Diff line change
@@ -0,0 +1,34 @@
// Author: Stephan Hageboeck, CERN 12 Sep 2019

/*****************************************************************************
* RooFit
* Authors: *
* WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
* DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
* *
* Copyright (c) 2000-2019, Regents of the University of California *
* and Stanford University. All rights reserved. *
* *
* Redistribution and use in source and binary forms, *
* with or without modification, are permitted according to the terms *
* listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
*****************************************************************************/

/**
* \class RooWrapperPdf
* The RooWrapperPdf is a class that can be used to convert a function into a PDF.
*
* During construction, a RooAbsReal has to be passed. When this function is evaluated, the wrapper pdf will
* in addition evaluate its integral, and normalise the returned value. It will further ensure that negative
* return values are clipped at zero.
*
* Functions calls such as analytical integral requests or plot sampling hints are simply forwarded to the RooAbsReal
* that was passed in the constructor.
*/


#include "RooWrapperPdf.h"

ClassImp(RooWrapperPdf)


1 change: 1 addition & 0 deletions roofit/roofitcore/test/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -9,3 +9,4 @@
ROOT_ADD_GTEST(simple simple.cxx LIBRARIES RooFitCore)
ROOT_ADD_GTEST(testWorkspace testWorkspace.cxx LIBRARIES RooFitCore RooFit RooStats)
ROOT_ADD_GTEST(testRooDataHist testRooDataHist.cxx LIBRARIES RooFitCore)
ROOT_ADD_GTEST(testRooWrapperPdf testRooWrapperPdf.cxx LIBRARIES Gpad RooFitCore)
119 changes: 119 additions & 0 deletions roofit/roofitcore/test/testRooWrapperPdf.cxx
Original file line number Diff line number Diff line change
@@ -0,0 +1,119 @@
// Tests for the RooWrapperPdf
// Author: Stephan Hageboeck, CERN 09/2019

#include "RooWrapperPdf.h"
#include "RooRealVar.h"
#include "RooProduct.h"
#include "RooRealSumFunc.h"
#include "RooConstVar.h"
#include "RooPlot.h"
#include "RooFitResult.h"
#include "RooDataSet.h"

#include "TCanvas.h"

#include "gtest/gtest.h"

TEST(RooWrapperPdf, Basics)
{
RooRealVar x("x", "x", -5., 5.);

// Implement a poor-man's polynomial. Value ranges are chosen to keep it positive.
RooRealVar a0("a0", "a0", 0.1, 0.1, 10.);
RooRealVar a1("a1", "a1", -0.01, -2.1, 0.);
RooRealVar a2("a2", "a2", 0.01, 0.01, 5.);
RooProduct xSq("xSq", "x^2", RooArgList(x, x));
RooConstVar one("one", "one", 1.);
RooRealSumFunc pol("pol", "pol", RooArgList(one, x, xSq), RooArgList(a0, a1, a2));

RooWrapperPdf polPdf("polPdf", "polynomial PDF", pol);

EXPECT_GT(pol.getVal(x)*1.05, polPdf.getVal(x)) << "Wrapper pdf normalises.";

RooArgSet intSet(x);
RooArgSet numSet;

EXPECT_NE(polPdf.getAnalyticalIntegralWN(intSet, numSet, &intSet, nullptr), 0)
<< "Test that PDF claims to have analytic integral with norm.";

// auto frame = x.frame();
// pol.plotOn(frame);
// polPdf.plotOn(frame);
// TCanvas canv;
// frame->Draw();
// canv.SaveAs("/tmp/testWrapperPdf.png");
}


TEST(RooWrapperPdf, GenerateAndFit) {
RooRealVar x("x", "x", -5., 5.);

// Implement a poor-man's polynomial. Value ranges are chosen to keep it positive.
RooRealVar a0("a0", "a0", 0.1);
RooRealVar a1("a1", "a1", -0.01);
RooRealVar a2("a2", "a2", 0.01, 0.001, 5.);
RooProduct xSq("xSq", "x^2", RooArgList(x, x));
RooConstVar one("one", "one", 1.);
RooRealSumFunc pol("pol", "pol", RooArgList(one, x, xSq), RooArgList(a0, a1, a2));

RooWrapperPdf polPdf("polPdf", "polynomial PDF", pol);

auto data = polPdf.generate(x, 50000);
a2.setVal(0.02);
auto result = polPdf.fitTo(*data, RooFit::Save());

EXPECT_EQ(result->status(), 0) << "Fit converged.";
EXPECT_LT(fabs(a2.getVal()-0.01), a2.getError());

// auto frame = x.frame();
// data->plotOn(frame);
// polPdf.plotOn(frame);
// TCanvas canv;
// frame->Draw();
// canv.SaveAs("/tmp/testWrapperPdf2.png");
}


TEST(RooWrapperPdf, DISABLED_FullAnalyticInt) {
RooRealVar x("x", "x", 4., 0., 10.);
RooRealVar y("y", "y", -0.5, -5., 5.);

RooProduct xy("xy", "x*y", RooArgList(x, y));
RooWrapperPdf prodPdf("prodPdf", "PDF(x*y)", xy);

RooArgSet intSet(x);
RooArgSet numSet;

EXPECT_NE(prodPdf.getAnalyticalIntegral(intSet, numSet, nullptr), 0)
<< "Test that PDF claims to have analytic integral.";

EXPECT_FLOAT_EQ(xy.getVal(), -1.5);
std::cout << "The following error is expected:\n----" << std::endl;
EXPECT_FLOAT_EQ(prodPdf.getVal(), 0.);
std::cout << "----" << std::endl;

constexpr double newY = 2.;
y.setVal(newY);
EXPECT_NEAR(prodPdf.getVal(), 3.*newY, 0.001);
EXPECT_NEAR(prodPdf.getVal(x), 0.3, 0.001);

// auto cdf = prodPdf.createCdf(x);
// for (unsigned int i=0; i<10; ++i) {
// x.setVal(i);
// std::cout << i << " " << std::setprecision(4)
// << std::setw(7) << xy.getVal(x)
// << std::setw(7) << prodPdf.getVal(x)
// << std::setw(7) << cdf->getVal() << std::endl;
// EXPECT_NEAR(xy.getVal(x), newY*i, 1.E-5);
// EXPECT_NEAR(prodPdf.getVal(x), newY*i/10., 1.E-5);
// EXPECT_NEAR(cdf->getVal(), newY*i, 1.E-5);
// }

// auto frame = x.frame();
// prodPdf.plotOn(frame);
// cdf->plotOn(frame, RooFit::LineColor(kRed));
// TCanvas canv;
// frame->Draw();
// canv.SaveAs("/tmp/testWrapperPdf3.png");
}

0 comments on commit 2712117

Please sign in to comment.