ROOT   Reference Guide
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_image
8## \macro_code
9## \macro_output
10##
11## \date February 2018
12## \authors Clemens Lange, Wouter Verkerke (C version)
13
14import ROOT
15from array import array
16
17
18def makeTH1():
19
20 # Create ROOT ROOT.TH1 filled with a Gaussian distribution
21
22 hh = ROOT.TH1D("hh", "hh", 25, -10, 10)
23 for i in range(100):
24 hh.Fill(ROOT.gRandom.Gaus(0, 3))
25 return hh
26
27
28def makeTTree():
29 # Create ROOT ROOT.TTree filled with a Gaussian distribution in x and a
30 # uniform distribution in y
31
32 tree = ROOT.TTree("tree", "tree")
33 px = array("d", [0])
34 py = array("d", [0])
35 tree.Branch("x", px, "x/D")
36 tree.Branch("y", py, "y/D")
37 for i in range(100):
38 px[0] = ROOT.gRandom.Gaus(0, 3)
39 py[0] = ROOT.gRandom.Uniform() * 30 - 15
40 tree.Fill()
41 return tree
42
43
44############################
45# Importing ROOT histograms
46############################
47# Import ROOT TH1 into a RooDataHist
48# ---------------------------------------------------------
49# Create a ROOT TH1 histogram
50hh = makeTH1()
51
52# Declare observable x
53x = ROOT.RooRealVar("x", "x", -10, 10)
54
55# Create a binned dataset that imports contents of ROOT.TH1 and associates
56# its contents to observable 'x'
57dh = ROOT.RooDataHist("dh", "dh", [x], Import=hh)
58
59# Plot and fit a RooDataHist
60# ---------------------------------------------------
61# Make plot of binned dataset showing Poisson error bars (RooFit 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, PrintLevel=-1)
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}, 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")
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
132
133print("\nOriginal data, line 20:")
134ds.get(20).Print("V")
135
136print("\nRead-back data, line 20:")
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)
162frame.GetYaxis().SetTitleOffset(1.4)
163frame.Draw()
164c.cd(2)
166frame2.GetYaxis().SetTitleOffset(1.4)
167frame2.Draw()
168c.cd(4)