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.
Author: Wouter Verkerke
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Tuesday, March 19, 2024 at 07:15 PM.
%%cpp -d
#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;
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);
Define convoluted pdf as function of psi: M=T(x)R = M(psi)
RooFFTConvPdf Mpsi("Mf", "Mf", psi, Tpsi, Rpsi);
[#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'.
Set the buffer fraction to zero to obtain a true cyclical convolution
Mpsi.setBufferFraction(0);
Generate some events in observable psi
std::unique_ptr<RooDataSet> data_psi{Mpsi.generate(psi, 10000)};
[#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 0x7efe94ebbd20 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[psi]_NORM_psi for nset (psi) with code 0 [#1] INFO:Caching -- RooAbsCachedPdf::getCache(Mf) creating new cache 0x7efe95b6ab90 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[psi]_NORM_psi for nset (psi) with code 0 from preexisting content.
Fit convoluted model as function of angle psi
Mpsi.fitTo(*data_psi, PrintLevel(-1));
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(Mf) creating new cache 0x7efe95c2e390 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[psi]_NORM_psi for nset (psi) with code 0 from preexisting content. [#1] INFO:Fitting -- RooAbsPdf::fitTo(Mf_over_Mf_Int[psi]) fixing normalization set for coefficient determination to observables in data [#1] INFO:Fitting -- using CPU computation library compiled with -mavx2 [#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_Mf_over_Mf_Int[psi]_MfData) Summation contains a RooNLLVar, using its error level [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization [#1] INFO:Caching -- RooAbsCachedPdf::getCache(Mf) creating new cache 0x7efe95de0a50 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[psi] for nset () with code 1 from preexisting content. [#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_Int[psi]) using numeric integrator RooIntegrator1D to calculate Int(psi) [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization
Plot cos(psi) frame with Mf(cpsi)
RooPlot *frame1 = psi.frame(Title("Cyclical convolution in angle psi"));
data_psi->plotOn(frame1);
Mpsi.plotOn(frame1);
[#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 0x7efe95dba100 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[psi]_NORM_psi for nset (psi) with code 0
Overlay comparison to unsmeared physics pdf T(psi)
Tpsi.plotOn(frame1, LineColor(kRed));
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_Int[psi]) using numeric integrator RooIntegrator1D to calculate Int(psi)
Define convoluted pdf as function of cos(psi): M=T(x)R = 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);
Generate some events
std::unique_ptr<RooDataSet> data_cpsi{Mcpsi.generate(cpsi, 10000)};
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(Mf) creating new cache 0x7efe957de7d0 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[cpsi]_NORM_cpsi for nset (cpsi) with code 0 from preexisting content. [#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:Caching -- RooAbsCachedPdf::getCache(Mf) creating new cache 0x7efe95e81440 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[cpsi]_NORM_cpsi for nset (cpsi) with code 0 from preexisting content. [#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_CONV_Rpsi_CACHE_Obs[cpsi]_NORM_cpsi_Int[cpsi]) using numeric integrator RooIntegrator1D to calculate Int(cpsi)
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));
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(Mf) creating new cache 0x7efe95f5b660 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[cpsi]_NORM_cpsi for nset (cpsi) with code 0 from preexisting content. [#1] INFO:Fitting -- RooAbsPdf::fitTo(Mf_over_Mf_Int[cpsi]) fixing normalization set for coefficient determination to observables in data [#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_Mf_over_Mf_Int[cpsi]_MfData) Summation contains a RooNLLVar, using its error level [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization [#1] INFO:Caching -- RooAbsCachedPdf::getCache(Mf) creating new cache 0x7efe95f65520 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[cpsi] for nset () with code 1 from preexisting content. [#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:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization
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);
[#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 0x7efe95e3f270 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)
Overlay comparison to unsmeared physics pdf Tf(cpsi)
Tcpsi.plotOn(frame2, LineColor(kRed));
[#1] INFO:NumericIntegration -- RooRealIntegral::init(T_Int[cpsi]) using numeric integrator RooIntegrator1D to calculate Int(cpsi)
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();
Draw all canvases
%jsroot on
gROOT->GetListOfCanvases()->Draw()