ROOT logo

From $ROOTSYS/tutorials/roofit/rf501_simultaneouspdf.C

//////////////////////////////////////////////////////////////////////////
//
// '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 "RooConstVar.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),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)) ;



  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() ;


}
 rf501_simultaneouspdf.C:1
 rf501_simultaneouspdf.C:2
 rf501_simultaneouspdf.C:3
 rf501_simultaneouspdf.C:4
 rf501_simultaneouspdf.C:5
 rf501_simultaneouspdf.C:6
 rf501_simultaneouspdf.C:7
 rf501_simultaneouspdf.C:8
 rf501_simultaneouspdf.C:9
 rf501_simultaneouspdf.C:10
 rf501_simultaneouspdf.C:11
 rf501_simultaneouspdf.C:12
 rf501_simultaneouspdf.C:13
 rf501_simultaneouspdf.C:14
 rf501_simultaneouspdf.C:15
 rf501_simultaneouspdf.C:16
 rf501_simultaneouspdf.C:17
 rf501_simultaneouspdf.C:18
 rf501_simultaneouspdf.C:19
 rf501_simultaneouspdf.C:20
 rf501_simultaneouspdf.C:21
 rf501_simultaneouspdf.C:22
 rf501_simultaneouspdf.C:23
 rf501_simultaneouspdf.C:24
 rf501_simultaneouspdf.C:25
 rf501_simultaneouspdf.C:26
 rf501_simultaneouspdf.C:27
 rf501_simultaneouspdf.C:28
 rf501_simultaneouspdf.C:29
 rf501_simultaneouspdf.C:30
 rf501_simultaneouspdf.C:31
 rf501_simultaneouspdf.C:32
 rf501_simultaneouspdf.C:33
 rf501_simultaneouspdf.C:34
 rf501_simultaneouspdf.C:35
 rf501_simultaneouspdf.C:36
 rf501_simultaneouspdf.C:37
 rf501_simultaneouspdf.C:38
 rf501_simultaneouspdf.C:39
 rf501_simultaneouspdf.C:40
 rf501_simultaneouspdf.C:41
 rf501_simultaneouspdf.C:42
 rf501_simultaneouspdf.C:43
 rf501_simultaneouspdf.C:44
 rf501_simultaneouspdf.C:45
 rf501_simultaneouspdf.C:46
 rf501_simultaneouspdf.C:47
 rf501_simultaneouspdf.C:48
 rf501_simultaneouspdf.C:49
 rf501_simultaneouspdf.C:50
 rf501_simultaneouspdf.C:51
 rf501_simultaneouspdf.C:52
 rf501_simultaneouspdf.C:53
 rf501_simultaneouspdf.C:54
 rf501_simultaneouspdf.C:55
 rf501_simultaneouspdf.C:56
 rf501_simultaneouspdf.C:57
 rf501_simultaneouspdf.C:58
 rf501_simultaneouspdf.C:59
 rf501_simultaneouspdf.C:60
 rf501_simultaneouspdf.C:61
 rf501_simultaneouspdf.C:62
 rf501_simultaneouspdf.C:63
 rf501_simultaneouspdf.C:64
 rf501_simultaneouspdf.C:65
 rf501_simultaneouspdf.C:66
 rf501_simultaneouspdf.C:67
 rf501_simultaneouspdf.C:68
 rf501_simultaneouspdf.C:69
 rf501_simultaneouspdf.C:70
 rf501_simultaneouspdf.C:71
 rf501_simultaneouspdf.C:72
 rf501_simultaneouspdf.C:73
 rf501_simultaneouspdf.C:74
 rf501_simultaneouspdf.C:75
 rf501_simultaneouspdf.C:76
 rf501_simultaneouspdf.C:77
 rf501_simultaneouspdf.C:78
 rf501_simultaneouspdf.C:79
 rf501_simultaneouspdf.C:80
 rf501_simultaneouspdf.C:81
 rf501_simultaneouspdf.C:82
 rf501_simultaneouspdf.C:83
 rf501_simultaneouspdf.C:84
 rf501_simultaneouspdf.C:85
 rf501_simultaneouspdf.C:86
 rf501_simultaneouspdf.C:87
 rf501_simultaneouspdf.C:88
 rf501_simultaneouspdf.C:89
 rf501_simultaneouspdf.C:90
 rf501_simultaneouspdf.C:91
 rf501_simultaneouspdf.C:92
 rf501_simultaneouspdf.C:93
 rf501_simultaneouspdf.C:94
 rf501_simultaneouspdf.C:95
 rf501_simultaneouspdf.C:96
 rf501_simultaneouspdf.C:97
 rf501_simultaneouspdf.C:98
 rf501_simultaneouspdf.C:99
 rf501_simultaneouspdf.C:100
 rf501_simultaneouspdf.C:101
 rf501_simultaneouspdf.C:102
 rf501_simultaneouspdf.C:103
 rf501_simultaneouspdf.C:104
 rf501_simultaneouspdf.C:105
 rf501_simultaneouspdf.C:106
 rf501_simultaneouspdf.C:107
 rf501_simultaneouspdf.C:108
 rf501_simultaneouspdf.C:109
 rf501_simultaneouspdf.C:110
 rf501_simultaneouspdf.C:111
 rf501_simultaneouspdf.C:112
 rf501_simultaneouspdf.C:113
 rf501_simultaneouspdf.C:114
 rf501_simultaneouspdf.C:115
 rf501_simultaneouspdf.C:116
 rf501_simultaneouspdf.C:117
 rf501_simultaneouspdf.C:118
 rf501_simultaneouspdf.C:119
 rf501_simultaneouspdf.C:120
 rf501_simultaneouspdf.C:121
 rf501_simultaneouspdf.C:122
 rf501_simultaneouspdf.C:123
 rf501_simultaneouspdf.C:124
 rf501_simultaneouspdf.C:125
 rf501_simultaneouspdf.C:126
 rf501_simultaneouspdf.C:127
 rf501_simultaneouspdf.C:128
 rf501_simultaneouspdf.C:129
 rf501_simultaneouspdf.C:130
 rf501_simultaneouspdf.C:131
 rf501_simultaneouspdf.C:132
 rf501_simultaneouspdf.C:133
 rf501_simultaneouspdf.C:134
 rf501_simultaneouspdf.C:135
 rf501_simultaneouspdf.C:136
 rf501_simultaneouspdf.C:137
 rf501_simultaneouspdf.C:138
 rf501_simultaneouspdf.C:139
 rf501_simultaneouspdf.C:140
 rf501_simultaneouspdf.C:141
 rf501_simultaneouspdf.C:142
 rf501_simultaneouspdf.C:143
 rf501_simultaneouspdf.C:144
 rf501_simultaneouspdf.C:145
 rf501_simultaneouspdf.C:146
 rf501_simultaneouspdf.C:147