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