Skip to content

Commit 36f0c3e

Browse files
taehyounparktaehyounpark
committed
docs
1 parent 5f19280 commit 36f0c3e

File tree

1 file changed

+200
-101
lines changed

1 file changed

+200
-101
lines changed

‎docs/pages/example.md

Lines changed: 200 additions & 101 deletions
Original file line numberDiff line numberDiff line change
@@ -54,114 +54,213 @@ int main() {
5454
}
5555
```
5656
57-
@section example-hww ROOT TTree
57+
@section example-hww ROOT TTree with systematic variations
5858
5959
- Simulated ggF HWW* events: [ATLAS open data](https://opendata.cern.ch/record/3825).
60-
- ROOT extensions for queryosity: [queryosity-hep](https://github.com/taehyounpark/AnalysisPlugins)
60+
- ROOT extensions for queryosity: [queryosity-hep](https://github.com/taehyounpark/queryosity-hep)
6161
6262
1. Apply the MC event weight.
6363
2. Select entries for which there are exactly two opposite-sign leptons in the event.
64-
3. Form a (2, 2) matrix of cut regions:
65-
- The leptons are same-/different-flavour.
66-
- The dilepton invariant mass is less/greater than 60 GeV.
67-
4. In each case, plot the distribution of the dilepton+MET transverse momentum.
68-
- Scale lepton tranvserse momenta by +/- 2% as systematic variations.
64+
3. Separate into different/same-flavour cases for electrons and muons.
65+
4. Apply signal region cuts: <i>E</i><sub>T</sub><sup>miss</sup> > 30 GeV, <i>m</i><sub>ll</sub> < 60 GeV.
66+
5. Merge back together to form flavour-inclusive opposite-sign signal region.
67+
6. In each case, plot the distribution of the dilepton+MET transverse momentum.
68+
- Scale lepton energy scale by +/- 2% as systematic variations.
6969
70-
(To-do: update example to newer API).
70+
(The order of selections is needlessly & purposefully convoluted to showcase the cutflow interface)
71+
72+
@cpp
73+
#include "qhep/Hist.h"
74+
#include "qhep/Tree.h"
75+
76+
#include "queryosity/queryosity.h"
77+
78+
namespace qty = queryosity;
79+
using dataflow = qty::dataflow;
80+
namespace multithread = qty::multithread;
81+
namespace dataset = qty::dataset;
82+
namespace column = qty::column;
83+
namespace query = qty::query;
84+
namespace systematic = qty::systematic;
85+
86+
#include "Math/Vector4D.h"
87+
#include "TCanvas.h"
88+
#include "TH1F.h"
89+
#include "TStyle.h"
90+
#include "TVector2.h"
91+
#include <ROOT/RVec.hxx>
92+
93+
using VecF = ROOT::RVec<float>;
94+
using VecD = ROOT::RVec<double>;
95+
using VecI = ROOT::RVec<int>;
96+
using VecUI = ROOT::RVec<unsigned int>;
97+
using P4 = ROOT::Math::PtEtaPhiEVector;
98+
99+
#include <chrono>
100+
#include <functional>
101+
#include <iostream>
102+
#include <memory>
103+
#include <sstream>
104+
105+
class NthP4 : public column::definition<P4(VecD, VecD, VecD, VecD)> {
106+
107+
public:
108+
NthP4(unsigned int index) : m_index(index) {}
109+
virtual ~NthP4() = default;
110+
111+
virtual P4 evaluate(column::observable<VecD> pt, column::observable<VecD> eta,
112+
column::observable<VecD> phi,
113+
column::observable<VecD> es) const override {
114+
return P4(pt->at(m_index), eta->at(m_index), phi->at(m_index),
115+
es->at(m_index));
116+
}
117+
118+
protected:
119+
const unsigned int m_index;
120+
};
121+
122+
int main() {
123+
124+
// the tree doesn't have enough events to multithread
125+
dataflow df(multithread::disable());
126+
127+
std::vector<std::string> tree_files{"hww.root"};
128+
std::string tree_name = "mini";
129+
auto ds = df.load(dataset::input<Tree>(tree_files, tree_name));
130+
131+
// weights
132+
auto mc_weight = ds.read(dataset::column<float>("mcWeight"));
133+
auto mu_sf = ds.read(dataset::column<float>("scaleFactor_MUON"));
134+
auto el_sf = ds.read(dataset::column<float>("scaleFactor_ELE"));
135+
136+
// leptons
137+
auto [lep_pt_MeV, lep_eta, lep_phi, lep_E_MeV, lep_Q, lep_type] = ds.read(
138+
dataset::column<VecF>("lep_pt"), dataset::column<VecF>("lep_eta"),
139+
dataset::column<VecF>("lep_phi"), dataset::column<VecF>("lep_E"),
140+
dataset::column<VecF>("lep_charge"), dataset::column<VecUI>("lep_type"));
141+
142+
// missing transverse energy
143+
auto [met_MeV, met_phi] = ds.read(dataset::column<float>("met_et"),
144+
dataset::column<float>("met_phi"));
145+
146+
// units
147+
auto MeV = df.define(column::constant(1000.0));
148+
auto lep_pt = lep_pt_MeV / MeV;
149+
auto lep_E = lep_E_MeV / MeV;
150+
auto met = met_MeV / MeV;
151+
152+
// vary the energy scale by +/-2%
153+
auto Escale =
154+
df.vary(column::expression([](VecD E) { return E; }),
155+
systematic::variation("eg_up", [](VecD E) { return E * 1.02; }),
156+
systematic::variation("eg_dn", [](VecD E) { return E * 0.98; }));
157+
158+
// apply the energy scale (uncertainties)
159+
// and select ones within acceptance
160+
auto lep_pt_min = df.define(column::constant(15.0));
161+
auto lep_eta_max = df.define(column::constant(2.4));
162+
auto lep_selection = (lep_eta < lep_eta_max) && (lep_eta > (-lep_eta_max)) &&
163+
(lep_pt > lep_pt_min);
164+
auto lep_pt_sel = Escale(lep_pt)[lep_selection];
165+
auto lep_E_sel = Escale(lep_E)[lep_selection];
166+
auto lep_eta_sel = lep_eta[lep_selection];
167+
auto lep_phi_sel = lep_phi[lep_selection];
168+
169+
// put (sub-)leading lepton into four-momentum
170+
auto l1p4 = df.define(column::definition<NthP4>(0), lep_pt_sel, lep_eta_sel,
171+
lep_phi_sel, lep_E_sel);
172+
auto l2p4 = df.define(column::definition<NthP4>(1), lep_pt_sel, lep_eta_sel,
173+
lep_phi_sel, lep_E_sel);
174+
175+
// compute dilepton invariant mass & higgs transverse momentum
176+
auto llp4 = l1p4 + l2p4;
177+
auto mll =
178+
df.define(column::expression([](const P4 &p4) { return p4.M(); }), llp4);
179+
auto higgs_pt =
180+
df.define(column::expression([](const P4 &p4, float q, float q_phi) {
181+
TVector2 p2;
182+
p2.SetMagPhi(p4.Pt(), p4.Phi());
183+
TVector2 q2;
184+
q2.SetMagPhi(q, q_phi);
185+
return (p2 + q2).Mod();
186+
}),
187+
llp4, met, met_phi);
188+
189+
// compute number of leptons
190+
auto nlep_req = df.define(column::constant<unsigned int>(2));
191+
auto nlep_sel =
192+
df.define(column::expression([](VecD const &lep) { return lep.size(); }),
193+
lep_pt_sel);
194+
195+
// apply cuts & weights
196+
auto weighted = df.weight(mc_weight * el_sf * mu_sf);
197+
auto cut_2l = weighted.filter(nlep_sel == nlep_req);
198+
auto cut_2los =
199+
cut_2l.filter(column::expression([](const VecF &lep_charge) {
200+
return lep_charge.at(0) + lep_charge.at(1) == 0;
201+
}),
202+
lep_Q);
203+
auto cut_2ldf =
204+
cut_2los.filter(column::expression([](const VecI &lep_type) {
205+
return lep_type.at(0) + lep_type.at(1) == 24;
206+
}),
207+
lep_type);
208+
auto cut_2lsf =
209+
cut_2los.filter(column::expression([](const VecI &lep_type) {
210+
return (lep_type.at(0) + lep_type.at(1) == 22) ||
211+
(lep_type.at(0) + lep_type.at(1) == 26);
212+
}),
213+
lep_type);
214+
auto mll_max = df.define(column::constant(80.0));
215+
auto met_min = df.define(column::constant(30.0));
216+
auto cut_2ldf_sr = cut_2ldf.filter((mll < mll_max) && (met > met_min));
217+
auto cut_2lsf_sr = cut_2lsf.filter((mll < mll_max) && (met > met_min));
218+
auto cut_2los_sr =
219+
df.filter(cut_2ldf_sr || cut_2lsf_sr).weight(mc_weight * el_sf * mu_sf);
220+
// once two selections are joined, they "forget" everything upstream
221+
// i.e. need to re-apply the event weight!
222+
223+
// histograms:
224+
// - at three regions
225+
// - nominal & eg_up & eg_dn
226+
auto [pth_2los_sr, pth_2ldf_sr, pth_2lsf_sr] =
227+
df.make(query::plan<Hist<1, float>>("pth", 30, 0, 150))
228+
.fill(higgs_pt)
229+
.book(cut_2los_sr, cut_2ldf_sr, cut_2lsf_sr);
230+
// all done at once :)
231+
232+
Double_t w = 1600;
233+
Double_t h = 800;
234+
TCanvas c("c", "c", w, h);
235+
c.SetWindowSize(w + (w - c.GetWw()), h + (h - c.GetWh()));
236+
c.Divide(3, 1);
237+
c.cd(1);
238+
pth_2los_sr.nominal()->SetLineColor(kBlack);
239+
pth_2los_sr.nominal()->Draw("ep");
240+
pth_2los_sr["eg_up"]->SetLineColor(kRed);
241+
pth_2los_sr["eg_up"]->Draw("same hist");
242+
pth_2los_sr["eg_dn"]->SetLineColor(kBlue);
243+
pth_2los_sr["eg_dn"]->Draw("same hist");
244+
c.cd(2);
245+
pth_2ldf_sr.nominal()->SetLineColor(kBlack);
246+
pth_2ldf_sr.nominal()->Draw("ep");
247+
pth_2ldf_sr["eg_up"]->SetLineColor(kRed);
248+
pth_2ldf_sr["eg_up"]->Draw("same hist");
249+
pth_2ldf_sr["eg_dn"]->SetLineColor(kBlue);
250+
pth_2ldf_sr["eg_dn"]->Draw("same hist");
251+
c.cd(3);
252+
pth_2lsf_sr.nominal()->SetLineColor(kBlack);
253+
pth_2lsf_sr.nominal()->Draw("ep");
254+
pth_2lsf_sr["eg_up"]->SetLineColor(kRed);
255+
pth_2lsf_sr["eg_up"]->Draw("same hist");
256+
pth_2lsf_sr["eg_dn"]->SetLineColor(kBlue);
257+
pth_2lsf_sr["eg_dn"]->Draw("same hist");
258+
c.SaveAs("pth.png");
259+
260+
return 0;
261+
}
262+
@endcpp
71263
72-
```cpp
73-
dataflow df;
74-
75-
auto tree_files = std::vector<std::string>{"hww.root"};
76-
auto tree_name = "mini";
77-
auto ds = df.load(dataset::input(tree_files, tree_name));
78-
79-
auto [mc_weight, el_sf, mu_sf] = ds.read(
80-
dataset::column<float>("mcWeight"),
81-
dataset::column<float>("scaleFactor_ELE"),
82-
dataset::column<float>("scaleFactor_MUON")
83-
);
84-
85-
auto [
86-
lep_pt_MeV,
87-
lep_eta,
88-
lep_phi,
89-
lep_E_MeV,
90-
lep_Q,
91-
lep_type
92-
] = ds.read(
93-
dataset::column<VecF>("lep_pt"),
94-
dataset::column<VecF>("lep_eta"),
95-
dataset::column<VecF>("lep_phi"),
96-
dataset::column<VecF>("lep_E"),
97-
dataset::column<VecF>("lep_E"),
98-
dataset::column<VecF>("lep_charge"),
99-
dataset::column<VecF>("lep_type")
100-
);
101-
102-
auto [met_MeV, met_phi] = ds.read(
103-
dataset::column<float>("met_et"),
104-
dataset::column<float>("met_phi")
105-
);
106-
107-
auto MeV = ana.constant(1000.0);
108-
auto lep_pt = lep_pt_MeV / MeV;
109-
auto lep_E = lep_pt_MeV / MeV;
110-
auto met = met_MeV / MeV;
111-
112-
auto lep_eta_max = df.constant(2.4);
113-
auto lep_pt_sel = lep_pt[ lep_eta < lep_eta_max && lep_eta > (-lep_eta_max) ];
114-
115-
auto p4l1 = df.define<NthP4>(0)(lep_pt_sel, lep_eta_sel, lep_phi_sel, lep_E_sel);
116-
auto p4l2 = df.define<NthP4>(1)(lep_pt_sel, lep_eta_sel, lep_phi_sel, lep_E_sel);
117-
118-
auto p4ll = p4l1+p4l2;
119-
120-
using P4 = TLorentzVector;
121-
auto mll = df.define([](const P4& p4){return p4.M();})(p4ll);
122-
123-
auto pth = df.define(
124-
[](const P4& p4, float q, float q_phi) {
125-
TVector2 p2; p2.SetMagPhi(p4.Pt(), p4.Phi());
126-
TVector2 q2; q2.SetMagPhi(q, q_phi);
127-
return (p2+q2).Mod();
128-
})(p4ll, met, met_phi);
129-
130-
auto n_lep_sel = df.define([](VecF const& lep){return lep.size();},lep_pt_sel);
131-
auto n_lep_req = df.define(column::constant<unsigned int>(2));
132-
133-
auto cut_2l = df.weight(mc_weight * el_sf * mu_sf)\
134-
.filter(n_lep_sel == n_lep_req);
135-
136-
auto cut_2los = cut_2l.filter([](const VecI& lep_charge){return lep_charge.book(0)+lep_charge.book(1)==0;},lep_Q);
137-
auto cut_2ldf = cut_2los.filter([](const VecUI& lep_type){return lep_type.book(0)+lep_type.book(1)==24;},lep_type);
138-
auto cut_2lsf = cut_2los.filter([](const VecUI& lep_type){return (lep_type.book(0)+lep_type.book(1)==22)||(lep_type.book(0)+lep_type.book(1)==26);},lep_type);
139-
140-
auto mll_cut = df.constant(60.0);
141-
auto cut_2ldf_sr = cut_2ldf.filter(mll < mll_cut);
142-
auto cut_2lsf_sr = cut_2lsf.filter(mll < mll_cut);
143-
auto cut_2ldf_wwcr = cut_2ldf.filter(mll > mll_cut);
144-
auto cut_2lsf_wwcr = cut_2lsf.filter(mll > mll_cut);
145-
146-
auto pth_hist = df.get<Hist<1,float>>("pth",100,0,400).fill(pth).book(cut_2los);
147-
148-
auto l1n2_pt_hists_wwcrs = l1n2_pt_hist.book(cut_2lsf_sr, cut_2lsf_wwcr);
149-
150-
auto Escale = df.define([](VecD E){return E;}).vary("lp4_up",[](VecD E){return E*1.02;}).vary("lp4_dn",[](VecD E){return E*0.98;});
151-
auto lep_pt_sel = Escale(lep_pt)[ lep_eta < lep_eta_max && lep_eta > (-lep_eta_max) ];
152-
auto lep_E_sel = Escale(lep_E)[ lep_eta < lep_eta_max && lep_eta > (-lep_eta_max) ];
153-
154-
auto p4l1 = df.define<NthP4>(0)(lep_pt, lep_eta, lep_phi, lep_E);
155-
auto p4l2 = df.define<NthP4>(1)(lep_pt, lep_eta, lep_phi, lep_E);
156-
157-
auto cut_2l = df.weight("weight")(mc_weight * el_sf * mu_sf)\
158-
.filter("2l")(n_lep_sel == n_lep_req);
159-
160-
auto mll_vars = df.get<Hist<1,float>>("mll",50,0,100).fill(mll).book(cut_2los);
161-
162-
mll_vars.nominal()->Draw();
163-
mll_vars["lp4_up"]->Draw("same");
164-
```
165264
@image html mll_varied.png
166265
167266
@section example-phys ATLAS DAOD_PHYS

0 commit comments

Comments
 (0)