rf210_angularconv.C File Reference

Detailed Description

Addition and convolution: convolution in cyclical angular observables theta

and construction of pdf in terms of transformed angular coordinates, e.g. cos(theta), where the convolution is performed in theta rather than cos(theta)

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

This tutorial requires FFT3 to be enabled.

 
#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 pdf 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 pdf 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 pdf 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
   std::unique_ptr<RooDataSet> data_psi{Mpsi.generate(psi, 10000)};
 
   // Fit convoluted model as function of angle psi
   Mpsi.fitTo(*data_psi, PrintLevel(-1));
 
   // 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 pdf 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 pdf 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
   std::unique_ptr<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, PrintLevel(-1));
 
   // 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 pdf 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();
}

[#1] INFO:Caching -- Changing internal binning of variable 'psi' in FFT 'Mf' from 100 to 930 to improve the precision of the numerical FFT. This can be done manually by setting an additional binning named 'cache'.
[#1] INFO:Eval -- RooRealVar::setRange(psi) new range named 'refrange_fft_Mf' created with bounds [0,3.14159]
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_Int[psi]) using numeric integrator RooIntegrator1D to calculate Int(psi)
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(Mf) creating new cache 0x37083a0 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[psi]_NORM_psi for nset (psi) with code 0
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_Int[psi]) using numeric integrator RooIntegrator1D to calculate Int(psi)
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(Mf) creating new cache 0x50c0c60 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[psi]_NORM_psi for nset (psi) with code 0
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_Int[psi]) using numeric integrator RooIntegrator1D to calculate Int(psi)
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(Mf) creating new cache 0x508c160 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[psi]_NORM_psi for nset (psi) with code 0
[#1] INFO:Fitting -- RooAbsPdf::fitTo(Mf_over_Mf_Int[psi]) fixing normalization set for coefficient determination to observables in data
[#1] INFO:Fitting -- using generic CPU library compiled with no vectorizations
[#1] INFO:Fitting -- Creation of NLL object took 1.1676 ms
[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_Mf_over_Mf_Int[psi]_MfData) Summation contains a RooNLLVar, using its error level
[#1] INFO:Minimization -- [fitFCN] No discrete parameters, performing continuous minimization only
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_Int[psi]) using numeric integrator RooIntegrator1D to calculate Int(psi)
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(Mf) creating new cache 0x51b45e0 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[psi]_NORM_psi for nset (psi) with code 0
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_Int[psi]) using numeric integrator RooIntegrator1D to calculate Int(psi)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_Int[psi]) using numeric integrator RooIntegrator1D to calculate Int(psi)
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(Mf) creating new cache 0x5057870 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[cpsi]_NORM_cpsi for nset (cpsi) with code 0
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_CONV_Rpsi_CACHE_Obs[cpsi]_NORM_cpsi_Int[cpsi]) using numeric integrator RooIntegrator1D to calculate Int(cpsi)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_Int[psi]) using numeric integrator RooIntegrator1D to calculate Int(psi)
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(Mf) creating new cache 0x24120a0 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[cpsi]_NORM_cpsi for nset (cpsi) with code 0
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_CONV_Rpsi_CACHE_Obs[cpsi]_NORM_cpsi_Int[cpsi]) using numeric integrator RooIntegrator1D to calculate Int(cpsi)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_Int[psi]) using numeric integrator RooIntegrator1D to calculate Int(psi)
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(Mf) creating new cache 0x5057870 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[cpsi]_NORM_cpsi for nset (cpsi) with code 0
[#1] INFO:Fitting -- RooAbsPdf::fitTo(Mf_over_Mf_Int[cpsi]) fixing normalization set for coefficient determination to observables in data
[#1] INFO:Fitting -- Creation of NLL object took 638.869 μs
[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_Mf_over_Mf_Int[cpsi]_MfData) Summation contains a RooNLLVar, using its error level
[#1] INFO:Minimization -- [fitFCN] No discrete parameters, performing continuous minimization only
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Mf_Int[cpsi]) using numeric integrator RooIntegrator1D to calculate Int(cpsi)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_Int[psi]) using numeric integrator RooIntegrator1D to calculate Int(psi)
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(Mf) creating new cache 0x526b380 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[cpsi] for nset () with code 1
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_Int[psi]) using numeric integrator RooIntegrator1D to calculate Int(psi)
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(Mf) creating new cache 0x4edf270 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[cpsi]_NORM_cpsi for nset (cpsi) with code 0
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_CONV_Rpsi_CACHE_Obs[cpsi]_NORM_cpsi_Int[cpsi]) using numeric integrator RooIntegrator1D to calculate Int(cpsi)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(T_Int[cpsi]) using numeric integrator RooIntegrator1D to calculate Int(cpsi)

Date

April 2009

Author

Wouter Verkerke

Definition in file rf210_angularconv.C.