doc/hackathon/rf501__simultaneouspdf_8C.html

#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 "TAxis.h"

#include "RooPlot.h"

using namespace RooFit;


void rf501_simultaneouspdf()

{

   // C r e a t e   m o d e l   f o r   p h y s i c s   s a m p l e

   // -------------------------------------------------------------


   // Create observables

   RooRealVar x("x", "x", -8, 8);


   // Construct signal pdf

   RooRealVar mean("mean", "mean", 0, -8, 8);

   RooRealVar sigma("sigma", "sigma", 0.3, 0.1, 10);

   RooGaussian gx("gx", "gx", x, mean, sigma);


   // Construct background pdf

   RooRealVar a0("a0", "a0", -0.1, -1, 1);

   RooRealVar a1("a1", "a1", 0.004, -1, 1);

   RooChebychev px("px", "px", x, RooArgSet(a0, a1));


   // Construct composite pdf

   RooRealVar f("f", "f", 0.2, 0., 1.);

   RooAddPdf model("model", "model", RooArgList(gx, px), f);


   // C r e a t e   m o d e l   f o r   c o n t r o l   s a m p l e

   // --------------------------------------------------------------


   // Construct signal pdf.

   // NOTE that sigma is shared with the signal sample model

   RooRealVar mean_ctl("mean_ctl", "mean_ctl", -3, -8, 8);

   RooGaussian gx_ctl("gx_ctl", "gx_ctl", x, mean_ctl, sigma);


   // Construct the background pdf

   RooRealVar a0_ctl("a0_ctl", "a0_ctl", -0.1, -1, 1);

   RooRealVar a1_ctl("a1_ctl", "a1_ctl", 0.5, -0.1, 1);

   RooChebychev px_ctl("px_ctl", "px_ctl", x, RooArgSet(a0_ctl, a1_ctl));


   // Construct the composite model

   RooRealVar f_ctl("f_ctl", "f_ctl", 0.5, 0., 1.);

   RooAddPdf model_ctl("model_ctl", "model_ctl", RooArgList(gx_ctl, px_ctl), f_ctl);


   // G e n e r a t e   e v e n t s   f o r   b o t h   s a m p l e s

   // ---------------------------------------------------------------


   // Generate 1000 events in x and y from model

   std::unique_ptr<RooDataSet> data{model.generate({x}, 1000)};

   std::unique_ptr<RooDataSet> data_ctl{model_ctl.generate({x}, 2000)};


   // C r e a t e   i n d e x   c a t e g o r y   a n d   j o i n   s a m p l e s

   // ---------------------------------------------------------------------------


   // Define category to distinguish physics and control samples events

   RooCategory sample("sample", "sample");

   sample.defineType("physics");

   sample.defineType("control");


   // Construct combined dataset in (x,sample)

   RooDataSet combData("combData", "combined data", x, Index(sample),

                       Import({{"physics", data.get()}, {"control", data_ctl.get()}}));


   // C o n s t r u c t   a   s i m u l t a n e o u s   p d f   i n   ( x , s a m p l e )

   // -----------------------------------------------------------------------------------


   // Construct a simultaneous pdf using category sample as index:

   // associate model with the physics state and model_ctl with the control state

   RooSimultaneous simPdf("simPdf", "simultaneous pdf", {{"physics", &model}, {"control", &model_ctl}}, sample);


   // 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

   std::unique_ptr<RooFitResult> fitResult{simPdf.fitTo(combData, Save(), PrintLevel(-1))};

   fitResult->Print();


   // 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(Title("Physics sample"));


   // Plot all data tagged as physics sample

   std::unique_ptr<RooAbsData> slicedData1{combData.reduce(Cut("sample==sample::physics"))};

   slicedData1->plotOn(frame1);


   // Plot "physics" slice of simultaneous pdf.

   simPdf.getPdf("physics")->plotOn(frame1);

   simPdf.getPdf("physics")->plotOn(frame1, Components("px"), LineStyle(kDashed));


   // The same plot for the control sample slice. We do this with a different

   // approach, using the data slice for the projection. This approach is more

   // general, because you can plot sums of slices by using logical or in the

   // Cut() command.

   // 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.

   // In other words: Since the PDF doesn't know the number of events in the different

   // category states, it doesn't know how much of each component it has to project out.

   // This information is read from the data.

   RooPlot *frame2 = x.frame(Bins(30), Title("Control sample"));

   std::unique_ptr<RooAbsData> slicedData2{combData.reduce(Cut("sample==sample::control"))};

   slicedData2->plotOn(frame2);

   simPdf.plotOn(frame2, ProjWData(sample, *slicedData2));

   simPdf.plotOn(frame2, Components("px_ctl"), ProjWData(sample, *slicedData2), LineStyle(kDashed));


   // The same plot for all the phase space. Here, we can just use the original

   // combined dataset.

   RooPlot *frame3 = x.frame(Title("Both samples"));

   combData.plotOn(frame3);

   simPdf.plotOn(frame3, ProjWData(sample, combData));

   simPdf.plotOn(frame3, Components("px,px_ctl"), ProjWData(sample, combData),

                 LineStyle(kDashed));


   TCanvas *c = new TCanvas("rf501_simultaneouspdf", "rf403_simultaneouspdf", 1200, 400);

   c->Divide(3);

   auto draw = [&](int i, RooPlot & frame) {

      c->cd(i);

      gPad->SetLeftMargin(0.15);

      frame.GetYaxis()->SetTitleOffset(1.4);

      frame.Draw();

   };

   draw(1, *frame1);

   draw(2, *frame2);

   draw(3, *frame3);

}

