# rf704_amplitudefit¶

Special pdf's: using a pdf defined by a sum of real-valued amplitude components

Author: Clemens Lange, Wouter Verkerke (C++ version)
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Wednesday, November 30, 2022 at 11:24 AM.

In [1]:
import ROOT

Welcome to JupyROOT 6.27/01


## Setup 2D amplitude functions¶

Observables

In [2]:
t = ROOT.RooRealVar("t", "time", -1.0, 15.0)
cosa = ROOT.RooRealVar("cosa", "cos(alpha)", -1.0, 1.0)


Use ROOT.RooTruthModel to obtain compiled implementation of sinh/cosh modulated decay functions

In [3]:
tau = ROOT.RooRealVar("tau", "#tau", 1.5)
deltaGamma = ROOT.RooRealVar("deltaGamma", "deltaGamma", 0.3)
tm = ROOT.RooTruthModel("tm", "tm", t)
coshGBasis = ROOT.RooFormulaVar("coshGBasis", "exp(-@0/ @1)*cosh(@0*@2/2)", [t, tau, deltaGamma])
sinhGBasis = ROOT.RooFormulaVar("sinhGBasis", "exp(-@0/ @1)*sinh(@0*@2/2)", [t, tau, deltaGamma])
coshGConv = tm.convolution(coshGBasis, t)
sinhGConv = tm.convolution(sinhGBasis, t)


Construct polynomial amplitudes in cos(a)

In [4]:
poly1 = ROOT.RooPolyVar("poly1", "poly1", cosa, [0.5, 0.2, 0.2], 0)
poly2 = ROOT.RooPolyVar("poly2", "poly2", cosa, [1.0, -0.2, 3.0], 0)


Construct 2D amplitude as uncorrelated product of amp(t)*amp(cosa)

In [5]:
ampl1 = ROOT.RooProduct("ampl1", "amplitude 1", [poly1, coshGConv])
ampl2 = ROOT.RooProduct("ampl2", "amplitude 2", [poly2, sinhGConv])


## Construct amplitude sum pdf¶

Amplitude strengths

In [6]:
f1 = ROOT.RooRealVar("f1", "f1", 1, 0, 2)
f2 = ROOT.RooRealVar("f2", "f2", 0.5, 0, 2)


Construct pdf

In [7]:
pdf = ROOT.RooRealSumPdf("pdf", "pdf", [ampl1, ampl2], [f1, f2])


Generate some toy data from pdf

In [8]:
data = pdf.generate({t, cosa}, 10000)

[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_coshGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_sinhGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_coshGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_sinhGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)


Fit pdf to toy data with only amplitude strength floating

In [9]:
pdf.fitTo(data)

