rf210_angularconv.C File Reference

Detailed Description

'ADDITION AND CONVOLUTION' RooFit tutorial macro #210

Convolution in cyclical angular observables theta, and construction of p.d.f in terms of transformed angular coordinates, e.g. cos(theta), where the convolution is performed in theta rather than cos(theta)

(require ROOT to be compiled with –enable-fftw3)

pdf(theta) = T(theta) (x) gauss(theta) pdf(cosTheta) = T(acos(cosTheta)) (x) gauss(acos(cosTheta))

Processing /mnt/build/workspace/root-makedoc-v608/rootspi/rdoc/src/v6-08-00-patches/tutorials/roofit/rf210_angularconv.C...

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



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

   // Define angle psi
   RooRealVar psi("psi","psi",0,3.14159268) ;

   // Define physics p.d.f T(psi)
   RooGenericPdf Tpsi("Tpsi","1+sin(2*@0)",psi) ;

   // Define resolution R(psi)
   RooRealVar gbias("gbias","gbias",0.2,0.,1) ;
   RooRealVar greso("greso","greso",0.3,0.1,1.0) ;
   RooGaussian Rpsi("Rpsi","Rpsi",psi,gbias,greso) ;

   // Define cos(psi) and function psif that calculates psi from cos(psi)
   RooRealVar cpsi("cpsi","cos(psi)",-1,1) ;
   RooFormulaVar psif("psif","acos(cpsi)",cpsi) ;

   // Define physics p.d.f. also as function of cos(psi): T(psif(cpsi)) = T(cpsi) ;
   RooGenericPdf Tcpsi("T","1+sin(2*@0)",psif) ;



   // C o n s t r u c t   c o n v o l u t i o n   p d f  i n   p s i
   // --------------------------------------------------------------

   // Define convoluted p.d.f. as function of psi: M=[T(x)R](psi) = M(psi)
   RooFFTConvPdf Mpsi("Mf","Mf",psi,Tpsi,Rpsi) ;

   // Set the buffer fraction to zero to obtain a true cyclical convolution
   Mpsi.setBufferFraction(0) ;



   // S a m p l e ,   f i t   a n d   p l o t   c o n v o l u t e d   p d f  ( p s i )
   // --------------------------------------------------------------------------------

   // Generate some events in observable psi
   RooDataSet* data_psi = Mpsi.generate(psi,10000) ;

   // Fit convoluted model as function of angle psi
   Mpsi.fitTo(*data_psi) ;

   // Plot cos(psi) frame with Mf(cpsi)
   RooPlot* frame1 = psi.frame(Title("Cyclical convolution in angle psi")) ;
   data_psi->plotOn(frame1) ;
   Mpsi.plotOn(frame1) ;

   // Overlay comparison to unsmeared physics p.d.f T(psi)
   Tpsi.plotOn(frame1,LineColor(kRed)) ;



   // C o n s t r u c t   c o n v o l u t i o n   p d f   i n   c o s ( p s i )
   // --------------------------------------------------------------------------


   // Define convoluted p.d.f. as function of cos(psi): M=[T(x)R](psif(cpsi)) = M(cpsi)
   //
   // Need to give both observable psi here (for definition of convolution)
   // and function psif here (for definition of observables, ultimately in cpsi)
   RooFFTConvPdf Mcpsi("Mf","Mf",psif,psi,Tpsi,Rpsi) ;

   // Set the buffer fraction to zero to obtain a true cyclical convolution
   Mcpsi.setBufferFraction(0) ;


   // S a m p l e ,   f i t   a n d   p l o t   c o n v o l u t e d   p d f  ( c o s p s i )
   // --------------------------------------------------------------------------------

   // Generate some events
   RooDataSet* data_cpsi = Mcpsi.generate(cpsi,10000) ;

   // set psi constant to exclude to be a parameter of the fit
   psi.setConstant(true);

   // Fit convoluted model as function of cos(psi)
   Mcpsi.fitTo(*data_cpsi) ;

   // Plot cos(psi) frame with Mf(cpsi)
   RooPlot* frame2 = cpsi.frame(Title("Same convolution in psi, expressed in cos(psi)")) ;
   data_cpsi->plotOn(frame2) ;
   Mcpsi.plotOn(frame2) ;

   // Overlay comparison to unsmeared physics p.d.f Tf(cpsi)
   Tcpsi.plotOn(frame2,LineColor(kRed)) ;



   // Draw frame on canvas
   TCanvas* c = new TCanvas("rf210_angularconv","rf210_angularconv",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() ;

}

Author

04/2009 - Wouter Verkerke

Definition in file rf210_angularconv.C.