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 
11 import ROOT
12 
13 
14 # Setup 2D amplitude functions
15 # -------------------------------------------------------
16 
17 # Observables
18 t = ROOT.RooRealVar("t", "time", -1., 15.)
19 cosa = ROOT.RooRealVar("cosa", "cos(alpha)", -1., 1.)
20 
21 # Use ROOT.RooTruthModel to obtain compiled implementation of sinh/cosh
22 # modulated decay functions
23 tau = ROOT.RooRealVar("tau", "#tau", 1.5)
24 deltaGamma = ROOT.RooRealVar("deltaGamma", "deltaGamma", 0.3)
25 tm = ROOT.RooTruthModel("tm", "tm", t)
26 coshGBasis = ROOT.RooFormulaVar(
27  "coshGBasis",
28  "exp(-@0/ @1)*cosh(@0*@2/2)",
29  ROOT.RooArgList(
30  t,
31  tau,
32  deltaGamma))
33 sinhGBasis = ROOT.RooFormulaVar(
34  "sinhGBasis",
35  "exp(-@0/ @1)*sinh(@0*@2/2)",
36  ROOT.RooArgList(
37  t,
38  tau,
39  deltaGamma))
40 coshGConv = tm.convolution(coshGBasis, t)
41 sinhGConv = tm.convolution(sinhGBasis, t)
42 
43 # Construct polynomial amplitudes in cos(a)
44 poly1 = 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)
53 poly2 = 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)
57 ampl1 = ROOT.RooProduct("ampl1", "amplitude 1",
58  ROOT.RooArgList(poly1, coshGConv))
59 ampl2 = ROOT.RooProduct("ampl2", "amplitude 2",
60  ROOT.RooArgList(poly2, sinhGConv))
61 
62 # Construct amplitude sum pdf
63 # -----------------------------------------------------
64 
65 # Amplitude strengths
66 f1 = ROOT.RooRealVar("f1", "f1", 1, 0, 2)
67 f2 = ROOT.RooRealVar("f2", "f2", 0.5, 0, 2)
68 
69 # Construct pdf
70 pdf = ROOT.RooRealSumPdf("pdf", "pdf", ROOT.RooArgList(
71  ampl1, ampl2), ROOT.RooArgList(f1, f2))
72 
73 # Generate some toy data from pdf
74 data = pdf.generate(ROOT.RooArgSet(t, cosa), 10000)
75 
76 # Fit pdf to toy data with only amplitude strength floating
77 pdf.fitTo(data)
78 
79 # Plot amplitude sum pdf
80 # -------------------------------------------
81 
82 # Make 2D plots of amplitudes
83 hh_cos = ampl1.createHistogram("hh_cos", t, ROOT.RooFit.Binning(
84  50), ROOT.RooFit.YVar(cosa, ROOT.RooFit.Binning(50)))
85 hh_sin = ampl2.createHistogram("hh_sin", t, ROOT.RooFit.Binning(
86  50), ROOT.RooFit.YVar(cosa, ROOT.RooFit.Binning(50)))
87 hh_cos.SetLineColor(ROOT.kBlue)
88 hh_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
93 frame1 = t.frame()
94 data.plotOn(frame1)
95 pdf.plotOn(frame1)
96 # workaround, see https://root.cern.ch/phpBB3/viewtopic.php?t=7764
97 ras_ampl1 = ROOT.RooArgSet(ampl1)
98 pdf.plotOn(frame1, ROOT.RooFit.Components(ras_ampl1),
99  ROOT.RooFit.LineStyle(ROOT.kDashed))
100 ras_ampl2 = ROOT.RooArgSet(ampl2)
101 pdf.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
107 frame2 = cosa.frame()
108 data.plotOn(frame2)
109 pdf.plotOn(frame2)
110 pdf.plotOn(frame2, ROOT.RooFit.Components(ras_ampl1),
111  ROOT.RooFit.LineStyle(ROOT.kDashed))
112 pdf.plotOn(frame2, ROOT.RooFit.Components(ras_ampl2), ROOT.RooFit.LineStyle(
113  ROOT.kDashed), ROOT.RooFit.LineColor(ROOT.kRed))
114 
115 c = ROOT.TCanvas("rf704_amplitudefit", "rf704_amplitudefit", 800, 800)
116 c.Divide(2, 2)
117 c.cd(1)
118 ROOT.gPad.SetLeftMargin(0.15)
119 frame1.GetYaxis().SetTitleOffset(1.4)
120 frame1.Draw()
121 c.cd(2)
122 ROOT.gPad.SetLeftMargin(0.15)
123 frame2.GetYaxis().SetTitleOffset(1.4)
124 frame2.Draw()
125 c.cd(3)
126 ROOT.gPad.SetLeftMargin(0.20)
127 hh_cos.GetZaxis().SetTitleOffset(2.3)
128 hh_cos.Draw("surf")
129 c.cd(4)
130 ROOT.gPad.SetLeftMargin(0.20)
131 hh_sin.GetZaxis().SetTitleOffset(2.3)
132 hh_sin.Draw("surf")
133 
134 c.SaveAs("rf704_amplitudefit.png")