#include <Roo1DTable.h>
#include <RooAbsReal.h>
#include <RooAddModel.h>
#include <RooAddPdf.h>
#include <RooAddition.h>
#include <RooArgSet.h>
#include <RooArgusBG.h>
#include <RooBCPEffDecay.h>
#include <RooBDecay.h>
#include <RooBMixDecay.h>
#include <RooBinning.h>
#include <RooBinningCategory.h>
#include <RooCFunction1Binding.h>
#include <RooCFunction3Binding.h>
#include <RooCategory.h>
#include <RooChebychev.h>
#include <RooChi2MCSModule.h>
#include <RooConstVar.h>
#include <RooDLLSignificanceMCSModule.h>
#include <RooDataHist.h>
#include <RooDataSet.h>
#include <RooDecay.h>
#include <RooEffProd.h>
#include <RooEfficiency.h>
#include <RooExponential.h>
#include <RooExtendPdf.h>
#include <RooFFTConvPdf.h>
#include <RooFitResult.h>
#include <RooFormulaVar.h>
#include <RooGaussModel.h>
#include <RooGaussian.h>
#include <RooGenericPdf.h>
#include <RooGlobalFunc.h>
#include <RooHelpers.h>
#include <RooHist.h>
#include <RooHistPdf.h>
#include <RooIntegralMorph.h>
#include <RooKeysPdf.h>
#include <RooLandau.h>
#include <RooMCStudy.h>
#include <RooMappedCategory.h>
#include <RooMinimizer.h>
#include <RooMultiCategory.h>
#include <RooNDKeysPdf.h>
#include <RooNumIntConfig.h>
#include <RooPlot.h>
#include <RooPolyVar.h>
#include <RooPolynomial.h>
#include <RooProdPdf.h>
#include <RooProduct.h>
#include <RooProfileLL.h>
#include <RooRandomizeParamMCSModule.h>
#include <RooRealConstant.h>
#include <RooRealSumPdf.h>
#include <RooRealVar.h>
#include <RooSimultaneous.h>
#include <RooSuperCategory.h>
#include <RooTFnBinding.h>
#include <RooThresholdCategory.h>
#include <RooTruthModel.h>
#include <RooUnitTest.h>
#include <RooWorkspace.h>
#include <Math/DistFunc.h>
#include <TCanvas.h>
#include <TDirectory.h>
#include <TF1.h>
#include <TFile.h>
#include <TH1.h>
#include <TH1D.h>
#include <TH2.h>
#include <TH2D.h>
#include <TMath.h>
#include <TPluginManager.h>
#include <TROOT.h>
#include <TRandom.h>
#include <TStyle.h>
#include <TTree.h>
#include <iomanip>
using namespace RooFit ;
//////////////////////////////////////////////////////////////////////////
//
// 'BASIC FUNCTIONALITY' RooFit tutorial macro #101
//
// Fitting, plotting, toy data generation on one-dimensional p.d.f
//
// pdf = gauss(x,m,s)
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
// Elementary operations on a gaussian PDF
class TestBasic101 : public RooUnitTest
{
public:
TestBasic101(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Fitting,plotting & event generation of basic p.d.f",refFile,writeRef,verbose) {} ;
bool testCode() {
// S e t u p m o d e l
// ---------------------
// Declare variables x,mean,sigma with associated name, title, initial value and allowed range
RooRealVar x("x","x",-10,10) ;
RooRealVar mean("mean","mean of gaussian",1,-10,10) ;
RooRealVar sigma("sigma","width of gaussian",1,0.1,10) ;
// Build gaussian p.d.f in terms of x,mean and sigma
RooGaussian gauss("gauss","gaussian PDF",x,mean,sigma) ;
// Construct plot frame in 'x'
RooPlot* xframe = x.frame(Title("Gaussian p.d.f.")) ;
// P l o t m o d e l a n d c h a n g e p a r a m e t e r v a l u e s
// ---------------------------------------------------------------------------
// Plot gauss in frame (i.e. in x)
gauss.plotOn(xframe) ;
// Change the value of sigma to 3
sigma.setVal(3) ;
// Plot gauss in frame (i.e. in x) and draw frame on canvas
gauss.plotOn(xframe,LineColor(kRed),Name("another")) ;
// G e n e r a t e e v e n t s
// -----------------------------
// Generate a dataset of 1000 events in x from gauss
std::unique_ptr<RooDataSet> data{gauss.generate(x,10000)};
// Make a second plot frame in x and draw both the
// data and the p.d.f in the frame
RooPlot* xframe2 = x.frame(Title("Gaussian p.d.f. with data")) ;
data->plotOn(xframe2) ;
gauss.plotOn(xframe2) ;
// F i t m o d e l t o d a t a
// -----------------------------
// Fit pdf to data
gauss.fitTo(*data) ;
// --- Post processing for stressRooFit ---
regPlot(xframe ,"rf101_plot1") ;
regPlot(xframe2,"rf101_plot2") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'BASIC FUNCTIONALITY' RooFit tutorial macro #102
//
// Importing data from ROOT TTrees and THx histograms
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic102 : public RooUnitTest
{
public:
TestBasic102(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Data import methods",refFile,writeRef,verbose) {} ;
std::unique_ptr<TH1> makeTH1()
{
// Create ROOT TH1 filled with a Gaussian distribution
auto hh = std::make_unique<TH1D>("hh","hh",25,-10,10);
for (int i=0 ; i<100 ; i++) {
hh->Fill(gRandom->Gaus(0,3)) ;
}
return hh ;
}
std::unique_ptr<TTree> makeTTree()
{
// Create ROOT TTree filled with a Gaussian distribution in x and a uniform distribution in y
auto tree = std::make_unique<TTree>("tree","tree");
double* px = new double ;
double* py = new double ;
tree->Branch("x",px,"x/D") ;
tree->Branch("y",py,"y/D") ;
for (int i=0 ; i<100 ; i++) {
*px = gRandom->Gaus(0,3) ;
*py = gRandom->Uniform()*30 - 15 ;
tree->Fill() ;
}
//delete px ;
//delete py ;
return tree ;
}
bool testCode() {
////////////////////////////////////////////////////////
// I m p o r t i n g R O O T h i s t o g r a m s //
////////////////////////////////////////////////////////
// I m p o r t T H 1 i n t o a R o o D a t a H i s t
// ---------------------------------------------------------
// Create a ROOT TH1 histogram
std::unique_ptr<TH1> hh{makeTH1()};
// Declare observable x
RooRealVar x("x","x",-10,10) ;
// Create a binned dataset that imports contents of TH1 and associates its contents to observable 'x'
RooDataHist dh("dh","dh",x,Import(*hh)) ;
// P l o t a n d f i t a R o o D a t a H i s t
// ---------------------------------------------------
// Make plot of binned dataset showing Poisson error bars (RooFit default)
RooPlot* frame = x.frame(Title("Imported TH1 with Poisson error bars")) ;
dh.plotOn(frame) ;
// Fit a Gaussian p.d.f to the data
RooRealVar mean("mean","mean",0,-10,10) ;
RooRealVar sigma("sigma","sigma",3,0.1,10) ;
RooGaussian gauss("gauss","gauss",x,mean,sigma) ;
gauss.fitTo(dh) ;
gauss.plotOn(frame) ;
// P l o t a n d f i t a R o o D a t a H i s t w i t h i n t e r n a l e r r o r s
// ---------------------------------------------------------------------------------------------
// If histogram has custom error (i.e. its contents is does not originate from a Poisson process
// but e.g. is a sum of weighted events) you can data with symmetric 'sum-of-weights' error instead
// (same error bars as shown by ROOT)
RooPlot* frame2 = x.frame(Title("Imported TH1 with internal errors")) ;
dh.plotOn(frame2,DataError(RooAbsData::SumW2)) ;
gauss.plotOn(frame2) ;
// Please note that error bars shown (Poisson or SumW2) are for visualization only, the are NOT used
// in a maximum likelihood fit
//
// A (binned) ML fit will ALWAYS assume the Poisson error interpretation of data (the mathematical definition
// of likelihood does not take any external definition of errors). Data with non-unit weights can only be correctly
// fitted with a chi^2 fit (see rf602_chi2fit.C)
////////////////////////////////////////////////
// I m p o r t i n g R O O T T T r e e s //
////////////////////////////////////////////////
// I m p o r t T T r e e i n t o a R o o D a t a S e t
// -----------------------------------------------------------
std::unique_ptr<TTree> tree{makeTTree()};
// Define 2nd observable y
RooRealVar y("y","y",-10,10) ;
// Construct unbinned dataset importing tree branches x and y matching between branches and RooRealVars
// is done by name of the branch/RRV
//
// Note that ONLY entries for which x,y have values within their allowed ranges as defined in
// RooRealVar x and y are imported. Since the y values in the import tree are in the range [-15,15]
// and RRV y defines a range [-10,10] this means that the RooDataSet below will have less entries than the TTree 'tree'
RooDataSet ds("ds","ds",RooArgSet(x,y),Import(*tree)) ;
// P l o t d a t a s e t w i t h m u l t i p l e b i n n i n g c h o i c e s
// ------------------------------------------------------------------------------------
// Print unbinned dataset with default frame binning (100 bins)
RooPlot* frame3 = y.frame(Title("Unbinned data shown in default frame binning")) ;
ds.plotOn(frame3) ;
// Print unbinned dataset with custom binning choice (20 bins)
RooPlot* frame4 = y.frame(Title("Unbinned data shown with custom binning")) ;
ds.plotOn(frame4,Binning(20)) ;
// Draw all frames on a canvas
regPlot(frame ,"rf102_plot1") ;
regPlot(frame2,"rf102_plot2") ;
regPlot(frame3,"rf102_plot3") ;
regPlot(frame4,"rf102_plot4") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'BASIC FUNCTIONALITY' RooFit tutorial macro #103
//
// Interpreted functions and p.d.f.s
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic103 : public RooUnitTest
{
public:
TestBasic103(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Interpreted expression p.d.f.",refFile,writeRef,verbose) {} ;
bool testCode() {
/////////////////////////////////////////////////////////
// G e n e r i c i n t e r p r e t e d p . d . f . //
/////////////////////////////////////////////////////////
// Declare observable x
RooRealVar x("x","x",-20,20) ;
// C o n s t r u c t g e n e r i c p d f f r o m i n t e r p r e t e d e x p r e s s i o n
// -------------------------------------------------------------------------------------------------
// To construct a proper p.d.f, the formula expression is explicitly normalized internally by dividing
// it by a numeric integral of the expresssion over x in the range [-20,20]
//
RooRealVar alpha("alpha","alpha",5,0.1,10) ;
RooGenericPdf genpdf("genpdf","genpdf","(1+0.1*abs(x)+sin(sqrt(abs(x*alpha+0.1))))",RooArgSet(x,alpha)) ;
// S a m p l e , f i t a n d p l o t g e n e r i c p d f
// ---------------------------------------------------------------
// Generate a toy dataset from the interpreted p.d.f
std::unique_ptr<RooDataSet> data{genpdf.generate(x,10000)};
// Fit the interpreted p.d.f to the generated data
genpdf.fitTo(*data) ;
// Make a plot of the data and the p.d.f overlaid
RooPlot* xframe = x.frame(Title("Interpreted expression pdf")) ;
data->plotOn(xframe) ;
genpdf.plotOn(xframe) ;
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// S t a n d a r d p . d . f a d j u s t w i t h i n t e r p r e t e d h e l p e r f u n c t i o n //
// //
// Make a gauss(x,sqrt(mean2),sigma) from a standard RooGaussian //
// //
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// C o n s t r u c t s t a n d a r d p d f w i t h f o r m u l a r e p l a c i n g p a r a m e t e r
// ------------------------------------------------------------------------------------------------------------
// Construct parameter mean2 and sigma
RooRealVar mean2("mean2","mean^2",10,0,200) ;
RooRealVar sigma("sigma","sigma",3,0.1,10) ;
// Construct interpreted function mean = sqrt(mean^2)
RooFormulaVar mean("mean","mean","sqrt(mean2)",mean2) ;
// Construct a gaussian g2(x,sqrt(mean2),sigma) ;
RooGaussian g2("g2","h2",x,mean,sigma) ;
// G e n e r a t e t o y d a t a
// ---------------------------------
// Construct a separate gaussian g1(x,10,3) to generate a toy Gaussian dataset with mean 10 and width 3
RooGaussian g1("g1","g1",x,10.0,3.0) ;
std::unique_ptr<RooDataSet> data2{g1.generate(x,1000)};
// F i t a n d p l o t t a i l o r e d s t a n d a r d p d f
// -------------------------------------------------------------------
// Fit g2 to data from g1
std::unique_ptr<RooFitResult> r{g2.fitTo(*data2,Save())};
// Plot data on frame and overlay projection of g2
RooPlot* xframe2 = x.frame(Title("Tailored Gaussian pdf")) ;
data2->plotOn(xframe2) ;
g2.plotOn(xframe2) ;
regPlot(xframe,"rf103_plot1") ;
regPlot(xframe2,"rf103_plot2") ;
regResult(std::move(r),"rf103_fit1") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'BASIC FUNCTIONALITY' RooFit tutorial macro #105
//
// Demonstration of binding ROOT Math functions as RooFit functions
// and pdfs
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic105 : public RooUnitTest
{
public:
TestBasic105(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("C++ function binding operator p.d.f",refFile,writeRef,verbose) {} ;
bool testCode() {
// B i n d T M a t h : : E r f C f u n c t i o n
// ---------------------------------------------------
// Bind one-dimensional TMath::Erf function as RooAbsReal function
RooRealVar x("x","x",-3,3) ;
std::unique_ptr<RooAbsReal> erf{bindFunction("erf",TMath::Erf,x)};
// Plot erf on frame
RooPlot* frame1 = x.frame(Title("TMath::Erf bound as RooFit function")) ;
erf->plotOn(frame1) ;
// B i n d R O O T : : M a t h : : b e t a _ p d f C f u n c t i o n
// -----------------------------------------------------------------------
// Bind pdf ROOT::Math::Beta with three variables as RooAbsPdf function
// exclude x=0 and x=1 points from range since beta_pdf diverges at x=0 and x=1
// when a < 1 and/or b <1
RooRealVar x2("x2","x2",0.001,0.999) ;
RooRealVar a("a","a",5,0,10) ;
RooRealVar b("b","b",2,0,10) ;
std::unique_ptr<RooAbsPdf> beta{bindPdf("beta",ROOT::Math::beta_pdf,x2,a,b)};
// Generate some events and fit
std::unique_ptr<RooDataSet> data{beta->generate(x2,10000)};
beta->fitTo(*data) ;
// Plot data and pdf on frame
RooPlot* frame2 = x2.frame(Title("ROOT::Math::Beta bound as RooFit pdf")) ;
data->plotOn(frame2) ;
beta->plotOn(frame2) ;
// B i n d R O O T T F 1 a s R o o F i t f u n c t i o n
// ---------------------------------------------------------------
// Create a ROOT TF1 function
auto fa1 = std::make_unique<TF1>("fa1","sin(x)/x",0,10);
// Create an observable
RooRealVar x3("x3","x3",0.01,20) ;
// Create binding of TF1 object to above observable
std::unique_ptr<RooAbsReal> rfa1{bindFunction(fa1.get(),x3)};
// Make plot frame in observable, plot TF1 binding function
RooPlot* frame3 = x3.frame(Title("TF1 bound as RooFit function")) ;
rfa1->plotOn(frame3) ;
regPlot(frame1,"rf105_plot1") ;
regPlot(frame2,"rf105_plot2") ;
regPlot(frame3,"rf105_plot3") ;
return true ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'BASIC FUNCTIONALITY' RooFit tutorial macro #108
//
// Plotting unbinned data with alternate and variable binnings
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic108 : public RooUnitTest
{
public:
TestBasic108(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Non-standard binning in counting and asymmetry plots",refFile,writeRef,verbose) {} ;
bool testCode() {
// S e t u p m o d e l
// ---------------------
// Build a B decay p.d.f with mixing
RooRealVar dt("dt","dt",-20,20) ;
RooRealVar dm("dm","dm",0.472) ;
RooRealVar tau("tau","tau",1.547) ;
RooRealVar w("w","mistag rate",0.1) ;
RooRealVar dw("dw","delta mistag rate",0.) ;
RooCategory mixState("mixState","B0/B0bar mixing state") ;
mixState.defineType("mixed",-1) ;
mixState.defineType("unmixed",1) ;
RooCategory tagFlav("tagFlav","Flavour of the tagged B0") ;
tagFlav.defineType("B0",1) ;
tagFlav.defineType("B0bar",-1) ;
// Build a gaussian resolution model
RooRealVar dterr("dterr","dterr",0.1,1.0) ;
RooRealVar bias1("bias1","bias1",0) ;
RooRealVar sigma1("sigma1","sigma1",0.1) ;
RooGaussModel gm1("gm1","gauss model 1",dt,bias1,sigma1) ;
// Construct Bdecay (x) gauss
RooBMixDecay bmix("bmix","decay",dt,mixState,tagFlav,tau,dm,w,dw,gm1,RooBMixDecay::DoubleSided) ;
// S a m p l e d a t a f r o m m o d e l
// --------------------------------------------
// Sample 2000 events in (dt,mixState,tagFlav) from bmix
std::unique_ptr<RooDataSet> data{bmix.generate(RooArgSet(dt,mixState,tagFlav),2000)};
// S h o w d t d i s t r i b u t i o n w i t h c u s t o m b i n n i n g
// -------------------------------------------------------------------------------
// Make plot of dt distribution of data in range (-15,15) with fine binning for dt>0 and coarse binning for dt<0
// Create binning object with range (-15,15)
RooBinning tbins(-15,15) ;
// Add 60 bins with uniform spacing in range (-15,0)
tbins.addUniform(60,-15,0) ;
// Add 15 bins with uniform spacing in range (0,15)
tbins.addUniform(15,0,15) ;
// Make plot with specified binning
RooPlot* dtframe = dt.frame(Range(-15,15),Title("dt distribution with custom binning")) ;
data->plotOn(dtframe,Binning(tbins)) ;
bmix.plotOn(dtframe) ;
// NB: Note that bin density for each bin is adjusted to that of default frame binning as shown
// in Y axis label (100 bins --> Events/0.4*Xaxis-dim) so that all bins represent a consistent density distribution
// S h o w m i x s t a t e a s y m m e t r y w i t h c u s t o m b i n n i n g
// ------------------------------------------------------------------------------------
// Make plot of dt distribution of data asymmetry in 'mixState' with variable binning
// Create binning object with range (-10,10)
RooBinning abins(-10,10) ;
// Add boundaries at 0, (-1,1), (-2,2), (-3,3), (-4,4) and (-6,6)
abins.addBoundary(0) ;
abins.addBoundaryPair(1) ;
abins.addBoundaryPair(2) ;
abins.addBoundaryPair(3) ;
abins.addBoundaryPair(4) ;
abins.addBoundaryPair(6) ;
// Create plot frame in dt
RooPlot* aframe = dt.frame(Range(-10,10),Title("mixState asymmetry distribution with custom binning")) ;
// Plot mixState asymmetry of data with specified customg binning
data->plotOn(aframe,Asymmetry(mixState),Binning(abins)) ;
// Plot corresponding property of p.d.f
bmix.plotOn(aframe,Asymmetry(mixState)) ;
// Adjust vertical range of plot to sensible values for an asymmetry
aframe->SetMinimum(-1.1) ;
aframe->SetMaximum(1.1) ;
// NB: For asymmetry distributions no density corrects are needed (and are thus not applied)
regPlot(dtframe,"rf108_plot1") ;
regPlot(aframe,"rf108_plot2") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'BASIC FUNCTIONALITY' RooFit tutorial macro #109
//
// Calculating chi^2 from histograms and curves in RooPlots,
// making histogram of residual and pull distributions
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic109 : public RooUnitTest
{
public:
TestBasic109(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Calculation of chi^2 and residuals in plots",refFile,writeRef,verbose) {} ;
bool testCode() {
// S e t u p m o d e l
// ---------------------
// Create observables
RooRealVar x("x","x",-10,10) ;
// Create Gaussian
RooRealVar sigma("sigma","sigma",3,0.1,10) ;
RooRealVar mean("mean","mean",0,-10,10) ;
RooGaussian gauss("gauss","gauss",x,0.0,sigma) ;
// Generate a sample of 1000 events with sigma=3
std::unique_ptr<RooDataSet> data{gauss.generate(x,10000)};
// Change sigma to 3.15
sigma=3.15 ;
// P l o t d a t a a n d s l i g h t l y d i s t o r t e d m o d e l
// ---------------------------------------------------------------------------
// Overlay projection of gauss with sigma=3.15 on data with sigma=3.0
RooPlot* frame1 = x.frame(Title("Data with distorted Gaussian pdf"),Bins(40)) ;
data->plotOn(frame1,DataError(RooAbsData::SumW2)) ;
gauss.plotOn(frame1) ;
// C a l c u l a t e c h i ^ 2
// ------------------------------
// Show the chi^2 of the curve w.r.t. the histogram
// If multiple curves or datasets live in the frame you can specify
// the name of the relevant curve and/or dataset in chiSquare()
regValue(frame1->chiSquare(),"rf109_chi2") ;
// S h o w r e s i d u a l a n d p u l l d i s t s
// -------------------------------------------------------
// Construct a histogram with the residuals of the data w.r.t. the curve
// we set `useAverage` to false for this test because this was done for the reference histogram
RooHist* hresid = frame1->residHist(nullptr, nullptr, false, false) ;
// Construct a histogram with the pulls of the data w.r.t the curve
// we set `useAverage` to false for this test because this was done for the reference histogram
RooHist* hpull = frame1->pullHist(nullptr, nullptr, false) ;
// Create a new frame to draw the residual distribution and add the distribution to the frame
RooPlot* frame2 = x.frame(Title("Residual Distribution")) ;
frame2->addPlotable(hresid,"P") ;
// Create a new frame to draw the pull distribution and add the distribution to the frame
RooPlot* frame3 = x.frame(Title("Pull Distribution")) ;
frame3->addPlotable(hpull,"P") ;
regPlot(frame1,"rf109_plot1") ;
regPlot(frame2,"rf109_plot2") ;
regPlot(frame3,"rf109_plot3") ;
//delete hresid ;
//delete hpull ;
return true ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'BASIC FUNCTIONALITY' RooFit tutorial macro #110
//
// Examples on normalization of p.d.f.s,
// integration of p.d.fs, construction
// of cumulative distribution functions from p.d.f.s
// in one dimension
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic110 : public RooUnitTest
{
public:
TestBasic110(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Normalization of p.d.f.s in 1D",refFile,writeRef,verbose) {} ;
bool testCode() {
// S e t u p m o d e l
// ---------------------
// Create observables x,y
RooRealVar x("x","x",-10,10) ;
// Create p.d.f. gaussx(x,-2,3)
RooGaussian gx("gx","gx",x,-2.0,3.0) ;
// R e t r i e v e r a w & n o r m a l i z e d v a l u e s o f R o o F i t p . d . f . s
// --------------------------------------------------------------------------------------------------
// Return 'raw' unnormalized value of gx
regValue(gx.getVal(),"rf110_gx") ;
// Return value of gx normalized over x in range [-10,10]
RooArgSet nset(x) ;
regValue(gx.getVal(&nset),"rf110_gx_Norm[x]") ;
// Create object representing integral over gx
// which is used to calculate gx_Norm[x] == gx / gx_Int[x]
std::unique_ptr<RooAbsReal> igx{gx.createIntegral(x)};
regValue(igx->getVal(),"rf110_gx_Int[x]") ;
// I n t e g r a t e n o r m a l i z e d p d f o v e r s u b r a n g e
// ----------------------------------------------------------------------------
// Define a range named "signal" in x from -5,5
x.setRange("signal",-5,5) ;
// Create an integral of gx_Norm[x] over x in range "signal"
// This is the fraction of of p.d.f. gx_Norm[x] which is in the
// range named "signal"
std::unique_ptr<RooAbsReal> igx_sig{gx.createIntegral(x,NormSet(x),Range("signal"))};
regValue(igx_sig->getVal(),"rf110_gx_Int[x|signal]_Norm[x]") ;
// C o n s t r u c t c u m u l a t i v e d i s t r i b u t i o n f u n c t i o n f r o m p d f
// -----------------------------------------------------------------------------------------------------
// Create the cumulative distribution function of gx
// i.e. calculate Int[-10,x] gx(x') dx'
std::unique_ptr<RooAbsReal> gx_cdf{gx.createCdf(x)};
// Plot cdf of gx versus x
RooPlot* frame = x.frame(Title("c.d.f of Gaussian p.d.f")) ;
gx_cdf->plotOn(frame) ;
regPlot(frame,"rf110_plot1") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'BASIC FUNCTIONALITY' RooFit tutorial macro #111
//
// Configuration and customization of how numeric (partial) integrals
// are executed
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic111 : public RooUnitTest
{
public:
TestBasic111(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Numeric integration configuration",refFile,writeRef,verbose) {} ;
bool testCode() {
// A d j u s t g l o b a l 1 D i n t e g r a t i o n p r e c i s i o n
// ----------------------------------------------------------------------------
// Example: Change global precision for 1D integrals from 1e-7 to 1e-6
//
// The relative epsilon (change as fraction of current best integral estimate) and
// absolute epsilon (absolute change w.r.t last best integral estimate) can be specified
// separately. For most p.d.f integrals the relative change criterium is the most important,
// however for certain non-p.d.f functions that integrate out to zero a separate absolute
// change criterium is necessary to declare convergence of the integral
//
// NB: This change is for illustration only. In general the precision should be at least 1e-7
// for normalization integrals for MINUIT to succeed.
//
RooAbsReal::defaultIntegratorConfig()->setEpsAbs(1e-6) ;
RooAbsReal::defaultIntegratorConfig()->setEpsRel(1e-6) ;
// N u m e r i c i n t e g r a t i o n o f l a n d a u p d f
// ------------------------------------------------------------------
RooRealVar x("x","x",-10,10) ;
RooLandau landau("landau","landau",x,0.0,0.1) ;
// Disable analytic integration from demonstration purposes
landau.forceNumInt(true);
// Calculate integral over landau with default choice of numeric integrator
std::unique_ptr<RooAbsReal> intLandau{landau.createIntegral(x)};
double val = intLandau->getVal() ;
regValue(val,"rf111_val1") ;
// S a m e w i t h c u s t o m c o n f i g u r a t i o n
// -----------------------------------------------------------
// Construct a custom configuration which uses the adaptive Gauss-Kronrod technique
// for closed 1D integrals
RooNumIntConfig customConfig(*RooAbsReal::defaultIntegratorConfig()) ;
#ifdef R__HAS_MATHMORE
customConfig.method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D") ;
#endif
// Calculate integral over landau with custom integral specification
std::unique_ptr<RooAbsReal> intLandau2{landau.createIntegral(x,NumIntConfig(customConfig))};
double val2 = intLandau2->getVal() ;
regValue(val2,"rf111_val2") ;
// A d j u s t i n g d e f a u l t c o n f i g f o r a s p e c i f i c p d f
// -------------------------------------------------------------------------------------
// Another possibility: associate custom numeric integration configuration as default for object 'landau'
landau.setIntegratorConfig(customConfig) ;
// Calculate integral over landau custom numeric integrator specified as object default
std::unique_ptr<RooAbsReal> intLandau3{landau.createIntegral(x)};
double val3 = intLandau3->getVal() ;
regValue(val3,"rf111_val3") ;
return true ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'ADDITION AND CONVOLUTION' RooFit tutorial macro #201
//
// Composite p.d.f with signal and background component
//
// pdf = f_bkg * bkg(x,a0,a1) + (1-fbkg) * (f_sig1 * sig1(x,m,s1 + (1-f_sig1) * sig2(x,m,s2)))
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic201 : public RooUnitTest
{
public:
TestBasic201(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Addition operator p.d.f.",refFile,writeRef,verbose) {} ;
bool testCode() {
// S e t u p c o m p o n e n t p d f s
// ---------------------------------------
// Declare observable x
RooRealVar x("x","x",0,10) ;
// Create two Gaussian PDFs g1(x,mean1,sigma) anf g2(x,mean2,sigma) and their parameters
RooRealVar mean("mean","mean of gaussians",5) ;
RooRealVar sigma1("sigma1","width of gaussians",0.5) ;
RooRealVar sigma2("sigma2","width of gaussians",1) ;
RooGaussian sig1("sig1","Signal component 1",x,mean,sigma1) ;
RooGaussian sig2("sig2","Signal component 2",x,mean,sigma2) ;
// Build Chebychev polynomial p.d.f.
RooRealVar a0("a0","a0",0.5,0.,1.) ;
RooRealVar a1("a1","a1",-0.2,-1.,1.) ;
RooChebychev bkg("bkg","Background",x,RooArgSet(a0,a1)) ;
////////////////////////////////////////////////////
// M E T H O D 1 - T w o R o o A d d P d f s //
////////////////////////////////////////////////////
// A d d s i g n a l c o m p o n e n t s
// ------------------------------------------
// Sum the signal components into a composite signal p.d.f.
RooRealVar sig1frac("sig1frac","fraction of component 1 in signal",0.8,0.,1.) ;
RooAddPdf sig("sig","Signal",RooArgList(sig1,sig2),sig1frac) ;
// A d d s i g n a l a n d b a c k g r o u n d
// ------------------------------------------------
// Sum the composite signal and background
RooRealVar bkgfrac("bkgfrac","fraction of background",0.5,0.,1.) ;
RooAddPdf model("model","g1+g2+a",RooArgList(bkg,sig),bkgfrac) ;
// S a m p l e , f i t a n d p l o t m o d e l
// ---------------------------------------------------
// Generate a data sample of 1000 events in x from model
std::unique_ptr<RooDataSet> data{model.generate(x,1000)};
// Fit model to data
model.fitTo(*data) ;
// Plot data and PDF overlaid
RooPlot* xframe = x.frame(Title("Example of composite pdf=(sig1+sig2)+bkg")) ;
data->plotOn(xframe) ;
model.plotOn(xframe) ;
// Overlay the background component of model with a dashed line
model.plotOn(xframe,Components(bkg),LineStyle(kDashed)) ;
// Overlay the background+sig2 components of model with a dotted line
model.plotOn(xframe,Components(RooArgSet(bkg,sig2)),LineStyle(kDotted)) ;
////////////////////////////////////////////////////////////////////////////////////////////////////
// M E T H O D 2 - O n e R o o A d d P d f w i t h r e c u r s i v e f r a c t i o n s //
////////////////////////////////////////////////////////////////////////////////////////////////////
// Construct sum of models on one go using recursive fraction interpretations
//
// model2 = bkg + (sig1 + sig2)
//
RooAddPdf model2("model","g1+g2+a",RooArgList(bkg,sig1,sig2),RooArgList(bkgfrac,sig1frac),true) ;
// NB: Each coefficient is interpreted as the fraction of the
// left-hand component of the i-th recursive sum, i.e.
//
// sum4 = A + ( B + ( C + D) with fraction fA, fB and fC expands to
//
// sum4 = fA*A + (1-fA)*(fB*B + (1-fB)*(fC*C + (1-fC)*D))
// P l o t r e c u r s i v e a d d i t i o n m o d e l
// ---------------------------------------------------------
model2.plotOn(xframe,LineColor(kRed),LineStyle(kDashed)) ;
model2.plotOn(xframe,Components(RooArgSet(bkg,sig2)),LineColor(kRed),LineStyle(kDashed)) ;
regPlot(xframe,"rf201_plot1") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'ADDITION AND CONVOLUTION' RooFit tutorial macro #202
//
// Setting up an extended maximum likelihood fit
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic202 : public RooUnitTest
{
public:
TestBasic202(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Extended ML fits to addition operator p.d.f.s",refFile,writeRef,verbose) {} ;
bool testCode() {
// S e t u p c o m p o n e n t p d f s
// ---------------------------------------
// Declare observable x
RooRealVar x("x","x",0,10) ;
// Create two Gaussian PDFs g1(x,mean1,sigma) anf g2(x,mean2,sigma) and their parameters
RooRealVar mean("mean","mean of gaussians",5) ;
RooRealVar sigma1("sigma1","width of gaussians",0.5) ;
RooRealVar sigma2("sigma2","width of gaussians",1) ;
RooGaussian sig1("sig1","Signal component 1",x,mean,sigma1) ;
RooGaussian sig2("sig2","Signal component 2",x,mean,sigma2) ;
// Build Chebychev polynomial p.d.f.
RooRealVar a0("a0","a0",0.5,0.,1.) ;
RooRealVar a1("a1","a1",-0.2,-1.,1.) ;
RooChebychev bkg("bkg","Background",x,RooArgSet(a0,a1)) ;
// Sum the signal components into a composite signal p.d.f.
RooRealVar sig1frac("sig1frac","fraction of component 1 in signal",0.8,0.,1.) ;
RooAddPdf sig("sig","Signal",RooArgList(sig1,sig2),sig1frac) ;
/////////////////////
// M E T H O D 1 //
/////////////////////
// C o n s t r u c t e x t e n d e d c o m p o s i t e m o d e l
// -------------------------------------------------------------------
// Sum the composite signal and background into an extended pdf nsig*sig+nbkg*bkg
RooRealVar nsig("nsig","number of signal events",500,0.,10000) ;
RooRealVar nbkg("nbkg","number of background events",500,0,10000) ;
RooAddPdf model("model","(g1+g2)+a",RooArgList(bkg,sig),RooArgList(nbkg,nsig)) ;
// S a m p l e , f i t a n d p l o t e x t e n d e d m o d e l
// ---------------------------------------------------------------------
// Generate a data sample of expected number events in x from model
// = model.expectedEvents() = nsig+nbkg
std::unique_ptr<RooDataSet> data{model.generate(x)};
// Fit model to data, extended ML term automatically included
model.fitTo(*data) ;
// Plot data and PDF overlaid, use expected number of events for p.d.f projection normalization
// rather than observed number of events (==data->numEntries())
RooPlot* xframe = x.frame(Title("extended ML fit example")) ;
data->plotOn(xframe) ;
model.plotOn(xframe,Normalization(1.0,RooAbsReal::RelativeExpected)) ;
// Overlay the background component of model with a dashed line
model.plotOn(xframe,Components(bkg),LineStyle(kDashed),Normalization(1.0,RooAbsReal::RelativeExpected)) ;
// Overlay the background+sig2 components of model with a dotted line
model.plotOn(xframe,Components(RooArgSet(bkg,sig2)),LineStyle(kDotted),Normalization(1.0,RooAbsReal::RelativeExpected)) ;
/////////////////////
// M E T H O D 2 //
/////////////////////
// C o n s t r u c t e x t e n d e d c o m p o n e n t s f i r s t
// ---------------------------------------------------------------------
// Associated nsig/nbkg as expected number of events with sig/bkg
RooExtendPdf esig("esig","extended signal p.d.f",sig,nsig) ;
RooExtendPdf ebkg("ebkg","extended background p.d.f",bkg,nbkg) ;
// S u m e x t e n d e d c o m p o n e n t s w i t h o u t c o e f s
// -------------------------------------------------------------------------
// Construct sum of two extended p.d.f. (no coefficients required)
RooAddPdf model2("model2","(g1+g2)+a",RooArgList(ebkg,esig)) ;
regPlot(xframe,"rf202_plot1") ;
return true ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'ADDITION AND CONVOLUTION' RooFit tutorial macro #203
//
// Fitting and plotting in sub ranges
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic203 : public RooUnitTest
{
public:
TestBasic203(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Basic fitting and plotting in ranges",refFile,writeRef,verbose) {} ;
bool testCode() {
// S e t u p m o d e l
// ---------------------
// Construct observables x
RooRealVar x("x","x",-10,10) ;
// Construct gaussx(x,mx,1)
RooRealVar mx("mx","mx",0,-10,10) ;
RooGaussian gx("gx","gx",x,mx,1.0) ;
// Construct px = 1 (flat in x)
RooPolynomial px("px","px",x) ;
// Construct model = f*gx + (1-f)px
RooRealVar f("f","f",0.,1.) ;
RooAddPdf model("model","model",RooArgList(gx,px),f) ;
// Generated 10000 events in (x,y) from p.d.f. model
std::unique_ptr<RooDataSet> modelData{model.generate(x,10000)};
// F i t f u l l r a n g e
// ---------------------------
// Fit p.d.f to all data
std::unique_ptr<RooFitResult> r_full{model.fitTo(*modelData,Save(true))};
// F i t p a r t i a l r a n g e
// ----------------------------------
// Define "signal" range in x as [-3,3]
x.setRange("signal",-3,3) ;
// Fit p.d.f only to data in "signal" range
std::unique_ptr<RooFitResult> r_sig{model.fitTo(*modelData,Save(true),Range("signal"))};
// P l o t / p r i n t r e s u l t s
// ---------------------------------------
// Make plot frame in x and add data and fitted model
RooPlot* frame = x.frame(Title("Fitting a sub range")) ;
modelData->plotOn(frame) ;
model.plotOn(frame,Range(""),LineStyle(kDashed),LineColor(kRed)) ; // Add shape in full ranged dashed
model.plotOn(frame) ; // By default only fitted range is shown
regPlot(frame,"rf203_plot") ;
regResult(std::move(r_full),"rf203_r_full") ;
regResult(std::move(r_sig),"rf203_r_sig") ;
return true;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'ADDITION AND CONVOLUTION' RooFit tutorial macro #204
//
// Extended maximum likelihood fit with alternate range definition
// for observed number of events.
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic204 : public RooUnitTest
{
public:
TestBasic204(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Extended ML fit in sub range",refFile,writeRef,verbose) {} ;
bool testCode() {
// S e t u p c o m p o n e n t p d f s
// ---------------------------------------
// Declare observable x
RooRealVar x("x","x",0,10) ;
// Create two Gaussian PDFs g1(x,mean1,sigma) anf g2(x,mean2,sigma) and their parameters
RooRealVar mean("mean","mean of gaussians",5) ;
RooRealVar sigma1("sigma1","width of gaussians",0.5) ;
RooRealVar sigma2("sigma2","width of gaussians",1) ;
RooGaussian sig1("sig1","Signal component 1",x,mean,sigma1) ;
RooGaussian sig2("sig2","Signal component 2",x,mean,sigma2) ;
// Build Chebychev polynomial p.d.f.
RooRealVar a0("a0","a0",0.5,0.,1.) ;
RooRealVar a1("a1","a1",-0.2,-1.,1.) ;
RooChebychev bkg("bkg","Background",x,RooArgSet(a0,a1)) ;
// Sum the signal components into a composite signal p.d.f.
RooRealVar sig1frac("sig1frac","fraction of component 1 in signal",0.8,0.,1.) ;
RooAddPdf sig("sig","Signal",RooArgList(sig1,sig2),sig1frac) ;
// C o n s t r u c t e x t e n d e d c o m p s wi t h r a n g e s p e c
// ------------------------------------------------------------------------------
// Define signal range in which events counts are to be defined
x.setRange("signalRange",4,6) ;
// Associated nsig/nbkg as expected number of events with sig/bkg _in_the_range_ "signalRange"
RooRealVar nsig("nsig","number of signal events in signalRange",500,0.,10000) ;
RooRealVar nbkg("nbkg","number of background events in signalRange",500,0,10000) ;
RooExtendPdf esig("esig","extended signal p.d.f",sig,nsig,"signalRange") ;
RooExtendPdf ebkg("ebkg","extended background p.d.f",bkg,nbkg,"signalRange") ;
// S u m e x t e n d e d c o m p o n e n t s
// ---------------------------------------------
// Construct sum of two extended p.d.f. (no coefficients required)
RooAddPdf model("model","(g1+g2)+a",RooArgList(ebkg,esig)) ;
// S a m p l e d a t a , f i t m o d e l
// -------------------------------------------
// Generate 1000 events from model so that nsig,nbkg come out to numbers <<500 in fit
std::unique_ptr<RooDataSet> data{model.generate(x,1000)};
// Perform unbinned extended ML fit to data
std::unique_ptr<RooFitResult> r{model.fitTo(*data,Extended(true),Save())};
regResult(std::move(r),"rf204_result") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'ADDITION AND CONVOLUTION' RooFit tutorial macro #205
//
// Options for plotting components of composite p.d.f.s.
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic205 : public RooUnitTest
{
public:
TestBasic205(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Component plotting variations",refFile,writeRef,verbose) {} ;
bool testCode() {
// S e t u p c o m p o s i t e p d f
// --------------------------------------
// Declare observable x
RooRealVar x("x","x",0,10) ;
// Create two Gaussian PDFs g1(x,mean1,sigma) anf g2(x,mean2,sigma) and their parameters
RooRealVar mean("mean","mean of gaussians",5) ;
RooRealVar sigma1("sigma1","width of gaussians",0.5) ;
RooRealVar sigma2("sigma2","width of gaussians",1) ;
RooGaussian sig1("sig1","Signal component 1",x,mean,sigma1) ;
RooGaussian sig2("sig2","Signal component 2",x,mean,sigma2) ;
// Sum the signal components into a composite signal p.d.f.
RooRealVar sig1frac("sig1frac","fraction of component 1 in signal",0.8,0.,1.) ;
RooAddPdf sig("sig","Signal",RooArgList(sig1,sig2),sig1frac) ;
// Build Chebychev polynomial p.d.f.
RooRealVar a0("a0","a0",0.5,0.,1.) ;
RooRealVar a1("a1","a1",-0.2,-1.,1.) ;
RooChebychev bkg1("bkg1","Background 1",x,RooArgSet(a0,a1)) ;
// Build expontential pdf
RooRealVar alpha("alpha","alpha", 1) ;
RooExponential bkg2("bkg2","Background 2",x,alpha, true);
// Sum the background components into a composite background p.d.f.
RooRealVar bkg1frac("bkg1frac","fraction of component 1 in background",0.8,0.,1.) ;
RooAddPdf bkg("bkg", "Total background", {bkg1,bkg2}, bkg1frac);
// Sum the composite signal and background
RooRealVar bkgfrac("bkgfrac","fraction of background",0.5,0.,1.) ;
RooAddPdf model("model", "g1+g2+a", {bkg,sig}, bkgfrac);
// S e t u p b a s i c p l o t w i t h d a t a a n d f u l l p d f
// ------------------------------------------------------------------------------
// Generate a data sample of 1000 events in x from model
std::unique_ptr<RooDataSet> data{model.generate(x,1000)};
// Plot data and complete PDF overlaid
RooPlot* xframe = x.frame(Title("Component plotting of pdf=(sig1+sig2)+(bkg1+bkg2)")) ;
data->plotOn(xframe) ;
model.plotOn(xframe) ;
// Clone xframe for use below
RooPlot* xframe2 = x.frame(Title("Component plotting of pdf=(sig1+sig2)+(bkg1+bkg2)"),Name("xframe2")) ;
// M a k e c o m p o n e n t b y o b j e c t r e f e r e n c e
// --------------------------------------------------------------------
// Plot single background component specified by object reference
model.plotOn(xframe,Components(bkg),LineColor(kRed)) ;
// Plot single background component specified by object reference
model.plotOn(xframe,Components(bkg2),LineStyle(kDashed),LineColor(kRed)) ;
// Plot multiple background components specified by object reference
// Note that specified components may occur at any level in object tree
// (e.g bkg is component of 'model' and 'sig2' is component 'sig')
model.plotOn(xframe,Components(RooArgSet(bkg,sig2)),LineStyle(kDotted)) ;
// M a k e c o m p o n e n t b y n a m e / r e g e x p
// ------------------------------------------------------------
// Plot single background component specified by name
model.plotOn(xframe2,Components("bkg"),LineColor(kCyan)) ;
// Plot multiple background components specified by name
model.plotOn(xframe2,Components("bkg1,sig2"),LineStyle(kDotted),LineColor(kCyan)) ;
// Plot multiple background components specified by regular expression on name
model.plotOn(xframe2,Components("sig*"),LineStyle(kDashed),LineColor(kCyan)) ;
// Plot multiple background components specified by multiple regular expressions on name
model.plotOn(xframe2,Components("bkg1,sig*"),LineStyle(kDashed),LineColor(kYellow),Invisible()) ;
regPlot(xframe,"rf205_plot1") ;
regPlot(xframe2,"rf205_plot2") ;
return true ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'ADDITION AND CONVOLUTION' RooFit tutorial macro #208
//
// One-dimensional numeric convolution
// (require ROOT to be compiled with --enable-fftw3)
//
// pdf = landau(t) (x) gauss(t)
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic208 : public RooUnitTest
{
public:
TestBasic208(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("FFT Convolution operator p.d.f.",refFile,writeRef,verbose) {} ;
bool isTestAvailable() {
// only if ROOT was build with fftw3 enabled
TString conffeatures = gROOT->GetConfigFeatures();
if(conffeatures.Contains("fftw3")) {
TPluginHandler *h;
if ((h = gROOT->GetPluginManager()->FindHandler("TVirtualFFT"))) {
if (h->LoadPlugin() == -1) {
gROOT->ProcessLine("new TNamed ;") ;
return false;
} else {
return true ;
}
}
}
return false ;
}
double ctol() { return 1e-2 ; } // Account for difficult shape of Landau distribution
bool testCode() {
// S e t u p c o m p o n e n t p d f s
// ---------------------------------------
// Construct observable
RooRealVar t("t","t",-10,30) ;
// Construct landau(t,ml,sl) ;
RooRealVar ml("ml","mean landau",5.,-20,20) ;
RooRealVar sl("sl","sigma landau",1,0.1,10) ;
RooLandau landau("lx","lx",t,ml,sl) ;
// Construct gauss(t,mg,sg)
RooRealVar mg("mg","mg",0) ;
RooRealVar sg("sg","sg",2,0.1,10) ;
RooGaussian gauss("gauss","gauss",t,mg,sg) ;
// C o n s t r u c t c o n v o l u t i o n p d f
// ---------------------------------------
// Set #bins to be used for FFT sampling to 10000
t.setBins(10000,"cache") ;
// Construct landau (x) gauss
RooFFTConvPdf lxg("lxg","landau (X) gauss",t,landau,gauss) ;
// S a m p l e , f i t a n d p l o t c o n v o l u t e d p d f
// ----------------------------------------------------------------------
// Sample 1000 events in x from gxlx
std::unique_ptr<RooDataSet> data{lxg.generate(t,10000)};
// Fit gxlx to data
{
// Get rid of the caching info prints
RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING, 0u, RooFit::Caching, true};
lxg.fitTo(*data) ;
}
// Plot data, landau pdf, landau (X) gauss pdf
RooPlot* frame = t.frame(Title("landau (x) gauss convolution")) ;
data->plotOn(frame) ;
lxg.plotOn(frame) ;
landau.plotOn(frame,LineStyle(kDashed)) ;
regPlot(frame,"rf208_plot1") ;
return true ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'ADDITION AND CONVOLUTION' RooFit tutorial macro #209
//
// Decay function p.d.fs with optional B physics
// effects (mixing and CP violation) that can be
// analytically convolved with e.g. Gaussian resolution
// functions
//
// pdf1 = decay(t,tau) (x) delta(t)
// pdf2 = decay(t,tau) (x) gauss(t,m,s)
// pdf3 = decay(t,tau) (x) (f*gauss1(t,m1,s1) + (1-f)*gauss2(t,m1,s1))
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic209 : public RooUnitTest
{
public:
TestBasic209(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Analytical convolution operator",refFile,writeRef,verbose) {} ;
bool testCode() {
// B - p h y s i c s p d f w i t h t r u t h r e s o l u t i o n
// ---------------------------------------------------------------------
// Variables of decay p.d.f.
RooRealVar dt("dt","dt",-10,10) ;
RooRealVar tau("tau","tau",1.548) ;
// Build a truth resolution model (delta function)
RooTruthModel tm("tm","truth model",dt) ;
// Construct decay(t) (x) delta(t)
RooDecay decay_tm("decay_tm","decay",dt,tau,tm,RooDecay::DoubleSided) ;
// Plot p.d.f. (dashed)
RooPlot* frame = dt.frame(Title("Bdecay (x) resolution")) ;
decay_tm.plotOn(frame,LineStyle(kDashed)) ;
// B - p h y s i c s p d f w i t h G a u s s i a n r e s o l u t i o n
// ----------------------------------------------------------------------------
// Build a gaussian resolution model
RooRealVar bias1("bias1","bias1",0) ;
RooRealVar sigma1("sigma1","sigma1",1) ;
RooGaussModel gm1("gm1","gauss model 1",dt,bias1,sigma1) ;
// Construct decay(t) (x) gauss1(t)
RooDecay decay_gm1("decay_gm1","decay",dt,tau,gm1,RooDecay::DoubleSided) ;
// Plot p.d.f.
decay_gm1.plotOn(frame) ;
// B - p h y s i c s p d f w i t h d o u b l e G a u s s i a n r e s o l u t i o n
// ------------------------------------------------------------------------------------------
// Build another gaussian resolution model
RooRealVar bias2("bias2","bias2",0) ;
RooRealVar sigma2("sigma2","sigma2",5) ;
RooGaussModel gm2("gm2","gauss model 2",dt,bias2,sigma2) ;
// Build a composite resolution model f*gm1+(1-f)*gm2
RooRealVar gm1frac("gm1frac","fraction of gm1",0.5) ;
RooAddModel gmsum("gmsum","sum of gm1 and gm2",RooArgList(gm1,gm2),gm1frac) ;
// Construct decay(t) (x) (f*gm1 + (1-f)*gm2)
RooDecay decay_gmsum("decay_gmsum","decay",dt,tau,gmsum,RooDecay::DoubleSided) ;
// Plot p.d.f. (red)
decay_gmsum.plotOn(frame,LineColor(kRed)) ;
regPlot(frame,"rf209_plot1") ;
return true ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #301
//
// Multi-dimensional p.d.f.s through composition, e.g. substituting a
// p.d.f parameter with a function that depends on other observables
//
// pdf = gauss(x,f(y),s) with f(y) = a0 + a1*y
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic301 : public RooUnitTest
{
public:
TestBasic301(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Composition extension of basic p.d.f",refFile,writeRef,verbose) {} ;
bool testCode() {
// S e t u p c o m p o s e d m o d e l g a u s s ( x , m ( y ) , s )
// -----------------------------------------------------------------------
// Create observables
RooRealVar x("x","x",-5,5) ;
RooRealVar y("y","y",-5,5) ;
// Create function f(y) = a0 + a1*y
RooRealVar a0("a0","a0",-0.5,-5,5) ;
RooRealVar a1("a1","a1",-0.5,-1,1) ;
RooPolyVar fy("fy","fy",y,RooArgSet(a0,a1)) ;
// Creat gauss(x,f(y),s)
RooRealVar sigma("sigma","width of gaussian",0.5) ;
RooGaussian model("model","Gaussian with shifting mean",x,fy,sigma) ;
// S a m p l e d a t a , p l o t d a t a a n d p d f o n x a n d y
// ---------------------------------------------------------------------------------
// Generate 10000 events in x and y from model
std::unique_ptr<RooDataSet> data{model.generate(RooArgSet(x,y),10000)};
// Plot x distribution of data and projection of model on x = Int(dy) model(x,y)
RooPlot* xframe = x.frame() ;
data->plotOn(xframe) ;
model.plotOn(xframe) ;
// Plot x distribution of data and projection of model on y = Int(dx) model(x,y)
RooPlot* yframe = y.frame() ;
data->plotOn(yframe) ;
model.plotOn(yframe) ;
// Make two-dimensional plot in x vs y
TH1* hh_model = model.createHistogram("hh_model",x,Binning(50),YVar(y,Binning(50))) ;
hh_model->SetLineColor(kBlue) ;
regPlot(xframe,"rf301_plot1") ;
regPlot(yframe,"rf302_plot2") ;
regTH(hh_model,"rf302_model2d") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #302
//
// Utility functions classes available for use in tailoring
// of composite (multidimensional) pdfs
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic302 : public RooUnitTest
{
public:
TestBasic302(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Sum and product utility functions",refFile,writeRef,verbose) {} ;
bool testCode() {
// C r e a t e o b s e r v a b l e s , p a r a m e t e r s
// -----------------------------------------------------------
// Create observables
RooRealVar x("x","x",-5,5) ;
RooRealVar y("y","y",-5,5) ;
// Create parameters
RooRealVar a0("a0","a0",-1.5,-5,5) ;
RooRealVar a1("a1","a1",-0.5,-1,1) ;
RooRealVar sigma("sigma","width of gaussian",0.5) ;
// U s i n g R o o F o r m u l a V a r t o t a i l o r p d f
// -----------------------------------------------------------------------
// Create interpreted function f(y) = a0 - a1*sqrt(10*abs(y))
RooFormulaVar fy_1("fy_1","a0-a1*sqrt(10*abs(y))",RooArgSet(y,a0,a1)) ;
// Create gauss(x,f(y),s)
RooGaussian model_1("model_1","Gaussian with shifting mean",x,fy_1,sigma) ;
// U s i n g R o o P o l y V a r t o t a i l o r p d f
// -----------------------------------------------------------------------
// Create polynomial function f(y) = a0 + a1*y
RooPolyVar fy_2("fy_2","fy_2",y,RooArgSet(a0,a1)) ;
// Create gauss(x,f(y),s)
RooGaussian model_2("model_2","Gaussian with shifting mean",x,fy_2,sigma) ;
// U s i n g R o o A d d i t i o n t o t a i l o r p d f
// -----------------------------------------------------------------------
// Create sum function f(y) = a0 + y
RooAddition fy_3("fy_3","a0+y",RooArgSet(a0,y)) ;
// Create gauss(x,f(y),s)
RooGaussian model_3("model_3","Gaussian with shifting mean",x,fy_3,sigma) ;
// U s i n g R o o P r o d u c t t o t a i l o r p d f
// -----------------------------------------------------------------------
// Create product function f(y) = a1*y
RooProduct fy_4("fy_4","a1*y",RooArgSet(a1,y)) ;
// Create gauss(x,f(y),s)
RooGaussian model_4("model_4","Gaussian with shifting mean",x,fy_4,sigma) ;
// P l o t a l l p d f s
// ----------------------------
// Make two-dimensional plots in x vs y
TH1* hh_model_1 = model_1.createHistogram("hh_model_1",x,Binning(50),YVar(y,Binning(50))) ;
TH1* hh_model_2 = model_2.createHistogram("hh_model_2",x,Binning(50),YVar(y,Binning(50))) ;
TH1* hh_model_3 = model_3.createHistogram("hh_model_3",x,Binning(50),YVar(y,Binning(50))) ;
TH1* hh_model_4 = model_4.createHistogram("hh_model_4",x,Binning(50),YVar(y,Binning(50))) ;
hh_model_1->SetLineColor(kBlue) ;
hh_model_2->SetLineColor(kBlue) ;
hh_model_3->SetLineColor(kBlue) ;
hh_model_4->SetLineColor(kBlue) ;
regTH(hh_model_1,"rf202_model2d_1") ;
regTH(hh_model_2,"rf202_model2d_2") ;
regTH(hh_model_3,"rf202_model2d_3") ;
regTH(hh_model_4,"rf202_model2d_4") ;
return true ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #303
//
// Use of tailored p.d.f as conditional p.d.fs.s
//
// pdf = gauss(x,f(y),sx | y ) with f(y) = a0 + a1*y
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic303 : public RooUnitTest
{
public:
RooDataSet* makeFakeDataXY()
{
RooRealVar x("x","x",-10,10) ;
RooRealVar y("y","y",-10,10) ;
RooArgSet coord(x,y) ;
RooDataSet* d = new RooDataSet("d","d",RooArgSet(x,y)) ;
for (int i=0 ; i<10000 ; i++) {
double tmpy = gRandom->Gaus(0,10) ;
double tmpx = gRandom->Gaus(0.5*tmpy,1) ;
if (fabs(tmpy)<10 && fabs(tmpx)<10) {
x = tmpx ;
y = tmpy ;
d->add(coord) ;
}
}
return d ;
}
TestBasic303(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Conditional use of F(x|y)",refFile,writeRef,verbose) {} ;
bool testCode() {
// S e t u p c o m p o s e d m o d e l g a u s s ( x , m ( y ) , s )
// -----------------------------------------------------------------------
// Create observables
RooRealVar x("x","x",-10,10) ;
RooRealVar y("y","y",-10,10) ;
// Create function f(y) = a0 + a1*y
RooRealVar a0("a0","a0",-0.5,-5,5) ;
RooRealVar a1("a1","a1",-0.5,-1,1) ;
RooPolyVar fy("fy","fy",y,RooArgSet(a0,a1)) ;
// Creat gauss(x,f(y),s)
RooRealVar sigma("sigma","width of gaussian",0.5,0.1,2.0) ;
RooGaussian model("model","Gaussian with shifting mean",x,fy,sigma) ;
// Obtain fake external experimental dataset with values for x and y
std::unique_ptr<RooDataSet> expDataXY{makeFakeDataXY()};
// G e n e r a t e d a t a f r o m c o n d i t i o n a l p . d . f m o d e l ( x | y )
// ---------------------------------------------------------------------------------------------
// Make subset of experimental data with only y values
std::unique_ptr<RooAbsData> expDataY{expDataXY->reduce(y)};
// Generate 10000 events in x obtained from _conditional_ model(x|y) with y values taken from experimental data
std::unique_ptr<RooDataSet> data{model.generate(x,ProtoData(static_cast<RooDataSet&>(*expDataY)))};
// F i t c o n d i t i o n a l p . d . f m o d e l ( x | y ) t o d a t a
// ---------------------------------------------------------------------------------------------
model.fitTo(*expDataXY,ConditionalObservables(y)) ;
// P r o j e c t c o n d i t i o n a l p . d . f o n x a n d y d i m e n s i o n s
// ---------------------------------------------------------------------------------------------
// Plot x distribution of data and projection of model on x = 1/Ndata sum(data(y_i)) model(x;y_i)
RooPlot* xframe = x.frame() ;
expDataXY->plotOn(xframe) ;
model.plotOn(xframe,ProjWData(*expDataY)) ;
// Speed up (and approximate) projection by using binned clone of data for projection
std::unique_ptr<RooAbsData> binnedDataY{static_cast<RooDataSet&>(*expDataY).binnedClone()};
model.plotOn(xframe,ProjWData(*binnedDataY),LineColor(kCyan),LineStyle(kDotted),Name("Alt1")) ;
// Show effect of projection with too coarse binning
static_cast<RooRealVar*>(expDataY->get()->find("y"))->setBins(5) ;
std::unique_ptr<RooAbsData> binnedDataY2{static_cast<RooDataSet&>(*expDataY).binnedClone()};
model.plotOn(xframe,ProjWData(*binnedDataY2),LineColor(kRed),Name("Alt2")) ;
regPlot(xframe,"rf303_plot1") ;
return true ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #304
//
// Simple uncorrelated multi-dimensional p.d.f.s
//
// pdf = gauss(x,mx,sx) * gauss(y,my,sy)
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic304 : public RooUnitTest
{
public:
TestBasic304(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Product operator p.d.f. with uncorrelated terms",refFile,writeRef,verbose) {} ;
bool testCode() {
// C r e a t e c o m p o n e n t p d f s i n x a n d y
// ----------------------------------------------------------------
// Create two p.d.f.s gaussx(x,meanx,sigmax) gaussy(y,meany,sigmay) and its variables
RooRealVar x("x","x",-5,5) ;
RooRealVar y("y","y",-5,5) ;
RooRealVar meanx("mean1","mean of gaussian x",2) ;
RooRealVar meany("mean2","mean of gaussian y",-2) ;
RooRealVar sigmax("sigmax","width of gaussian x",1) ;
RooRealVar sigmay("sigmay","width of gaussian y",5) ;
RooGaussian gaussx("gaussx","gaussian PDF",x,meanx,sigmax) ;
RooGaussian gaussy("gaussy","gaussian PDF",y,meany,sigmay) ;
// C o n s t r u c t u n c o r r e l a t e d p r o d u c t p d f
// -------------------------------------------------------------------
// Multiply gaussx and gaussy into a two-dimensional p.d.f. gaussxy
RooProdPdf gaussxy("gaussxy","gaussx*gaussy",RooArgList(gaussx,gaussy)) ;
// S a m p l e p d f , p l o t p r o j e c t i o n o n x a n d y
// ---------------------------------------------------------------------------
// Generate 10000 events in x and y from gaussxy
std::unique_ptr<RooDataSet> data{gaussxy.generate(RooArgSet(x,y),10000)};
// Plot x distribution of data and projection of gaussxy on x = Int(dy) gaussxy(x,y)
RooPlot* xframe = x.frame(Title("X projection of gauss(x)*gauss(y)")) ;
data->plotOn(xframe) ;
gaussxy.plotOn(xframe) ;
// Plot x distribution of data and projection of gaussxy on y = Int(dx) gaussxy(x,y)
RooPlot* yframe = y.frame(Title("Y projection of gauss(x)*gauss(y)")) ;
data->plotOn(yframe) ;
gaussxy.plotOn(yframe) ;
regPlot(xframe,"rf304_plot1") ;
regPlot(yframe,"rf304_plot2") ;
return true ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #305
//
// Multi-dimensional p.d.f.s with conditional p.d.fs in product
//
// pdf = gauss(x,f(y),sx | y ) * gauss(y,ms,sx) with f(y) = a0 + a1*y
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic305 : public RooUnitTest
{
public:
TestBasic305(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Product operator p.d.f. with conditional term",refFile,writeRef,verbose) {} ;
bool testCode() {
// C r e a t e c o n d i t i o n a l p d f g x ( x | y )
// -----------------------------------------------------------
// Create observables
RooRealVar x("x","x",-5,5) ;
RooRealVar y("y","y",-5,5) ;
// Create function f(y) = a0 + a1*y
RooRealVar a0("a0","a0",-0.5,-5,5) ;
RooRealVar a1("a1","a1",-0.5,-1,1) ;
RooPolyVar fy("fy","fy",y,RooArgSet(a0,a1)) ;
// Create gaussx(x,f(y),sx)
RooRealVar sigmax("sigma","width of gaussian",0.5) ;
RooGaussian gaussx("gaussx","Gaussian in x with shifting mean in y",x,fy,sigmax) ;
// C r e a t e p d f g y ( y )
// -----------------------------------------------------------
// Create gaussy(y,0,5)
RooGaussian gaussy("gaussy","Gaussian in y",y,0.0,3.0) ;
// C r e a t e p r o d u c t g x ( x | y ) * g y ( y )
// -------------------------------------------------------
// Create gaussx(x,sx|y) * gaussy(y)
RooProdPdf model("model","gaussx(x|y)*gaussy(y)",gaussy,Conditional(gaussx,x)) ;
// S a m p l e , f i t a n d p l o t p r o d u c t p d f
// ---------------------------------------------------------------
// Generate 1000 events in x and y from model
std::unique_ptr<RooDataSet> data{model.generate(RooArgSet(x,y),10000)};
// Plot x distribution of data and projection of model on x = Int(dy) model(x,y)
RooPlot* xframe = x.frame() ;
data->plotOn(xframe) ;
model.plotOn(xframe) ;
// Plot x distribution of data and projection of model on y = Int(dx) model(x,y)
RooPlot* yframe = y.frame() ;
data->plotOn(yframe) ;
model.plotOn(yframe) ;
// Make two-dimensional plot in x vs y
TH1* hh_model = model.createHistogram("hh_model_rf305",x,Binning(50),YVar(y,Binning(50))) ;
hh_model->SetLineColor(kBlue) ;
regTH(hh_model,"rf305_model2d") ;
regPlot(xframe,"rf305_plot1") ;
regPlot(yframe,"rf305_plot2") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #306
//
// Complete example with use of conditional p.d.f. with per-event errors
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic306 : public RooUnitTest
{
public:
TestBasic306(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Conditional use of per-event error p.d.f. F(t|dt)",refFile,writeRef,verbose) {} ;
bool testCode() {
// B - p h y s i c s p d f w i t h p e r - e v e n t G a u s s i a n r e s o l u t i o n
// ----------------------------------------------------------------------------------------------
// Observables
RooRealVar dt("dt","dt",-10,10) ;
RooRealVar dterr("dterr","per-event error on dt",0.01,10) ;
// Build a gaussian resolution model scaled by the per-event error = gauss(dt,bias,sigma*dterr)
RooRealVar bias("bias","bias",0,-10,10) ;
RooRealVar sigma("sigma","per-event error scale factor",1,0.1,10) ;
RooGaussModel gm("gm1","gauss model scaled bt per-event error",dt,bias,sigma,dterr) ;
// Construct decay(dt) (x) gauss1(dt|dterr)
RooRealVar tau("tau","tau",1.548) ;
RooDecay decay_gm("decay_gm","decay",dt,tau,gm,RooDecay::DoubleSided) ;
// C o n s t r u c t f a k e ' e x t e r n a l ' d a t a w i t h p e r - e v e n t e r r o r
// ------------------------------------------------------------------------------------------------------
// Use landau p.d.f to get somewhat realistic distribution with long tail
RooLandau pdfDtErr("pdfDtErr","pdfDtErr",dterr,1.0,0.25) ;
std::unique_ptr<RooDataSet> expDataDterr{pdfDtErr.generate(dterr,10000)};
// S a m p l e d a t a f r o m c o n d i t i o n a l d e c a y _ g m ( d t | d t e r r )
// ---------------------------------------------------------------------------------------------
// Specify external dataset with dterr values to use decay_dm as conditional p.d.f.
std::unique_ptr<RooDataSet> data{decay_gm.generate(dt,ProtoData(*expDataDterr))};
// F i t c o n d i t i o n a l d e c a y _ d m ( d t | d t e r r )
// ---------------------------------------------------------------------
// Specify dterr as conditional observable
decay_gm.fitTo(*data,ConditionalObservables(dterr)) ;
// P l o t c o n d i t i o n a l d e c a y _ d m ( d t | d t e r r )
// ---------------------------------------------------------------------
// Make two-dimensional plot of conditional p.d.f in (dt,dterr)
TH1* hh_decay = decay_gm.createHistogram("hh_decay",dt,Binning(50),YVar(dterr,Binning(50))) ;
hh_decay->SetLineColor(kBlue) ;
// Plot decay_gm(dt|dterr) at various values of dterr
RooPlot* frame = dt.frame(Title("Slices of decay(dt|dterr) at various dterr")) ;
for (int ibin=0 ; ibin<100 ; ibin+=20) {
dterr.setBin(ibin) ;
decay_gm.plotOn(frame,Normalization(5.),Name(Form("curve_slice_%d",ibin))) ;
}
// Make projection of data an dt
RooPlot* frame2 = dt.frame(Title("Projection of decay(dt|dterr) on dt")) ;
data->plotOn(frame2) ;
// Make projection of decay(dt|dterr) on dt.
//
// Instead of integrating out dterr, make a weighted average of curves
// at values dterr_i as given in the external dataset.
// (The true argument bins the data before projection to speed up the process)
decay_gm.plotOn(frame2,ProjWData(*expDataDterr,true)) ;
regTH(hh_decay,"rf306_model2d") ;
regPlot(frame,"rf306_plot1") ;
regPlot(frame2,"rf306_plot2") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #307
//
// Complete example with use of full p.d.f. with per-event errors
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic307 : public RooUnitTest
{
public:
TestBasic307(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Full per-event error p.d.f. F(t|dt)G(dt)",refFile,writeRef,verbose) {} ;
bool testCode() {
// B - p h y s i c s p d f w i t h p e r - e v e n t G a u s s i a n r e s o l u t i o n
// ----------------------------------------------------------------------------------------------
// Observables
RooRealVar dt("dt","dt",-10,10) ;
RooRealVar dterr("dterr","per-event error on dt",0.1,10) ;
// Build a gaussian resolution model scaled by the per-event error = gauss(dt,bias,sigma*dterr)
RooRealVar bias("bias","bias",0,-10,10) ;
RooRealVar sigma("sigma","per-event error scale factor",1,0.1,10) ;
RooGaussModel gm("gm1","gauss model scaled bt per-event error",dt,bias,sigma,dterr) ;
// Construct decay(dt) (x) gauss1(dt|dterr)
RooRealVar tau("tau","tau",1.548) ;
RooDecay decay_gm("decay_gm","decay",dt,tau,gm,RooDecay::DoubleSided) ;
// C o n s t r u c t e m p i r i c a l p d f f o r p e r - e v e n t e r r o r
// -----------------------------------------------------------------
// Use landau p.d.f to get empirical distribution with long tail
RooLandau pdfDtErr("pdfDtErr","pdfDtErr",dterr,1.0,0.25) ;
std::unique_ptr<RooDataSet> expDataDterr{pdfDtErr.generate(dterr,10000)};
// Construct a histogram pdf to describe the shape of the dtErr distribution
std::unique_ptr<RooDataHist> expHistDterr{expDataDterr->binnedClone()};
RooHistPdf pdfErr("pdfErr","pdfErr",dterr,*expHistDterr) ;
// C o n s t r u c t c o n d i t i o n a l p r o d u c t d e c a y _ d m ( d t | d t e r r ) * p d f ( d t e r r )
// ----------------------------------------------------------------------------------------------------------------------
// Construct production of conditional decay_dm(dt|dterr) with empirical pdfErr(dterr)
RooProdPdf model("model","model",pdfErr,Conditional(decay_gm,dt)) ;
// (Alternatively you could also use the landau shape pdfDtErr)
//RooProdPdf model("model","model",pdfDtErr,Conditional(decay_gm,dt)) ;
// S a m p l e, f i t a n d p l o t p r o d u c t m o d e l
// ------------------------------------------------------------------
// Specify external dataset with dterr values to use model_dm as conditional p.d.f.
std::unique_ptr<RooDataSet> data{model.generate(RooArgSet(dt,dterr),10000)};
// F i t c o n d i t i o n a l d e c a y _ d m ( d t | d t e r r )
// ---------------------------------------------------------------------
// Specify dterr as conditional observable
model.fitTo(*data) ;
// P l o t c o n d i t i o n a l d e c a y _ d m ( d t | d t e r r )
// ---------------------------------------------------------------------
// Make projection of data an dt
RooPlot* frame = dt.frame(Title("Projection of model(dt|dterr) on dt")) ;
data->plotOn(frame) ;
model.plotOn(frame) ;
regPlot(frame,"rf307_plot1") ;
return true ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #308
//
// Examples on normalization of p.d.f.s,
// integration of p.d.fs, construction
// of cumulative distribution functions from p.d.f.s
// in two dimensions
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic308 : public RooUnitTest
{
public:
TestBasic308(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Normalization of p.d.f.s in 2D",refFile,writeRef,verbose) {} ;
bool testCode() {
// S e t u p m o d e l
// ---------------------
// Create observables x,y
RooRealVar x("x","x",-10,10) ;
RooRealVar y("y","y",-10,10) ;
// Create p.d.f. gaussx(x,-2,3), gaussy(y,2,2)
RooGaussian gx("gx","gx",x,-2.0,3.0) ;
RooGaussian gy("gy","gy",y,+2.0,2.0) ;
// Create gxy = gx(x)*gy(y)
RooProdPdf gxy("gxy","gxy",RooArgSet(gx,gy)) ;
// R e t r i e v e r a w & n o r m a l i z e d v a l u e s o f R o o F i t p . d . f . s
// --------------------------------------------------------------------------------------------------
// Return 'raw' unnormalized value of gx
regValue(gxy.getVal(),"rf308_gxy") ;
// Return value of gxy normalized over x _and_ y in range [-10,10]
RooArgSet nset_xy(x,y) ;
regValue(gxy.getVal(&nset_xy),"rf308_gx_Norm[x,y]") ;
// Create object representing integral over gx
// which is used to calculate gx_Norm[x,y] == gx / gx_Int[x,y]
std::unique_ptr<RooAbsReal> igxy{gxy.createIntegral(RooArgSet(x,y))};
regValue(igxy->getVal(),"rf308_gx_Int[x,y]") ;
// NB: it is also possible to do the following
// Return value of gxy normalized over x in range [-10,10] (i.e. treating y as parameter)
RooArgSet nset_x(x) ;
regValue(gxy.getVal(&nset_x),"rf308_gx_Norm[x]") ;
// Return value of gxy normalized over y in range [-10,10] (i.e. treating x as parameter)
RooArgSet nset_y(y) ;
regValue(gxy.getVal(&nset_y),"rf308_gx_Norm[y]") ;
// I n t e g r a t e n o r m a l i z e d p d f o v e r s u b r a n g e
// ----------------------------------------------------------------------------
// Define a range named "signal" in x from -5,5
x.setRange("signal",-5,5) ;
y.setRange("signal",-3,3) ;
// Create an integral of gxy_Norm[x,y] over x and y in range "signal"
// This is the fraction of of p.d.f. gxy_Norm[x,y] which is in the
// range named "signal"
std::unique_ptr<RooAbsReal> igxy_sig{gxy.createIntegral(RooArgSet(x,y),NormSet(RooArgSet(x,y)),Range("signal"))};
regValue(igxy_sig->getVal(),"rf308_gx_Int[x,y|signal]_Norm[x,y]") ;
// C o n s t r u c t c u m u l a t i v e d i s t r i b u t i o n f u n c t i o n f r o m p d f
// -----------------------------------------------------------------------------------------------------
// Create the cumulative distribution function of gx
// i.e. calculate Int[-10,x] gx(x') dx'
std::unique_ptr<RooAbsReal> gxy_cdf{gxy.createCdf(RooArgSet(x,y))};
// Plot cdf of gx versus x
TH1* hh_cdf = gxy_cdf->createHistogram("hh_cdf",x,Binning(40),YVar(y,Binning(40))) ;
regTH(hh_cdf,"rf308_cdf") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #309
//
// Projecting p.d.f and data slices in discrete observables
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic310 : public RooUnitTest
{
public:
TestBasic310(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Data and p.d.f projection in category slice",refFile,writeRef,verbose) {} ;
bool testCode() {
// C r e a t e B d e c a y p d f w it h m i x i n g
// ----------------------------------------------------------
// Decay time observables
RooRealVar dt("dt","dt",-20,20) ;
// Discrete observables mixState (B0tag==B0reco?) and tagFlav (B0tag==B0(bar)?)
RooCategory mixState("mixState","B0/B0bar mixing state") ;
RooCategory tagFlav("tagFlav","Flavour of the tagged B0") ;
// Define state labels of discrete observables
mixState.defineType("mixed",-1) ;
mixState.defineType("unmixed",1) ;
tagFlav.defineType("B0",1) ;
tagFlav.defineType("B0bar",-1) ;
// Model parameters
RooRealVar dm("dm","delta m(B)",0.472,0.,1.0) ;
RooRealVar tau("tau","B0 decay time",1.547,1.0,2.0) ;
RooRealVar w("w","Flavor Mistag rate",0.03,0.0,1.0) ;
RooRealVar dw("dw","Flavor Mistag rate difference between B0 and B0bar",0.01) ;
// Build a gaussian resolution model
RooRealVar bias1("bias1","bias1",0) ;
RooRealVar sigma1("sigma1","sigma1",0.01) ;
RooGaussModel gm1("gm1","gauss model 1",dt,bias1,sigma1) ;
// Construct a decay pdf, smeared with single gaussian resolution model
RooBMixDecay bmix_gm1("bmix","decay",dt,mixState,tagFlav,tau,dm,w,dw,gm1,RooBMixDecay::DoubleSided) ;
// Generate BMixing data with above set of event errors
std::unique_ptr<RooDataSet> data{bmix_gm1.generate(RooArgSet(dt,tagFlav,mixState),20000)};
// P l o t f u l l d e c a y d i s t r i b u t i o n
// ----------------------------------------------------------
// Create frame, plot data and pdf projection (integrated over tagFlav and mixState)
RooPlot* frame = dt.frame(Title("Inclusive decay distribution")) ;
data->plotOn(frame) ;
bmix_gm1.plotOn(frame) ;
// P l o t d e c a y d i s t r . f o r m i x e d a n d u n m i x e d s l i c e o f m i x S t a t e
// ------------------------------------------------------------------------------------------------------------------
// Create frame, plot data (mixed only)
RooPlot* frame2 = dt.frame(Title("Decay distribution of mixed events")) ;
data->plotOn(frame2,Cut("mixState==mixState::mixed")) ;
// Position slice in mixState at "mixed" and plot slice of pdf in mixstate over data (integrated over tagFlav)
bmix_gm1.plotOn(frame2,Slice(mixState,"mixed")) ;
// Create frame, plot data (unmixed only)
RooPlot* frame3 = dt.frame(Title("Decay distribution of unmixed events")) ;
data->plotOn(frame3,Cut("mixState==mixState::unmixed")) ;
// Position slice in mixState at "unmixed" and plot slice of pdf in mixstate over data (integrated over tagFlav)
bmix_gm1.plotOn(frame3,Slice(mixState,"unmixed")) ;
regPlot(frame,"rf310_plot1") ;
regPlot(frame2,"rf310_plot2") ;
regPlot(frame3,"rf310_plot3") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #310
//
// Projecting p.d.f and data ranges in continuous observables
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic311 : public RooUnitTest
{
public:
TestBasic311(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Data and p.d.f projection in sub range",refFile,writeRef,verbose) {} ;
bool testCode() {
// C r e a t e 3 D p d f a n d d a t a
// -------------------------------------------
// Create observables
RooRealVar x("x","x",-5,5) ;
RooRealVar y("y","y",-5,5) ;
RooRealVar z("z","z",-5,5) ;
// Create signal pdf gauss(x)*gauss(y)*gauss(z)
RooGaussian gx("gx","gx",x,0.0,1.0) ;
RooGaussian gy("gy","gy",y,0.0,1.0) ;
RooGaussian gz("gz","gz",z,0.0,1.0) ;
RooProdPdf sig("sig","sig",RooArgSet(gx,gy,gz)) ;
// Create background pdf poly(x)*poly(y)*poly(z)
RooPolynomial px("px","px",x,RooArgSet(-0.1,0.004)) ;
RooPolynomial py("py","py",y,RooArgSet(0.1,-0.004)) ;
RooPolynomial pz("pz","pz",z) ;
RooProdPdf bkg("bkg","bkg",RooArgSet(px,py,pz)) ;
// Create composite pdf sig+bkg
RooRealVar fsig("fsig","signal fraction",0.1,0.,1.) ;
RooAddPdf model("model","model",RooArgList(sig,bkg),fsig) ;
std::unique_ptr<RooDataSet> data{model.generate(RooArgSet(x,y,z),20000)};
// P r o j e c t p d f a n d d a t a o n x
// -------------------------------------------------
// Make plain projection of data and pdf on x observable
RooPlot* frame = x.frame(Title("Projection of 3D data and pdf on X"),Bins(40)) ;
data->plotOn(frame) ;
model.plotOn(frame) ;
// P r o j e c t p d f a n d d a t a o n x i n s i g n a l r a n g e
// ----------------------------------------------------------------------------------
// Define signal region in y and z observables
y.setRange("sigRegion",-1,1) ;
z.setRange("sigRegion",-1,1) ;
// Make plot frame
RooPlot* frame2 = x.frame(Title("Same projection on X in signal range of (Y,Z)"),Bins(40)) ;
// Plot subset of data in which all observables are inside "sigRegion"
// For observables that do not have an explicit "sigRegion" range defined (e.g. observable)
// an implicit definition is used that is identical to the full range (i.e. [-5,5] for x)
data->plotOn(frame2,CutRange("sigRegion")) ;
// Project model on x, integrating projected observables (y,z) only in "sigRegion"
model.plotOn(frame2,ProjectionRange("sigRegion")) ;
regPlot(frame,"rf311_plot1") ;
regPlot(frame2,"rf312_plot2") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #312
//
// Performing fits in multiple (disjoint) ranges in one or more dimensions
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic312 : public RooUnitTest
{
public:
TestBasic312(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Fit in multiple rectangular ranges",refFile,writeRef,verbose) {} ;
bool testCode() {
// C r e a t e 2 D p d f a n d d a t a
// -------------------------------------------
// Define observables x,y
RooRealVar x("x","x",-10,10) ;
RooRealVar y("y","y",-10,10) ;
// Construct the signal pdf gauss(x)*gauss(y)
RooRealVar mx("mx","mx",1,-10,10) ;
RooRealVar my("my","my",1,-10,10) ;
RooGaussian gx("gx","gx",x,mx,1.0) ;
RooGaussian gy("gy","gy",y,my,1.0) ;
RooProdPdf sig("sig","sig",gx,gy) ;
// Construct the background pdf (flat in x,y)
RooPolynomial px("px","px",x) ;
RooPolynomial py("py","py",y) ;
RooProdPdf bkg("bkg","bkg",px,py) ;
// Construct the composite model sig+bkg
RooRealVar f("f","f",0.,1.) ;
RooAddPdf model("model","model",RooArgList(sig,bkg),f) ;
// Sample 10000 events in (x,y) from the model
std::unique_ptr<RooDataSet> modelData{model.generate(RooArgSet(x,y),10000)};
// D e f i n e s i g n a l a n d s i d e b a n d r e g i o n s
// -------------------------------------------------------------------
// Construct the SideBand1,SideBand2,Signal regions
//
// |
// +-------------+-----------+
// | | |
// | Side | Sig |
// | Band1 | nal |
// | | |
// --+-------------+-----------+--
// | |
// | Side |
// | Band2 |
// | |
// +-------------+-----------+
// |
x.setRange("SB1",-10,+10) ;
y.setRange("SB1",-10,0) ;
x.setRange("SB2",-10,0) ;
y.setRange("SB2",0,+10) ;
x.setRange("SIG",0,+10) ;
y.setRange("SIG",0,+10) ;
x.setRange("FULL",-10,+10) ;
y.setRange("FULL",-10,+10) ;
// P e r f o r m f i t s i n i n d i v i d u a l s i d e b a n d r e g i o n s
// -------------------------------------------------------------------------------------
// Perform fit in SideBand1 region (RooAddPdf coefficients will be interpreted in full range)
std::unique_ptr<RooFitResult> r_sb1{model.fitTo(*modelData,Range("SB1"),Save())};
// Perform fit in SideBand2 region (RooAddPdf coefficients will be interpreted in full range)
std::unique_ptr<RooFitResult> r_sb2{model.fitTo(*modelData,Range("SB2"),Save())};
// P e r f o r m f i t s i n j o i n t s i d e b a n d r e g i o n s
// -----------------------------------------------------------------------------
// Now perform fit to joint 'L-shaped' sideband region 'SB1|SB2'
// (RooAddPdf coefficients will be interpreted in full range)
std::unique_ptr<RooFitResult> r_sb12{model.fitTo(*modelData,Range("SB1,SB2"),Save())};
regResult(std::move(r_sb1),"rf312_fit_sb1") ;
regResult(std::move(r_sb2),"rf312_fit_sb2") ;
regResult(std::move(r_sb12),"rf312_fit_sb12") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #313
//
// Working with parameterized ranges to define non-rectangular regions
// for fitting and integration
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic313 : public RooUnitTest
{
public:
TestBasic313(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Integration over non-rectangular regions",refFile,writeRef,verbose) {} ;
bool testCode() {
// C r e a t e 3 D p d f
// -------------------------
// Define observable (x,y,z)
RooRealVar x("x","x",0,10) ;
RooRealVar y("y","y",0,10) ;
RooRealVar z("z","z",0,10) ;
// Define 3 dimensional pdf
RooRealVar z0("z0","z0",-0.1,1) ;
RooPolynomial px("px","px",x,RooConst(0)) ;
RooPolynomial py("py","py",y,RooConst(0)) ;
RooPolynomial pz("pz","pz",z,z0) ;
RooProdPdf pxyz("pxyz","pxyz",RooArgSet(px,py,pz)) ;
// D e f i n e d n o n - r e c t a n g u l a r r e g i o n R i n ( x , y , z )
// -------------------------------------------------------------------------------------
//
// R = Z[0 - 0.1*Y^2] * Y[0.1*X - 0.9*X] * X[0 - 10]
//
// Construct range parameterized in "R" in y [ 0.1*x, 0.9*x ]
RooFormulaVar ylo("ylo","0.1*x",x) ;
RooFormulaVar yhi("yhi","0.9*x",x) ;
y.setRange("R",ylo,yhi) ;
// Construct parameterized ranged "R" in z [ 0, 0.1*y^2 ]
RooFormulaVar zlo("zlo","0.0*y",y) ;
RooFormulaVar zhi("zhi","0.1*y*y",y) ;
z.setRange("R",zlo,zhi) ;
// C a l c u l a t e i n t e g r a l o f n o r m a l i z e d p d f i n R
// ----------------------------------------------------------------------------------
{
// To remove the INFO:NumericIntegration ouput from the stressRooFit output,
// change the message level locally.
RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::INFO, 0u, RooFit::NumIntegration, false};
// Create integral over normalized pdf model over x,y,z in "R" region
std::unique_ptr<RooAbsReal> intPdf{pxyz.createIntegral(RooArgSet(x,y,z),RooArgSet(x,y,z),"R")};
// Plot value of integral as function of pdf parameter z0
RooPlot* frame = z0.frame(Title("Integral of pxyz over x,y,z in region R")) ;
intPdf->plotOn(frame) ;
regPlot(frame,"rf313_plot1") ;
}
return true;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #314
//
// Working with parameterized ranges in a fit. This an example of a
// fit with an acceptance that changes per-event
//
// pdf = exp(-t/tau) with t[tmin,5]
//
// where t and tmin are both observables in the dataset
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic314 : public RooUnitTest
{
public:
TestBasic314(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Fit with non-rectangular observable boundaries",refFile,writeRef,verbose) {} ;
bool testCode() {
// D e f i n e o b s e r v a b l e s a n d d e c a y p d f
// ---------------------------------------------------------------
// Declare observables
RooRealVar t("t","t",0,5) ;
RooRealVar tmin("tmin","tmin",0,0,5) ;
// Make parameterized range in t : [tmin,5]
t.setRange(tmin,RooConst(t.getMax())) ;
// Make pdf
RooRealVar tau("tau","tau",-1.54,-10,-0.1) ;
RooExponential model("model","model",t,tau) ;
// C r e a t e i n p u t d a t a
// ------------------------------------
// Generate complete dataset without acceptance cuts (for reference)
std::unique_ptr<RooDataSet> dall{model.generate(t,10000)};
// Generate a (fake) prototype dataset for acceptance limit values
std::unique_ptr<RooDataSet> tmp{RooGaussian("gmin","gmin",tmin,0.0,0.5).generate(tmin,5000)};
// Generate dataset with t values that observe (t>tmin)
std::unique_ptr<RooDataSet> dacc{model.generate(t,ProtoData(*tmp))};
// F i t p d f t o d a t a i n a c c e p t a n c e r e g i o n
// -----------------------------------------------------------------------
std::unique_ptr<RooFitResult> r{model.fitTo(*dacc,Save())};
// P l o t f i t t e d p d f o n f u l l a n d a c c e p t e d d a t a
// ---------------------------------------------------------------------------------
// Make plot frame, add datasets and overlay model
RooPlot* frame = t.frame(Title("Fit to data with per-event acceptance")) ;
dall->plotOn(frame,MarkerColor(kRed),LineColor(kRed)) ;
model.plotOn(frame) ;
dacc->plotOn(frame,Name("dacc")) ;
// Print fit results to demonstrate absence of bias
regResult(std::move(r),"rf314_fit") ;
regPlot(frame,"rf314_plot1") ;
return true;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #315
//
// Marginizalization of multi-dimensional p.d.f.s through integration
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic315 : public RooUnitTest
{
public:
TestBasic315(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("P.d.f. marginalization through integration",refFile,writeRef,verbose) {} ;
bool testCode() {
// C r e a t e p d f m ( x , y ) = g x ( x | y ) * g ( y )
// --------------------------------------------------------------
// Increase default precision of numeric integration
// as this exercise has high sensitivity to numeric integration precision
RooAbsPdf::defaultIntegratorConfig()->setEpsRel(1e-8) ;
RooAbsPdf::defaultIntegratorConfig()->setEpsAbs(1e-8) ;
// Create observables
RooRealVar x("x","x",-5,5) ;
RooRealVar y("y","y",-2,2) ;
// Create function f(y) = a0 + a1*y
RooRealVar a0("a0","a0",0) ;
RooRealVar a1("a1","a1",-1.5,-3,1) ;
RooPolyVar fy("fy","fy",y,RooArgSet(a0,a1)) ;
// Create gaussx(x,f(y),sx)
RooRealVar sigmax("sigmax","width of gaussian",0.5) ;
RooGaussian gaussx("gaussx","Gaussian in x with shifting mean in y",x,fy,sigmax) ;
// Create gaussy(y,0,2)
RooGaussian gaussy("gaussy","Gaussian in y",y,0.0,2.0) ;
// Create gaussx(x,sx|y) * gaussy(y)
RooProdPdf model("model","gaussx(x|y)*gaussy(y)",gaussy,Conditional(gaussx,x)) ;
// M a r g i n a l i z e m ( x , y ) t o m ( x )
// ----------------------------------------------------
// modelx(x) = Int model(x,y) dy
std::unique_ptr<RooAbsPdf> modelx{model.createProjection(y)};
// U s e m a r g i n a l i z e d p . d . f . a s r e g u l a r 1 - D p . d . f .
// ------------------------------------------------------------------------------------------
// Sample 1000 events from modelx
std::unique_ptr<RooAbsData> data{modelx->generateBinned(x,1000)};
// Fit modelx to toy data
modelx->fitTo(*data) ;
// Plot modelx over data
RooPlot* frame = x.frame(40) ;
data->plotOn(frame) ;
modelx->plotOn(frame) ;
regPlot(frame,"rf315_frame") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #316
//
// Using the likelihood ratio techique to construct a signal enhanced
// one-dimensional projection of a multi-dimensional p.d.f.
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic316 : public RooUnitTest
{
public:
TestBasic316(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Likelihood ratio projection plot",refFile,writeRef,verbose) {} ;
bool testCode() {
// C r e a t e 3 D p d f a n d d a t a
// -------------------------------------------
// Create observables
RooRealVar x("x","x",-5,5) ;
RooRealVar y("y","y",-5,5) ;
RooRealVar z("z","z",-5,5) ;
// Create signal pdf gauss(x)*gauss(y)*gauss(z)
RooGaussian gx("gx","gx",x,0.0,1.0) ;
RooGaussian gy("gy","gy",y,0.0,1.0) ;
RooGaussian gz("gz","gz",z,0.0,1.0) ;
RooProdPdf sig("sig","sig",RooArgSet(gx,gy,gz)) ;
// Create background pdf poly(x)*poly(y)*poly(z)
RooPolynomial px("px","px",x,RooArgSet(-0.1,0.004)) ;
RooPolynomial py("py","py",y,RooArgSet(0.1,-0.004)) ;
RooPolynomial pz("pz","pz",z) ;
RooProdPdf bkg("bkg","bkg",RooArgSet(px,py,pz)) ;
// Create composite pdf sig+bkg
RooRealVar fsig("fsig","signal fraction",0.1,0.,1.) ;
RooAddPdf model("model","model",RooArgList(sig,bkg),fsig) ;
std::unique_ptr<RooDataSet> data{model.generate(RooArgSet(x,y,z),20000)};
// P r o j e c t p d f a n d d a t a o n x
// -------------------------------------------------
// Make plain projection of data and pdf on x observable
RooPlot* frame = x.frame(Title("Projection of 3D data and pdf on X"),Bins(40)) ;
data->plotOn(frame) ;
model.plotOn(frame) ;
// D e f i n e p r o j e c t e d s i g n a l l i k e l i h o o d r a t i o
// ----------------------------------------------------------------------------------
// Calculate projection of signal and total likelihood on (y,z) observables
// i.e. integrate signal and composite model over x
std::unique_ptr<RooAbsPdf> sigyz{sig.createProjection(x)};
std::unique_ptr<RooAbsPdf> totyz{model.createProjection(x)};
// Construct the log of the signal / signal+background probability
RooFormulaVar llratio_func("llratio","log10(@0)-log10(@1)",RooArgList(*sigyz,*totyz)) ;
// P l o t d a t a w i t h a L L r a t i o c u t
// -------------------------------------------------------
// Calculate the llratio value for each event in the dataset
data->addColumn(llratio_func) ;
// Extract the subset of data with large signal likelihood
std::unique_ptr<RooAbsData> dataSel{data->reduce(Cut("llratio>0.7"))};
// Make plot frame
RooPlot* frame2 = x.frame(Title("Same projection on X with LLratio(y,z)>0.7"),Bins(40)) ;
// Plot select data on frame
dataSel->plotOn(frame2) ;
// M a k e M C p r o j e c t i o n o f p d f w i t h s a m e L L r a t i o c u t
// ---------------------------------------------------------------------------------------------
// Generate large number of events for MC integration of pdf projection
std::unique_ptr<RooDataSet> mcprojData{model.generate(RooArgSet(x,y,z),10000)};
// Calculate LL ratio for each generated event and select MC events with llratio)0.7
mcprojData->addColumn(llratio_func) ;
std::unique_ptr<RooAbsData> mcprojDataSel{mcprojData->reduce(Cut("llratio>0.7"))};
// Project model on x, integrating projected observables (y,z) with Monte Carlo technique
// on set of events with the same llratio cut as was applied to data
model.plotOn(frame2,ProjWData(static_cast<RooDataSet&>(*mcprojDataSel)));
regPlot(frame,"rf316_plot1") ;
regPlot(frame2,"rf316_plot2") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'DATA AND CATEGORIES' RooFit tutorial macro #402
//
// Tools for manipulation of (un)binned datasets
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic402 : public RooUnitTest
{
public:
TestBasic402(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Basic operations on datasets",refFile,writeRef,verbose) {} ;
bool testCode() {
// Binned (RooDataHist) and unbinned datasets (RooDataSet) share
// many properties and inherit from a common abstract base class
// (RooAbsData), that provides an interface for all operations
// that can be performed regardless of the data format
RooRealVar x("x","x",-10,10) ;
RooRealVar y("y","y", 0, 40) ;
RooCategory c("c","c") ;
c.defineType("Plus",+1) ;
c.defineType("Minus",-1) ;
// B a s i c O p e r a t i o n s o n u n b i n n e d d a t a s e t s
// --------------------------------------------------------------
// RooDataSet is an unbinned dataset (a collection of points in N-dimensional space)
RooDataSet d("d","d",RooArgSet(x,y,c)) ;
// Unlike RooAbsArgs (RooAbsPdf,RooFormulaVar,....) datasets are not attached to
// the variables they are constructed from. Instead they are attached to an internal
// clone of the supplied set of arguments
// Fill d with dummy values
int i ;
for (i=0 ; i<1000 ; i++) {
x = i/50 - 10 ;
y = sqrt(1.0*i) ;
c.setLabel((i%2)?"Plus":"Minus") ;
// We must explicitly refer to x,y,c here to pass the values because
// d is not linked to them (as explained above)
d.add(RooArgSet(x,y,c)) ;
}
// R e d u c i n g , A p p e n d i n g a n d M e r g i n g
// -------------------------------------------------------------
// The reduce() function returns a new dataset which is a subset of the original
std::unique_ptr<RooAbsData> d1{d.reduce(RooArgSet(x,c))};
std::unique_ptr<RooAbsData> d2{d.reduce(RooArgSet(y))};
std::unique_ptr<RooAbsData> d3{d.reduce("y>5.17")};
std::unique_ptr<RooAbsData> d4{d.reduce(RooArgSet(x,c),"y>5.17")};
regValue(d3->numEntries(),"rf403_nd3") ;
regValue(d4->numEntries(),"rf403_nd4") ;
// The merge() function adds two data set column-wise
static_cast<RooDataSet&>(*d1).merge(static_cast<RooDataSet*>(d2.get()));
// The append() function addes two datasets row-wise
static_cast<RooDataSet&>(*d1).append(static_cast<RooDataSet&>(*d3));
regValue(d1->numEntries(),"rf403_nd1") ;
// O p e r a t i o n s o n b i n n e d d a t a s e t s
// ---------------------------------------------------------
// A binned dataset can be constructed empty, from an unbinned dataset, or
// from a ROOT native histogram (TH1,2,3)
// The binning of real variables (like x,y) is done using their fit range
// 'get/setRange()' and number of specified fit bins 'get/setBins()'.
// Category dimensions of binned datasets get one bin per defined category state
x.setBins(10) ;
y.setBins(10) ;
RooDataHist dh("dh","binned version of d",RooArgSet(x,y),d) ;
RooPlot* yframe = y.frame(Bins(10),Title("Operations on binned datasets")) ;
dh.plotOn(yframe) ; // plot projection of 2D binned data on y
// Reduce the 2-dimensional binned dataset to a 1-dimensional binned dataset
//
// All reduce() methods are interfaced in RooAbsData. All reduction techniques
// demonstrated on unbinned datasets can be applied to binned datasets as well.
std::unique_ptr<RooAbsData> dh2{dh.reduce(y,"x>0")};
// Add dh2 to yframe and redraw
dh2->plotOn(yframe,LineColor(kRed),MarkerColor(kRed),Name("dh2")) ;
regPlot(yframe,"rf402_plot1") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'DATA AND CATEGORIES' RooFit tutorial macro #403
//
// Using weights in unbinned datasets
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic403 : public RooUnitTest
{
public:
TestBasic403(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Fits with weighted datasets",refFile,writeRef,verbose) {} ;
bool testCode() {
// C r e a t e o b s e r v a b l e a n d u n w e i g h t e d d a t a s e t
// -------------------------------------------------------------------------------
// Declare observable
RooRealVar x("x","x",-10,10) ;
x.setBins(40) ;
// Construction a uniform pdf
RooPolynomial p0("px","px",x) ;
// Sample 1000 events from pdf
std::unique_ptr<RooDataSet> data{p0.generate(x,1000)};
// C a l c u l a t e w e i g h t a n d m a k e d a t a s e t w e i g h t e d
// -----------------------------------------------------------------------------------
// Construct formula to calculate (fake) weight for events
RooFormulaVar wFunc("w","event weight","(x*x+10)",x) ;
// Add column with variable w to previously generated dataset
RooRealVar* w = (RooRealVar*) data->addColumn(wFunc) ;
// Instruct dataset d in interpret w as event weight rather than as observable
RooDataSet dataw(data->GetName(),data->GetTitle(),data.get(),*data->get(),nullptr,w->GetName()) ;
//data->setWeightVar(*w) ;
// U n b i n n e d M L f i t t o w e i g h t e d d a t a
// ---------------------------------------------------------------
// Construction quadratic polynomial pdf for fitting
RooRealVar a0("a0","a0",1) ;
RooRealVar a1("a1","a1",0,-1,1) ;
RooRealVar a2("a2","a2",1,0,10) ;
RooPolynomial p2("p2","p2",x,RooArgList(a0,a1,a2),0) ;
// Fit quadratic polynomial to weighted data
// NOTE: Maximum likelihood fit to weighted data does in general
// NOT result in correct error estimates, unless individual
// event weights represent Poisson statistics themselves.
// In general, parameter error reflect precision of SumOfWeights
// events rather than NumEvents events. See comparisons below
std::unique_ptr<RooFitResult> r_ml_wgt{p2.fitTo(dataw,Save())};
// P l o t w e i g h e d d a t a a n d f i t r e s u l t
// ---------------------------------------------------------------
// Construct plot frame
RooPlot* frame = x.frame(Title("Unbinned ML fit, binned chi^2 fit to weighted data")) ;
// Plot data using sum-of-weights-squared error rather than Poisson errors
dataw.plotOn(frame,DataError(RooAbsData::SumW2)) ;
// Overlay result of 2nd order polynomial fit to weighted data
p2.plotOn(frame) ;
// M L F i t o f p d f t o e q u i v a l e n t u n w e i g h t e d d a t a s e t
// -----------------------------------------------------------------------------------------
// Construct a pdf with the same shape as p0 after weighting
RooGenericPdf genPdf("genPdf","x*x+10",x) ;
// Sample a dataset with the same number of events as data
std::unique_ptr<RooDataSet> data2{genPdf.generate(x,1000)};
// Sample a dataset with the same number of weights as data
std::unique_ptr<RooDataSet> data3{genPdf.generate(x,43000)};
// Fit the 2nd order polynomial to both unweighted datasets and save the results for comparison
std::unique_ptr<RooFitResult> r_ml_unw10{p2.fitTo(*data2,Save())};
std::unique_ptr<RooFitResult> r_ml_unw43{p2.fitTo(*data3,Save())};
// C h i 2 f i t o f p d f t o b i n n e d w e i g h t e d d a t a s e t
// ------------------------------------------------------------------------------------
// Construct binned clone of unbinned weighted dataset
std::unique_ptr<RooDataHist> binnedData{dataw.binnedClone()};
// Perform chi2 fit to binned weighted dataset using sum-of-weights errors
//
// NB: Within the usual approximations of a chi2 fit, a chi2 fit to weighted
// data using sum-of-weights-squared errors does give correct error
// estimates
std::unique_ptr<RooAbsReal> chi2{p2.createChi2(*binnedData)};
RooMinimizer m(*chi2) ;
m.migrad() ;
m.hesse() ;
// Plot chi^2 fit result on frame as well
std::unique_ptr<RooFitResult> r_chi2_wgt{m.save()};
p2.plotOn(frame,LineStyle(kDashed),LineColor(kRed),Name("p2_alt")) ;
// C o m p a r e f i t r e s u l t s o f c h i 2 , M L f i t s t o ( u n ) w e i g h t e d d a t a
// ---------------------------------------------------------------------------------------------------------------
// Note that ML fit on 1Kevt of weighted data is closer to result of ML fit on 43Kevt of unweighted data
// than to 1Kevt of unweighted data, whereas the reference chi^2 fit with SumW2 error gives a result closer to
// that of an unbinned ML fit to 1Kevt of unweighted data.
regResult(std::move(r_ml_unw10),"rf403_ml_unw10") ;
regResult(std::move(r_ml_unw43),"rf403_ml_unw43") ;
regResult(std::move(r_ml_wgt) ,"rf403_ml_wgt") ;
regResult(std::move(r_chi2_wgt) ,"rf403_ml_chi2") ;
regPlot(frame,"rf403_plot1") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'DATA AND CATEGORIES' RooFit tutorial macro #404
//
// Working with RooCategory objects to describe discrete variables
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic404 : public RooUnitTest
{
public:
TestBasic404(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Categories basic functionality",refFile,writeRef,verbose) {} ;
bool testCode() {
// C o n s t r u c t a c a t e g o r y w i t h l a b e l s
// ----------------------------------------------------------------
// Define a category with labels only
RooCategory tagCat("tagCat","Tagging category") ;
tagCat.defineType("Lepton") ;
tagCat.defineType("Kaon") ;
tagCat.defineType("NetTagger-1") ;
tagCat.defineType("NetTagger-2") ;
// C o n s t r u c t a c a t e g o r y w i t h l a b e l s a n d i n d e c e s
// ----------------------------------------------------------------------------------------
// Define a category with explicitly numbered states
RooCategory b0flav("b0flav","B0 flavour eigenstate") ;
b0flav.defineType("B0",-1) ;
b0flav.defineType("B0bar",1) ;
// G e n e r a t e d u m m y d a t a f o r t a b u l a t i o n d e m o
// ----------------------------------------------------------------------------
// Generate a dummy dataset
RooRealVar x("x","x",0,10) ;
std::unique_ptr<RooDataSet> data{RooPolynomial("p","p",x).generate(RooArgSet(x,b0flav,tagCat),10000)};
// P r i n t t a b l e s o f c a t e g o r y c o n t e n t s o f d a t a s e t s
// ------------------------------------------------------------------------------------------
// Tables are equivalent of plots for categories
Roo1DTable* btable = data->table(b0flav) ;
// Create table for subset of events matching cut expression
Roo1DTable* ttable = data->table(tagCat,"x>8.23") ;
// Create table for all (tagCat x b0flav) state combinations
Roo1DTable* bttable = data->table(RooArgSet(tagCat,b0flav)) ;
// Retrieve number of events from table
// Number can be non-integer if source dataset has weighed events
double nb0 = btable->get("B0") ;
regValue(nb0,"rf404_nb0") ;
// Retrieve fraction of events with "Lepton" tag
double fracLep = ttable->getFrac("Lepton") ;
regValue(fracLep,"rf404_fracLep") ;
// D e f i n i n g r a n g e s f o r p l o t t i n g , f i t t i n g o n c a t e g o r i e s
// ------------------------------------------------------------------------------------------------------
// Define named range as comma separated list of labels
tagCat.setRange("good","Lepton,Kaon") ;
// Or add state names one by one
tagCat.addToRange("soso","NetTagger-1") ;
tagCat.addToRange("soso","NetTagger-2") ;
// Use category range in dataset reduction specification
std::unique_ptr<RooAbsData> goodData{data->reduce(CutRange("good"))};
Roo1DTable* gtable = goodData->table(tagCat) ;
regTable(btable,"rf404_btable") ;
regTable(ttable,"rf404_ttable") ;
regTable(bttable,"rf404_bttable") ;
regTable(gtable,"rf404_gtable") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'DATA AND CATEGORIES' RooFit tutorial macro #405
//
// Demonstration of real-->discrete mapping functions
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic405 : public RooUnitTest
{
public:
TestBasic405(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Real-to-category functions",refFile,writeRef,verbose) {} ;
bool testCode() {
// D e f i n e p d f i n x , s a m p l e d a t a s e t i n x
// ------------------------------------------------------------------------
// Define a dummy PDF in x
RooRealVar x("x","x",0,10) ;
RooArgusBG a("a","argus(x)",x,RooRealConstant::value(10),RooRealConstant::value(-1)) ;
// Generate a dummy dataset
std::unique_ptr<RooDataSet> data{a.generate(x,10000)};
// C r e a t e a t h r e s h o l d r e a l - > c a t f u n c t i o n
// --------------------------------------------------------------------------
// A RooThresholdCategory is a category function that maps regions in a real-valued
// input observable observables to state names. At construction time a 'default'
// state name must be specified to which all values of x are mapped that are not
// otherwise assigned
RooThresholdCategory xRegion("xRegion","region of x",x,"Background") ;
// Specify thresholds and state assignments one-by-one.
// Each statement specifies that all values _below_ the given value
// (and above any lower specified threshold) are mapped to the
// category state with the given name
//
// Background | SideBand | Signal | SideBand | Background
// 4.23 5.23 8.23 9.23
xRegion.addThreshold(4.23,"Background") ;
xRegion.addThreshold(5.23,"SideBand") ;
xRegion.addThreshold(8.23,"Signal") ;
xRegion.addThreshold(9.23,"SideBand") ;
// U s e t h r e s h o l d f u n c t i o n t o p l o t d a t a r e g i o n s
// -------------------------------------------------------------------------------------
// Add values of threshold function to dataset so that it can be used as observable
data->addColumn(xRegion) ;
// Make plot of data in x
RooPlot* xframe = x.frame(Title("Demo of threshold and binning mapping functions")) ;
data->plotOn(xframe) ;
// Use calculated category to select sideband data
data->plotOn(xframe,Cut("xRegion==xRegion::SideBand"),MarkerColor(kRed),LineColor(kRed),Name("data_cut")) ;
// C r e a t e a b i n n i n g r e a l - > c a t f u n c t i o n
// ----------------------------------------------------------------------
// A RooBinningCategory is a category function that maps bins of a (named) binning definition
// in a real-valued input observable observables to state names. The state names are automatically
// constructed from the variable name, the binning name and the bin number. If no binning name
// is specified the default binning is mapped
x.setBins(10,"coarse") ;
RooBinningCategory xBins("xBins","coarse bins in x",x,"coarse") ;
// U s e b i n n i n g f u n c t i o n f o r t a b u l a t i o n a n d p l o t t i n g
// -----------------------------------------------------------------------------------------------
// Print table of xBins state multiplicity. Note that xBins does not need to be an observable in data
// it can be a function of observables in data as well
Roo1DTable* xbtable = data->table(xBins) ;
// Add values of xBins function to dataset so that it can be used as observable
RooCategory* xb = (RooCategory*) data->addColumn(xBins) ;
// Define range "alt" as including bins 1,3,5,7,9
xb->setRange("alt","x_coarse_bin1,x_coarse_bin3,x_coarse_bin5,x_coarse_bin7,x_coarse_bin9") ;
// Construct subset of data matching range "alt" but only for the first 5000 events and plot it on the fram
std::unique_ptr<RooAbsData> dataSel{data->reduce(CutRange("alt"),EventRange(0,5000))};
// dataSel->plotOn(xframe,MarkerColor(kGreen),LineColor(kGreen),Name("data_sel")) ;
regTable(xbtable,"rf405_xbtable") ;
regPlot(xframe,"rf405_plot1") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'DATA AND CATEGORIES' RooFit tutorial macro #406
//
// Demonstration of discrete-->discrete (invertable) functions
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic406 : public RooUnitTest
{
public:
TestBasic406(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Category-to-category functions",refFile,writeRef,verbose) {} ;
bool testCode() {
// C o n s t r u c t t w o c a t e g o r i e s
// ----------------------------------------------
// Define a category with labels only
RooCategory tagCat("tagCat","Tagging category") ;
tagCat.defineType("Lepton") ;
tagCat.defineType("Kaon") ;
tagCat.defineType("NetTagger-1") ;
tagCat.defineType("NetTagger-2") ;
// Define a category with explicitly numbered states
RooCategory b0flav("b0flav","B0 flavour eigenstate") ;
b0flav.defineType("B0",-1) ;
b0flav.defineType("B0bar",1) ;
// Construct a dummy dataset with random values of tagCat and b0flav
RooRealVar x("x","x",0,10) ;
RooPolynomial p("p","p",x) ;
std::unique_ptr<RooDataSet> data{p.generate(RooArgSet(x,b0flav,tagCat),10000)};
// C r e a t e a c a t - > c a t m a p p i n g c a t e g o r y
// ---------------------------------------------------------------------
// A RooMappedCategory is category->category mapping function based on string expression
// The constructor takes an input category an a default state name to which unassigned
// states are mapped
RooMappedCategory tcatType("tcatType","tagCat type",tagCat,"Cut based") ;
// Enter fully specified state mappings
tcatType.map("Lepton","Cut based") ;
tcatType.map("Kaon","Cut based") ;
// Enter a wilcard expression mapping
tcatType.map("NetTagger*","Neural Network") ;
// Make a table of the mapped category state multiplicit in data
Roo1DTable* mtable = data->table(tcatType) ;
// C r e a t e a c a t X c a t p r o d u c t c a t e g o r y
// ----------------------------------------------------------------------
// A SUPER-category is 'product' of _lvalue_ categories. The state names of a super
// category is a composite of the state labels of the input categories
RooSuperCategory b0Xtcat("b0Xtcat","b0flav X tagCat",RooArgSet(b0flav,tagCat)) ;
// Make a table of the product category state multiplicity in data
Roo1DTable* stable = data->table(b0Xtcat) ;
// Since the super category is an lvalue, assignment is explicitly possible
b0Xtcat.setLabel("{B0bar;Lepton}") ;
// A MULTI-category is a 'product' of any category (function). The state names of a super
// category is a composite of the state labels of the input categories
RooMultiCategory b0Xttype("b0Xttype","b0flav X tagType",RooArgSet(b0flav,tcatType)) ;
// Make a table of the product category state multiplicity in data
Roo1DTable* xtable = data->table(b0Xttype) ;
regTable(mtable,"rf406_mtable") ;
regTable(stable,"rf406_stable") ;
regTable(xtable,"rf406_xtable") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'ORGANIZATION AND SIMULTANEOUS FITS' RooFit tutorial macro #501
//
// Using simultaneous p.d.f.s to describe simultaneous fits to multiple
// datasets
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic501 : public RooUnitTest
{
public:
TestBasic501(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Simultaneous p.d.f. operator",refFile,writeRef,verbose) {} ;
bool testCode() {
// C r e a t e m o d e l f o r p h y s i c s s a m p l e
// -------------------------------------------------------------
// Create observables
RooRealVar x("x","x",-8,8) ;
// Construct signal pdf
RooRealVar mean("mean","mean",0,-8,8) ;
RooRealVar sigma("sigma","sigma",0.3,0.1,10) ;
RooGaussian gx("gx","gx",x,mean,sigma) ;
// Construct background pdf
RooRealVar a0("a0","a0",-0.1,-1,1) ;
RooRealVar a1("a1","a1",0.004,-1,1) ;
RooChebychev px("px","px",x,RooArgSet(a0,a1)) ;
// Construct composite pdf
RooRealVar f("f","f",0.2,0.,1.) ;
RooAddPdf model("model","model",RooArgList(gx,px),f) ;
// C r e a t e m o d e l f o r c o n t r o l s a m p l e
// --------------------------------------------------------------
// Construct signal pdf.
// NOTE that sigma is shared with the signal sample model
RooRealVar mean_ctl("mean_ctl","mean_ctl",-3,-8,8) ;
RooGaussian gx_ctl("gx_ctl","gx_ctl",x,mean_ctl,sigma) ;
// Construct the background pdf
RooRealVar a0_ctl("a0_ctl","a0_ctl",-0.1,-1,1) ;
RooRealVar a1_ctl("a1_ctl","a1_ctl",0.5,-0.1,1) ;
RooChebychev px_ctl("px_ctl","px_ctl",x,RooArgSet(a0_ctl,a1_ctl)) ;
// Construct the composite model
RooRealVar f_ctl("f_ctl","f_ctl",0.5,0.,1.) ;
RooAddPdf model_ctl("model_ctl","model_ctl",RooArgList(gx_ctl,px_ctl),f_ctl) ;
// G e n e r a t e e v e n t s f o r b o t h s a m p l e s
// ---------------------------------------------------------------
// Generate 1000 events in x and y from model
std::unique_ptr<RooDataSet> data{model.generate(RooArgSet(x),100)};
std::unique_ptr<RooDataSet> data_ctl{model_ctl.generate(RooArgSet(x),2000)};
// C r e a t e i n d e x c a t e g o r y a n d j o i n s a m p l e s
// ---------------------------------------------------------------------------
// Define category to distinguish physics and control samples events
RooCategory sample("sample","sample") ;
sample.defineType("physics") ;
sample.defineType("control") ;
// Construct combined dataset in (x,sample)
RooDataSet combData("combData","combined data",x,Index(sample),Import({{"physics",data.get()}, {"control",data_ctl.get()}})) ;
// C o n s t r u c t a s i m u l t a n e o u s p d f i n ( x , s a m p l e )
// -----------------------------------------------------------------------------------
// Construct a simultaneous pdf using category sample as index
RooSimultaneous simPdf("simPdf","simultaneous pdf",sample) ;
// Associate model with the physics state and model_ctl with the control state
simPdf.addPdf(model,"physics") ;
simPdf.addPdf(model_ctl,"control") ;
// P e r f o r m a s i m u l t a n e o u s f i t
// ---------------------------------------------------
// Perform simultaneous fit of model to data and model_ctl to data_ctl
simPdf.fitTo(combData) ;
// P l o t m o d e l s l i c e s o n d a t a s l i c e s
// ----------------------------------------------------------------
// Make a frame for the physics sample
RooPlot* frame1 = x.frame(Bins(30),Title("Physics sample")) ;
// Plot all data tagged as physics sample
combData.plotOn(frame1,Cut("sample==sample::physics")) ;
// Plot "physics" slice of simultaneous pdf.
// NBL You _must_ project the sample index category with data using ProjWData
// as a RooSimultaneous makes no prediction on the shape in the index category
// and can thus not be integrated
simPdf.plotOn(frame1,Slice(sample,"physics"),ProjWData(sample,combData)) ;
simPdf.plotOn(frame1,Slice(sample,"physics"),Components("px"),ProjWData(sample,combData),LineStyle(kDashed)) ;
// The same plot for the control sample slice
RooPlot* frame2 = x.frame(Bins(30),Title("Control sample")) ;
combData.plotOn(frame2,Cut("sample==sample::control")) ;
simPdf.plotOn(frame2,Slice(sample,"control"),ProjWData(sample,combData)) ;
simPdf.plotOn(frame2,Slice(sample,"control"),Components("px_ctl"),ProjWData(sample,combData),LineStyle(kDashed)) ;
regPlot(frame1,"rf501_plot1") ;
regPlot(frame2,"rf501_plot2") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'ORGANIZATION AND SIMULTANEOUS FITS' RooFit tutorial macro #501
//
// Using simultaneous p.d.f.s to describe simultaneous fits to multiple
// datasets
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic599 : public RooUnitTest
{
public:
TestBasic599(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Workspace and p.d.f. persistence",refFile,writeRef,verbose) {} ;
bool testCode() {
if (_write) {
RooWorkspace *w = new RooWorkspace("TestBasic11_ws") ;
regWS(w,"Basic11_ws") ;
// Build Gaussian PDF in X
RooRealVar x("x","x",-10,10) ;
RooRealVar meanx("meanx","mean of gaussian",-1) ;
RooRealVar sigmax("sigmax","width of gaussian",3) ;
RooGaussian gaussx("gaussx","gaussian PDF",x,meanx,sigmax) ;
// Build Gaussian PDF in Y
RooRealVar y("y","y",-10,10) ;
RooRealVar meany("meany","mean of gaussian",-1) ;
RooRealVar sigmay("sigmay","width of gaussian",3) ;
RooGaussian gaussy("gaussy","gaussian PDF",y,meany,sigmay) ;
// Make product of X and Y
RooProdPdf gaussxy("gaussxy","gaussx*gaussy",RooArgSet(gaussx,gaussy)) ;
// Make flat bkg in X and Y
RooPolynomial flatx("flatx","flatx",x) ;
RooPolynomial flaty("flaty","flaty",x) ;
RooProdPdf flatxy("flatxy","flatx*flaty",RooArgSet(flatx,flaty)) ;
// Make sum of gaussxy and flatxy
RooRealVar frac("frac","frac",0.5,0.,1.) ;
RooAddPdf sumxy("sumxy","sumxy",RooArgList(gaussxy,flatxy),frac) ;
// Store p.d.f in workspace
w->import(gaussx) ;
w->import(gaussxy,RenameConflictNodes("set2")) ;
w->import(sumxy,RenameConflictNodes("set3")) ;
// Make reference plot of GaussX
RooPlot* frame1 = x.frame() ;
gaussx.plotOn(frame1) ;
regPlot(frame1,"Basic11_gaussx_framex") ;
// Make reference plots for GaussXY
RooPlot* frame2 = x.frame() ;
gaussxy.plotOn(frame2) ;
regPlot(frame2,"Basic11_gaussxy_framex") ;
RooPlot* frame3 = y.frame() ;
gaussxy.plotOn(frame3) ;
regPlot(frame3,"Basic11_gaussxy_framey") ;
// Make reference plots for SumXY
RooPlot* frame4 = x.frame() ;
sumxy.plotOn(frame4) ;
regPlot(frame4,"Basic11_sumxy_framex") ;
RooPlot* frame5 = y.frame() ;
sumxy.plotOn(frame5) ;
regPlot(frame5,"Basic11_sumxy_framey") ;
// Analytically convolved p.d.f.s
// Build a simple decay PDF
RooRealVar dt("dt","dt",-20,20) ;
RooRealVar tau("tau","tau",1.548) ;
// Build a gaussian resolution model
RooRealVar bias1("bias1","bias1",0) ;
RooRealVar sigma1("sigma1","sigma1",1) ;
RooGaussModel gm1("gm1","gauss model 1",dt,bias1,sigma1) ;
// Construct a decay PDF, smeared with single gaussian resolution model
RooDecay decay_gm1("decay_gm1","decay",dt,tau,gm1,RooDecay::DoubleSided) ;
// Build another gaussian resolution model
RooRealVar bias2("bias2","bias2",0) ;
RooRealVar sigma2("sigma2","sigma2",5) ;
RooGaussModel gm2("gm2","gauss model 2",dt,bias2,sigma2) ;
// Build a composite resolution model
RooRealVar gm1frac("gm1frac","fraction of gm1",0.5) ;
RooAddModel gmsum("gmsum","sum of gm1 and gm2",RooArgList(gm1,gm2),gm1frac) ;
// Construct a decay PDF, smeared with double gaussian resolution model
RooDecay decay_gmsum("decay_gmsum","decay",dt,tau,gmsum,RooDecay::DoubleSided) ;
w->import(decay_gm1) ;
w->import(decay_gmsum,RenameConflictNodes("set3")) ;
RooPlot* frame6 = dt.frame() ;
decay_gm1.plotOn(frame6) ;
regPlot(frame6,"Basic11_decay_gm1_framedt") ;
RooPlot* frame7 = dt.frame() ;
decay_gmsum.plotOn(frame7) ;
regPlot(frame7,"Basic11_decay_gmsum_framedt") ;
// Construct simultaneous p.d.f
RooCategory cat("cat","cat") ;
cat.defineType("A") ;
cat.defineType("B") ;
RooSimultaneous sim("sim","sim",cat) ;
sim.addPdf(gaussxy,"A") ;
sim.addPdf(flatxy,"B") ;
w->import(sim,RenameConflictNodes("set4")) ;
// Make plot with dummy dataset for index projection
RooDataHist dh("dh","dh",cat) ;
cat.setLabel("A") ;
dh.add(cat) ;
cat.setLabel("B") ;
dh.add(cat) ;
RooPlot* frame8 = x.frame() ;
sim.plotOn(frame8,ProjWData(cat,dh),Project(cat)) ;
regPlot(frame8,"Basic11_sim_framex") ;
} else {
RooWorkspace* w = getWS("Basic11_ws") ;
if (!w) return false ;
// Retrieve p.d.f from workspace
RooAbsPdf* gaussx = w->pdf("gaussx") ;
// Make test plot and offer for comparison against ref plot
RooPlot* frame1 = w->var("x")->frame() ;
gaussx->plotOn(frame1) ;
regPlot(frame1,"Basic11_gaussx_framex") ;
// Retrieve p.d.f from workspace
RooAbsPdf* gaussxy = w->pdf("gaussxy") ;
// Make test plot and offer for comparison against ref plot
RooPlot* frame2 = w->var("x")->frame() ;
gaussxy->plotOn(frame2) ;
regPlot(frame2,"Basic11_gaussxy_framex") ;
// Make test plot and offer for comparison against ref plot
RooPlot* frame3 = w->var("y")->frame() ;
gaussxy->plotOn(frame3) ;
regPlot(frame3,"Basic11_gaussxy_framey") ;
// Retrieve p.d.f from workspace
RooAbsPdf* sumxy = w->pdf("sumxy") ;
// Make test plot and offer for comparison against ref plot
RooPlot* frame4 = w->var("x")->frame() ;
sumxy->plotOn(frame4) ;
regPlot(frame4,"Basic11_sumxy_framex") ;
// Make test plot and offer for comparison against ref plot
RooPlot* frame5 = w->var("y")->frame() ;
sumxy->plotOn(frame5) ;
regPlot(frame5,"Basic11_sumxy_framey") ;
// Retrieve p.d.f from workspace
RooAbsPdf* decay_gm1 = w->pdf("decay_gm1") ;
// Make test plot and offer for comparison against ref plot
RooPlot* frame6 = w->var("dt")->frame() ;
decay_gm1->plotOn(frame6) ;
regPlot(frame6,"Basic11_decay_gm1_framedt") ;
// Retrieve p.d.f from workspace
RooAbsPdf* decay_gmsum = w->pdf("decay_gmsum") ;
// Make test plot and offer for comparison against ref plot
RooPlot* frame7 = w->var("dt")->frame() ;
decay_gmsum->plotOn(frame7) ;
regPlot(frame7,"Basic11_decay_gmsum_framedt") ;
// Retrieve p.d.f. from workspace
RooAbsPdf* sim = w->pdf("sim") ;
RooCategory* cat = w->cat("cat") ;
// Make plot with dummy dataset for index projection
RooPlot* frame8 = w->var("x")->frame() ;
RooDataHist dh("dh","dh",*cat) ;
cat->setLabel("A") ;
dh.add(*cat) ;
cat->setLabel("B") ;
dh.add(*cat) ;
sim->plotOn(frame8,ProjWData(*cat,dh),Project(*cat)) ;
regPlot(frame8,"Basic11_sim_framex") ;
}
// "Workspace persistence"
return true ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #601
//
// Interactive minimization with MINUIT
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic601 : public RooUnitTest
{
public:
TestBasic601(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Interactive Minuit",refFile,writeRef,verbose) {} ;
bool testCode() {
// S e t u p p d f a n d l i k e l i h o o d
// -----------------------------------------------
// Observable
RooRealVar x("x","x",-20,20) ;
// Model (intentional strong correlations)
RooRealVar mean("mean","mean of g1 and g2",0) ;
RooRealVar sigma_g1("sigma_g1","width of g1",3) ;
RooGaussian g1("g1","g1",x,mean,sigma_g1) ;
RooRealVar sigma_g2("sigma_g2","width of g2",4,3.0,6.0) ;
RooGaussian g2("g2","g2",x,mean,sigma_g2) ;
RooRealVar frac("frac","frac",0.5,0.0,1.0) ;
RooAddPdf model("model","model",RooArgList(g1,g2),frac) ;
// Generate 1000 events
std::unique_ptr<RooDataSet> data{model.generate(x,1000)};
// Construct unbinned likelihood
std::unique_ptr<RooAbsReal> nll{model.createNLL(*data)};
// I n t e r a c t i v e m i n i m i z a t i o n , e r r o r a n a l y s i s
// -------------------------------------------------------------------------------
// Create MINUIT interface object
RooMinimizer m(*nll) ;
// Call MIGRAD to minimize the likelihood
m.migrad() ;
// Run HESSE to calculate errors from d2L/dp2
m.hesse() ;
// Run MINOS on sigma_g2 parameter only
m.minos(sigma_g2) ;
// S a v i n g r e s u l t s , c o n t o u r p l o t s
// ---------------------------------------------------------
// Save a snapshot of the fit result. This object contains the initial
// fit parameters, the final fit parameters, the complete correlation
// matrix, the EDM, the minimized FCN , the last MINUIT status code and
// the number of times the RooFit function object has indicated evaluation
// problems (e.g. zero probabilities during likelihood evaluation)
std::unique_ptr<RooFitResult> r{m.save()};
// C h a n g e p a r a m e t e r v a l u e s , f l o a t i n g
// -----------------------------------------------------------------
// At any moment you can manually change the value of a (constant)
// parameter
mean = 0.3 ;
// Rerun MIGRAD,HESSE
m.migrad() ;
m.hesse() ;
// Now fix sigma_g2
sigma_g2.setConstant(true) ;
// Rerun MIGRAD,HESSE
m.migrad() ;
m.hesse() ;
std::unique_ptr<RooFitResult> r2{m.save()};
regResult(std::move(r),"rf601_r") ;
regResult(std::move(r2),"rf601_r2") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #602
//
// Setting up a binning chi^2 fit
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic602 : public RooUnitTest
{
public:
TestBasic602(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Chi2 minimization",refFile,writeRef,verbose) {} ;
bool testCode() {
// S e t u p m o d e l
// ---------------------
// Declare observable x
RooRealVar x("x","x",0,10) ;
// Create two Gaussian PDFs g1(x,mean1,sigma) anf g2(x,mean2,sigma) and their parameters
RooRealVar mean("mean","mean of gaussians",5) ;
RooRealVar sigma1("sigma1","width of gaussians",0.5) ;
RooRealVar sigma2("sigma2","width of gaussians",1) ;
RooGaussian sig1("sig1","Signal component 1",x,mean,sigma1) ;
RooGaussian sig2("sig2","Signal component 2",x,mean,sigma2) ;
// Build Chebychev polynomial p.d.f.
RooRealVar a0("a0","a0",0.5,0.,1.) ;
RooRealVar a1("a1","a1",-0.2,-1.,1.) ;
RooChebychev bkg("bkg","Background",x,RooArgSet(a0,a1)) ;
// Sum the signal components into a composite signal p.d.f.
RooRealVar sig1frac("sig1frac","fraction of component 1 in signal",0.8,0.,1.) ;
RooAddPdf sig("sig","Signal",RooArgList(sig1,sig2),sig1frac) ;
// Sum the composite signal and background
RooRealVar bkgfrac("bkgfrac","fraction of background",0.5,0.,1.) ;
RooAddPdf model("model","g1+g2+a",RooArgList(bkg,sig),bkgfrac) ;
// C r e a t e b i n n e d d a t a s e t
// -----------------------------------------
std::unique_ptr<RooDataSet> d{model.generate(x,10000)};
std::unique_ptr<RooDataHist> dh{d->binnedClone()};
// Construct a chi^2 of the data and the model,
// which is the input probability density scaled
// by the number of events in the dataset
std::unique_ptr<RooAbsReal> chi2{model.createChi2(*dh, DataError(RooAbsData::SumW2), Verbose(false))};
// Use RooMinimizer interface to minimize chi^2
RooMinimizer m(*chi2) ;
m.migrad() ;
m.hesse() ;
std::unique_ptr<RooFitResult> r{ m.save()};
regResult(std::move(r),"rf602_r");
return true ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #604
//
// Fitting with constraints
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic604 : public RooUnitTest
{
public:
TestBasic604(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Auxiliary observable constraints",refFile,writeRef,verbose) {} ;
bool testCode() {
// C r e a t e m o d e l a n d d a t a s e t
// ----------------------------------------------
// Construct a Gaussian p.d.f
RooRealVar x("x","x",-10,10) ;
RooRealVar m("m","m",0,-10,10) ;
RooRealVar s("s","s",2,0.1,10) ;
RooGaussian gauss("gauss","gauss(x,m,s)",x,m,s) ;
// Construct a flat p.d.f (polynomial of 0th order)
RooPolynomial poly("poly","poly(x)",x) ;
// Construct model = f*gauss + (1-f)*poly
RooRealVar f("f","f",0.5,0.,1.) ;
RooAddPdf model("model","model",RooArgSet(gauss,poly),f) ;
// Generate small dataset for use in fitting below
std::unique_ptr<RooDataSet> d{model.generate(x,50)};
// C r e a t e c o n s t r a i n t p d f
// -----------------------------------------
// Construct Gaussian constraint p.d.f on parameter f at 0.8 with resolution of 0.1
RooGaussian fconstraint("fconstraint","fconstraint",f,0.8,0.1) ;
// M E T H O D 1 - A d d i n t e r n a l c o n s t r a i n t t o m o d e l
// -------------------------------------------------------------------------------------
// Multiply constraint term with regular p.d.f using RooProdPdf
// Specify in fitTo() that internal constraints on parameter f should be used
// Multiply constraint with p.d.f
RooProdPdf modelc("modelc","model with constraint",RooArgSet(model,fconstraint)) ;
// Fit modelc without use of constraint term
std::unique_ptr<RooFitResult> r1{modelc.fitTo(*d,Save())};
// Fit modelc with constraint term on parameter f
std::unique_ptr<RooFitResult> r2{modelc.fitTo(*d,Constrain(f),Save())};
// M E T H O D 2 - S p e c i f y e x t e r n a l c o n s t r a i n t w h e n f i t t i n g
// -------------------------------------------------------------------------------------------------------
// Construct another Gaussian constraint p.d.f on parameter f at 0.8 with resolution of 0.1
RooGaussian fconstext("fconstext","fconstext",f,0.2,0.1) ;
// Fit with external constraint
std::unique_ptr<RooFitResult> r3{model.fitTo(*d,ExternalConstraints(fconstext),Save())};
regResult(std::move(r1),"rf604_r1") ;
regResult(std::move(r2),"rf604_r2") ;
regResult(std::move(r3),"rf604_r3") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #605
//
// Working with the profile likelihood estimator
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic605 : public RooUnitTest
{
public:
TestBasic605(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Profile Likelihood operator",refFile,writeRef,verbose) {} ;
bool testCode() {
// C r e a t e m o d e l a n d d a t a s e t
// -----------------------------------------------
// Observable
RooRealVar x("x","x",-20,20) ;
// Model (intentional strong correlations)
RooRealVar mean("mean","mean of g1 and g2",0,-10,10) ;
RooRealVar sigma_g1("sigma_g1","width of g1",3) ;
RooGaussian g1("g1","g1",x,mean,sigma_g1) ;
RooRealVar sigma_g2("sigma_g2","width of g2",4,3.0,6.0) ;
RooGaussian g2("g2","g2",x,mean,sigma_g2) ;
RooRealVar frac("frac","frac",0.5,0.0,1.0) ;
RooAddPdf model("model","model",RooArgList(g1,g2),frac) ;
// Generate 1000 events
std::unique_ptr<RooDataSet> data{model.generate(x,1000)};
// C o n s t r u c t p l a i n l i k e l i h o o d
// ---------------------------------------------------
// Construct unbinned likelihood
std::unique_ptr<RooAbsReal> nll{model.createNLL(*data)};
nll->SetName("nll");
// Minimize likelihood w.r.t all parameters before making plots
RooMinimizer(*nll).migrad() ;
// Plot likelihood scan frac
RooPlot* frame1 = frac.frame(Bins(10),Range(0.01,0.95),Title("LL and profileLL in frac")) ;
nll->plotOn(frame1,ShiftToZero()) ;
// Plot likelihood scan in sigma_g2
RooPlot* frame2 = sigma_g2.frame(Bins(10),Range(3.3,5.0),Title("LL and profileLL in sigma_g2")) ;
nll->plotOn(frame2,ShiftToZero()) ;
// C o n s t r u c t p r o f i l e l i k e l i h o o d i n f r a c
// -----------------------------------------------------------------------
// The profile likelihood estimator on nll for frac will minimize nll w.r.t
// all floating parameters except frac for each evaluation
RooProfileLL pll_frac("pll_frac","pll_frac",*nll,frac) ;
// Plot the profile likelihood in frac
pll_frac.plotOn(frame1,LineColor(kRed)) ;
// Adjust frame maximum for visual clarity
frame1->SetMinimum(0) ;
frame1->SetMaximum(3) ;
// C o n s t r u c t p r o f i l e l i k e l i h o o d i n s i g m a _ g 2
// -------------------------------------------------------------------------------
// The profile likelihood estimator on nll for sigma_g2 will minimize nll
// w.r.t all floating parameters except sigma_g2 for each evaluation
RooProfileLL pll_sigmag2("pll_sigmag2","pll_sigmag2",*nll,sigma_g2) ;
// Plot the profile likelihood in sigma_g2
pll_sigmag2.plotOn(frame2,LineColor(kRed)) ;
// Adjust frame maximum for visual clarity
frame2->SetMinimum(0) ;
frame2->SetMaximum(3) ;
regPlot(frame1,"rf605_plot1") ;
regPlot(frame2,"rf605_plot2") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #606
//
// Understanding and customizing error handling in likelihood evaluations
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic606 : public RooUnitTest
{
public:
TestBasic606(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("NLL error handling",refFile,writeRef,verbose) {} ;
bool testCode() {
// C r e a t e m o d e l a n d d a t a s e t
// ----------------------------------------------
// Observable
RooRealVar m("m","m",5.20,5.30) ;
// Parameters
RooRealVar m0("m0","m0",5.291,5.20,5.30) ;
RooRealVar k("k","k",-30,-50,-10) ;
// Pdf
RooArgusBG argus("argus","argus",m,m0,k) ;
// Sample 1000 events in m from argus
std::unique_ptr<RooDataSet> data{argus.generate(m,1000)};
// P l o t m o d e l a n d d a t a
// --------------------------------------
RooPlot* frame1 = m.frame(Bins(40),Title("Argus model and data")) ;
data->plotOn(frame1) ;
argus.plotOn(frame1) ;
// F i t m o d e l t o d a t a
// ---------------------------------
argus.fitTo(*data,PrintEvalErrors(10),Warnings(false)) ;
m0.setError(0.1) ;
argus.fitTo(*data,PrintEvalErrors(0),EvalErrorWall(false),Warnings(false)) ;
// P l o t l i k e l i h o o d a s f u n c t i o n o f m 0
// ------------------------------------------------------------------
// Construct likelihood function of model and data
std::unique_ptr<RooAbsReal> nll{argus.createNLL(*data)};
nll->SetName("nll");
// Plot likelihood in m0 in range that includes problematic values
// In this configuration no messages are printed for likelihood evaluation errors,
// but if an likelihood value evaluates with error, the corresponding value
// on the curve will be set to the value given in EvalErrorValue().
RooPlot* frame2 = m0.frame(Range(5.288,5.293),Title("-log(L) scan vs m0, problematic regions masked")) ;
nll->plotOn(frame2,PrintEvalErrors(-1),ShiftToZero(),EvalErrorValue(nll->getVal()+10),LineColor(kRed)) ;
regPlot(frame1,"rf606_plot1") ;
regPlot(frame2,"rf606_plot3") ; // 3 is the reference of the plot
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #607
//
// Demonstration of options of the RooFitResult class
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic607 : public RooUnitTest
{
public:
TestBasic607(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Fit Result functionality",refFile,writeRef,verbose) {} ;
bool testCode() {
// C r e a t e p d f , d a t a
// --------------------------------
// Declare observable x
RooRealVar x("x","x",0,10) ;
// Create two Gaussian PDFs g1(x,mean1,sigma) anf g2(x,mean2,sigma) and their parameters
RooRealVar mean("mean","mean of gaussians",5,-10,10) ;
RooRealVar sigma1("sigma1","width of gaussians",0.5,0.1,10) ;
RooRealVar sigma2("sigma2","width of gaussians",1,0.1,10) ;
RooGaussian sig1("sig1","Signal component 1",x,mean,sigma1) ;
RooGaussian sig2("sig2","Signal component 2",x,mean,sigma2) ;
// Build Chebychev polynomial p.d.f.
RooRealVar a0("a0","a0",0.5,0.,1.) ;
RooRealVar a1("a1","a1",-0.2,-1.,1.) ;
RooChebychev bkg("bkg","Background",x,RooArgSet(a0,a1)) ;
// Sum the signal components into a composite signal p.d.f.
RooRealVar sig1frac("sig1frac","fraction of component 1 in signal",0.8,0.,1.) ;
RooAddPdf sig("sig","Signal",RooArgList(sig1,sig2),sig1frac) ;
// Sum the composite signal and background
RooRealVar bkgfrac("bkgfrac","fraction of background",0.5,0.,1.) ;
RooAddPdf model("model","g1+g2+a",RooArgList(bkg,sig),bkgfrac) ;
// Generate 1000 events
std::unique_ptr<RooDataSet> data{model.generate(x,1000)};
// F i t p d f t o d a t a , s a v e f i t r e s u l t
// -------------------------------------------------------------
// Perform fit and save result
std::unique_ptr<RooFitResult> r{model.fitTo(*data,Save())};
// V i s u a l i z e c o r r e l a t i o n m a t r i x
// -------------------------------------------------------
// Construct 2D color plot of correlation matrix
gStyle->SetOptStat(0) ;
gStyle->SetPalette(1) ;
TH2* hcorr = r->correlationHist() ;
// Sample dataset with parameter values according to distribution
// of covariance matrix of fit result
RooDataSet randPars("randPars","randPars",r->floatParsFinal()) ;
for (int i=0 ; i<10000 ; i++) {
randPars.add(r->randomizePars()) ;
}
// make histogram of 2D distribution in sigma1 vs sig1frac
TH1* hhrand = randPars.createHistogram("hhrand",sigma1,Binning(35,0.25,0.65),YVar(sig1frac,Binning(35,0.3,1.1))) ;
regTH(hcorr,"rf607_hcorr") ;
regTH(hhrand,"rf607_hhand") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #609
//
// Setting up a chi^2 fit to an unbinned dataset with X,Y,err(Y)
// values (and optionally err(X) values)
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic609 : public RooUnitTest
{
public:
TestBasic609(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Chi^2 fit to X-Y dataset",refFile,writeRef,verbose) {} ;
bool testCode() {
// C r e a t e d a t a s e t w i t h X a n d Y v a l u e s
// -------------------------------------------------------------------
// Make weighted XY dataset with asymmetric errors stored
// The StoreError() argument is essential as it makes
// the dataset store the error in addition to the values
// of the observables. If errors on one or more observables
// are asymmetric, one can store the asymmetric error
// using the StoreAsymError() argument
RooRealVar x("x","x",-11,11) ;
RooRealVar y("y","y",-10,200) ;
RooDataSet dxy("dxy","dxy",RooArgSet(x,y),StoreError(RooArgSet(x,y))) ;
// Fill an example dataset with X,err(X),Y,err(Y) values
for (int i=0 ; i<=10 ; i++) {
// Set X value and error
x = -10 + 2*i;
x.setError( i<5 ? 0.5/1. : 1.0/1. ) ;
// Set Y value and error
y = x.getVal() * x.getVal() + 4*fabs(RooRandom::randomGenerator()->Gaus()) ;
y.setError(sqrt(y.getVal())) ;
dxy.add(RooArgSet(x,y)) ;
}
// P e r f o r m c h i 2 f i t t o X + / - d x a n d Y + / - d Y v a l u e s
// ---------------------------------------------------------------------------------------
// Make fit function
RooRealVar a("a","a",0.0,-10,10) ;
RooRealVar b("b","b",0.0,-100,100) ;
RooRealVar c("c", "c", 0.0, -100, 100);
RooPolyVar f("f", "f", x, RooArgList(b, a, c));
// Plot dataset in X-Y interpretation
RooPlot* frame = x.frame(Title("Chi^2 fit of function set of (X#pmdX,Y#pmdY) values")) ;
dxy.plotOnXY(frame,YVar(y)) ;
// Fit chi^2 using X and Y errors
std::unique_ptr<RooFitResult> fit1{f.chi2FitTo(dxy, YVar(y), Save(), PrintLevel(-1))};
// Overlay fitted function
f.plotOn(frame) ;
// Alternative: fit chi^2 integrating f(x) over ranges defined by X errors, rather
// than taking point at center of bin
std::unique_ptr<RooFitResult> fit2{f.chi2FitTo(dxy, YVar(y), Save(), PrintLevel(-1), Integrate(true))};
// Overlay alternate fit result
f.plotOn(frame,LineStyle(kDashed),LineColor(kRed),Name("alternate")) ;
regPlot(frame,"rf609_frame") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'SPECIAL PDFS' RooFit tutorial macro #701
//
// Unbinned maximum likelihood fit of an efficiency eff(x) function to
// a dataset D(x,cut), where cut is a category encoding a selection, of which
// the efficiency as function of x should be described by eff(x)
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic701 : public RooUnitTest
{
public:
TestBasic701(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Efficiency operator p.d.f. 1D",refFile,writeRef,verbose) {} ;
bool testCode() {
// C o n s t r u c t e f f i c i e n c y f u n c t i o n e ( x )
// -------------------------------------------------------------------
// Declare variables x,mean,sigma with associated name, title, initial value and allowed range
RooRealVar x("x","x",-10,10) ;
// Efficiency function eff(x;a,b)
RooRealVar a("a","a",0.4,0,1) ;
RooRealVar b("b","b",5) ;
RooRealVar c("c","c",-1,-10,10) ;
RooFormulaVar effFunc("effFunc","(1-a)+a*cos((x-c)/b)",RooArgList(a,b,c,x)) ;
// C o n s t r u c t c o n d i t i o n a l e f f i c i e n c y p d f E ( c u t | x )
// ------------------------------------------------------------------------------------------
// Acceptance state cut (1 or 0)
RooCategory cut("cut","cutr") ;
cut.defineType("accept",1) ;
cut.defineType("reject",0) ;
// Construct efficiency p.d.f eff(cut|x)
RooEfficiency effPdf("effPdf","effPdf",effFunc,cut,"accept") ;
// G e n e r a t e d a t a ( x , c u t ) f r o m a t o y m o d e l
// -----------------------------------------------------------------------------
// Construct global shape p.d.f shape(x) and product model(x,cut) = eff(cut|x)*shape(x)
// (These are _only_ needed to generate some toy MC here to be used later)
RooPolynomial shapePdf("shapePdf","shapePdf",x,RooConst(-0.095)) ;
RooProdPdf model("model","model",shapePdf,Conditional(effPdf,cut)) ;
// Generate some toy data from model
std::unique_ptr<RooDataSet> data{model.generate(RooArgSet(x,cut),10000)};
// F i t c o n d i t i o n a l e f f i c i e n c y p d f t o d a t a
// --------------------------------------------------------------------------
// Fit conditional efficiency p.d.f to data
effPdf.fitTo(*data,ConditionalObservables(x)) ;
// P l o t f i t t e d , d a t a e f f i c i e n c y
// --------------------------------------------------------
// Plot distribution of all events and accepted fraction of events on frame
RooPlot* frame1 = x.frame(Bins(20),Title("Data (all, accepted)")) ;
data->plotOn(frame1) ;
data->plotOn(frame1,Cut("cut==cut::accept"),MarkerColor(kRed),LineColor(kRed)) ;
// Plot accept/reject efficiency on data overlay fitted efficiency curve
RooPlot* frame2 = x.frame(Bins(20),Title("Fitted efficiency")) ;
data->plotOn(frame2,Efficiency(cut)) ; // needs ROOT version >= 5.21
effFunc.plotOn(frame2,LineColor(kRed)) ;
regPlot(frame1,"rf701_plot1") ;
regPlot(frame2,"rf701_plot2") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'SPECIAL PDFS' RooFit tutorial macro #702
//
// Unbinned maximum likelihood fit of an efficiency eff(x) function to
// a dataset D(x,cut), where cut is a category encoding a selection whose
// efficiency as function of x should be described by eff(x)
//
//
/////////////////////////////////////////////////////////////////////////
class TestBasic702 : public RooUnitTest
{
public:
TestBasic702(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Efficiency operator p.d.f. 2D",refFile,writeRef,verbose) {} ;
bool testCode() {
bool flat=false ;
// C o n s t r u c t e f f i c i e n c y f u n c t i o n e ( x , y )
// -----------------------------------------------------------------------
// Declare variables x,mean,sigma with associated name, title, initial value and allowed range
RooRealVar x("x","x",-10,10) ;
RooRealVar y("y","y",-10,10) ;
// Efficiency function eff(x;a,b)
RooRealVar ax("ax","ay",0.6,0,1) ;
RooRealVar bx("bx","by",5) ;
RooRealVar cx("cx","cy",-1,-10,10) ;
RooRealVar ay("ay","ay",0.2,0,1) ;
RooRealVar by("by","by",5) ;
RooRealVar cy("cy","cy",-1,-10,10) ;
RooFormulaVar effFunc("effFunc","((1-ax)+ax*cos((x-cx)/bx))*((1-ay)+ay*cos((y-cy)/by))",RooArgList(ax,bx,cx,x,ay,by,cy,y)) ;
// Acceptance state cut (1 or 0)
RooCategory cut("cut","cutr") ;
cut.defineType("accept",1) ;
cut.defineType("reject",0) ;
// C o n s t r u c t c o n d i t i o n a l e f f i c i e n c y p d f E ( c u t | x , y )
// ---------------------------------------------------------------------------------------------
// Construct efficiency p.d.f eff(cut|x)
RooEfficiency effPdf("effPdf","effPdf",effFunc,cut,"accept") ;
// G e n e r a t e d a t a ( x , y , c u t ) f r o m a t o y m o d e l
// -------------------------------------------------------------------------------
// Construct global shape p.d.f shape(x) and product model(x,cut) = eff(cut|x)*shape(x)
// (These are _only_ needed to generate some toy MC here to be used later)
RooPolynomial shapePdfX("shapePdfX","shapePdfX",x,RooConst(flat?0:-0.095)) ;
RooPolynomial shapePdfY("shapePdfY","shapePdfY",y,RooConst(flat?0:+0.095)) ;
RooProdPdf shapePdf("shapePdf","shapePdf",RooArgSet(shapePdfX,shapePdfY)) ;
RooProdPdf model("model","model",shapePdf,Conditional(effPdf,cut)) ;
// Generate some toy data from model
std::unique_ptr<RooDataSet> data{model.generate(RooArgSet(x,y,cut),10000)};
// F i t c o n d i t i o n a l e f f i c i e n c y p d f t o d a t a
// --------------------------------------------------------------------------
// Fit conditional efficiency p.d.f to data
effPdf.fitTo(*data,ConditionalObservables(RooArgSet(x,y))) ;
// P l o t f i t t e d , d a t a e f f i c i e n c y
// --------------------------------------------------------
// Make 2D histograms of all data, selected data and efficiency function
TH1* hh_data_all = data->createHistogram("hh_data_all",x,Binning(8),YVar(y,Binning(8))) ;
TH1* hh_data_sel = data->createHistogram("hh_data_sel",x,Binning(8),YVar(y,Binning(8)),Cut("cut==cut::accept")) ;
TH1* hh_eff = effFunc.createHistogram("hh_eff",x,Binning(50),YVar(y,Binning(50))) ;
// Some adjustsment for good visualization
hh_data_all->SetMinimum(0) ;
hh_data_sel->SetMinimum(0) ;
hh_eff->SetMinimum(0) ;
hh_eff->SetLineColor(kBlue) ;
regTH(hh_data_all,"rf702_hh_data_all") ;
regTH(hh_data_sel,"rf702_hh_data_sel") ;
regTH(hh_eff,"rf702_hh_eff") ;
return true;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'SPECIAL PDFS' RooFit tutorial macro #703
//
// Using a product of an (acceptance) efficiency and a p.d.f as p.d.f.
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic703 : public RooUnitTest
{
public:
TestBasic703(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Efficiency product operator p.d.f",refFile,writeRef,verbose) {} ;
bool testCode() {
// D e f i n e o b s e r v a b l e s a n d d e c a y p d f
// ---------------------------------------------------------------
// Declare observables
RooRealVar t("t","t",0,5) ;
// Make pdf
RooRealVar tau("tau","tau",-1.54,-4,-0.1) ;
RooExponential model("model","model",t,tau) ;
// D e f i n e e f f i c i e n c y f u n c t i o n
// ---------------------------------------------------
// Use error function to simulate turn-on slope
RooFormulaVar eff("eff","0.5*(TMath::Erf((t-1)/0.5)+1)",t) ;
// D e f i n e d e c a y p d f w i t h e f f i c i e n c y
// ---------------------------------------------------------------
// Multiply pdf(t) with efficiency in t
RooEffProd modelEff("modelEff","model with efficiency",model,eff) ;
// P l o t e f f i c i e n c y , p d f
// ----------------------------------------
RooPlot* frame1 = t.frame(Title("Efficiency")) ;
eff.plotOn(frame1,LineColor(kRed)) ;
RooPlot* frame2 = t.frame(Title("Pdf with and without efficiency")) ;
model.plotOn(frame2,LineStyle(kDashed)) ;
modelEff.plotOn(frame2) ;
// G e n e r a t e t o y d a t a , f i t m o d e l E f f t o d a t a
// ------------------------------------------------------------------------------
// Generate events. If the input pdf has an internal generator, the internal generator
// is used and an accept/reject sampling on the efficiency is applied.
std::unique_ptr<RooDataSet> data{modelEff.generate(t,10000)};
// Fit pdf. The normalization integral is calculated numerically.
modelEff.fitTo(*data) ;
// Plot generated data and overlay fitted pdf
RooPlot* frame3 = t.frame(Title("Fitted pdf with efficiency")) ;
data->plotOn(frame3) ;
modelEff.plotOn(frame3) ;
regPlot(frame1,"rf703_plot1") ;
regPlot(frame2,"rf703_plot2") ;
regPlot(frame3,"rf703_plot3") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'SPECIAL PDFS' RooFit tutorial macro #704
//
// Using a p.d.f defined by a sum of real-valued amplitude components
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic704 : public RooUnitTest
{
public:
TestBasic704(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Amplitude sum operator p.d.f",refFile,writeRef,verbose) {} ;
bool testCode() {
// S e t u p 2 D a m p l i t u d e f u n c t i o n s
// -------------------------------------------------------
// Observables
RooRealVar t("t","time",-1.,15.);
RooRealVar cosa("cosa","cos(alpha)",-1.,1.);
RooRealVar tau("tau", "#tau", 1.5);
RooRealVar deltaGamma("deltaGamma","deltaGamma", 0.3);
RooFormulaVar coshG("coshGBasis","exp(-@0/ @1)*cosh(@0*@2/2)", {t,tau,deltaGamma});
RooFormulaVar sinhG("sinhGBasis","exp(-@0/ @1)*sinh(@0*@2/2)", {t,tau,deltaGamma});
// Construct polynomial amplitudes in cos(a)
RooPolyVar poly1("poly1","poly1",cosa,RooArgList(0.5, 0.2, 0.2),0);
RooPolyVar poly2("poly2","poly2",cosa,RooArgList(1.0, -0.2, 3.0),0);
// Construct 2D amplitude as uncorrelated product of amp(t)*amp(cosa)
RooProduct ampl1("ampl1","amplitude 1", {poly1, coshG});
RooProduct ampl2("ampl2","amplitude 2", {poly2, sinhG});
// C o n s t r u c t a m p l i t u d e s u m p d f
// -----------------------------------------------------
// Amplitude strengths
RooRealVar f1("f1","f1",1,0,2) ;
RooRealVar f2("f2","f2",0.5,0,2) ;
// Construct pdf
RooRealSumPdf pdf("pdf","pdf",RooArgList(ampl1,ampl2),RooArgList(f1,f2)) ;
// Generate some toy data from pdf
std::unique_ptr<RooDataSet> data{pdf.generate(RooArgSet(t,cosa),10000)};
// Fit pdf to toy data with only amplitude strength floating
pdf.fitTo(*data) ;
// P l o t a m p l i t u d e s u m p d f
// -------------------------------------------
// Make 2D plots of amplitudes
TH1* hh_cos = ampl1.createHistogram("hh_cos",t,Binning(50),YVar(cosa,Binning(50))) ;
TH1* hh_sin = ampl2.createHistogram("hh_sin",t,Binning(50),YVar(cosa,Binning(50))) ;
hh_cos->SetLineColor(kBlue) ;
hh_sin->SetLineColor(kBlue) ;
// Make projection on t, plot data, pdf and its components
// Note component projections may be larger than sum because amplitudes can be negative
RooPlot* frame1 = t.frame();
data->plotOn(frame1);
pdf.plotOn(frame1);
pdf.plotOn(frame1,Components(ampl1),LineStyle(kDashed));
pdf.plotOn(frame1,Components(ampl2),LineStyle(kDashed),LineColor(kRed));
// Make projection on cosa, plot data, pdf and its components
// Note that components projection may be larger than sum because amplitudes can be negative
RooPlot* frame2 = cosa.frame();
data->plotOn(frame2);
pdf.plotOn(frame2);
pdf.plotOn(frame2,Components(ampl1),LineStyle(kDashed));
pdf.plotOn(frame2,Components(ampl2),LineStyle(kDashed),LineColor(kRed));
regPlot(frame1,"rf704_plot1") ;
regPlot(frame2,"rf704_plot2") ;
regTH(hh_cos,"rf704_hh_cos") ;
regTH(hh_sin,"rf704_hh_sin") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'SPECIAL PDFS' RooFit tutorial macro #705
//
// Linear interpolation between p.d.f shapes using the 'Alex Read' algorithm
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic705 : public RooUnitTest
{
public:
double ctol() { return 5e-2 ; } // very conservative, this is a numerically difficult test
TestBasic705(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Linear morph operator p.d.f.",refFile,writeRef,verbose) {} ;
bool testCode() {
// C r e a t e e n d p o i n t p d f s h a p e s
// ------------------------------------------------------
// Observable
RooRealVar x("x","x",-20,20) ;
// Lower end point shape: a Gaussian
RooRealVar g1mean("g1mean","g1mean",-10) ;
RooGaussian g1("g1","g1",x,g1mean,2.0) ;
// Upper end point shape: a Polynomial
RooPolynomial g2("g2","g2",x,RooArgSet(-0.03,-0.001)) ;
// C r e a t e i n t e r p o l a t i n g p d f
// -----------------------------------------------
// Create interpolation variable
RooRealVar alpha("alpha","alpha",0,1.0) ;
// Specify sampling density on observable and interpolation variable
x.setBins(1000,"cache") ;
alpha.setBins(50,"cache") ;
// Construct interpolating pdf in (x,a) represent g1(x) at a=a_min
// and g2(x) at a=a_max
RooIntegralMorph lmorph("lmorph","lmorph",g1,g2,x,alpha) ;
// P l o t i n t e r p o l a t i n g p d f a t v a r i o u s a l p h a
// -----------------------------------------------------------------------------
// Show end points as blue curves
RooPlot* frame1 = x.frame() ;
g1.plotOn(frame1) ;
g2.plotOn(frame1) ;
// Show interpolated shapes in red
alpha.setVal(0.125) ;
lmorph.plotOn(frame1,LineColor(kRed),Name("alt1")) ;
alpha.setVal(0.25) ;
lmorph.plotOn(frame1,LineColor(kRed),Name("alt2")) ;
alpha.setVal(0.375) ;
lmorph.plotOn(frame1,LineColor(kRed),Name("alt3")) ;
alpha.setVal(0.50) ;
lmorph.plotOn(frame1,LineColor(kRed),Name("alt4")) ;
alpha.setVal(0.625) ;
lmorph.plotOn(frame1,LineColor(kRed),Name("alt5")) ;
alpha.setVal(0.75) ;
lmorph.plotOn(frame1,LineColor(kRed),Name("alt6")) ;
alpha.setVal(0.875) ;
lmorph.plotOn(frame1,LineColor(kRed),Name("alt7")) ;
alpha.setVal(0.95) ;
lmorph.plotOn(frame1,LineColor(kRed),Name("alt8")) ;
// S h o w 2 D d i s t r i b u t i o n o f p d f ( x , a l p h a )
// -----------------------------------------------------------------------
// Create 2D histogram
TH1* hh = lmorph.createHistogram("hh",x,Binning(40),YVar(alpha,Binning(40))) ;
hh->SetLineColor(kBlue) ;
// F i t p d f t o d a t a s e t w i t h a l p h a = 0 . 8
// -----------------------------------------------------------------
// Generate a toy dataset at alpha = 0.8
alpha=0.8 ;
std::unique_ptr<RooDataSet> data{lmorph.generate(x,1000)};
// Fit pdf to toy data
lmorph.setCacheAlpha(true) ;
{
// Get rid of the caching info prints
RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING, 0u, RooFit::Caching, true};
lmorph.fitTo(*data) ;
}
// Plot fitted pdf and data overlaid
RooPlot* frame2 = x.frame(Bins(100)) ;
data->plotOn(frame2) ;
lmorph.plotOn(frame2) ;
// S c a n - l o g ( L ) v s a l p h a
// -----------------------------------------
// Show scan -log(L) of dataset w.r.t alpha
RooPlot* frame3 = alpha.frame(Bins(100),Range(0.5,0.9)) ;
// Make 2D pdf of histogram
{
// Get rid of the caching info prints
RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING, 0u, RooFit::Caching, true};
std::unique_ptr<RooAbsReal> nll{lmorph.createNLL(*data)};
nll->SetName("nll");
nll->plotOn(frame3,ShiftToZero()) ;
}
lmorph.setCacheAlpha(false) ;
regPlot(frame1,"rf705_plot1") ;
regPlot(frame2,"rf705_plot2") ;
regPlot(frame3,"rf705_plot3") ;
regTH(hh,"rf705_hh") ;
return true;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'SPECIAL PDFS' RooFit tutorial macro #706
//
// Histogram based p.d.f.s and functions
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic706 : public RooUnitTest
{
public:
TestBasic706(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Histogram based p.d.f.s",refFile,writeRef,verbose) {} ;
bool testCode() {
// C r e a t e p d f f o r s a m p l i n g
// ---------------------------------------------
RooRealVar x("x","x",0,20) ;
RooPolynomial p("p","p",x,RooArgList(0.01,-0.01,0.0004)) ;
// C r e a t e l o w s t a t s h i s t o g r a m
// ---------------------------------------------------
// Sample 500 events from p
x.setBins(20) ;
std::unique_ptr<RooDataSet> data1{p.generate(x,500)};
// Create a binned dataset with 20 bins and 500 events
std::unique_ptr<RooDataHist> hist1{data1->binnedClone()};
// Represent data in dh as pdf in x
RooHistPdf histpdf1("histpdf1","histpdf1",x,*hist1,0) ;
// Plot unbinned data and histogram pdf overlaid
RooPlot* frame1 = x.frame(Title("Low statistics histogram pdf"),Bins(100)) ;
data1->plotOn(frame1) ;
histpdf1.plotOn(frame1) ;
// C r e a t e h i g h s t a t s h i s t o g r a m
// -----------------------------------------------------
// Sample 100000 events from p
x.setBins(10) ;
std::unique_ptr<RooDataSet> data2{p.generate(x,100000)};
// Create a binned dataset with 10 bins and 100K events
std::unique_ptr<RooDataHist> hist2{data2->binnedClone()};
// Represent data in dh as pdf in x, apply 2nd order interpolation
RooHistPdf histpdf2("histpdf2","histpdf2",x,*hist2,2) ;
// Plot unbinned data and histogram pdf overlaid
RooPlot* frame2 = x.frame(Title("High stats histogram pdf with interpolation"),Bins(100)) ;
data2->plotOn(frame2) ;
histpdf2.plotOn(frame2) ;
regPlot(frame1,"rf607_plot1") ;
regPlot(frame2,"rf607_plot2") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'SPECIAL PDFS' RooFit tutorial macro #707
//
// Using non-parametric (multi-dimensional) kernel estimation p.d.f.s
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic707 : public RooUnitTest
{
public:
TestBasic707(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Kernel estimation p.d.f.s",refFile,writeRef,verbose) {} ;
bool testCode() {
// C r e a t e l o w s t a t s 1 - D d a t a s e t
// -------------------------------------------------------
// Create a toy pdf for sampling
RooRealVar x("x","x",0,20) ;
RooPolynomial p("p","p",x,RooArgList(0.01,-0.01,0.0004)) ;
// Sample 500 events from p
std::unique_ptr<RooDataSet> data1{p.generate(x,200)};
// C r e a t e 1 - D k e r n e l e s t i m a t i o n p d f
// ---------------------------------------------------------------
// Create adaptive kernel estimation pdf. In this configuration the input data
// is mirrored over the boundaries to minimize edge effects in distribution
// that do not fall to zero towards the edges
RooKeysPdf kest1("kest1","kest1",x,*data1,RooKeysPdf::MirrorBoth) ;
// An adaptive kernel estimation pdf on the same data without mirroring option
// for comparison
RooKeysPdf kest2("kest2","kest2",x,*data1,RooKeysPdf::NoMirror) ;
// Adaptive kernel estimation pdf with increased bandwidth scale factor
// (promotes smoothness over detail preservation)
RooKeysPdf kest3("kest3","kest3",x,*data1,RooKeysPdf::MirrorBoth,2) ;
// Plot kernel estimation pdfs with and without mirroring over data
RooPlot* frame = x.frame(Title("Adaptive kernel estimation pdf with and w/o mirroring"),Bins(20)) ;
data1->plotOn(frame) ;
kest1.plotOn(frame) ;
kest2.plotOn(frame,LineStyle(kDashed),LineColor(kRed)) ;
// Plot kernel estimation pdfs with regular and increased bandwidth
RooPlot* frame2 = x.frame(Title("Adaptive kernel estimation pdf with regular, increased bandwidth")) ;
kest1.plotOn(frame2) ;
kest3.plotOn(frame2,LineColor(kMagenta)) ;
// C r e a t e l o w s t a t s 2 - D d a t a s e t
// -------------------------------------------------------
// Construct a 2D toy pdf for sampleing
RooRealVar y("y","y",0,20) ;
RooPolynomial py("py","py",y,RooArgList(0.01,0.01,-0.0004)) ;
RooProdPdf pxy("pxy","pxy",RooArgSet(p,py)) ;
std::unique_ptr<RooDataSet> data2{pxy.generate(RooArgSet(x,y),1000)};
// C r e a t e 2 - D k e r n e l e s t i m a t i o n p d f
// ---------------------------------------------------------------
// Create 2D adaptive kernel estimation pdf with mirroring
RooNDKeysPdf kest4("kest4","kest4",RooArgSet(x,y),*data2,"am") ;
// Create 2D adaptive kernel estimation pdf with mirroring and double bandwidth
RooNDKeysPdf kest5("kest5","kest5",RooArgSet(x,y),*data2,"am",2) ;
// Create a histogram of the data
TH1* hh_data = data2->createHistogram("hh_data",x,Binning(10),YVar(y,Binning(10))) ;
// Create histogram of the 2d kernel estimation pdfs
TH1* hh_pdf = kest4.createHistogram("hh_pdf",x,Binning(25),YVar(y,Binning(25))) ;
TH1* hh_pdf2 = kest5.createHistogram("hh_pdf2",x,Binning(25),YVar(y,Binning(25))) ;
hh_pdf->SetLineColor(kBlue) ;
hh_pdf2->SetLineColor(kMagenta) ;
regPlot(frame,"rf707_plot1") ;
regPlot(frame2,"rf707_plot2") ;
regTH(hh_data,"rf707_hhdata") ;
regTH(hh_pdf,"rf707_hhpdf") ;
regTH(hh_pdf2,"rf707_hhpdf2") ;
return true ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'SPECIAL PDFS' RooFit tutorial macro #708
//
// Special decay pdf for B physics with mixing and/or CP violation
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic708 : public RooUnitTest
{
public:
TestBasic708(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("B Physics p.d.f.s",refFile,writeRef,verbose) {} ;
bool testCode() {
////////////////////////////////////////////////////
// B - D e c a y w i t h m i x i n g //
////////////////////////////////////////////////////
// C o n s t r u c t p d f
// -------------------------
// Observable
RooRealVar dt("dt","dt",-10,10) ;
dt.setBins(40) ;
// Parameters
RooRealVar dm("dm","delta m(B0)",0.472) ;
RooRealVar tau("tau","tau (B0)",1.547) ;
RooRealVar w("w","flavour mistag rate",0.1) ;
RooRealVar dw("dw","delta mistag rate for B0/B0bar",0.1) ;
RooCategory mixState("mixState","B0/B0bar mixing state") ;
mixState.defineType("mixed",-1) ;
mixState.defineType("unmixed",1) ;
RooCategory tagFlav("tagFlav","Flavour of the tagged B0") ;
tagFlav.defineType("B0",1) ;
tagFlav.defineType("B0bar",-1) ;
// Use delta function resolution model
RooTruthModel tm("tm","truth model",dt) ;
// Construct Bdecay with mixing
RooBMixDecay bmix("bmix","decay",dt,mixState,tagFlav,tau,dm,w,dw,tm,RooBMixDecay::DoubleSided) ;
// P l o t p d f i n v a r i o u s s l i c e s
// ---------------------------------------------------
// Generate some data
std::unique_ptr<RooDataSet> data{bmix.generate(RooArgSet(dt,mixState,tagFlav),10000)};
// Plot B0 and B0bar tagged data separately
// For all plots below B0 and B0 tagged data will look somewhat differently
// if the flavor tagging mistag rate for B0 and B0 is different (i.e. dw!=0)
RooPlot* frame1 = dt.frame(Title("B decay distribution with mixing (B0/B0bar)")) ;
data->plotOn(frame1,Cut("tagFlav==tagFlav::B0")) ;
bmix.plotOn(frame1,Slice(tagFlav,"B0")) ;
data->plotOn(frame1,Cut("tagFlav==tagFlav::B0bar"),MarkerColor(kCyan)) ;
bmix.plotOn(frame1,Slice(tagFlav,"B0bar"),LineColor(kCyan),Name("alt")) ;
// Plot mixed slice for B0 and B0bar tagged data separately
RooPlot* frame2 = dt.frame(Title("B decay distribution of mixed events (B0/B0bar)")) ;
data->plotOn(frame2,Cut("mixState==mixState::mixed&&tagFlav==tagFlav::B0")) ;
bmix.plotOn(frame2,Slice({{&tagFlav,"B0"}, {&mixState,"mixed"}})) ;
data->plotOn(frame2,Cut("mixState==mixState::mixed&&tagFlav==tagFlav::B0bar"),MarkerColor(kCyan)) ;
bmix.plotOn(frame2,Slice({{&tagFlav,"B0bar"}, {&mixState,"mixed"}}),LineColor(kCyan),Name("alt")) ;
// Plot unmixed slice for B0 and B0bar tagged data separately
RooPlot* frame3 = dt.frame(Title("B decay distribution of unmixed events (B0/B0bar)")) ;
data->plotOn(frame3,Cut("mixState==mixState::unmixed&&tagFlav==tagFlav::B0")) ;
bmix.plotOn(frame3,Slice({{&tagFlav,"B0"}, {&mixState,"unmixed"}})) ;
data->plotOn(frame3,Cut("mixState==mixState::unmixed&&tagFlav==tagFlav::B0bar"),MarkerColor(kCyan)) ;
bmix.plotOn(frame3,Slice({{&tagFlav,"B0bar"}, {&mixState,"unmixed"}}),LineColor(kCyan),Name("alt")) ;
///////////////////////////////////////////////////////
// B - D e c a y w i t h C P v i o l a t i o n //
///////////////////////////////////////////////////////
// C o n s t r u c t p d f
// -------------------------
// Additional parameters needed for B decay with CPV
RooRealVar CPeigen("CPeigen","CP eigen value",-1) ;
RooRealVar absLambda("absLambda","|lambda|",1,0,2) ;
RooRealVar argLambda("absLambda","|lambda|",0.7,-1,1) ;
RooRealVar effR("effR","B0/B0bar reco efficiency ratio",1) ;
// Construct Bdecay with CP violation
RooBCPEffDecay bcp("bcp","bcp", dt, tagFlav, tau, dm, w, CPeigen, absLambda, argLambda, effR, dw, tm, RooBCPEffDecay::DoubleSided) ;
// P l o t s c e n a r i o 1 - s i n ( 2 b ) = 0 . 7 , | l | = 1
// ---------------------------------------------------------------------------
// Generate some data
std::unique_ptr<RooDataSet> data2{bcp.generate(RooArgSet(dt,tagFlav),10000)};
// Plot B0 and B0bar tagged data separately
RooPlot* frame4 = dt.frame(Title("B decay distribution with CPV(|l|=1,Im(l)=0.7) (B0/B0bar)")) ;
data2->plotOn(frame4,Cut("tagFlav==tagFlav::B0")) ;
bcp.plotOn(frame4,Slice(tagFlav,"B0")) ;
data2->plotOn(frame4,Cut("tagFlav==tagFlav::B0bar"),MarkerColor(kCyan)) ;
bcp.plotOn(frame4,Slice(tagFlav,"B0bar"),LineColor(kCyan),Name("alt")) ;
// P l o t s c e n a r i o 2 - s i n ( 2 b ) = 0 . 7 , | l | = 0 . 7
// -------------------------------------------------------------------------------
absLambda=0.7 ;
// Generate some data
std::unique_ptr<RooDataSet> data3{bcp.generate(RooArgSet(dt,tagFlav),10000)};
// Plot B0 and B0bar tagged data separately (sin2b = 0.7 plus direct CPV |l|=0.5)
RooPlot* frame5 = dt.frame(Title("B decay distribution with CPV(|l|=0.7,Im(l)=0.7) (B0/B0bar)")) ;
data3->plotOn(frame5,Cut("tagFlav==tagFlav::B0")) ;
bcp.plotOn(frame5,Slice(tagFlav,"B0")) ;
data3->plotOn(frame5,Cut("tagFlav==tagFlav::B0bar"),MarkerColor(kCyan)) ;
bcp.plotOn(frame5,Slice(tagFlav,"B0bar"),LineColor(kCyan),Name("alt")) ;
//////////////////////////////////////////////////////////////////////////////////
// G e n e r i c B d e c a y w i t h u s e r c o e f f i c i e n t s //
//////////////////////////////////////////////////////////////////////////////////
// C o n s t r u c t p d f
// -------------------------
// Model parameters
RooRealVar DGbG("DGbG","DGamma/GammaAvg",0.5,-1,1);
RooRealVar Adir("Adir","-[1-abs(l)**2]/[1+abs(l)**2]",0);
RooRealVar Amix("Amix","2Im(l)/[1+abs(l)**2]",0.7);
RooRealVar Adel("Adel","2Re(l)/[1+abs(l)**2]",0.7);
// Derived input parameters for pdf
RooFormulaVar DG("DG","Delta Gamma","@1/@0",RooArgList(tau,DGbG));
// Construct coefficient functions for sin,cos,sinh modulations of decay distribution
RooFormulaVar fsin("fsin","fsin","@0*@1*(1-2*@2)",RooArgList(Amix,tagFlav,w));
RooFormulaVar fcos("fcos","fcos","@0*@1*(1-2*@2)",RooArgList(Adir,tagFlav,w));
RooFormulaVar fsinh("fsinh","fsinh","@0",RooArgList(Adel));
// Construct generic B decay pdf using above user coefficients
RooBDecay bcpg("bcpg","bcpg",dt,tau,DG,RooConst(1),fsinh,fcos,fsin,dm,tm,RooBDecay::DoubleSided);
// P l o t - I m ( l ) = 0 . 7 , R e ( l ) = 0 . 7 | l | = 1, d G / G = 0 . 5
// -------------------------------------------------------------------------------------
// Generate some data
std::unique_ptr<RooDataSet> data4{bcpg.generate(RooArgSet(dt,tagFlav),10000)};
// Plot B0 and B0bar tagged data separately
RooPlot* frame6 = dt.frame(Title("B decay distribution with CPV(Im(l)=0.7,Re(l)=0.7,|l|=1,dG/G=0.5) (B0/B0bar)")) ;
data4->plotOn(frame6,Cut("tagFlav==tagFlav::B0")) ;
bcpg.plotOn(frame6,Slice(tagFlav,"B0")) ;
data4->plotOn(frame6,Cut("tagFlav==tagFlav::B0bar"),MarkerColor(kCyan)) ;
bcpg.plotOn(frame6,Slice(tagFlav,"B0bar"),LineColor(kCyan),Name("alt")) ;
regPlot(frame1,"rf708_plot1") ;
regPlot(frame2,"rf708_plot2") ;
regPlot(frame3,"rf708_plot3") ;
regPlot(frame4,"rf708_plot4") ;
regPlot(frame5,"rf708_plot5") ;
regPlot(frame6,"rf708_plot6") ;
return true ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'VALIDATION AND MC STUDIES' RooFit tutorial macro #801
//
// A Toy Monte Carlo study that perform cycles of
// event generation and fittting
//
//
/////////////////////////////////////////////////////////////////////////
class TestBasic801 : public RooUnitTest
{
public:
TestBasic801(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("Automated MC studies",refFile,writeRef,verbose) {} ;
bool testCode() {
// C r e a t e m o d e l
// -----------------------
// Declare observable x
RooRealVar x("x","x",0,10) ;
x.setBins(40) ;
// Create two Gaussian PDFs g1(x,mean1,sigma) anf g2(x,mean2,sigma) and their parameters
RooRealVar mean("mean","mean of gaussians",5,0,10) ;
RooRealVar sigma1("sigma1","width of gaussians",0.5) ;
RooRealVar sigma2("sigma2","width of gaussians",1) ;
RooGaussian sig1("sig1","Signal component 1",x,mean,sigma1) ;
RooGaussian sig2("sig2","Signal component 2",x,mean,sigma2) ;
// Build Chebychev polynomial p.d.f.
RooRealVar a0("a0","a0",0.5,0.,1.) ;
RooRealVar a1("a1","a1",-0.2,-1.,1.) ;
RooChebychev bkg("bkg","Background",x,RooArgSet(a0,a1)) ;
// Sum the signal components into a composite signal p.d.f.
RooRealVar sig1frac("sig1frac","fraction of component 1 in signal",0.8,0.,1.) ;
RooAddPdf sig("sig","Signal",RooArgList(sig1,sig2),sig1frac) ;
// Sum the composite signal and background
RooRealVar nbkg("nbkg","number of background events,",150,0,1000) ;
RooRealVar nsig("nsig","number of signal events",150,0,1000) ;
RooAddPdf model("model","g1+g2+a",RooArgList(bkg,sig),RooArgList(nbkg,nsig)) ;
// C r e a t e m a n a g e r
// ---------------------------
// Instantiate RooMCStudy manager on model with x as observable and given choice of fit options
//
// The Silence() option kills all messages below the PROGRESS level, leaving only a single message
// per sample executed, and any error message that occur during fitting
//
// The Extended() option has two effects:
// 1) The extended ML term is included in the likelihood and
// 2) A poisson fluctuation is introduced on the number of generated events
//
// The FitOptions() given here are passed to the fitting stage of each toy experiment.
// If Save() is specified, the fit result of each experiment is saved by the manager
//
// A Binned() option is added in this example to bin the data between generation and fitting
// to speed up the study at the expemse of some precision
RooMCStudy mcstudy{model,x,Binned(true),Silence(),Extended(),
FitOptions(Save(true),PrintEvalErrors(0))};
// G e n e r a t e a n d f i t e v e n t s
// ---------------------------------------------
// Generate and fit 100 samples of Poisson(nExpected) events
mcstudy.generateAndFit(100) ;
// E x p l o r e r e s u l t s o f s t u d y
// ------------------------------------------------
// Make plots of the distributions of mean, the error on mean and the pull of mean
RooPlot* frame1 = mcstudy.plotParam(mean,Bins(40)) ;
RooPlot* frame2 = mcstudy.plotError(mean,Bins(40)) ;
RooPlot* frame3 = mcstudy.plotPull(mean,Bins(40),FitGauss(true)) ;
// Plot distribution of minimized likelihood
RooPlot* frame4 = mcstudy.plotNLL(Bins(40)) ;
regPlot(frame1,"rf801_plot1") ;
regPlot(frame2,"rf801_plot2") ;
regPlot(frame3,"rf801_plot3") ;
regPlot(frame4,"rf801_plot4") ;
return true ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'VALIDATION AND MC STUDIES' RooFit tutorial macro #802
//
// RooMCStudy: using separate fit and generator models, using the chi^2 calculator model
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic802 : public RooUnitTest
{
public:
TestBasic802(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("MC Study with chi^2 calculator",refFile,writeRef,verbose) {} ;
bool testCode() {
// C r e a t e m o d e l
// -----------------------
// Observables, parameters
RooRealVar x("x","x",-10,10) ;
x.setBins(10) ;
RooRealVar mean("mean","mean of gaussian",0) ;
RooRealVar sigma("sigma","width of gaussian",5,1,10) ;
// Create Gaussian pdf
RooGaussian gauss("gauss","gaussian PDF",x,mean,sigma) ;
// C r e a t e m a n a g e r w i t h c h i ^ 2 a d d - o n m o d u l e
// ----------------------------------------------------------------------------
// Create study manager for binned likelihood fits of a Gaussian pdf in 10 bins
RooMCStudy mcs{gauss,x,Silence(),Binned()};
// Add chi^2 calculator module to mcs
RooChi2MCSModule chi2mod ;
mcs.addModule(chi2mod) ;
// Generate 200 samples of 1000 events
mcs.generateAndFit(200,1000) ;
// Fill histograms with distributions chi2 and prob(chi2,ndf) that
// are calculated by RooChiMCSModule
RooRealVar* chi2 = (RooRealVar*) mcs.fitParDataSet().get()->find("chi2") ;
RooRealVar* prob = (RooRealVar*) mcs.fitParDataSet().get()->find("prob") ;
TH1* h_chi2 = new TH1F("h_chi2","",40,0,20) ;
TH1* h_prob = new TH1F("h_prob","",40,0,1) ;
mcs.fitParDataSet().fillHistogram(h_chi2,*chi2) ;
mcs.fitParDataSet().fillHistogram(h_prob,*prob) ;
// C r e a t e m a n a g e r w i t h s e p a r a t e f i t m o d e l
// ----------------------------------------------------------------------------
// Create alternate pdf with shifted mean
RooRealVar mean2("mean2","mean of gaussian 2",0.5) ;
RooGaussian gauss2("gauss2","gaussian PDF2",x,mean2,sigma) ;
// Create study manager with separate generation and fit model. This configuration
// is set up to generate bad fits as the fit and generator model have different means
// and the mean parameter is not floating in the fit
RooMCStudy mcs2{gauss2,x,FitModel(gauss),Silence(),Binned()};
// Add chi^2 calculator module to mcs
RooChi2MCSModule chi2mod2 ;
mcs2.addModule(chi2mod2) ;
// Generate 200 samples of 1000 events
mcs2.generateAndFit(200,1000) ;
// Fill histograms with distributions chi2 and prob(chi2,ndf) that
// are calculated by RooChiMCSModule
TH1* h2_chi2 = new TH1F("h2_chi2","",40,0,20) ;
TH1* h2_prob = new TH1F("h2_prob","",40,0,1) ;
mcs2.fitParDataSet().fillHistogram(h2_chi2,*chi2) ;
mcs2.fitParDataSet().fillHistogram(h2_prob,*prob) ;
h_chi2->SetLineColor(kRed) ;
h_prob->SetLineColor(kRed) ;
regTH(h_chi2,"rf802_hist_chi2") ;
regTH(h2_chi2,"rf802_hist2_chi2") ;
regTH(h_prob,"rf802_hist_prob") ;
regTH(h2_prob,"rf802_hist2_prob") ;
return true ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'VALIDATION AND MC STUDIES' RooFit tutorial macro #803
//
// RooMCStudy: Using the randomizer and profile likelihood add-on models
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic803 : public RooUnitTest
{
public:
TestBasic803(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("MC Study with param rand. and Z calc",refFile,writeRef,verbose) {} ;
bool testCode() {
// C r e a t e m o d e l
// -----------------------
// Simulation of signal and background of top quark decaying into
// 3 jets with background
// Observable
RooRealVar mjjj("mjjj","m(3jet) (GeV)",100,85.,350.) ;
// Signal component (Gaussian)
RooRealVar mtop("mtop","m(top)",162) ;
RooRealVar wtop("wtop","m(top) resolution",15.2) ;
RooGaussian sig("sig","top signal",mjjj,mtop,wtop) ;
// Background component (Chebychev)
RooRealVar c0("c0","Chebychev coefficient 0",-0.846,-1.,1.) ;
RooRealVar c1("c1","Chebychev coefficient 1", 0.112, 0.,1.) ;
RooRealVar c2("c2","Chebychev coefficient 2", 0.076, 0.,1.) ;
RooChebychev bkg("bkg","combinatorial background",mjjj,RooArgList(c0,c1,c2)) ;
// Composite model
RooRealVar nsig("nsig","number of signal events",53,0,1e3) ;
RooRealVar nbkg("nbkg","number of background events",103,0,5e3) ;
RooAddPdf model("model","model",RooArgList(sig,bkg),RooArgList(nsig,nbkg)) ;
// C r e a t e m a n a g e r
// ---------------------------
// Configure manager to perform binned extended likelihood fits (Binned(),Extended()) on data generated
// with a Poisson fluctuation on Nobs (Extended())
auto mcs = std::make_unique<RooMCStudy>(model,mjjj,Binned(),Silence(),Extended(true),
FitOptions(Extended(true),PrintEvalErrors(-1))) ;
// C u s t o m i z e m a n a g e r
// ---------------------------------
// Add module that randomizes the summed value of nsig+nbkg
// sampling from a uniform distribution between 0 and 1000
//
// In general one can randomize a single parameter, or a
// sum of N parameters, using either a uniform or a Gaussian
// distribution. Multiple randomization can be executed
// by a single randomizer module
RooRandomizeParamMCSModule randModule ;
randModule.sampleSumUniform(RooArgSet(nsig,nbkg),50,500) ;
mcs->addModule(randModule) ;
// Add profile likelihood calculation of significance. Redo each
// fit while keeping parameter nsig fixed to zero. For each toy,
// the difference in -log(L) of both fits is stored, as well
// a simple significance interpretation of the delta(-logL)
// using Dnll = 0.5 sigma^2
RooDLLSignificanceMCSModule sigModule(nsig,0) ;
mcs->addModule(sigModule) ;
// R u n m a n a g e r , m a k e p l o t s
// ---------------------------------------------
mcs->generateAndFit(50) ;
// Make some plots
RooRealVar* ngen = (RooRealVar*) mcs->fitParDataSet().get()->find("ngen") ;
RooRealVar* dll = (RooRealVar*) mcs->fitParDataSet().get()->find("dll_nullhypo_nsig") ;
RooRealVar* z = (RooRealVar*) mcs->fitParDataSet().get()->find("significance_nullhypo_nsig") ;
RooRealVar* nsigerr = (RooRealVar*) mcs->fitParDataSet().get()->find("nsigerr") ;
TH1* dll_vs_ngen = new TH2F("h_dll_vs_ngen" ,"",40,0,500,40,0,50) ;
TH1* z_vs_ngen = new TH2F("h_z_vs_ngen" ,"",40,0,500,40,0,10) ;
TH1* errnsig_vs_ngen = new TH2F("h_nsigerr_vs_ngen","",40,0,500,40,0,30) ;
TH1* errnsig_vs_nsig = new TH2F("h_nsigerr_vs_nsig","",40,0,200,40,0,30) ;
mcs->fitParDataSet().fillHistogram(dll_vs_ngen,RooArgList(*ngen,*dll)) ;
mcs->fitParDataSet().fillHistogram(z_vs_ngen,RooArgList(*ngen,*z)) ;
mcs->fitParDataSet().fillHistogram(errnsig_vs_ngen,RooArgList(*ngen,*nsigerr)) ;
mcs->fitParDataSet().fillHistogram(errnsig_vs_nsig,RooArgList(nsig,*nsigerr)) ;
regTH(dll_vs_ngen,"rf803_dll_vs_ngen") ;
regTH(z_vs_ngen,"rf803_z_vs_ngen") ;
regTH(errnsig_vs_ngen,"rf803_errnsig_vs_ngen") ;
regTH(errnsig_vs_nsig,"rf803_errnsig_vs_nsig") ;
return true ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'VALIDATION AND MC STUDIES' RooFit tutorial macro #804
//
// Using RooMCStudy on models with constrains
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
class TestBasic804 : public RooUnitTest
{
public:
TestBasic804(TFile* refFile, bool writeRef, int verbose) : RooUnitTest("MC Studies with aux. obs. constraints",refFile,writeRef,verbose) {} ;
double htol() { return 0.1 ; } // numerically very difficult test
bool testCode() {
// C r e a t e m o d e l w i t h p a r a m e t e r c o n s t r a i n t
// ---------------------------------------------------------------------------
// Observable
RooRealVar x("x","x",-10,10) ;
// Signal component
RooRealVar m("m","m",0,-10,10) ;
RooRealVar s("s","s",2,0.1,10) ;
RooGaussian g("g","g",x,m,s) ;
// Background component
RooPolynomial p("p","p",x) ;
// Composite model
RooRealVar f("f","f",0.4,0.,1.) ;
RooAddPdf sum("sum","sum",RooArgSet(g,p),f) ;
// Construct constraint on parameter f
RooGaussian fconstraint("fconstraint","fconstraint",f,0.7,0.1) ;
// Multiply constraint with p.d.f
RooProdPdf sumc("sumc","sum with constraint",RooArgSet(sum,fconstraint)) ;
// S e t u p t o y s t u d y w i t h m o d e l
// ---------------------------------------------------
// Perform toy study with internal constraint on f
//RooMCStudy mcs(sumc,x,Constrain(f),Silence(),Binned(),FitOptions(PrintLevel(-1))) ;
RooMCStudy mcs(sumc,x,Constrain(f),Binned());
// Run 50 toys of 2000 events.
// Before each toy is generated, a value for the f is sampled from the constraint pdf and
// that value is used for the generation of that toy.
mcs.generateAndFit(50,2000) ;
// Make plot of distribution of generated value of f parameter
RooRealVar* f_gen = (RooRealVar*) mcs.fitParDataSet().get()->find("f_gen") ;
TH1* h_f_gen = new TH1F("h_f_gen","",40,0,1) ;
mcs.fitParDataSet().fillHistogram(h_f_gen,*f_gen) ;
// Make plot of distribution of fitted value of f parameter
RooPlot* frame1 = mcs.plotParam(f,Bins(40),Range(0.4,1)) ;
frame1->SetTitle("Distribution of fitted f values") ;
// Make plot of pull distribution on f
RooPlot* frame2 = mcs.plotPull(f,Bins(40),Range(-3,3)) ;
frame1->SetTitle("Distribution of f pull values") ;
regTH(h_f_gen,"rf804_h_f_gen") ;
regPlot(frame1,"rf804_plot1") ;
regPlot(frame2,"rf804_plot2") ;
return true ;
}
} ;