ROOT logo

From $ROOTSYS/tutorials/roofit/rf303_conditional.C

/////////////////////////////////////////////////////////////////////////
//
// '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 "RooDataHist.h"
#include "RooGaussian.h"
#include "RooPolyVar.h"
#include "RooProdPdf.h"
#include "RooPlot.h"
#include "TRandom.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "TH1.h"
using namespace RooFit ;


RooDataSet* makeFakeDataXY() ;

void rf303_conditional()
{
  // 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)) ;
  data->Print() ;



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


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


  // Make canvas and draw RooPlots
  new TCanvas("rf303_conditional","rf303_conditional",600, 460);
  gPad->SetLeftMargin(0.15) ; xframe->GetYaxis()->SetTitleOffset(1.2) ; xframe->Draw() ;
}




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 ;
}

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