rf211_paramconv.C

Go to the documentation of this file.

/// \file
/// \ingroup tutorial_roofit
/// \notebook
///  'ADDITION AND CONVOLUTION' RooFit tutorial macro #211
///
///  Working a with a p.d.f. with a convolution operator in terms
///  of a parameter
///
///  (require ROOT to be compiled with --enable-fftw3)
///
/// \macro_image
/// \macro_output
/// \macro_code
/// \author 04/2009 - Wouter Verkerke 


#include "RooRealVar.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "RooGenericPdf.h"
#include "RooFormulaVar.h"
#include "RooFFTConvPdf.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "TH2.h"
using namespace RooFit ;


void rf211_paramconv()
{
   // S e t u p   c o m p o n e n t   p d f s 
   // ---------------------------------------

   // Gaussian g(x ; mean,sigma) 
   RooRealVar x("x","x",-10,10) ;
   RooRealVar mean("mean","mean",-3,3) ;
   RooRealVar sigma("sigma","sigma",0.5,0.1,10) ;
   RooGaussian modelx("gx","gx",x,mean,sigma) ;

   // Block function in mean
   RooRealVar a("a","a",2,1,10) ;
   RooGenericPdf model_mean("model_mean","abs(mean)<a",RooArgList(mean,a)) ;

   // Convolution in mean parameter model = g(x,mean,sigma) (x) block(mean)
   x.setBins(1000,"cache") ;
   mean.setBins(50,"cache") ;
   RooFFTConvPdf model("model","model",mean,modelx,model_mean) ;

   // Configure convolution to construct a 2-D cache in (x,mean)
   // rather than a 1-d cache in mean that needs to be recalculated
   // for each value of x
   model.setCacheObservables(x) ;
   model.setBufferFraction(1.0) ;

   // Integrate model over mean projModel = Int model dmean
   RooAbsPdf* projModel = model.createProjection(mean) ;

   // Generate 1000 toy events
   RooDataHist* d = projModel->generateBinned(x,1000) ;

   // Fit p.d.f. to toy data
   projModel->fitTo(*d,Verbose()) ;

   // Plot data and fitted p.d.f.
   RooPlot* frame = x.frame(Bins(25)) ;
   d->plotOn(frame) ;
   projModel->plotOn(frame) ;

   // Make 2d histogram of model(x;mean)
   TH1* hh = model.createHistogram("hh",x,Binning(50),YVar(mean,Binning(50)),ConditionalObservables(mean)) ;
   hh->SetTitle("histogram of model(x|mean)") ;
   hh->SetLineColor(kBlue) ;

   // Draw frame on canvas
   TCanvas* c = new TCanvas("rf211_paramconv","rf211_paramconv",800,400) ;
   c->Divide(2) ;
   c->cd(1) ; gPad->SetLeftMargin(0.15) ; frame->GetYaxis()->SetTitleOffset(1.4) ; frame->Draw() ;
   c->cd(2) ; gPad->SetLeftMargin(0.20) ; hh->GetZaxis()->SetTitleOffset(2.5) ; hh->Draw("surf") ;

}