Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
rf801_mcstudy.py
Go to the documentation of this file.
1 ## \file
2 ## \ingroup tutorial_roofit
3 ## \notebook
4 ##
5 ## Validation and MC studies: toy Monte Carlo study that perform cycles of event generation and fitting
6 ##
7 ## \macro_code
8 ##
9 ## \date February 2018
10 ## \author Clemens Lange, Wouter Verkerke (C++ version)
11 
12 import ROOT
13 
14 
15 # Create model
16 # -----------------------
17 
18 # Declare observable x
19 x = ROOT.RooRealVar("x", "x", 0, 10)
20 x.setBins(40)
21 
22 # Create two Gaussian PDFs g1(x,mean1,sigma) anf g2(x,mean2,sigma) and
23 # their parameters
24 mean = ROOT.RooRealVar("mean", "mean of gaussians", 5, 0, 10)
25 sigma1 = ROOT.RooRealVar("sigma1", "width of gaussians", 0.5)
26 sigma2 = ROOT.RooRealVar("sigma2", "width of gaussians", 1)
27 
28 sig1 = ROOT.RooGaussian("sig1", "Signal component 1", x, mean, sigma1)
29 sig2 = ROOT.RooGaussian("sig2", "Signal component 2", x, mean, sigma2)
30 
31 # Build Chebychev polynomial p.d.f.
32 a0 = ROOT.RooRealVar("a0", "a0", 0.5, 0., 1.)
33 a1 = ROOT.RooRealVar("a1", "a1", -0.2, -1, 1.)
34 bkg = ROOT.RooChebychev("bkg", "Background", x, ROOT.RooArgList(a0, a1))
35 
36 # Sum the signal components into a composite signal p.d.f.
37 sig1frac = ROOT.RooRealVar(
38  "sig1frac", "fraction of component 1 in signal", 0.8, 0., 1.)
39 sig = ROOT.RooAddPdf(
40  "sig", "Signal", ROOT.RooArgList(sig1, sig2), ROOT.RooArgList(sig1frac))
41 
42 # Sum the composite signal and background
43 nbkg = ROOT.RooRealVar(
44  "nbkg", "number of background events, ", 150, 0, 1000)
45 nsig = ROOT.RooRealVar("nsig", "number of signal events", 150, 0, 1000)
46 model = ROOT.RooAddPdf(
47  "model", "g1+g2+a", ROOT.RooArgList(bkg, sig), ROOT.RooArgList(nbkg, nsig))
48 
49 # Create manager
50 # ---------------------------
51 
52 # Instantiate ROOT.RooMCStudy manager on model with x as observable and given choice of fit options
53 #
54 # The Silence() option kills all messages below the PROGRESS level, only a single message
55 # per sample executed, any error message that occur during fitting
56 #
57 # The Extended() option has two effects:
58 # 1) The extended ML term is included in the likelihood and
59 # 2) A poisson fluctuation is introduced on the number of generated events
60 #
61 # The FitOptions() given here are passed to the fitting stage of each toy experiment.
62 # If Save() is specified, fit result of each experiment is saved by the manager
63 #
64 # A Binned() option is added in self example to bin the data between generation and fitting
65 # to speed up the study at the expemse of some precision
66 
67 mcstudy = ROOT.RooMCStudy(
68  model,
69  ROOT.RooArgSet(x),
70  ROOT.RooFit.Binned(
71  ROOT.kTRUE),
72  ROOT.RooFit.Silence(),
73  ROOT.RooFit.Extended(),
74  ROOT.RooFit.FitOptions(
75  ROOT.RooFit.Save(
76  ROOT.kTRUE),
77  ROOT.RooFit.PrintEvalErrors(0)))
78 
79 # Generate and fit events
80 # ---------------------------------------------
81 
82 # Generate and fit 1000 samples of Poisson(nExpected) events
83 mcstudy.generateAndFit(1000)
84 
85 # Explore results of study
86 # ------------------------------------------------
87 
88 # Make plots of the distributions of mean, error on mean and the pull of
89 # mean
90 frame1 = mcstudy.plotParam(mean, ROOT.RooFit.Bins(40))
91 frame2 = mcstudy.plotError(mean, ROOT.RooFit.Bins(40))
92 frame3 = mcstudy.plotPull(mean, ROOT.RooFit.Bins(
93  40), ROOT.RooFit.FitGauss(ROOT.kTRUE))
94 
95 # Plot distribution of minimized likelihood
96 frame4 = mcstudy.plotNLL(ROOT.RooFit.Bins(40))
97 
98 # Make some histograms from the parameter dataset
99 hh_cor_a0_s1f = ROOT.RooAbsData.createHistogram(
100  mcstudy.fitParDataSet(), "hh", a1, ROOT.RooFit.YVar(sig1frac))
101 hh_cor_a0_a1 = ROOT.RooAbsData.createHistogram(mcstudy.fitParDataSet(),
102  "hh", a0, ROOT.RooFit.YVar(a1))
103 
104 # Access some of the saved fit results from individual toys
105 corrHist000 = mcstudy.fitResult(0).correlationHist("c000")
106 corrHist127 = mcstudy.fitResult(127).correlationHist("c127")
107 corrHist953 = mcstudy.fitResult(953).correlationHist("c953")
108 
109 # Draw all plots on a canvas
110 ROOT.gStyle.SetPalette(1)
111 ROOT.gStyle.SetOptStat(0)
112 c = ROOT.TCanvas("rf801_mcstudy", "rf801_mcstudy", 900, 900)
113 c.Divide(3, 3)
114 c.cd(1)
115 ROOT.gPad.SetLeftMargin(0.15)
116 frame1.GetYaxis().SetTitleOffset(1.4)
117 frame1.Draw()
118 c.cd(2)
119 ROOT.gPad.SetLeftMargin(0.15)
120 frame2.GetYaxis().SetTitleOffset(1.4)
121 frame2.Draw()
122 c.cd(3)
123 ROOT.gPad.SetLeftMargin(0.15)
124 frame3.GetYaxis().SetTitleOffset(1.4)
125 frame3.Draw()
126 c.cd(4)
127 ROOT.gPad.SetLeftMargin(0.15)
128 frame4.GetYaxis().SetTitleOffset(1.4)
129 frame4.Draw()
130 c.cd(5)
131 ROOT.gPad.SetLeftMargin(0.15)
132 hh_cor_a0_s1f.GetYaxis().SetTitleOffset(1.4)
133 hh_cor_a0_s1f.Draw("box")
134 c.cd(6)
135 ROOT.gPad.SetLeftMargin(0.15)
136 hh_cor_a0_a1.GetYaxis().SetTitleOffset(1.4)
137 hh_cor_a0_a1.Draw("box")
138 c.cd(7)
139 ROOT.gPad.SetLeftMargin(0.15)
140 corrHist000.GetYaxis().SetTitleOffset(1.4)
141 corrHist000.Draw("colz")
142 c.cd(8)
143 ROOT.gPad.SetLeftMargin(0.15)
144 corrHist127.GetYaxis().SetTitleOffset(1.4)
145 corrHist127.Draw("colz")
146 c.cd(9)
147 ROOT.gPad.SetLeftMargin(0.15)
148 corrHist953.GetYaxis().SetTitleOffset(1.4)
149 corrHist953.Draw("colz")
150 
151 c.SaveAs("rf801_mcstudy.png")
152 
153 # Make ROOT.RooMCStudy object available on command line after
154 # macro finishes
155 ROOT.gDirectory.Add(mcstudy)