rf210_angularconv.C File Reference

'ADDITION AND CONVOLUTION' RooFit tutorial macro #210 More...

Go to the source code of this file.

Detailed Description

'ADDITION AND CONVOLUTION' RooFit tutorial macro #210

Convolution in cyclical angular observables theta, and construction of p.d.f in terms of trasnformed 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))

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#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.