ROOT   Reference Guide
Searching...
No Matches
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.0, 15.0)
19cosa = ROOT.RooRealVar("cosa", "cos(alpha)", -1.0, 1.0)
20
21tau = ROOT.RooRealVar("tau", "#tau", 1.5)
22deltaGamma = ROOT.RooRealVar("deltaGamma", "deltaGamma", 0.3)
23coshG = ROOT.RooFormulaVar("coshGBasis", "exp(-@0/ @1)*cosh(@0*@2/2)", [t, tau, deltaGamma])
24sinhG = ROOT.RooFormulaVar("sinhGBasis", "exp(-@0/ @1)*sinh(@0*@2/2)", [t, tau, deltaGamma])
25
26# Construct polynomial amplitudes in cos(a)
27poly1 = ROOT.RooPolyVar("poly1", "poly1", cosa, [0.5, 0.2, 0.2], 0)
28poly2 = ROOT.RooPolyVar("poly2", "poly2", cosa, [1.0, -0.2, 3.0], 0)
29
30# Construct 2D amplitude as uncorrelated product of amp(t)*amp(cosa)
31ampl1 = ROOT.RooProduct("ampl1", "amplitude 1", [poly1, coshG])
32ampl2 = ROOT.RooProduct("ampl2", "amplitude 2", [poly2, sinhG])
33
34# Construct amplitude sum pdf
35# -----------------------------------------------------
36
37# Amplitude strengths
38f1 = ROOT.RooRealVar("f1", "f1", 1, 0, 2)
39f2 = ROOT.RooRealVar("f2", "f2", 0.5, 0, 2)
40
41# Construct pdf
42pdf = ROOT.RooRealSumPdf("pdf", "pdf", [ampl1, ampl2], [f1, f2])
43
44# Generate some toy data from pdf
45data = pdf.generate({t, cosa}, 10000)
46
47# Fit pdf to toy data with only amplitude strength floating
48pdf.fitTo(data)
49
50# Plot amplitude sum pdf
51# -------------------------------------------
52
53# Make 2D plots of amplitudes
54hh_cos = ampl1.createHistogram("hh_cos", t, Binning=50, YVar=dict(var=cosa, Binning=50))
55hh_sin = ampl2.createHistogram("hh_sin", t, Binning=50, YVar=dict(var=cosa, Binning=50))
56hh_cos.SetLineColor(ROOT.kBlue)
57hh_sin.SetLineColor(ROOT.kRed)
58
59# Make projection on t, data, and its components
60# Note component projections may be larger than sum because amplitudes can
61# be negative
62frame1 = t.frame()
63data.plotOn(frame1)
64pdf.plotOn(frame1)
65pdf.plotOn(frame1, Components=ampl1, LineStyle="--")
66pdf.plotOn(frame1, Components=ampl2, LineStyle="--", LineColor="r")
67
68# Make projection on cosa, data, and its components
69# Note that components projection may be larger than sum because
70# amplitudes can be negative
71frame2 = cosa.frame()
72data.plotOn(frame2)
73pdf.plotOn(frame2)
74pdf.plotOn(frame2, Components=ampl1, LineStyle="--")
75pdf.plotOn(frame2, Components=ampl2, LineStyle="--", LineColor="r")
76
77c = ROOT.TCanvas("rf704_amplitudefit", "rf704_amplitudefit", 800, 800)
78c.Divide(2, 2)
79c.cd(1)
81frame1.GetYaxis().SetTitleOffset(1.4)
82frame1.Draw()
83c.cd(2)