rf501_simultaneouspdf.C File Reference

Detailed Description

Organisation and simultaneous fits: using simultaneous pdfs to describe simultaneous fits to multiple datasets

 
#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
   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}, {"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:
   // 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, PrintLevel(-1), 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(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.
   // 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.
   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));
 
   TCanvas *c = new TCanvas("rf501_simultaneouspdf", "rf403_simultaneouspdf", 800, 400);
   c->Divide(2);
   c->cd(1);
   gPad->SetLeftMargin(0.15);
   frame1->GetYaxis()->SetTitleOffset(1.4);
   frame1->Draw();
   c->cd(2);
   gPad->SetLeftMargin(0.15);
   frame2->GetYaxis()->SetTitleOffset(1.4);
   frame2->Draw();
}

 
RooAbsTestStatistic::initSimMode: creating slave calculator #0 for state control (2000 dataset entries)
RooAbsTestStatistic::initSimMode: creating slave calculator #1 for state physics (100 dataset entries)
[#1] INFO:Fitting -- RooAbsTestStatistic::initSimMode: created 2 slave calculators.
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
[#1] INFO:Minimization --  The following expressions will be evaluated in cache-and-track mode: (gx_ctl,px_ctl)
[#1] INFO:Minimization --  The following expressions will be evaluated in cache-and-track mode: (gx,px)
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization
 
  RooFitResult: minimized FCN value: 5773.95, estimated distance to minimum: 9.09933e-05
                covariance matrix quality: Full, accurate covariance matrix
                Status : MINIMIZE=0 HESSE=0 
 
    Floating Parameter    FinalValue +/-  Error   
  --------------------  --------------------------
                    a0    2.4369e-01 +/-  1.75e-01
                a0_ctl   -4.3957e-03 +/-  5.29e-02
                    a1    5.2692e-02 +/-  1.78e-01
                a1_ctl    5.4474e-01 +/-  3.70e-02
                     f    6.8391e-02 +/-  3.86e-02
                 f_ctl    5.0279e-01 +/-  1.24e-02
                  mean   -4.6480e-01 +/-  2.34e-01
              mean_ctl   -3.0263e+00 +/-  1.08e-02
                 sigma    3.0786e-01 +/-  8.77e-03
 
[#1] INFO:Plotting -- RooTreeData::plotOn: plotting 100 events out of 2100 total events
[#1] INFO:Plotting -- RooSimultaneous::plotOn(simPdf) plot on x represents a slice in the index category (sample)
[#1] INFO:Plotting -- RooAbsReal::plotOn(model) slice variable sample was not projected anyway
[#1] INFO:Plotting -- RooSimultaneous::plotOn(simPdf) plot on x represents a slice in the index category (sample)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (px)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()
[#1] INFO:Plotting -- RooAbsReal::plotOn(model) slice variable sample was not projected anyway
[#1] INFO:Plotting -- RooTreeData::plotOn: plotting 2000 events out of 2100 total events
[#1] INFO:Plotting -- RooSimultaneous::plotOn(simPdf) plot on x represents a slice in the index category (sample)
[#1] INFO:Plotting -- RooAbsReal::plotOn(model_ctl) slice variable sample was not projected anyway
[#1] INFO:Plotting -- RooSimultaneous::plotOn(simPdf) plot on x represents a slice in the index category (sample)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model_ctl) directly selected PDF components: (px_ctl)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model_ctl) indirectly selected PDF components: ()
[#1] INFO:Plotting -- RooAbsReal::plotOn(model_ctl) slice variable sample was not projected anyway

Date

July 2008

Author

Wouter Verkerke

Definition in file rf501_simultaneouspdf.C.