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