Logo ROOT   6.14/05
Reference Guide
rf210_angularconv.C File Reference
Tutorials » RooFit Tutorials

Detailed Description

View in nbviewer Open in SWAN '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))

pict1_rf210_angularconv.C.png
Processing /mnt/build/workspace/root-makedoc-v614/rootspi/rdoc/src/v6-14-00-patches/tutorials/roofit/rf210_angularconv.C...
RooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby
Copyright (C) 2000-2013 NIKHEF, University of California & Stanford University
All rights reserved, please read http://roofit.sourceforge.net/license.txt
[#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 0x2010ed0 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[psi] for nset (psi) with code 0
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(Mf) creating new cache 0x26e46b0 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[psi] for nset (psi) with code 0 from preexisting content.
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(Mf) creating new cache 0x2573420 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[psi] for nset (psi) with code 0 from preexisting content.
[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_Int[psi]) using numeric integrator RooIntegrator1D to calculate Int(psi)
[#1] INFO:Minization -- The following expressions have been identified as constant and will be precalculated and cached: (Tpsi)
**********
** 1 **SET PRINT 1
**********
**********
** 2 **SET NOGRAD
**********
PARAMETER DEFINITIONS:
NO. NAME VALUE STEP SIZE LIMITS
1 gbias 2.00000e-01 1.00000e-01 0.00000e+00 1.00000e+00
2 greso 3.00000e-01 9.00000e-02 1.00000e-01 1.00000e+00
**********
** 3 **SET ERR 0.5
**********
**********
** 4 **SET PRINT 1
**********
**********
** 5 **SET STR 1
**********
NOW USING STRATEGY 1: TRY TO BALANCE SPEED AGAINST RELIABILITY
**********
** 6 **MIGRAD 1000 1
**********
FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4.
START MIGRAD MINIMIZATION. STRATEGY 1. CONVERGENCE WHEN EDM .LT. 1.00e-03
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_Int[psi]) using numeric integrator RooIntegrator1D to calculate Int(psi)
FCN=9461.22 FROM MIGRAD STATUS=INITIATE 8 CALLS 9 TOTAL
EDM= unknown STRATEGY= 1 NO ERROR MATRIX
EXT PARAMETER CURRENT GUESS STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 gbias 2.00000e-01 1.00000e-01 2.57889e-01 4.71903e+01
2 greso 3.00000e-01 9.00000e-02 2.46497e-01 1.58737e+01
ERR DEF= 0.5
MIGRAD MINIMIZATION HAS CONVERGED.
MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
COVARIANCE MATRIX CALCULATED SUCCESSFULLY
FCN=9460.75 FROM MIGRAD STATUS=CONVERGED 32 CALLS 33 TOTAL
EDM=1.54674e-08 STRATEGY= 1 ERROR MATRIX ACCURATE
EXT PARAMETER STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 gbias 1.93440e-01 7.42287e-03 1.26201e-03 5.36276e-03
2 greso 2.96325e-01 9.18155e-03 1.65893e-03 -2.96198e-03
ERR DEF= 0.5
EXTERNAL ERROR MATRIX. NDIM= 25 NPAR= 2 ERR DEF=0.5
5.511e-05 2.076e-07
2.076e-07 8.432e-05
PARAMETER CORRELATION COEFFICIENTS
NO. GLOBAL 1 2
1 0.00305 1.000 0.003
2 0.00305 0.003 1.000
**********
** 7 **SET ERR 0.5
**********
**********
** 8 **SET PRINT 1
**********
**********
** 9 **HESSE 1000
**********
COVARIANCE MATRIX CALCULATED SUCCESSFULLY
FCN=9460.75 FROM HESSE STATUS=OK 10 CALLS 43 TOTAL
EDM=1.54739e-08 STRATEGY= 1 ERROR MATRIX ACCURATE
EXT PARAMETER INTERNAL INTERNAL
NO. NAME VALUE ERROR STEP SIZE VALUE
1 gbias 1.93440e-01 7.42286e-03 5.04803e-05 -6.60005e-01
2 greso 2.96325e-01 9.18153e-03 6.63572e-05 -5.98886e-01
ERR DEF= 0.5
EXTERNAL ERROR MATRIX. NDIM= 25 NPAR= 2 ERR DEF=0.5
5.511e-05 1.636e-07
1.636e-07 8.432e-05
PARAMETER CORRELATION COEFFICIENTS
NO. GLOBAL 1 2
1 0.00240 1.000 0.002
2 0.00240 0.002 1.000
[#1] INFO:Minization -- RooMinimizer::optimizeConst: deactivating const optimization
[#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 0x1869270 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[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 0x2010ed0 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[cpsi] for nset (cpsi) with code 0 from preexisting content.
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_CONV_Rpsi_CACHE_Obs[cpsi]_Int[cpsi]) using numeric integrator RooIntegrator1D to calculate Int(cpsi)
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(Mf) creating new cache 0x25cf3f0 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[cpsi] for nset (cpsi) with code 0 from preexisting content.
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_CONV_Rpsi_CACHE_Obs[cpsi]_Int[cpsi]) using numeric integrator RooIntegrator1D to calculate Int(cpsi)
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(Mf) creating new cache 0x27b40f0 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[cpsi] for nset (cpsi) with code 0 from preexisting content.
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_CONV_Rpsi_CACHE_Obs[cpsi]_Int[cpsi]) using numeric integrator RooIntegrator1D to calculate Int(cpsi)
[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
[#1] INFO:Minization -- The following expressions have been identified as constant and will be precalculated and cached: (psif)
**********
** 10 **SET PRINT 1
**********
**********
** 11 **SET NOGRAD
**********
PARAMETER DEFINITIONS:
NO. NAME VALUE STEP SIZE LIMITS
1 gbias 1.93440e-01 7.42286e-03 0.00000e+00 1.00000e+00
2 greso 2.96325e-01 9.18153e-03 1.00000e-01 1.00000e+00
**********
** 12 **SET ERR 0.5
**********
**********
** 13 **SET PRINT 1
**********
**********
** 14 **SET STR 1
**********
NOW USING STRATEGY 1: TRY TO BALANCE SPEED AGAINST RELIABILITY
**********
** 15 **MIGRAD 1000 1
**********
FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4.
START MIGRAD MINIMIZATION. STRATEGY 1. CONVERGENCE WHEN EDM .LT. 1.00e-03
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_Int[psi]) using numeric integrator RooIntegrator1D to calculate Int(psi)
FCN=5106.07 FROM MIGRAD STATUS=INITIATE 4 CALLS 5 TOTAL
EDM= unknown STRATEGY= 1 NO ERROR MATRIX
EXT PARAMETER CURRENT GUESS STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 gbias 1.93440e-01 7.42286e-03 1.87954e-02 2.56568e+01
2 greso 2.96325e-01 9.18153e-03 2.47085e-02 2.78281e+01
ERR DEF= 0.5
MIGRAD MINIMIZATION HAS CONVERGED.
MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
COVARIANCE MATRIX CALCULATED SUCCESSFULLY
FCN=5105.7 FROM MIGRAD STATUS=CONVERGED 30 CALLS 31 TOTAL
EDM=2.64459e-09 STRATEGY= 1 ERROR MATRIX ACCURATE
EXT PARAMETER STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 gbias 1.89592e-01 9.16795e-03 1.12282e-03 2.12532e-03
2 greso 2.89839e-01 1.02617e-02 1.33840e-03 1.05386e-03
ERR DEF= 0.5
EXTERNAL ERROR MATRIX. NDIM= 25 NPAR= 2 ERR DEF=0.5
8.407e-05 -2.232e-05
-2.232e-05 1.053e-04
PARAMETER CORRELATION COEFFICIENTS
NO. GLOBAL 1 2
1 0.23715 1.000 -0.237
2 0.23715 -0.237 1.000
**********
** 16 **SET ERR 0.5
**********
**********
** 17 **SET PRINT 1
**********
**********
** 18 **HESSE 1000
**********
COVARIANCE MATRIX CALCULATED SUCCESSFULLY
FCN=5105.7 FROM HESSE STATUS=OK 10 CALLS 41 TOTAL
EDM=2.64028e-09 STRATEGY= 1 ERROR MATRIX ACCURATE
EXT PARAMETER INTERNAL INTERNAL
NO. NAME VALUE ERROR STEP SIZE VALUE
1 gbias 1.89592e-01 9.16907e-03 4.49130e-05 -6.69783e-01
2 greso 2.89839e-01 1.02630e-02 5.35358e-05 -6.16443e-01
ERR DEF= 0.5
EXTERNAL ERROR MATRIX. NDIM= 25 NPAR= 2 ERR DEF=0.5
8.409e-05 -2.237e-05
-2.237e-05 1.054e-04
PARAMETER CORRELATION COEFFICIENTS
NO. GLOBAL 1 2
1 0.23764 1.000 -0.238
2 0.23764 -0.238 1.000
[#1] INFO:Minization -- RooMinimizer::optimizeConst: deactivating const optimization
[#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 0x25c0950 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[cpsi] for nset (cpsi) with code 0
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_CONV_Rpsi_CACHE_Obs[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)
#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.