Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
rf102_dataimport.py
Go to the documentation of this file.
1## \file
2## \ingroup tutorial_roofit
3## \notebook
4## 'BASIC FUNCTIONALITY' RooFit tutorial macro #102
5## Importing data from ROOT TTrees and THx histograms
6##
7## \macro_code
8##
9## \date February 2018
10## \author Clemens Lange
11## \author Wouter Verkerke (C version)
12
13import ROOT
14from array import array
15
16
17def makeTH1():
18
19 # Create ROOT ROOT.TH1 filled with a Gaussian distribution
20
21 hh = ROOT.TH1D("hh", "hh", 25, -10, 10)
22 for i in range(100):
23 hh.Fill(ROOT.gRandom.Gaus(0, 3))
24 return hh
25
26
27def makeTTree():
28 # Create ROOT ROOT.TTree filled with a Gaussian distribution in x and a
29 # uniform distribution in y
30
31 tree = ROOT.TTree("tree", "tree")
32 px = array("d", [0])
33 py = array("d", [0])
34 tree.Branch("x", px, "x/D")
35 tree.Branch("y", py, "y/D")
36 for i in range(100):
37 px[0] = ROOT.gRandom.Gaus(0, 3)
38 py[0] = ROOT.gRandom.Uniform() * 30 - 15
39 tree.Fill()
40 return tree
41
42
43############################
44# Importing ROOT histograms
45############################
46# Import ROOT TH1 into a RooDataHist
47# ---------------------------------------------------------
48# Create a ROOT TH1 histogram
49hh = makeTH1()
50
51# Declare observable x
52x = ROOT.RooRealVar("x", "x", -10, 10)
53
54# Create a binned dataset that imports contents of ROOT.TH1 and associates
55# its contents to observable 'x'
56dh = ROOT.RooDataHist("dh", "dh", [x], Import=hh)
57
58# Plot and fit a RooDataHist
59# ---------------------------------------------------
60# Make plot of binned dataset showing Poisson error bars (ROOT.RooFit
61# default)
62frame = x.frame(Title="Imported ROOT.TH1 with Poisson error bars")
63dh.plotOn(frame)
64
65# Fit a Gaussian p.d.f to the data
66mean = ROOT.RooRealVar("mean", "mean", 0, -10, 10)
67sigma = ROOT.RooRealVar("sigma", "sigma", 3, 0.1, 10)
68gauss = ROOT.RooGaussian("gauss", "gauss", x, mean, sigma)
69gauss.fitTo(dh)
70gauss.plotOn(frame)
71
72# Plot and fit a RooDataHist with internal errors
73# ---------------------------------------------------------------------------------------------
74
75# If histogram has custom error (i.e. its contents is does not originate from a Poisson process
76# but e.g. is a sum of weighted events) you can data with symmetric 'sum-of-weights' error instead
77# (same error bars as shown by ROOT)
78frame2 = x.frame(Title="Imported ROOT.TH1 with internal errors")
79dh.plotOn(frame2, DataError="SumW2")
80gauss.plotOn(frame2)
81
82# Please note that error bars shown (Poisson or SumW2) are for visualization only, the are NOT used
83# in a maximum likelihood fit
84#
85# A (binned) ML fit will ALWAYS assume the Poisson error interpretation of data (the mathematical definition
86# of likelihood does not take any external definition of errors). Data with non-unit weights can only be correctly
87# fitted with a chi^2 fit (see rf602_chi2fit.py)
88#
89# Importing ROOT TTrees
90# -----------------------------------------------------------
91# Import ROOT TTree into a RooDataSet
92
93tree = makeTTree()
94
95# Define 2nd observable y
96y = ROOT.RooRealVar("y", "y", -10, 10)
97
98# Construct unbinned dataset importing tree branches x and y matching between branches and ROOT.RooRealVars
99# is done by name of the branch/RRV
100#
101# Note that ONLY entries for which x,y have values within their allowed ranges as defined in
102# ROOT.RooRealVar x and y are imported. Since the y values in the import tree are in the range [-15,15]
103# and RRV y defines a range [-10,10] this means that the ROOT.RooDataSet
104# below will have less entries than the ROOT.TTree 'tree'
105
106ds = ROOT.RooDataSet("ds", "ds", {x, y}, ROOT.RooFit.Import(tree))
107
108# Use ascii import/export for datasets
109# ------------------------------------------------------------------------------------
110
111
112def write_dataset(ds, filename):
113 # Write data to output stream
114 outstream = ROOT.std.ofstream(filename)
115 # Optionally, adjust the stream here (e.g. std::setprecision)
116 ds.write(outstream)
117 outstream.close()
118
119
120write_dataset(ds, "rf102_testData.txt")
121
122# Read data from input stream. The variables of the dataset need to be supplied
123# to the RooDataSet::read() function.
124print("\n-----------------------\nReading data from ASCII")
125dataReadBack = ROOT.RooDataSet.read(
126 "rf102_testData.txt",
127 [x, y], # variables to be read. If the file has more fields, these are ignored.
128 "D", # Prints if a RooFit message stream listens for debug messages. Use Q for quiet.
129)
130
131dataReadBack.Print("V")
132
133print("\nOriginal data, line 20:")
134ds.get(20).Print("V")
135
136print("\nRead-back data, line 20:")
137dataReadBack.get(20).Print("V")
138
139
140# Plot data set with multiple binning choices
141# ------------------------------------------------------------------------------------
142# Print number of events in dataset
143ds.Print()
144
145# Print unbinned dataset with default frame binning (100 bins)
146frame3 = y.frame(Title="Unbinned data shown in default frame binning")
147ds.plotOn(frame3)
148
149# Print unbinned dataset with custom binning choice (20 bins)
150frame4 = y.frame(Title="Unbinned data shown with custom binning")
151ds.plotOn(frame4, Binning=20)
152
153frame5 = y.frame(Title="Unbinned data read back from ASCII file")
154ds.plotOn(frame5, Binning=20)
155dataReadBack.plotOn(frame5, Binning=20, MarkerColor="r", MarkerStyle=5)
156
157# Draw all frames on a canvas
158c = ROOT.TCanvas("rf102_dataimport", "rf102_dataimport", 800, 800)
159c.Divide(3, 2)
160c.cd(1)
161ROOT.gPad.SetLeftMargin(0.15)
162frame.GetYaxis().SetTitleOffset(1.4)
163frame.Draw()
164c.cd(2)
165ROOT.gPad.SetLeftMargin(0.15)
166frame2.GetYaxis().SetTitleOffset(1.4)
167frame2.Draw()
168c.cd(4)
169ROOT.gPad.SetLeftMargin(0.15)
170frame3.GetYaxis().SetTitleOffset(1.4)
171frame3.Draw()
172c.cd(5)
173ROOT.gPad.SetLeftMargin(0.15)
174frame4.GetYaxis().SetTitleOffset(1.4)
175frame4.Draw()
176c.cd(6)
177ROOT.gPad.SetLeftMargin(0.15)
178frame4.GetYaxis().SetTitleOffset(1.4)
179frame5.Draw()
180
181c.SaveAs("rf102_dataimport.png")