Logo ROOT  
Reference Guide
df104.py
Go to the documentation of this file.
1## \file
2## \ingroup tutorial_rcanvas
3## The Higgs to two photons analysis from the ATLAS Open Data 2020 release, with RDataFrame.
4##
5## This tutorial is the Higgs to two photons analysis from the ATLAS Open Data release in 2020
6## (http://opendata.atlas.cern/release/2020/documentation/). The data was taken with the ATLAS detector
7## during 2016 at a center-of-mass energy of 13 TeV. Although the Higgs to two photons decay is very rare,
8## the contribution of the Higgs can be seen as a narrow peak around 125 GeV because of the excellent
9## reconstruction and identification efficiency of photons at the ATLAS experiment.
10##
11## The analysis is translated to a RDataFrame workflow processing 1.7 GB of simulated events and data.
12##
13## This macro is replica of tutorials/dataframe/df104_HiggsToTwoPhotons.py, but with usage of ROOT7 graphics
14## Run macro with python3 -i df104.py command to get interactive canvas
15##
16## \macro_image (rcanvas_js)
17## \macro_code
18##
19## \date 2021-06-15
20## \authors Stefan Wunsch (KIT, CERN) Sergey Linev (GSI)
21
22import ROOT
23import os
24from ROOT.Experimental import RCanvas, RText, RAttrText, RAttrFont, RAttrFill, RAttrLine, RLegend, RPadPos, RPadExtent, TObjectDrawable
25
26# Enable multi-threading
27ROOT.ROOT.EnableImplicitMT()
28
29# Create a ROOT dataframe for each dataset
30path = "root://eospublic.cern.ch//eos/opendata/atlas/OutreachDatasets/2020-01-22"
31df = {}
32df["data"] = ROOT.RDataFrame("mini", (os.path.join(path, "GamGam/Data/data_{}.GamGam.root".format(x)) for x in ("A", "B", "C", "D")))
33df["ggH"] = ROOT.RDataFrame("mini", os.path.join(path, "GamGam/MC/mc_343981.ggH125_gamgam.GamGam.root"))
34df["VBF"] = ROOT.RDataFrame("mini", os.path.join(path, "GamGam/MC/mc_345041.VBFH125_gamgam.GamGam.root"))
35processes = list(df.keys())
36
37# Apply scale factors and MC weight for simulated events and a weight of 1 for the data
38for p in ["ggH", "VBF"]:
39 df[p] = df[p].Define("weight",
40 "scaleFactor_PHOTON * scaleFactor_PhotonTRIGGER * scaleFactor_PILEUP * mcWeight");
41df["data"] = df["data"].Define("weight", "1.0")
42
43# Select the events for the analysis
44for p in processes:
45 # Apply preselection cut on photon trigger
46 df[p] = df[p].Filter("trigP")
47
48 # Find two good muons with tight ID, pt > 25 GeV and not in the transition region between barrel and encap
49 df[p] = df[p].Define("goodphotons", "photon_isTightID && (photon_pt > 25000) && (abs(photon_eta) < 2.37) && ((abs(photon_eta) < 1.37) || (abs(photon_eta) > 1.52))")\
50 .Filter("Sum(goodphotons) == 2")
51
52 # Take only isolated photons
53 df[p] = df[p].Filter("Sum(photon_ptcone30[goodphotons] / photon_pt[goodphotons] < 0.065) == 2")\
54 .Filter("Sum(photon_etcone20[goodphotons] / photon_pt[goodphotons] < 0.065) == 2")
55
56# Compile a function to compute the invariant mass of the diphoton system
57ROOT.gInterpreter.Declare(
58"""
59using Vec_t = const ROOT::VecOps::RVec<float>;
60float ComputeInvariantMass(Vec_t& pt, Vec_t& eta, Vec_t& phi, Vec_t& e) {
61 ROOT::Math::PtEtaPhiEVector p1(pt[0], eta[0], phi[0], e[0]);
62 ROOT::Math::PtEtaPhiEVector p2(pt[1], eta[1], phi[1], e[1]);
63 return (p1 + p2).mass() / 1000.0;
64}
65""")
66
67# Define a new column with the invariant mass and perform final event selection
68hists = {}
69for p in processes:
70 # Make four vectors and compute invariant mass
71 df[p] = df[p].Define("m_yy", "ComputeInvariantMass(photon_pt[goodphotons], photon_eta[goodphotons], photon_phi[goodphotons], photon_E[goodphotons])")
72
73 # Make additional kinematic cuts and select mass window
74 df[p] = df[p].Filter("photon_pt[goodphotons][0] / 1000.0 / m_yy > 0.35")\
75 .Filter("photon_pt[goodphotons][1] / 1000.0 / m_yy > 0.25")\
76 .Filter("m_yy > 105 && m_yy < 160")
77
78 # Book histogram of the invariant mass with this selection
79 hists[p] = df[p].Histo1D(
80 ROOT.RDF.TH1DModel(p, "Diphoton invariant mass; m_{#gamma#gamma} [GeV];Events", 30, 105, 160),
81 "m_yy", "weight")
82
83# Run the event loop
84
85# RunGraphs allows to run the event loops of the separate RDataFrame graphs
86# concurrently. This results in an improved usage of the available resources
87# if each separate RDataFrame can not utilize all available resources, e.g.,
88# because not enough data is available.
89ROOT.RDF.RunGraphs([hists[s] for s in ["ggH", "VBF", "data"]])
90
91ggh = hists["ggH"].GetValue()
92vbf = hists["VBF"].GetValue()
93data = hists["data"].GetValue()
94
95# Create the plot
96
97# Set styles - not yet available for v7
98# ROOT.gROOT.SetStyle("ATLAS")
99
100# Create canvas with pads for main plot and data/MC ratio
101c = RCanvas.Create("df104_HiggsToTwoPhotons")
102
103lower_pad = c.AddPad(RPadPos(0,0.65), RPadExtent(1, 0.35))
104upper_pad = c.AddPad(RPadPos(0,0), RPadExtent(1, 0.65))
105
106upper_frame = upper_pad.AddFrame()
107upper_frame.margins.bottom = 0
108upper_frame.margins.left = 0.14
109upper_frame.margins.right = 0.05
110upper_frame.x.labels.hide = True
111
112lower_frame = lower_pad.AddFrame()
113lower_frame.margins.top = 0
114lower_frame.margins.left = 0.14
115lower_frame.margins.right = 0.05
116lower_frame.margins.bottom = 0.3
117
118# Fit signal + background model to data
119fit = ROOT.TF1("fit", "([0]+[1]*x+[2]*x^2+[3]*x^3)+[4]*exp(-0.5*((x-[5])/[6])^2)", 105, 160)
120fit.FixParameter(5, 125.0)
121fit.FixParameter(4, 119.1)
122fit.FixParameter(6, 2.39)
123fit.SetLineColor(2)
124fit.SetLineStyle(1)
125fit.SetLineWidth(2)
126# do not draw fit function
127data.Fit("fit", "0", "", 105, 160)
128
129# Draw data
130data.SetMarkerStyle(20)
131data.SetMarkerSize(1.2)
132data.SetLineWidth(2)
133data.SetLineColor(ROOT.kBlack)
134data.SetMinimum(1e-3)
135data.SetMaximum(8e3)
136data.GetYaxis().SetLabelSize(0.045)
137data.GetYaxis().SetTitleSize(0.05)
138data.SetStats(0)
139data.SetTitle("")
140
141data_drawable = upper_pad.Add[TObjectDrawable]()
142data_drawable.Set(data, "E")
143
144# Draw fit
145fit_drawable = upper_pad.Add[TObjectDrawable]()
146fit_drawable.Set(fit, "SAME")
147
148# Draw background
149bkg = ROOT.TF1("bkg", "([0]+[1]*x+[2]*x^2+[3]*x^3)", 105, 160)
150for i in range(4):
151 bkg.SetParameter(i, fit.GetParameter(i))
152bkg.SetLineColor(4)
153bkg.SetLineStyle(2)
154bkg.SetLineWidth(2)
155bkg_drawable = upper_pad.Add[TObjectDrawable]()
156bkg_drawable.Set(bkg, "SAME")
157
158# Scale simulated events with luminosity * cross-section / sum of weights
159# and merge to single Higgs signal
160lumi = 10064.0
161ggh.Scale(lumi * 0.102 / 55922617.6297)
162vbf.Scale(lumi * 0.008518764 / 3441426.13711)
163higgs = ggh.Clone()
164higgs.Add(vbf)
165higgs_drawable = upper_pad.Add[TObjectDrawable]()
166higgs_drawable.Set(higgs, "HIST SAME")
167
168# Draw ratio
169
170ratiobkg = ROOT.TH1I("zero", "", 10, 105, 160)
171ratiobkg.SetLineColor(4)
172ratiobkg.SetLineStyle(2)
173ratiobkg.SetLineWidth(2)
174ratiobkg.SetMinimum(-125)
175ratiobkg.SetMaximum(250)
176ratiobkg.GetXaxis().SetLabelSize(0.08)
177ratiobkg.GetXaxis().SetTitleSize(0.12)
178ratiobkg.GetXaxis().SetTitleOffset(1.0)
179ratiobkg.GetYaxis().SetLabelSize(0.08)
180ratiobkg.GetYaxis().SetTitleSize(0.09)
181ratiobkg.GetYaxis().SetTitle("Data - Bkg.")
182ratiobkg.GetYaxis().CenterTitle()
183ratiobkg.GetYaxis().SetNdivisions(503, False)
184ratiobkg.GetYaxis().ChangeLabel(-1, -1, 0)
185ratiobkg.GetXaxis().SetTitle("m_{#gamma#gamma} [GeV]")
186lower_pad.Add[TObjectDrawable]().Set(ratiobkg, "AXIS")
187
188ratiosig = ROOT.TH1F("ratiosig", "ratiosig", 5500, 105, 160)
189ratiosig.Eval(fit)
190ratiosig.SetLineColor(2)
191ratiosig.SetLineStyle(1)
192ratiosig.SetLineWidth(2)
193ratiosig.Add(bkg, -1)
194lower_pad.Add[TObjectDrawable]().Set(ratiosig, "SAME")
195
196ratiodata = data.Clone()
197ratiodata.Add(bkg, -1)
198for i in range(1, data.GetNbinsX()):
199 ratiodata.SetBinError(i, data.GetBinError(i))
200
201lower_pad.Add[TObjectDrawable]().Set(ratiodata, "E SAME")
202
203# Add RLegend
204legend = upper_pad.Draw[RLegend](RPadPos(0.01, 0.01), RPadExtent(0.3, 0.4))
205legend.text.size = 0.05
206legend.text.align = RAttrText.kRightCenter
207legend.text.font = RAttrFont.kArial
208legend.border.style = RAttrLine.kNone
209legend.border.width = 0
210legend.fill.style = RAttrFill.kHollow
211
212legend.AddEntry(data_drawable, "Data", "lme")
213legend.AddEntry(bkg_drawable, "Background", "l")
214legend.AddEntry(fit_drawable, "Signal + Bkg.", "l")
215legend.AddEntry(higgs_drawable, "Signal", "l")
216
217# Add ATLAS labels
218lbl1 = upper_pad.Draw[RText](RPadPos(0.05, 0.88), "ATLAS")
219lbl1.onFrame = True
220lbl1.text.font = RAttrFont.kArialBoldOblique
221lbl1.text.size = 0.04
222lbl1.text.align = RAttrText.kLeftBottom
223lbl2 = upper_pad.Draw[RText](RPadPos(0.05 + 0.16, 0.88), "Open Data")
224lbl2.onFrame = True
225lbl2.text.font = RAttrFont.kArial
226lbl2.text.size = 0.04
227lbl2.text.align = RAttrText.kLeftBottom
228lbl3 = upper_pad.Draw[RText](RPadPos(0.05, 0.82), "#sqrt{s} = 13 TeV, 10 fb^{-1}")
229lbl3.onFrame = True
230lbl3.text.font = RAttrFont.kArial
231lbl3.text.size = 0.03
232lbl3.text.align = RAttrText.kLeftBottom
233
234# show canvas finally
235c.SetSize(700, 780)
236c.Show()
237
238# Save plot in PNG file
239c.SaveAs("df104.png")
240print("Saved figure to df104.png")
ROOT's RDataFrame offers a high level interface for analyses of data stored in TTree,...
Definition: RDataFrame.hxx:40
RVec< T > Filter(const RVec< T > &v, F &&f)
Create a new collection with the elements passing the filter expressed by the predicate.
Definition: RVec.hxx:2052
void RunGraphs(std::vector< RResultHandle > handles)
Trigger the event loop of multiple RDataFrames concurrently.
Definition: RDFHelpers.cxx:23
A struct which stores the parameters of a TH1D.
Definition: HistoModels.hxx:27