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 0 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 0x10f5de40 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 0x14284670 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 0x142bb880 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.25032 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 0x142eab40 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 0x14384ab0 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 0x143d7990 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 0x14253910 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 633.46 μ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 0x143d48b0 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 0x144b4b70 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.