f
#define f(i)
Definition RSha256.hxx:104

c
#define c(i)
Definition RSha256.hxx:101

RooAddPdf.h

RooCategory.h

RooChebychev.h

RooDataSet.h

RooGaussian.h

RooPlot.h

RooRealVar.h

RooSimultaneous.h

kDashed
@ kDashed
Definition TAttLine.h:54

TAxis.h

TCanvas.h

gPad
#define gPad
Definition TVirtualPad.h:322

RooAbsArg::Print
void Print(Option_t *options=nullptr) const override
Print the object to the defaultPrintStream().
Definition RooAbsArg.h:238

RooAddPdf
Efficient implementation of a sum of PDFs of the form.
Definition RooAddPdf.h:33

RooArgList
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22

RooArgSet
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:24

RooCategory
Object to represent discrete states.
Definition RooCategory.h:28

RooChebychev
Chebychev polynomial p.d.f.
Definition RooChebychev.h:25

RooDataSet
Container class to hold unbinned data.
Definition RooDataSet.h:32

RooGaussian
Plain Gaussian p.d.f.
Definition RooGaussian.h:24

RooPlot
Plot frame and a container for graphics objects within that frame.
Definition RooPlot.h:43

RooRealVar
Variable that can be changed from the outside.
Definition RooRealVar.h:37

RooSimultaneous
Facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset.
Definition RooSimultaneous.h:41

TCanvas
The Canvas class.
Definition TCanvas.h:23

RooFit::Title
RooCmdArg Title(const char *name)
Definition RooGlobalFunc.cxx:532

RooFit::Bins
RooCmdArg Bins(Int_t nbin)
Definition RooGlobalFunc.cxx:536

RooFit::Import
RooCmdArg Import(const char *state, TH1 &histo)
Definition RooGlobalFunc.cxx:375

RooFit::Save
RooCmdArg Save(bool flag=true)
Definition RooGlobalFunc.cxx:651

RooFit::PrintLevel
RooCmdArg PrintLevel(Int_t code)
Definition RooGlobalFunc.cxx:659

RooFit::Components
RooCmdArg Components(Args_t &&... argsOrArgSet)
Definition RooGlobalFunc.h:141

RooFit::ProjWData
RooCmdArg ProjWData(const RooAbsData &projData, bool binData=false)
Definition RooGlobalFunc.cxx:193

RooFit::Cut
RooCmdArg Cut(const char *cutSpec)
Definition RooGlobalFunc.cxx:309

RooFit::LineStyle
RooCmdArg LineStyle(Style_t style)
Definition RooGlobalFunc.cxx:233

sigma
const Double_t sigma
Definition h1analysisProxy.h:11

x
Double_t x[n]
Definition legend1.C:17

RooFit
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition CodegenImpl.h:72

rf501_simultaneouspdf
Definition rf501_simultaneouspdf.py:1

rf503_wspaceread.data
data
Definition rf503_wspaceread.py:37

rf503_wspaceread.model
model
Definition rf503_wspaceread.py:36