Logo ROOT  
Reference Guide
rf704_amplitudefit.py
Go to the documentation of this file.
1## \file
2## \ingroup tutorial_roofit
3## \notebook
4## Special pdf's: using a pdf defined by a sum of real-valued amplitude components
5##
6## \macro_code
7##
8## \date February 2018
9## \authors Clemens Lange, Wouter Verkerke (C++ version)
10
11import ROOT
12
13
14# Setup 2D amplitude functions
15# -------------------------------------------------------
16
17# Observables
18t = ROOT.RooRealVar("t", "time", -1., 15.)
19cosa = ROOT.RooRealVar("cosa", "cos(alpha)", -1., 1.)
20
21# Use ROOT.RooTruthModel to obtain compiled implementation of sinh/cosh
22# modulated decay functions
23tau = ROOT.RooRealVar("tau", "#tau", 1.5)
24deltaGamma = ROOT.RooRealVar("deltaGamma", "deltaGamma", 0.3)
25tm = ROOT.RooTruthModel("tm", "tm", t)
26coshGBasis = ROOT.RooFormulaVar(
27 "coshGBasis",
28 "exp(-@0/ @1)*cosh(@0*@2/2)",
29 ROOT.RooArgList(
30 t,
31 tau,
32 deltaGamma))
33sinhGBasis = ROOT.RooFormulaVar(
34 "sinhGBasis",
35 "exp(-@0/ @1)*sinh(@0*@2/2)",
36 ROOT.RooArgList(
37 t,
38 tau,
39 deltaGamma))
40coshGConv = tm.convolution(coshGBasis, t)
41sinhGConv = tm.convolution(sinhGBasis, t)
42
43# Construct polynomial amplitudes in cos(a)
44poly1 = ROOT.RooPolyVar(
45 "poly1",
46 "poly1",
47 cosa,
48 ROOT.RooArgList(
49 ROOT.RooFit.RooConst(0.5),
50 ROOT.RooFit.RooConst(0.2),
51 ROOT.RooFit.RooConst(0.2)),
52 0)
53poly2 = ROOT.RooPolyVar("poly2", "poly2", cosa, ROOT.RooArgList(
54 ROOT.RooFit.RooConst(1), ROOT.RooFit.RooConst(-0.2), ROOT.RooFit.RooConst(3)), 0)
55
56# Construct 2D amplitude as uncorrelated product of amp(t)*amp(cosa)
57ampl1 = ROOT.RooProduct("ampl1", "amplitude 1",
58 ROOT.RooArgList(poly1, coshGConv))
59ampl2 = ROOT.RooProduct("ampl2", "amplitude 2",
60 ROOT.RooArgList(poly2, sinhGConv))
61
62# Construct amplitude sum pdf
63# -----------------------------------------------------
64
65# Amplitude strengths
66f1 = ROOT.RooRealVar("f1", "f1", 1, 0, 2)
67f2 = ROOT.RooRealVar("f2", "f2", 0.5, 0, 2)
68
69# Construct pdf
70pdf = ROOT.RooRealSumPdf("pdf", "pdf", ROOT.RooArgList(
71 ampl1, ampl2), ROOT.RooArgList(f1, f2))
72
73# Generate some toy data from pdf
74data = pdf.generate(ROOT.RooArgSet(t, cosa), 10000)
75
76# Fit pdf to toy data with only amplitude strength floating
77pdf.fitTo(data)
78
79# Plot amplitude sum pdf
80# -------------------------------------------
81
82# Make 2D plots of amplitudes
83hh_cos = ampl1.createHistogram("hh_cos", t, ROOT.RooFit.Binning(
84 50), ROOT.RooFit.YVar(cosa, ROOT.RooFit.Binning(50)))
85hh_sin = ampl2.createHistogram("hh_sin", t, ROOT.RooFit.Binning(
86 50), ROOT.RooFit.YVar(cosa, ROOT.RooFit.Binning(50)))
87hh_cos.SetLineColor(ROOT.kBlue)
88hh_sin.SetLineColor(ROOT.kRed)
89
90# Make projection on t, data, and its components
91# Note component projections may be larger than sum because amplitudes can
92# be negative
93frame1 = t.frame()
94data.plotOn(frame1)
95pdf.plotOn(frame1)
96# workaround, see https://root.cern.ch/phpBB3/viewtopic.php?t=7764
97ras_ampl1 = ROOT.RooArgSet(ampl1)
98pdf.plotOn(frame1, ROOT.RooFit.Components(ras_ampl1),
99 ROOT.RooFit.LineStyle(ROOT.kDashed))
100ras_ampl2 = ROOT.RooArgSet(ampl2)
101pdf.plotOn(frame1, ROOT.RooFit.Components(ras_ampl2), ROOT.RooFit.LineStyle(
102 ROOT.kDashed), ROOT.RooFit.LineColor(ROOT.kRed))
103
104# Make projection on cosa, data, and its components
105# Note that components projection may be larger than sum because
106# amplitudes can be negative
107frame2 = cosa.frame()
108data.plotOn(frame2)
109pdf.plotOn(frame2)
110pdf.plotOn(frame2, ROOT.RooFit.Components(ras_ampl1),
111 ROOT.RooFit.LineStyle(ROOT.kDashed))
112pdf.plotOn(frame2, ROOT.RooFit.Components(ras_ampl2), ROOT.RooFit.LineStyle(
113 ROOT.kDashed), ROOT.RooFit.LineColor(ROOT.kRed))
114
115c = ROOT.TCanvas("rf704_amplitudefit", "rf704_amplitudefit", 800, 800)
116c.Divide(2, 2)
117c.cd(1)
118ROOT.gPad.SetLeftMargin(0.15)
119frame1.GetYaxis().SetTitleOffset(1.4)
120frame1.Draw()
121c.cd(2)
122ROOT.gPad.SetLeftMargin(0.15)
123frame2.GetYaxis().SetTitleOffset(1.4)
124frame2.Draw()
125c.cd(3)
126ROOT.gPad.SetLeftMargin(0.20)
127hh_cos.GetZaxis().SetTitleOffset(2.3)
128hh_cos.Draw("surf")
129c.cd(4)
130ROOT.gPad.SetLeftMargin(0.20)
131hh_sin.GetZaxis().SetTitleOffset(2.3)
132hh_sin.Draw("surf")
133
134c.SaveAs("rf704_amplitudefit.png")