Out[9]:
<cppyy.gbl.RooFitResult object at 0x(nil)>
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_coshGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_sinhGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
[#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: (ampl1,ampl2)
**********
**    1 **SET PRINT           1
**********
**********
**********
PARAMETER DEFINITIONS:
NO.   NAME         VALUE      STEP SIZE      LIMITS
1 f1           1.00000e+00  2.00000e-01    0.00000e+00  2.00000e+00
2 f2           5.00000e-01  2.00000e-01    0.00000e+00  2.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
**********
**********
FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4.
START MIGRAD MINIMIZATION.  STRATEGY  1.  CONVERGENCE WHEN EDM .LT. 1.00e-03
FCN=24323.4 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  f1           1.00000e+00   2.00000e-01   2.01358e-01   2.36561e+01
2  f2           5.00000e-01   2.00000e-01   2.35352e-01  -4.09909e+01
ERR DEF= 0.5
MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
COVARIANCE MATRIX CALCULATED SUCCESSFULLY
FCN=24323.1 FROM MIGRAD    STATUS=CONVERGED      33 CALLS          34 TOTAL
EDM=1.04489e-07    STRATEGY= 1      ERROR MATRIX ACCURATE
EXT PARAMETER                                   STEP         FIRST
NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE
1  f1           9.85261e-01   4.18474e-01   3.64059e-03  -9.57093e-03
2  f2           5.07031e-01   2.19792e-01   2.15145e-03   1.61567e-02
ERR DEF= 0.5
EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  2    ERR DEF=0.5
1.864e-01  9.565e-02
9.565e-02  4.937e-02
PARAMETER  CORRELATION COEFFICIENTS
NO.  GLOBAL      1      2
1  0.99694   1.000  0.997
2  0.99694   0.997  1.000
**********
**    7 **SET ERR         0.5
**********
**********
**    8 **SET PRINT           1
**********
**********
**    9 **HESSE        1000
**********
COVARIANCE MATRIX CALCULATED SUCCESSFULLY
FCN=24323.1 FROM HESSE     STATUS=OK             10 CALLS          44 TOTAL
EDM=1.04101e-07    STRATEGY= 1      ERROR MATRIX ACCURATE
EXT PARAMETER                                INTERNAL      INTERNAL
NO.   NAME      VALUE            ERROR       STEP SIZE       VALUE
1  f1           9.85261e-01   1.40907e+00   1.45624e-04  -1.47397e-02
2  f2           5.07031e-01   1.24234e+00   8.60579e-05  -5.15499e-01
ERR DEF= 0.5
EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  2    ERR DEF=0.5
4.590e+00  2.362e+00
2.362e+00  1.215e+00
PARAMETER  CORRELATION COEFFICIENTS
NO.  GLOBAL      1      2
1  0.99988   1.000  1.000
2  0.99988   1.000  1.000
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization


## Plot amplitude sum pdf¶

Make 2D plots of amplitudes

In [10]:
hh_cos = ampl1.createHistogram("hh_cos", t, Binning=50, YVar=dict(var=cosa, Binning=50))
hh_sin = ampl2.createHistogram("hh_sin", t, Binning=50, YVar=dict(var=cosa, Binning=50))
hh_cos.SetLineColor(ROOT.kBlue)
hh_sin.SetLineColor(ROOT.kRed)

[#0] WARNING:InputArguments -- RooAbsReal::createHistogram(ampl1) WARNING extended mode requested for a non-pdf object, ignored
[#0] WARNING:InputArguments -- RooAbsReal::createHistogram(ampl2) WARNING extended mode requested for a non-pdf object, ignored


Make projection on t, data, and its components Note component projections may be larger than sum because amplitudes can be negative

In [11]:
frame1 = t.frame()
data.plotOn(frame1)
pdf.plotOn(frame1)
pdf.plotOn(frame1, Components=ampl1, LineStyle="--")
pdf.plotOn(frame1, Components=ampl2, LineStyle="--", LineColor="r")

Out[11]:
<cppyy.gbl.RooPlot object at 0x9438ef0>
[#1] INFO:Plotting -- RooAbsReal::plotOn(pdf) plot on t integrates over variables (cosa)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_coshGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_sinhGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(pdf) directly selected PDF components: (ampl1)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(pdf) indirectly selected PDF components: (poly1,tm_conv_coshGBasis_[t],coshGBasis)
[#1] INFO:Plotting -- RooAbsReal::plotOn(pdf) plot on t integrates over variables (cosa)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_coshGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_sinhGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(pdf) directly selected PDF components: (ampl2)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(pdf) indirectly selected PDF components: (poly2,tm_conv_sinhGBasis_[t],sinhGBasis)
[#1] INFO:Plotting -- RooAbsReal::plotOn(pdf) plot on t integrates over variables (cosa)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_coshGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_sinhGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)


Make projection on cosa, data, and its components Note that components projection may be larger than sum because amplitudes can be negative

In [12]:
frame2 = cosa.frame()
data.plotOn(frame2)
pdf.plotOn(frame2)
pdf.plotOn(frame2, Components=ampl1, LineStyle="--")
pdf.plotOn(frame2, Components=ampl2, LineStyle="--", LineColor="r")

c = ROOT.TCanvas("rf704_amplitudefit", "rf704_amplitudefit", 800, 800)
c.Divide(2, 2)
c.cd(1)
frame1.GetYaxis().SetTitleOffset(1.4)
frame1.Draw()
c.cd(2)
frame2.GetYaxis().SetTitleOffset(1.4)
frame2.Draw()
c.cd(3)
hh_cos.GetZaxis().SetTitleOffset(2.3)
hh_cos.Draw("surf")
c.cd(4)
hh_sin.GetZaxis().SetTitleOffset(2.3)
hh_sin.Draw("surf")

c.SaveAs("rf704_amplitudefit.png")

[#1] INFO:Plotting -- RooAbsReal::plotOn(pdf) plot on cosa integrates over variables (t)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_coshGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_sinhGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_coshGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_sinhGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(pdf) directly selected PDF components: (ampl1)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(pdf) indirectly selected PDF components: (poly1,tm_conv_coshGBasis_[t],coshGBasis)
[#1] INFO:Plotting -- RooAbsReal::plotOn(pdf) plot on cosa integrates over variables (t)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_coshGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_sinhGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_coshGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(pdf) directly selected PDF components: (ampl2)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(pdf) indirectly selected PDF components: (poly2,tm_conv_sinhGBasis_[t],sinhGBasis)
[#1] INFO:Plotting -- RooAbsReal::plotOn(pdf) plot on cosa integrates over variables (t)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_coshGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_sinhGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(tm_conv_sinhGBasis_[t]_Int[t]) using numeric integrator RooIntegrator1D to calculate Int(t)

Info in <TCanvas::Print>: png file rf704_amplitudefit.png has been created


Draw all canvases

In [13]:
from ROOT import gROOT
gROOT.GetListOfCanvases().Draw()