Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
df104_HiggsToTwoPhotons.py
Go to the documentation of this file.
1## \file
2## \ingroup tutorial_dataframe
3## \notebook -draw
4## The Higgs to two photons analysis from the ATLAS Open Data 2020 release, with RDataFrame.
5##
6## This tutorial is the Higgs to two photons analysis from the ATLAS Open Data release in 2020
7## (http://opendata.atlas.cern/release/2020/documentation/). The data was taken with the ATLAS detector
8## during 2016 at a center-of-mass energy of 13 TeV. Although the Higgs to two photons decay is very rare,
9## the contribution of the Higgs can be seen as a narrow peak around 125 GeV because of the excellent
10## reconstruction and identification efficiency of photons at the ATLAS experiment.
11##
12## The analysis is translated to a RDataFrame workflow processing 1.7 GB of simulated events and data.
13##
14## \macro_image
15## \macro_code
16## \macro_output
17##
18## \date February 2020
19## \author Stefan Wunsch (KIT, CERN)
20
21import ROOT
22import os
23
24# Enable multi-threading
25ROOT.ROOT.EnableImplicitMT()
26
27# Create a ROOT dataframe for each dataset
28path = "root://eospublic.cern.ch//eos/opendata/atlas/OutreachDatasets/2020-01-22"
29df = {}
30df["data"] = ROOT.RDataFrame("mini", (os.path.join(path, "GamGam/Data/data_{}.GamGam.root".format(x)) for x in ("A", "B", "C", "D")))
31df["ggH"] = ROOT.RDataFrame("mini", os.path.join(path, "GamGam/MC/mc_343981.ggH125_gamgam.GamGam.root"))
32df["VBF"] = ROOT.RDataFrame("mini", os.path.join(path, "GamGam/MC/mc_345041.VBFH125_gamgam.GamGam.root"))
33processes = list(df.keys())
34
35# Apply scale factors and MC weight for simulated events and a weight of 1 for the data
36for p in ["ggH", "VBF"]:
37 df[p] = df[p].Define("weight",
38 "scaleFactor_PHOTON * scaleFactor_PhotonTRIGGER * scaleFactor_PILEUP * mcWeight");
39df["data"] = df["data"].Define("weight", "1.0")
40
41# Select the events for the analysis
42for p in processes:
43 # Apply preselection cut on photon trigger
44 df[p] = df[p].Filter("trigP")
45
46 # Find two good muons with tight ID, pt > 25 GeV and not in the transition region between barrel and encap
47 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))")\
48 .Filter("Sum(goodphotons) == 2")
49
50 # Take only isolated photons
51 df[p] = df[p].Filter("Sum(photon_ptcone30[goodphotons] / photon_pt[goodphotons] < 0.065) == 2")\
52 .Filter("Sum(photon_etcone20[goodphotons] / photon_pt[goodphotons] < 0.065) == 2")
53
54# Compile a function to compute the invariant mass of the diphoton system
55ROOT.gInterpreter.Declare(
56"""
57using Vec_t = const ROOT::VecOps::RVec<float>;
58float ComputeInvariantMass(Vec_t& pt, Vec_t& eta, Vec_t& phi, Vec_t& e) {
59 ROOT::Math::PtEtaPhiEVector p1(pt[0], eta[0], phi[0], e[0]);
60 ROOT::Math::PtEtaPhiEVector p2(pt[1], eta[1], phi[1], e[1]);
61 return (p1 + p2).mass() / 1000.0;
62}
63""")
64
65# Define a new column with the invariant mass and perform final event selection
66hists = {}
67for p in processes:
68 # Make four vectors and compute invariant mass
69 df[p] = df[p].Define("m_yy", "ComputeInvariantMass(photon_pt[goodphotons], photon_eta[goodphotons], photon_phi[goodphotons], photon_E[goodphotons])")
70
71 # Make additional kinematic cuts and select mass window
72 df[p] = df[p].Filter("photon_pt[goodphotons][0] / 1000.0 / m_yy > 0.35")\
73 .Filter("photon_pt[goodphotons][1] / 1000.0 / m_yy > 0.25")\
74 .Filter("m_yy > 105 && m_yy < 160")
75
76 # Book histogram of the invariant mass with this selection
77 hists[p] = df[p].Histo1D(
78 ROOT.RDF.TH1DModel(p, "Diphoton invariant mass; m_{#gamma#gamma} [GeV];Events", 30, 105, 160),
79 "m_yy", "weight")
80
81# Run the event loop
82
83# RunGraphs allows to run the event loops of the separate RDataFrame graphs
84# concurrently. This results in an improved usage of the available resources
85# if each separate RDataFrame can not utilize all available resources, e.g.,
86# because not enough data is available.
87ROOT.RDF.RunGraphs([hists[s] for s in ["ggH", "VBF", "data"]])
88
89ggh = hists["ggH"].GetValue()
90vbf = hists["VBF"].GetValue()
91data = hists["data"].GetValue()
92
93# Create the plot
94
95# Set styles
96ROOT.gROOT.SetStyle("ATLAS")
97
98# Create canvas with pads for main plot and data/MC ratio
99c = ROOT.TCanvas("c", "", 700, 750)
100
101upper_pad = ROOT.TPad("upper_pad", "", 0, 0.35, 1, 1)
102lower_pad = ROOT.TPad("lower_pad", "", 0, 0, 1, 0.35)
103for p in [upper_pad, lower_pad]:
104 p.SetLeftMargin(0.14)
105 p.SetRightMargin(0.05)
106 p.SetTickx(False)
107 p.SetTicky(False)
108upper_pad.SetBottomMargin(0)
109lower_pad.SetTopMargin(0)
110lower_pad.SetBottomMargin(0.3)
111
112upper_pad.Draw()
113lower_pad.Draw()
114
115# Fit signal + background model to data
116upper_pad.cd()
117fit = ROOT.TF1("fit", "([0]+[1]*x+[2]*x^2+[3]*x^3)+[4]*exp(-0.5*((x-[5])/[6])^2)", 105, 160)
118fit.FixParameter(5, 125.0)
119fit.FixParameter(4, 119.1)
120fit.FixParameter(6, 2.39)
121fit.SetLineColor(2)
122fit.SetLineStyle(1)
123fit.SetLineWidth(2)
124data.Fit("fit", "", "E SAME", 105, 160)
125fit.Draw("SAME")
126
127# Draw background
128bkg = ROOT.TF1("bkg", "([0]+[1]*x+[2]*x^2+[3]*x^3)", 105, 160)
129for i in range(4):
130 bkg.SetParameter(i, fit.GetParameter(i))
131bkg.SetLineColor(4)
132bkg.SetLineStyle(2)
133bkg.SetLineWidth(2)
134bkg.Draw("SAME")
135
136# Draw data
137data.SetMarkerStyle(20)
138data.SetMarkerSize(1.2)
139data.SetLineWidth(2)
140data.SetLineColor(ROOT.kBlack)
141data.Draw("E SAME")
142data.SetMinimum(1e-3)
143data.SetMaximum(8e3)
144data.GetYaxis().SetLabelSize(0.045)
145data.GetYaxis().SetTitleSize(0.05)
146data.SetStats(0)
147data.SetTitle("")
148
149# Scale simulated events with luminosity * cross-section / sum of weights
150# and merge to single Higgs signal
151lumi = 10064.0
152ggh.Scale(lumi * 0.102 / 55922617.6297)
153vbf.Scale(lumi * 0.008518764 / 3441426.13711)
154higgs = ggh.Clone()
155higgs.Add(vbf)
156higgs.Draw("HIST SAME")
157
158# Draw ratio
159lower_pad.cd()
160
161ratiobkg = ROOT.TF1("zero", "0", 105, 160)
162ratiobkg.SetLineColor(4)
163ratiobkg.SetLineStyle(2)
164ratiobkg.SetLineWidth(2)
165ratiobkg.SetMinimum(-125)
166ratiobkg.SetMaximum(250)
167ratiobkg.GetXaxis().SetLabelSize(0.08)
168ratiobkg.GetXaxis().SetTitleSize(0.12)
169ratiobkg.GetXaxis().SetTitleOffset(1.0)
170ratiobkg.GetYaxis().SetLabelSize(0.08)
171ratiobkg.GetYaxis().SetTitleSize(0.09)
172ratiobkg.GetYaxis().SetTitle("Data - Bkg.")
173ratiobkg.GetYaxis().CenterTitle()
174ratiobkg.GetYaxis().SetTitleOffset(0.7)
175ratiobkg.GetYaxis().SetNdivisions(503, False)
176ratiobkg.GetYaxis().ChangeLabel(-1, -1, 0)
177ratiobkg.GetXaxis().SetTitle("m_{#gamma#gamma} [GeV]")
178ratiobkg.Draw()
179
180ratiosig = ROOT.TH1F("ratiosig", "ratiosig", 5500, 105, 160)
181ratiosig.Eval(fit)
182ratiosig.SetLineColor(2)
183ratiosig.SetLineStyle(1)
184ratiosig.SetLineWidth(2)
185ratiosig.Add(bkg, -1)
186ratiosig.Draw("SAME")
187
188ratiodata = data.Clone()
189ratiodata.Add(bkg, -1)
190ratiodata.Draw("E SAME")
191for i in range(1, data.GetNbinsX()):
192 ratiodata.SetBinError(i, data.GetBinError(i))
193
194# Add legend
195upper_pad.cd()
196legend = ROOT.TLegend(0.55, 0.55, 0.89, 0.85)
197legend.SetTextFont(42)
198legend.SetFillStyle(0)
199legend.SetBorderSize(0)
200legend.SetTextSize(0.05)
201legend.SetTextAlign(32)
202legend.AddEntry(data, "Data" ,"lep")
203legend.AddEntry(bkg, "Background", "l")
204legend.AddEntry(fit, "Signal + Bkg.", "l")
205legend.AddEntry(higgs, "Signal", "l")
206legend.Draw("SAME")
207
208# Add ATLAS label
209text = ROOT.TLatex()
210text.SetNDC()
211text.SetTextFont(72)
212text.SetTextSize(0.05)
213text.DrawLatex(0.18, 0.84, "ATLAS")
214text.SetTextFont(42)
215text.DrawLatex(0.18 + 0.13, 0.84, "Open Data")
216text.SetTextSize(0.04)
217text.DrawLatex(0.18, 0.78, "#sqrt{s} = 13 TeV, 10 fb^{-1}")
218
219# Save the plot
220c.SaveAs("df104_HiggsToTwoPhotons.png")
221print("Saved figure to df104_HiggsToTwoPhotons.png")
ROOT's RDataFrame offers a high level interface for analyses of data stored in TTrees,...
void RunGraphs(std::vector< RResultHandle > handles)
Trigger the event loop of multiple RDataFrames concurrently.
A struct which stores the parameters of a TH1D.