-
Lorenzo Moneta authored
Fix a bug in the definition of the range for a coefficient a1 in the Roofit tutorial and stressRooFit. The parameter was considered constant instead should have been varying. Need to update stressRooFit file since the fit now are different (parameter a1 is varying in the fit instead of being constant )
Lorenzo Moneta authoredFix a bug in the definition of the range for a coefficient a1 in the Roofit tutorial and stressRooFit. The parameter was considered constant instead should have been varying. Need to update stressRooFit file since the fit now are different (parameter a1 is varying in the fit instead of being constant )
stressRooFit_tests.cxx 198.96 KiB
//////////////////////////////////////////////////////////////////////////
//
// '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
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooUnitTest.h"
using namespace RooFit ;
// Elementary operations on a gaussian PDF
class TestBasic101 : public RooUnitTest
{
public:
TestBasic101(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Fitting,plotting & event generation of basic p.d.f",refFile,writeRef,verbose) {} ;
Bool_t 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
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") ;
delete data ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'BASIC FUNCTIONALITY' RooFit tutorial macro #102
//
// Importing data from ROOT TTrees and THx histograms
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TTree.h"
#include "TH1D.h"
#include "TRandom.h"
using namespace RooFit ;
class TestBasic102 : public RooUnitTest
{
public:
TestBasic102(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Data import methods",refFile,writeRef,verbose) {} ;
TH1* makeTH1()
{
// Create ROOT TH1 filled with a Gaussian distribution
TH1D* hh = new TH1D("hh","hh",25,-10,10) ;
for (int i=0 ; i<100 ; i++) {
hh->Fill(gRandom->Gaus(0,3)) ;
}
return hh ;
}
TTree* makeTTree()
{
// Create ROOT TTree filled with a Gaussian distribution in x and a uniform distribution in y
TTree* tree = new TTree("tree","tree") ;
Double_t* px = new Double_t ;
Double_t* py = new Double_t ;
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_t 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
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
// -----------------------------------------------------------
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") ;
delete hh ;
delete tree ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'BASIC FUNCTIONALITY' RooFit tutorial macro #103
//
// Interpreted functions and p.d.f.s
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"#include "TCanvas.h"
#include "RooConstVar.h"
#include "RooPlot.h"
#include "RooFitResult.h"
#include "RooGenericPdf.h"
using namespace RooFit ;
class TestBasic103 : public RooUnitTest
{
public:
TestBasic103(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Interpreted expression p.d.f.",refFile,writeRef,verbose) {} ;
Bool_t 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
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,RooConst(10),RooConst(3)) ;
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
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(r,"rf103_fit1") ;
delete data ;
delete data2 ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'BASIC FUNCTIONALITY' RooFit tutorial macro #105
//
// Demonstration of binding ROOT Math functions as RooFit functions
// and pdfs
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TMath.h"
#include "TF1.h"
#include "Math/DistFunc.h"
#include "RooCFunction1Binding.h"
#include "RooCFunction3Binding.h"
#include "RooTFnBinding.h"
using namespace RooFit ;
class TestBasic105 : public RooUnitTest
{
public:
TestBasic105(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("C++ function binding operator p.d.f",refFile,writeRef,verbose) {} ;
Bool_t 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) ;
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
RooRealVar x2("x2","x2",0,0.999) ;
RooRealVar a("a","a",5,0,10) ;
RooRealVar b("b","b",2,0,10) ;
RooAbsPdf* beta = bindPdf("beta",ROOT::Math::beta_pdf,x2,a,b) ;
// Generate some events and fit
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
TF1 *fa1 = new 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
RooAbsReal* rfa1 = bindFunction(fa1,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") ;
delete erf ;
delete beta ;
delete fa1 ;
delete rfa1 ;
delete data ;
return kTRUE ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'BASIC FUNCTIONALITY' RooFit tutorial macro #108
//
// Plotting unbinned data with alternate and variable binnings
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"#include "RooGaussModel.h"
#include "RooDecay.h"
#include "RooBMixDecay.h"
#include "RooCategory.h"
#include "RooBinning.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;
class TestBasic108 : public RooUnitTest
{
public:
TestBasic108(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Non-standard binning in counting and asymmetry plots",refFile,writeRef,verbose) {} ;
Bool_t 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
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") ;
delete data ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// '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
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooHist.h"
using namespace RooFit ;
class TestBasic109 : public RooUnitTest{
public:
TestBasic109(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Calculation of chi^2 and residuals in plots",refFile,writeRef,verbose) {} ;
Bool_t 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,RooConst(0),sigma) ;
// Generate a sample of 1000 events with sigma=3
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
RooHist* hresid = frame1->residHist() ;
// Construct a histogram with the pulls of the data w.r.t the curve
RooHist* hpull = frame1->pullHist() ;
// 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 data ;
//delete hresid ;
//delete hpull ;
return kTRUE ;
}
} ;
///////////////////////////////////////////////////////////////////////////
// '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
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooGaussian.h"
#include "RooAbsReal.h"
#include "RooPlot.h"
#include "TCanvas.h"
using namespace RooFit ;
class TestBasic110 : public RooUnitTest
{
public:
TestBasic110(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Normalization of p.d.f.s in 1D",refFile,writeRef,verbose) {} ;
Bool_t 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,RooConst(-2),RooConst(3)) ;
// 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]
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"
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'
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") ;
delete igx ;
delete igx_sig ;
delete gx_cdf ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'BASIC FUNCTIONALITY' RooFit tutorial macro #111
//
// Configuration and customization of how numeric (partial) integrals
// are executed
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooNumIntConfig.h"
#include "RooLandau.h"
#include "RooArgSet.h"
#include <iomanip>
using namespace RooFit ;
class TestBasic111 : public RooUnitTest
{
public:
TestBasic111(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Numeric integration configuration",refFile,writeRef,verbose) {} ;
Bool_t 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
// ------------------------------------------------------------------
// Construct p.d.f without support for analytical integrator for demonstration purposes
RooRealVar x("x","x",-10,10) ;
RooLandau landau("landau","landau",x,RooConst(0),RooConst(0.1)) ;
// Calculate integral over landau with default choice of numeric integrator
RooAbsReal* intLandau = landau.createIntegral(x) ;
Double_t 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()) ;
customConfig.method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D") ;
// Calculate integral over landau with custom integral specification
RooAbsReal* intLandau2 = landau.createIntegral(x,NumIntConfig(customConfig)) ;
Double_t 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
RooAbsReal* intLandau3 = landau.createIntegral(x) ;
Double_t val3 = intLandau3->getVal() ;
regValue(val3,"rf111_val3") ;
delete intLandau ;
delete intLandau2 ;
delete intLandau3 ;
return kTRUE ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// '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
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;
class TestBasic201 : public RooUnitTest
{
public:
TestBasic201(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Addition operator p.d.f.",refFile,writeRef,verbose) {} ;
Bool_t 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
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),kTRUE) ;
// 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") ;
delete data ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'ADDITION AND CONVOLUTION' RooFit tutorial macro #202
//
// Setting up an extended maximum likelihood fit
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooExtendPdf.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;
class TestBasic202 : public RooUnitTest
{
public:
TestBasic202(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Extended ML fits to addition operator p.d.f.s",refFile,writeRef,verbose) {} ;
Bool_t 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
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") ;
delete data ;
return kTRUE ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'ADDITION AND CONVOLUTION' RooFit tutorial macro #203
//
// Fitting and plotting in sub ranges
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolynomial.h"
#include "RooAddPdf.h"
#include "RooFitResult.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;
class TestBasic203 : public RooUnitTest
{
public:
TestBasic203(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Basic fitting and plotting in ranges",refFile,writeRef,verbose) {} ;
Bool_t 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,RooConst(1)) ;
// 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
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
RooFitResult* r_full = model.fitTo(*modelData,Save(kTRUE)) ;
// 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
RooFitResult* r_sig = model.fitTo(*modelData,Save(kTRUE),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("Full"),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(r_full,"rf203_r_full") ;
regResult(r_sig,"rf203_r_sig") ;
delete modelData ;
return kTRUE;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'ADDITION AND CONVOLUTION' RooFit tutorial macro #204
//
// Extended maximum likelihood fit with alternate range definition
// for observed number of events.
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooExtendPdf.h"
#include "RooFitResult.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;
class TestBasic204 : public RooUnitTest
{
public:
TestBasic204(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Extended ML fit in sub range",refFile,writeRef,verbose) {} ;
Bool_t 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
RooDataSet *data = model.generate(x,1000) ;
// Perform unbinned extended ML fit to data
RooFitResult* r = model.fitTo(*data,Extended(kTRUE),Save()) ;
regResult(r,"rf204_result") ;
delete data ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'ADDITION AND CONVOLUTION' RooFit tutorial macro #205
//
// Options for plotting components of composite p.d.f.s.
//
//
//
// 07/2008 - Wouter Verkerke
///////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooAddPdf.h"
#include "RooChebychev.h"
#include "RooExponential.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;
class TestBasic205 : public RooUnitTest
{
public:
TestBasic205(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Component plotting variations",refFile,writeRef,verbose) {} ;
Bool_t 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) ;
// Sum the background components into a composite background p.d.f.
RooRealVar bkg1frac("sig1frac","fraction of component 1 in background",0.2,0.,1.) ;
RooAddPdf bkg("bkg","Signal",RooArgList(bkg1,bkg2),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) ;
// 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
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") ;
delete data ;
return kTRUE ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// '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
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooLandau.h"
#include "RooFFTConvPdf.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TPluginManager.h"#include "TROOT.h"
using namespace RooFit ;
class TestBasic208 : public RooUnitTest
{
public:
TestBasic208(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("FFT Convolution operator p.d.f.",refFile,writeRef,verbose) {} ;
Bool_t 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 kFALSE;
} else {
return kTRUE ;
}
}
}
return kFALSE ;
}
Double_t ctol() { return 1e-2 ; } // Account for difficult shape of Landau distribution
Bool_t 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
RooDataSet* data = lxg.generate(t,10000) ;
// Fit gxlx to data
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") ;
delete data ;
return kTRUE ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// '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
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussModel.h"
#include "RooAddModel.h"
#include "RooTruthModel.h"
#include "RooDecay.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;
class TestBasic209 : public RooUnitTest
{
public:
TestBasic209(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Analytical convolution operator",refFile,writeRef,verbose) {} ;
Bool_t 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 kTRUE ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// '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
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolyVar.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;
class TestBasic301 : public RooUnitTest
{public:
TestBasic301(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Composition extension of basic p.d.f",refFile,writeRef,verbose) {} ;
Bool_t 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
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") ;
delete data ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #302
//
// Utility functions classes available for use in tailoring
// of composite (multidimensional) pdfs
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"#include "TCanvas.h"
#include "RooPlot.h"
#include "RooFormulaVar.h"
#include "RooAddition.h"
#include "RooProduct.h"
#include "RooPolyVar.h"
#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;
class TestBasic302 : public RooUnitTest
{
public:
TestBasic302(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Sum and product utility functions",refFile,writeRef,verbose) {} ;
Bool_t 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 kTRUE ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// '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
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolyVar.h"
#include "RooProdPdf.h"
#include "RooPlot.h"
#include "TRandom.h"
#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;
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_t tmpy = gRandom->Gaus(0,10) ;
Double_t 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_t writeRef, Int_t verbose) : RooUnitTest("Conditional use of F(x|y)",refFile,writeRef,verbose) {} ;
Bool_t 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
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
RooDataSet* expDataY= (RooDataSet*) expDataXY->reduce(y) ;
// Generate 10000 events in x obtained from _conditional_ model(x|y) with y values taken from experimental data
RooDataSet *data = model.generate(x,ProtoData(*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
RooAbsData* binnedDataY = expDataY->binnedClone() ;
model.plotOn(xframe,ProjWData(*binnedDataY),LineColor(kCyan),LineStyle(kDotted),Name("Alt1")) ;
// Show effect of projection with too coarse binning ((RooRealVar*)expDataY->get()->find("y"))->setBins(5) ;
RooAbsData* binnedDataY2 = expDataY->binnedClone() ;
model.plotOn(xframe,ProjWData(*binnedDataY2),LineColor(kRed),Name("Alt2")) ;
regPlot(xframe,"rf303_plot1") ;
delete binnedDataY ;
delete binnedDataY2 ;
delete expDataXY ;
delete expDataY ;
delete data ;
return kTRUE ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// '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
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooProdPdf.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;
class TestBasic304 : public RooUnitTest
{
public:
TestBasic304(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Product operator p.d.f. with uncorrelated terms",refFile,writeRef,verbose) {} ;
Bool_t 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
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") ;
delete data ;
return kTRUE ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// '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
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolyVar.h"
#include "RooProdPdf.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;
class TestBasic305 : public RooUnitTest
{
public:
TestBasic305(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Product operator p.d.f. with conditional term",refFile,writeRef,verbose) {} ;
Bool_t 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,RooConst(0),RooConst(3)) ;
// 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
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") ;
delete data ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #306
//
// Complete example with use of conditional p.d.f. with per-event errors//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooGaussModel.h"
#include "RooDecay.h"
#include "RooLandau.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH2D.h"
using namespace RooFit ;
class TestBasic306 : public RooUnitTest
{
public:
TestBasic306(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Conditional use of per-event error p.d.f. F(t|dt)",refFile,writeRef,verbose) {} ;
Bool_t 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,RooConst(1),RooConst(0.25)) ;
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.
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_t 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 kTRUE argument bins the data before projection to speed up the process)
decay_gm.plotOn(frame2,ProjWData(*expDataDterr,kTRUE)) ;
regTH(hh_decay,"rf306_model2d") ;
regPlot(frame,"rf306_plot1") ;
regPlot(frame2,"rf306_plot2") ;
delete expDataDterr ;
delete data ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #307
//
// Complete example with use of full p.d.f. with per-event errors
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooGaussModel.h"
#include "RooDecay.h"
#include "RooLandau.h"
#include "RooProdPdf.h"
#include "RooHistPdf.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;
class TestBasic307 : public RooUnitTest
{
public:
TestBasic307(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Full per-event error p.d.f. F(t|dt)G(dt)",refFile,writeRef,verbose) {} ;
Bool_t 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,RooConst(1),RooConst(0.25)) ;
RooDataSet* expDataDterr = pdfDtErr.generate(dterr,10000) ;
// Construct a histogram pdf to describe the shape of the dtErr distribution
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.
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") ;
delete expDataDterr ;
delete expHistDterr ;
delete data ;
return kTRUE ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// '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
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooGaussian.h"
#include "RooProdPdf.h"
#include "RooAbsReal.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;
class TestBasic308 : public RooUnitTest
{
public:
TestBasic308(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Normalization of p.d.f.s in 2D",refFile,writeRef,verbose) {} ;
Bool_t 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,RooConst(-2),RooConst(3)) ;
RooGaussian gy("gy","gy",y,RooConst(+2),RooConst(2)) ;
// 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]
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"
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'
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") ;
delete igxy_sig ;
delete igxy ;
delete gxy_cdf ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #309
//
// Projecting p.d.f and data slices in discrete observables
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"#include "RooGaussModel.h"
#include "RooDecay.h"
#include "RooBMixDecay.h"
#include "RooCategory.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;
class TestBasic310 : public RooUnitTest
{
public:
TestBasic310(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Data and p.d.f projection in category slice",refFile,writeRef,verbose) {} ;
Bool_t 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
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") ;
delete data ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #310
//
// Projecting p.d.f and data ranges in continuous observables
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooProdPdf.h"
#include "RooAddPdf.h"
#include "RooPolynomial.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;
class TestBasic311 : public RooUnitTest
{
public:
TestBasic311(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Data and p.d.f projection in sub range",refFile,writeRef,verbose) {} ;
Bool_t 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,RooConst(0),RooConst(1)) ;
RooGaussian gy("gy","gy",y,RooConst(0),RooConst(1)) ;
RooGaussian gz("gz","gz",z,RooConst(0),RooConst(1)) ;
RooProdPdf sig("sig","sig",RooArgSet(gx,gy,gz)) ;
// Create background pdf poly(x)*poly(y)*poly(z)
RooPolynomial px("px","px",x,RooArgSet(RooConst(-0.1),RooConst(0.004))) ;
RooPolynomial py("py","py",y,RooArgSet(RooConst(0.1),RooConst(-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) ;
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") ;
delete data ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #312
//
// Performing fits in multiple (disjoint) ranges in one or more dimensions
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooProdPdf.h"
#include "RooAddPdf.h"
#include "RooPolynomial.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooFitResult.h"
using namespace RooFit ;
class TestBasic312 : public RooUnitTest
{
public:
TestBasic312(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Fit in multiple rectangular ranges",refFile,writeRef,verbose) {} ;
Bool_t 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,RooConst(1)) ;
RooGaussian gy("gy","gy",y,my,RooConst(1)) ;
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
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)
RooFitResult* r_sb1 = model.fitTo(*modelData,Range("SB1"),Save()) ;
// Perform fit in SideBand2 region (RooAddPdf coefficients will be interpreted in full range)
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)
RooFitResult* r_sb12 = model.fitTo(*modelData,Range("SB1,SB2"),Save()) ;
regResult(r_sb1,"rf312_fit_sb1") ;
regResult(r_sb2,"rf312_fit_sb2") ;
regResult(r_sb12,"rf312_fit_sb12") ;
delete modelData ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #313
//
// Working with parameterized ranges to define non-rectangular regions
// for fitting and integration
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolynomial.h"
#include "RooProdPdf.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;
class TestBasic313 : public RooUnitTest
{
public:
TestBasic313(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Integration over non-rectangular regions",refFile,writeRef,verbose) {} ;
Bool_t 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
// ----------------------------------------------------------------------------------
// Create integral over normalized pdf model over x,y,z in "R" region
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") ;
delete intPdf ;
return kTRUE;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// '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
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooExponential.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooFitResult.h"
using namespace RooFit ;
class TestBasic314 : public RooUnitTest
{
public:
TestBasic314(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Fit with non-rectangular observable boundaries",refFile,writeRef,verbose) {} ;
Bool_t 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)
RooDataSet* dall = model.generate(t,10000) ;
// Generate a (fake) prototype dataset for acceptance limit values
RooDataSet* tmp = RooGaussian("gmin","gmin",tmin,RooConst(0),RooConst(0.5)).generate(tmin,5000) ;
// Generate dataset with t values that observe (t>tmin)
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
// -----------------------------------------------------------------------
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(r,"rf314_fit") ;
regPlot(frame,"rf314_plot1") ;
delete tmp ;
delete dacc ;
delete dall ;
return kTRUE;
}
} ;
//////////////////////////////////////////////////////////////////////////
//// 'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #315
//
// Marginizalization of multi-dimensional p.d.f.s through integration
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooProdPdf.h"
#include "RooPolyVar.h"
#include "TH1.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooNumIntConfig.h"
#include "RooConstVar.h"
using namespace RooFit ;
class TestBasic315 : public RooUnitTest
{
public:
TestBasic315(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("P.d.f. marginalization through integration",refFile,writeRef,verbose) {} ;
Bool_t 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,RooConst(0),RooConst(2)) ;
// 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
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
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") ;
delete data ;
delete modelx ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// '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
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolynomial.h"
#include "RooAddPdf.h"
#include "RooProdPdf.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;
class TestBasic316 : public RooUnitTest
{
public:
TestBasic316(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Likelihood ratio projection plot",refFile,writeRef,verbose) {} ;
Bool_t 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,RooConst(0),RooConst(1)) ; RooGaussian gy("gy","gy",y,RooConst(0),RooConst(1)) ;
RooGaussian gz("gz","gz",z,RooConst(0),RooConst(1)) ;
RooProdPdf sig("sig","sig",RooArgSet(gx,gy,gz)) ;
// Create background pdf poly(x)*poly(y)*poly(z)
RooPolynomial px("px","px",x,RooArgSet(RooConst(-0.1),RooConst(0.004))) ;
RooPolynomial py("py","py",y,RooArgSet(RooConst(0.1),RooConst(-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) ;
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
RooAbsPdf* sigyz = sig.createProjection(x) ;
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
RooDataSet* dataSel = (RooDataSet*) 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
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) ;
RooDataSet* mcprojDataSel = (RooDataSet*) 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(*mcprojDataSel)) ;
regPlot(frame,"rf316_plot1") ;
regPlot(frame2,"rf316_plot2") ;
delete data ;
delete dataSel ;
delete mcprojData ;
delete mcprojDataSel ;
delete sigyz ;
delete totyz ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'DATA AND CATEGORIES' RooFit tutorial macro #402
//
// Tools for manipulation of (un)binned datasets
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "RooCategory.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TFile.h"
using namespace RooFit ;
class TestBasic402 : public RooUnitTest
{
public:
TestBasic402(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Basic operations on datasets",refFile,writeRef,verbose) {} ;
Bool_t 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_t 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
RooDataSet* d1 = (RooDataSet*) d.reduce(RooArgSet(x,c)) ;
RooDataSet* d2 = (RooDataSet*) d.reduce(RooArgSet(y)) ;
RooDataSet* d3 = (RooDataSet*) d.reduce("y>5.17") ;
RooDataSet* d4 = (RooDataSet*) 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
d1->merge(d2) ;
// The append() function addes two datasets row-wise
d1->append(*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.
RooDataHist* dh2 = (RooDataHist*) dh.reduce(y,"x>0") ;
// Add dh2 to yframe and redraw
dh2->plotOn(yframe,LineColor(kRed),MarkerColor(kRed),Name("dh2")) ;
regPlot(yframe,"rf402_plot1") ;
delete d1 ;
delete d2 ;
delete d3 ;
delete d4 ;
delete dh2 ;
return kTRUE ; }
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'DATA AND CATEGORIES' RooFit tutorial macro #403
//
// Using weights in unbinned datasets
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "RooFormulaVar.h"
#include "RooGenericPdf.h"
#include "RooPolynomial.h"
#include "RooChi2Var.h"
#include "RooMinuit.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooFitResult.h"
using namespace RooFit ;
class TestBasic403 : public RooUnitTest
{
public:
TestBasic403(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Fits with weighted datasets",refFile,writeRef,verbose) {} ;
Bool_t 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
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,*data->get(),0,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
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
RooDataSet* data2 = genPdf.generate(x,1000) ;
// Sample a dataset with the same number of weights as data
RooDataSet* data3 = genPdf.generate(x,43000) ;
// Fit the 2nd order polynomial to both unweighted datasets and save the results for comparison
RooFitResult* r_ml_unw10 = p2.fitTo(*data2,Save()) ;
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
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
RooChi2Var chi2("chi2","chi2",p2,*binnedData) ;
RooMinuit m(chi2) ;
m.migrad() ;
m.hesse() ;
// Plot chi^2 fit result on frame as well
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(r_ml_unw10,"rf403_ml_unw10") ;
regResult(r_ml_unw43,"rf403_ml_unw43") ;
regResult(r_ml_wgt ,"rf403_ml_wgt") ;
regResult(r_chi2_wgt ,"rf403_ml_chi2") ;
regPlot(frame,"rf403_plot1") ;
delete binnedData ;
delete data ;
delete data2 ;
delete data3 ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'DATA AND CATEGORIES' RooFit tutorial macro #404
//
// Working with RooCategory objects to describe discrete variables
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooPolynomial.h"
#include "RooCategory.h"
#include "Roo1DTable.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;
class TestBasic404 : public RooUnitTest
{
public:
TestBasic404(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Categories basic functionality",refFile,writeRef,verbose) {} ;
Bool_t 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) ;
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_t nb0 = btable->get("B0") ;
regValue(nb0,"rf404_nb0") ;
// Retrieve fraction of events with "Lepton" tag
Double_t 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
RooDataSet* goodData = (RooDataSet*) 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") ;
delete goodData ;
delete data ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'DATA AND CATEGORIES' RooFit tutorial macro #405
//
// Demonstration of real-->discrete mapping functions
////
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooCategory.h"
#include "RooThresholdCategory.h"
#include "RooBinningCategory.h"
#include "Roo1DTable.h"
#include "RooArgusBG.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooRealConstant.h"
using namespace RooFit ;
class TestBasic405 : public RooUnitTest
{
public:
TestBasic405(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Real-to-category functions",refFile,writeRef,verbose) {} ;
Bool_t 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
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
RooDataSet* dataSel = (RooDataSet*) 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") ;
delete data ;
delete dataSel ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'DATA AND CATEGORIES' RooFit tutorial macro #406
//
// Demonstration of discrete-->discrete (invertable) functions
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"#include "RooPolynomial.h"
#include "RooCategory.h"
#include "RooMappedCategory.h"
#include "RooMultiCategory.h"
#include "RooSuperCategory.h"
#include "Roo1DTable.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;
class TestBasic406 : public RooUnitTest
{
public:
TestBasic406(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Category-to-category functions",refFile,writeRef,verbose) {} ;
Bool_t 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) ;
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") ;
delete data ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// '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
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooSimultaneous.h"
#include "RooCategory.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;
class TestBasic501 : public RooUnitTest
{
public:
TestBasic501(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Simultaneous p.d.f. operator",refFile,writeRef,verbose) {} ;
Bool_t 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
RooDataSet *data = model.generate(RooArgSet(x),100) ;
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),Import("control",*data_ctl)) ;
// 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") ;
delete data ;
delete data_ctl ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// '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
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooWorkspace.h"
#include "RooProdPdf.h"
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooGaussModel.h"
#include "RooAddModel.h"
#include "RooDecay.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooSimultaneous.h"
#include "RooCategory.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;
class TestBasic599 : public RooUnitTest
{
public:
TestBasic599(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Workspace and p.d.f. persistence",refFile,writeRef,verbose) {} ;
Bool_t 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 kFALSE ;
// 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 kTRUE ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #601
//
// Interactive minimization with MINUIT
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooProdPdf.h"
#include "RooAddPdf.h"
#include "RooMinuit.h"
#include "RooNLLVar.h"
#include "RooFitResult.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;
class TestBasic601 : public RooUnitTest
{
public:
TestBasic601(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Interactive Minuit",refFile,writeRef,verbose) {} ;
Bool_t 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
RooDataSet* data = model.generate(x,1000) ;
// Construct unbinned likelihood
RooNLLVar nll("nll","nll",model,*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
RooMinuit 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)
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(kTRUE) ;
// Rerun MIGRAD,HESSE
m.migrad() ;
m.hesse() ;
RooFitResult* r2 = m.save() ;
regResult(r,"rf601_r") ;
regResult(r2,"rf601_r2") ;
delete data ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #602
//
// Setting up a binning chi^2 fit
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooChi2Var.h"
#include "RooMinuit.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;
class TestBasic602 : public RooUnitTest
{
public:
TestBasic602(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Chi2 minimization",refFile,writeRef,verbose) {} ;
Bool_t 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
// -----------------------------------------
RooDataSet* d = model.generate(x,10000) ;
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
RooChi2Var chi2("chi2","chi2",model,*dh) ;
// Use RooMinuit interface to minimize chi^2
RooMinuit m(chi2) ;
m.migrad() ;
m.hesse() ;
RooFitResult* r = m.save() ;
regResult(r,"rf602_r") ;
delete d ;
delete dh ;
return kTRUE ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #604
//
// Fitting with constraints
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolynomial.h"
#include "RooAddPdf.h"
#include "RooProdPdf.h"
#include "RooFitResult.h"
#include "RooPlot.h"#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;
class TestBasic604 : public RooUnitTest
{
public:
TestBasic604(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Auxiliary observable constraints",refFile,writeRef,verbose) {} ;
Bool_t 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
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,RooConst(0.8),RooConst(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
RooFitResult* r1 = modelc.fitTo(*d,Save()) ;
// Fit modelc with constraint term on parameter f
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,RooConst(0.2),RooConst(0.1)) ;
// Fit with external constraint
RooFitResult* r3 = model.fitTo(*d,ExternalConstraints(fconstext),Save()) ;
regResult(r1,"rf604_r1") ;
regResult(r2,"rf604_r2") ; regResult(r3,"rf604_r3") ;
delete d ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #605
//
// Working with the profile likelihood estimator
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooAddPdf.h"
#include "RooNLLVar.h"
#include "RooProfileLL.h"
#include "RooMinuit.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;
class TestBasic605 : public RooUnitTest
{
public:
TestBasic605(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Profile Likelihood operator",refFile,writeRef,verbose) {} ;
Bool_t 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
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
RooNLLVar nll("nll","nll",model,*data) ;
// Minimize likelihood w.r.t all parameters before making plots
RooMinuit(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") ;
delete data ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #606
//
// Understanding and customizing error handling in likelihood evaluations
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooArgusBG.h"
#include "RooNLLVar.h"#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;
class TestBasic606 : public RooUnitTest
{
public:
TestBasic606(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("NLL error handling",refFile,writeRef,verbose) {} ;
Bool_t 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
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(kFALSE)) ;
m0.setError(0.1) ;
argus.fitTo(*data,PrintEvalErrors(0),EvalErrorWall(kFALSE),Warnings(kFALSE)) ;
// 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
RooNLLVar nll("nll","nll",argus,*data) ;
// 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
delete data ;
return kTRUE ;
}
} ;//////////////////////////////////////////////////////////////////////////
//
// 'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #607
//
// Demonstration of options of the RooFitResult class
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooAddPdf.h"
#include "RooChebychev.h"
#include "RooFitResult.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TFile.h"
#include "TStyle.h"
#include "TH2.h"
using namespace RooFit ;
class TestBasic607 : public RooUnitTest
{
public:
TestBasic607(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Fit Result functionality",refFile,writeRef,verbose) {} ;
Bool_t 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
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
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_t 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") ;
delete data ;
delete r ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// '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
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooPolyVar.h"
#include "RooChi2Var.h"
#include "RooMinuit.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TRandom.h"
using namespace RooFit ;
class TestBasic609 : public RooUnitTest
{
public:
TestBasic609(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Chi^2 fit to X-Y dataset",refFile,writeRef,verbose) {} ;
Bool_t 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) ;
RooPolyVar f("f","f",x,RooArgList(b,a,RooConst(1))) ;
// 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
f.chi2FitTo(dxy,YVar(y)) ;
// 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
f.chi2FitTo(dxy,YVar(y),Integrate(kTRUE)) ;
// Overlay alternate fit result
f.plotOn(frame,LineStyle(kDashed),LineColor(kRed),Name("alternate")) ;
regPlot(frame,"rf609_frame") ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// '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
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooFormulaVar.h"
#include "RooProdPdf.h"
#include "RooEfficiency.h"
#include "RooPolynomial.h"
#include "RooCategory.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;
// Elementary operations on a gaussian PDF
class TestBasic701 : public RooUnitTest
{
public:
TestBasic701(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Efficiency operator p.d.f. 1D",refFile,writeRef,verbose) {} ;
Bool_t 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
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") ;
delete data ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// '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)
//
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooCategory.h"
#include "RooEfficiency.h"
#include "RooPolynomial.h"
#include "RooProdPdf.h"
#include "RooFormulaVar.h"
#include "TCanvas.h"
#include "TH1.h"
#include "RooPlot.h"
using namespace RooFit ;
// Elementary operations on a gaussian PDF
class TestBasic702 : public RooUnitTest
{
public:
TestBasic702(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Efficiency operator p.d.f. 2D",refFile,writeRef,verbose) {} ;
Bool_t testCode() {
Bool_t flat=kFALSE ;
// 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
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") ;
delete data ;
return kTRUE;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// '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
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooExponential.h"
#include "RooEffProd.h"
#include "RooFormulaVar.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;
// Elementary operations on a gaussian PDF
class TestBasic703 : public RooUnitTest
{
public:
TestBasic703(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Efficiency product operator p.d.f",refFile,writeRef,verbose) {} ;
Bool_t 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.
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") ;
delete data ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'SPECIAL PDFS' RooFit tutorial macro #704
//
// Using a p.d.f defined by a sum of real-valued amplitude components
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooTruthModel.h"
#include "RooFormulaVar.h"
#include "RooRealSumPdf.h"
#include "RooPolyVar.h"
#include "RooProduct.h"
#include "TH1.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;
// Elementary operations on a gaussian PDF
class TestBasic704 : public RooUnitTest
{
public:
TestBasic704(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Amplitude sum operator p.d.f",refFile,writeRef,verbose) {} ;
Bool_t 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.);
// Use RooTruthModel to obtain compiled implementation of sinh/cosh modulated decay functions
RooRealVar tau("tau","#tau",1.5);
RooRealVar deltaGamma("deltaGamma","deltaGamma", 0.3);
RooTruthModel tm("tm","tm",t) ;
RooFormulaVar coshGBasis("coshGBasis","exp(-@0/ @1)*cosh(@0*@2/2)",RooArgList(t,tau,deltaGamma));
RooFormulaVar sinhGBasis("sinhGBasis","exp(-@0/ @1)*sinh(@0*@2/2)",RooArgList(t,tau,deltaGamma));
RooAbsReal* coshGConv = tm.convolution(&coshGBasis,&t);
RooAbsReal* sinhGConv = tm.convolution(&sinhGBasis,&t);
// Construct polynomial amplitudes in cos(a)
RooPolyVar poly1("poly1","poly1",cosa,RooArgList(RooConst(0.5),RooConst(0.2),RooConst(0.2)),0);
RooPolyVar poly2("poly2","poly2",cosa,RooArgList(RooConst(1),RooConst(-0.2),RooConst(3)),0);
// Construct 2D amplitude as uncorrelated product of amp(t)*amp(cosa)
RooProduct ampl1("ampl1","amplitude 1",RooArgSet(poly1,*coshGConv));
RooProduct ampl2("ampl2","amplitude 2",RooArgSet(poly2,*sinhGConv));
// 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
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") ;
delete data ;
delete coshGConv ;
delete sinhGConv ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'SPECIAL PDFS' RooFit tutorial macro #705
//
// Linear interpolation between p.d.f shapes using the 'Alex Read' algorithm
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolynomial.h"
#include "RooIntegralMorph.h"
#include "RooNLLVar.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TH1.h"
using namespace RooFit ;
// Elementary operations on a gaussian PDF
class TestBasic705 : public RooUnitTest
{
public:
Double_t ctol() { return 5e-2 ; } // very conservative, this is a numerically difficult test
TestBasic705(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Linear morph operator p.d.f.",refFile,writeRef,verbose) {} ;
Bool_t 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,RooConst(2)) ;
// Upper end point shape: a Polynomial
RooPolynomial g2("g2","g2",x,RooArgSet(RooConst(-0.03),RooConst(-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 ;
RooDataSet* data = lmorph.generate(x,1000) ;
// Fit pdf to toy data
lmorph.setCacheAlpha(kTRUE) ;
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
RooNLLVar nll("nll","nll",lmorph,*data) ;
nll.plotOn(frame3,ShiftToZero()) ;
lmorph.setCacheAlpha(kFALSE) ;
regPlot(frame1,"rf705_plot1") ;
regPlot(frame2,"rf705_plot2") ;
regPlot(frame3,"rf705_plot3") ;
regTH(hh,"rf705_hh") ;
delete data ;
return kTRUE;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'SPECIAL PDFS' RooFit tutorial macro #706
//
// Histogram based p.d.f.s and functions
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolynomial.h"
#include "RooHistPdf.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;
// Elementary operations on a gaussian PDF
class TestBasic706 : public RooUnitTest
{
public:
TestBasic706(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Histogram based p.d.f.s",refFile,writeRef,verbose) {} ;
Bool_t 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(RooConst(0.01),RooConst(-0.01),RooConst(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) ;
RooDataSet* data1 = p.generate(x,500) ;
// Create a binned dataset with 20 bins and 500 events
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) ;
RooDataSet* data2 = p.generate(x,100000) ;
// Create a binned dataset with 10 bins and 100K events
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") ;
delete data1 ;
delete hist1 ;
delete data2 ;
delete hist2 ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'SPECIAL PDFS' RooFit tutorial macro #707
//
// Using non-parametric (multi-dimensional) kernel estimation p.d.f.s
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolynomial.h"
#include "RooKeysPdf.h"
#include "RooNDKeysPdf.h"
#include "RooProdPdf.h"
#include "TCanvas.h"
#include "TH1.h"
#include "RooPlot.h"
using namespace RooFit ;
// Elementary operations on a gaussian PDF
class TestBasic707 : public RooUnitTest
{
public:
TestBasic707(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Kernel estimation p.d.f.s",refFile,writeRef,verbose) {} ; Bool_t 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(RooConst(0.01),RooConst(-0.01),RooConst(0.0004))) ;
// Sample 500 events from p
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(RooConst(0.01),RooConst(0.01),RooConst(-0.0004))) ;
RooProdPdf pxy("pxy","pxy",RooArgSet(p,py)) ;
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") ;
delete data1 ;
delete data2 ;
return kTRUE ;
}
} ;
//////////////////////////////////////////////////////////////////////////
//
// 'SPECIAL PDFS' RooFit tutorial macro #708
//
// Special decay pdf for B physics with mixing and/or CP violation
//
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooCategory.h"
#include "RooBMixDecay.h"
#include "RooBCPEffDecay.h"
#include "RooBDecay.h"
#include "RooFormulaVar.h"
#include "RooTruthModel.h"
#include "TCanvas.h"
#include "RooPlot.h"
using namespace RooFit ;
// Elementary operations on a gaussian PDF
class TestBasic708 : public RooUnitTest
{
public:
TestBasic708(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("B Physics p.d.f.s",refFile,writeRef,verbose) {} ;
Bool_t 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
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"),Slice(mixState,"mixed")) ;
data->plotOn(frame2,Cut("mixState==mixState::mixed&&tagFlav==tagFlav::B0bar"),MarkerColor(kCyan)) ;
bmix.plotOn(frame2,Slice(tagFlav,"B0bar"),Slice(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"),Slice(mixState,"unmixed")) ;
data->plotOn(frame3,Cut("mixState==mixState::unmixed&&tagFlav==tagFlav::B0bar"),MarkerColor(kCyan)) ;
bmix.plotOn(frame3,Slice(tagFlav,"B0bar"),Slice(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
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
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 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") ;
delete data ;
delete data2 ;
delete data3 ;
delete data4 ;
return kTRUE ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'VALIDATION AND MC STUDIES' RooFit tutorial macro #801
//
// A Toy Monte Carlo study that perform cycles of
// event generation and fittting
//
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooMCStudy.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH2.h"
#include "RooFitResult.h"
#include "TStyle.h"
#include "TDirectory.h"
using namespace RooFit ;
// Elementary operations on a gaussian PDF
class TestBasic801 : public RooUnitTest
{
public:
TestBasic801(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("Automated MC studies",refFile,writeRef,verbose) {} ;
Bool_t 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 = new RooMCStudy(model,x,Binned(kTRUE),Silence(),Extended(),
FitOptions(Save(kTRUE),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(kTRUE)) ;
// 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") ;
delete mcstudy ;
return kTRUE ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// '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
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooMCStudy.h"
#include "RooChi2MCSModule.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TDirectory.h"
using namespace RooFit ;
// Elementary operations on a gaussian PDF
class TestBasic802 : public RooUnitTest
{
public:
TestBasic802(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("MC Study with chi^2 calculator",refFile,writeRef,verbose) {} ;
Bool_t 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 = new RooMCStudy(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 = new RooMCStudy(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") ;
delete mcs ;
delete mcs2 ;
return kTRUE ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'VALIDATION AND MC STUDIES' RooFit tutorial macro #803
//
// RooMCStudy: Using the randomizer and profile likelihood add-on models
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooMCStudy.h"
#include "RooRandomizeParamMCSModule.h"
#include "RooDLLSignificanceMCSModule.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TDirectory.h"
using namespace RooFit ;
class TestBasic803 : public RooUnitTest
{
public:
TestBasic803(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("MC Study with param rand. and Z calc",refFile,writeRef,verbose) {} ;
Bool_t 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())
RooMCStudy* mcs = new RooMCStudy(model,mjjj,Binned(),Silence(),Extended(kTRUE),
FitOptions(Extended(kTRUE),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") ;
delete mcs ;
return kTRUE ;
}
} ;
/////////////////////////////////////////////////////////////////////////
//
// 'VALIDATION AND MC STUDIES' RooFit tutorial macro #804
//
// Using RooMCStudy on models with constrains
//
//
// 07/2008 - Wouter Verkerke
//
/////////////////////////////////////////////////////////////////////////
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooPolynomial.h"
#include "RooAddPdf.h"
#include "RooProdPdf.h"#include "RooMCStudy.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;
class TestBasic804 : public RooUnitTest
{
public:
TestBasic804(TFile* refFile, Bool_t writeRef, Int_t verbose) : RooUnitTest("MC Studies with aux. obs. constraints",refFile,writeRef,verbose) {} ;
Double_t htol() { return 0.1 ; } // numerically very difficult test
Bool_t 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,RooConst(0.7),RooConst(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 kTRUE ;
}
} ;