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.

 
 
␛[1mRooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby␛[0m 
                Copyright (C) 2000-2013 NIKHEF, University of California & Stanford University
                All rights reserved, please read http://roofit.sourceforge.net/license.txt
 
[#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 0x562873dffef0 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 0x56287415f780 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[psi]_NORM_psi for nset (psi) with code 0 from preexisting content.
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(Mf) creating new cache 0x56287415f780 with pdf Tpsi_CONV_Rpsi_CACHE_Obs[psi]_NORM_psi for nset (psi) with code 0 from preexisting content.
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Tpsi_Int[psi]) using numeric integrator RooIntegrator1D to calculate Int(psi)
[#1] INFO:Minimization --  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=9480.16 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   5.40962e+01
   2  greso        3.00000e-01   9.00000e-02   2.46497e-01   5.76705e+00
                               ERR DEF= 0.5
 MIGRAD MINIMIZATION HAS CONVERGED.
 MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
 COVARIANCE MATRIX CALCULATED SUCCESSFULLY
 FCN=9479.64 FROM MIGRAD    STATUS=CONVERGED      32 CALLS          33 TOTAL
                     EDM=1.07647e-07    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  gbias        1.92485e-01   7.44189e-03   1.26877e-03   1.36618e-02
   2  greso        2.98582e-01   9.19693e-03   1.65656e-03  -8.25561e-03
                               ERR DEF= 0.5
 EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  2    ERR DEF=0.5
  5.539e-05  1.690e-07 
  1.690e-07  8.460e-05 
 PARAMETER  CORRELATION COEFFICIENTS  
       NO.  GLOBAL      1      2
        1  0.00247   1.000  0.002
        2  0.00247   0.002  1.000
 **********
 **    7 **SET ERR         0.5
 **********
 **********
 **    8 **SET PRINT           1
 **********
 **********
 **    9 **HESSE        1000
 **********
 COVARIANCE MATRIX CALCULATED SUCCESSFULLY
 FCN=9479.64 FROM HESSE     STATUS=OK             10 CALLS          43 TOTAL
                     EDM=1.07707e-07    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                INTERNAL      INTERNAL  
  NO.   NAME      VALUE            ERROR       STEP SIZE       VALUE   
   1  gbias        1.92485e-01   7.44188e-03   5.07510e-05  -6.62425e-01
   2  greso        2.98582e-01   9.19691e-03   6.62625e-05  -5.92825e-01
                               ERR DEF= 0.5
 EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  2    ERR DEF=0.5
  5.539e-05  1.248e-07 
  1.248e-07  8.460e-05 
 PARAMETER  CORRELATION COEFFICIENTS  
       NO.  GLOBAL      1      2
        1  0.00182   1.000  0.002
        2  0.00182   0.002  1.000
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: 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 0x56287415f4c0 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 0x562874347640 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 0x5628742f6df0 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 0x5628742f6df0 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:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
[#1] INFO:Minimization --  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.92485e-01  7.44188e-03    0.00000e+00  1.00000e+00
     2 greso        2.98582e-01  9.19691e-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=5178.91 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.92485e-01   7.44188e-03   1.88791e-02   2.79436e+01
   2  greso        2.98582e-01   9.19691e-03   2.46483e-02  -8.41273e+00
                               ERR DEF= 0.5
 MIGRAD MINIMIZATION HAS CONVERGED.
 MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
 COVARIANCE MATRIX CALCULATED SUCCESSFULLY
 FCN=5178.63 FROM MIGRAD    STATUS=CONVERGED      28 CALLS          29 TOTAL
                     EDM=3.37796e-07    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  gbias        1.85814e-01   9.36461e-03   1.16493e-03   2.44691e-02
   2  greso        3.02469e-01   1.02714e-02   1.32671e-03   1.52048e-03
                               ERR DEF= 0.5
 EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  2    ERR DEF=0.5
  8.771e-05 -2.174e-05 
 -2.174e-05  1.055e-04 
 PARAMETER  CORRELATION COEFFICIENTS  
       NO.  GLOBAL      1      2
        1  0.22594   1.000 -0.226
        2  0.22594  -0.226  1.000
 **********
 **   16 **SET ERR         0.5
 **********
 **********
 **   17 **SET PRINT           1
 **********
 **********
 **   18 **HESSE        1000
 **********
 COVARIANCE MATRIX CALCULATED SUCCESSFULLY
 FCN=5178.63 FROM HESSE     STATUS=OK             10 CALLS          39 TOTAL
                     EDM=3.36905e-07    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                INTERNAL      INTERNAL  
  NO.   NAME      VALUE            ERROR       STEP SIZE       VALUE   
   1  gbias        1.85814e-01   9.36561e-03   2.32986e-04  -6.79459e-01
   2  greso        3.02469e-01   1.02725e-02   5.30684e-05  -5.82446e-01
                               ERR DEF= 0.5
 EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  2    ERR DEF=0.5
  8.773e-05 -2.179e-05 
 -2.179e-05  1.056e-04 
 PARAMETER  CORRELATION COEFFICIENTS  
       NO.  GLOBAL      1      2
        1  0.22639   1.000 -0.226
        2  0.22639  -0.226  1.000
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: 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 0x5628742f6df0 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)

 
#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
   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 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
   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 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();
}

Date

April 2009

Author

Wouter Verkerke

Definition in file rf210_angularconv.C.