Rivet analyses referenceATLAS_2017_I1625109Measurement of $ZZ -> 4\ell$ production at 13 TeVExperiment: ATLAS (LHC) Inspire ID: 1625109 Status: VALIDATED Authors:
Beam energies: (6500.0, 6500.0) GeV Run details:
Measurements of $ZZ$ production in the $\ell^+\ell^-\ell^{\prime +}\ell^{\prime -}$-channel in proton--proton collisions at 13 TeV center-of-mass energy at the Large Hadron Collider are presented. The data correspond to 36.1 fb${}^{-1}$ of collisions collected by the ATLAS experiment in 2015 and 2016. Here $\ell$ and $\ell^\prime$ stand for electrons or muons. Integrated and differential $ZZ\to \ell^+\ell^-\ell^{\prime +}\ell^{\prime -}$-cross sections with $Z\to\ell^+\ell^-$-candidate masses in the range of 66 GeV to 116 GeV are measured in a fiducial phase space corresponding to the detector acceptance and corrected for detector effects. The differential cross sections are presented in bins of twenty observables, including several that describe the jet activity. The integrated cross section is also extrapolated to a total phase space and to all standard model decays of $Z$ bosons with mass between 66 GeV and 116 GeV, resulting in a value of 17.3$\pm$0.9[$\pm$0.6(stat)$\pm$0.5(syst)$\pm$0.6(lumi)] pb. The measurements are found to be in good agreement with the standard model. A search for neutral triple gauge couplings is performed using the transverse momentum distribution of the leading $Z$ boson candidate. No evidence for such couplings is found and exclusion limits are set on their parameters. Source code: ATLAS_2017_I1625109.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/PromptFinalState.hh"
5#include "Rivet/Projections/VetoedFinalState.hh"
6#include "Rivet/Projections/LeptonFinder.hh"
7#include "Rivet/Projections/FastJets.hh"
8
9namespace Rivet {
10
11
12 /// @brief measurement of on-shell ZZ at 13 TeV
13 class ATLAS_2017_I1625109 : public Analysis {
14 public:
15
16 /// Constructor
17 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2017_I1625109);
18
19 /// @name Analysis methods
20 /// @{
21
22 struct Dilepton {
23 Dilepton() {};
24 Dilepton(const ParticlePair & _leptons) : leptons(_leptons) {}
25
26 FourMomentum momentum() const {
27 return leptons.first.mom() + leptons.second.mom();
28 }
29
30 ParticlePair leptons;
31 };
32
33
34 struct Quadruplet {
35
36 DressedLeptons getLeptonsSortedByPt() const {
37 DressedLeptons out = { leadingDilepton.leptons.first, leadingDilepton.leptons.second,
38 subleadingDilepton.leptons.first, subleadingDilepton.leptons.second };
39 std::sort(out.begin(), out.end(), cmpMomByPt);
40 return out;
41 }
42
43 Quadruplet(const Dilepton& dilepton1, const Dilepton& dilepton2) {
44 if (dilepton1.momentum().pt() > dilepton2.momentum().pt()) {
45 leadingDilepton = dilepton1;
46 subleadingDilepton = dilepton2;
47 }
48 else {
49 leadingDilepton = dilepton2;
50 subleadingDilepton = dilepton1;
51 }
52 leptonsSortedByPt = getLeptonsSortedByPt();
53 }
54
55 FourMomentum momentum() const {
56 return leadingDilepton.momentum() + subleadingDilepton.momentum();
57 }
58
59 double distanceFromZMass() const {
60 return abs(leadingDilepton.momentum().mass() - Z_mass) + abs(subleadingDilepton.momentum().mass() - Z_mass);
61 }
62
63 Dilepton leadingDilepton;
64 Dilepton subleadingDilepton;
65 DressedLeptons leptonsSortedByPt;
66 };
67
68 typedef vector<Quadruplet> Quadruplets;
69
70 typedef std::pair<size_t, size_t> IndexPair;
71
72
73 vector<IndexPair> getOppositeChargePairsIndices(const DressedLeptons& leptons) {
74 vector<IndexPair> indices = {};
75 if (leptons.size() < 2) return indices;
76 for (size_t i = 0; i < leptons.size(); ++i) {
77 for (size_t k = i+1; k < leptons.size(); ++k) {
78 const auto charge_i = leptons.at(i).charge();
79 const auto charge_k = leptons.at(k).charge();
80 if (charge_i == -charge_k) {
81 indices.push_back(std::make_pair(i, k));
82 }
83 }
84 }
85 return indices;
86 }
87
88 bool indicesOverlap(IndexPair a, IndexPair b) {
89 return (a.first == b.first || a.first == b.second || a.second == b.first || a.second == b.second);
90 }
91
92
93 bool passesHierarchicalPtRequirements(const Quadruplet& quadruplet) {
94 const auto& sorted_leptons = quadruplet.leptonsSortedByPt;
95 if (sorted_leptons.at(0).pt() < 20*GeV) return false;
96 if (sorted_leptons.at(1).pt() < 15*GeV) return false;
97 if (sorted_leptons.at(2).pt() < 10*GeV) return false;
98 return true;
99 }
100
101 bool passesDileptonMinimumMassRequirement(const Quadruplet& quadruplet) {
102 const auto& leptons = quadruplet.leptonsSortedByPt;
103 for (const Particle& l1 : leptons) {
104 for (const Particle& l2 : leptons) {
105 if (isSame(l1, l2)) continue;
106 if ((l1.pid() + l2.pid() == 0) && ((l1.mom() + l2.mom()).mass() < 5.0*GeV)) return false;
107 }
108 }
109 return true;
110 }
111
112 bool passesLeptonDeltaRRequirements(const Quadruplet& quadruplet) {
113 const auto& leptons = quadruplet.leptonsSortedByPt;
114 for (const Particle& l1 : leptons) {
115 for (const Particle& l2 : leptons) {
116 if (isSame(l1, l2)) continue;
117 // Any lepton flavour:
118 if (deltaR(l1.mom(), l2.mom()) < 0.1) return false;
119 // Different lepton flavour:
120 if ((l1.abspid() - l2.abspid() != 0) && (deltaR(l1.mom(), l2.mom()) < 0.2)) return false;
121 }
122 }
123 return true;
124 }
125
126 Quadruplets formQuadrupletsByChannel(const DressedLeptons& same_flavour_leptons, vector<IndexPair> indices) {
127 Quadruplets quadruplets = {};
128 for (size_t i = 0; i < indices.size(); ++i) {
129 for (size_t k = i+1; k < indices.size(); ++k) {
130 const auto& pair_i = indices.at(i);
131 const auto& pair_k = indices.at(k);
132 if (indicesOverlap(pair_i, pair_k)) continue;
133 const auto d1 = Dilepton({same_flavour_leptons.at(pair_i.first), same_flavour_leptons.at(pair_i.second)});
134 const auto d2 = Dilepton({same_flavour_leptons.at(pair_k.first), same_flavour_leptons.at(pair_k.second)});
135 const auto quadruplet = Quadruplet(d1, d2);
136 if (passesHierarchicalPtRequirements(quadruplet)) quadruplets.push_back(quadruplet);
137 }
138 }
139 return quadruplets;
140 }
141
142 Quadruplets formQuadrupletsByChannel(const DressedLeptons& electrons, vector<IndexPair> e_indices,
143 const DressedLeptons& muons, vector<IndexPair> m_indices) {
144 Quadruplets quadruplets = {};
145 for (const auto& pair_e : e_indices) {
146 for (const auto& pair_m : m_indices) {
147 const auto d1 = Dilepton({electrons.at(pair_e.first), electrons.at(pair_e.second)});
148 const auto d2 = Dilepton({muons.at(pair_m.first), muons.at(pair_m.second)});
149 const auto quadruplet = Quadruplet(d1, d2);
150 if (passesHierarchicalPtRequirements(quadruplet)) quadruplets.push_back(quadruplet);
151 }
152 }
153 return quadruplets;
154 }
155
156
157 Quadruplets getQuadruplets(const DressedLeptons& electrons, const DressedLeptons& muons) {
158 const auto oc_electrons_indices = getOppositeChargePairsIndices(electrons);
159 const auto oc_muons_indices = getOppositeChargePairsIndices(muons);
160
161 const auto electron_quadruplets = formQuadrupletsByChannel(electrons, oc_electrons_indices);
162 const auto muon_quadruplets = formQuadrupletsByChannel(muons, oc_muons_indices);
163 const auto mixed_quadruplets = formQuadrupletsByChannel(electrons, oc_electrons_indices, muons, oc_muons_indices);
164
165 auto quadruplets = electron_quadruplets;
166 quadruplets.insert(quadruplets.end(), muon_quadruplets.begin(), muon_quadruplets.end());
167 quadruplets.insert(quadruplets.end(), mixed_quadruplets.begin(), mixed_quadruplets.end());
168
169 return quadruplets;
170 }
171
172
173 Quadruplet selectQuadruplet(const Quadruplets& quadruplets) {
174 if (quadruplets.empty()) throw std::logic_error("Expect at least one quadruplet! The user should veto events without quadruplets.");
175 Quadruplets sortedQuadruplets = quadruplets;
176 std::sort(sortedQuadruplets.begin(), sortedQuadruplets.end(), [](const Quadruplet& a, const Quadruplet& b) {
177 return a.distanceFromZMass() < b.distanceFromZMass();
178 });
179 return sortedQuadruplets.at(0);
180 }
181
182 /// Book histograms and initialise projections before the run
183 void init() {
184 const Cut presel = Cuts::abseta < 5 && Cuts::pT > 100*MeV;
185 const FinalState fs(presel);
186
187 // Prompt leptons, photons, neutrinos
188 // Excluding those from tau decay
189 const PromptFinalState photons(presel && Cuts::abspid == PID::PHOTON, TauDecaysAs::NONPROMPT);
190 const PromptFinalState bare_elecs(presel && Cuts::abspid == PID::ELECTRON, TauDecaysAs::NONPROMPT);
191 const PromptFinalState bare_muons(presel && Cuts::abspid == PID::MUON, TauDecaysAs::NONPROMPT);
192
193 // Baseline lepton and jet declaration
194 const Cut lepton_baseline_cuts = Cuts::abseta < 2.7 && Cuts::pT > 5*GeV;
195 const LeptonFinder elecs = LeptonFinder(bare_elecs, photons, 0.1, lepton_baseline_cuts);
196 const LeptonFinder muons = LeptonFinder(bare_muons, photons, 0.1, lepton_baseline_cuts);
197 declare(elecs, "electrons");
198 declare(muons, "muons");
199
200 VetoedFinalState jet_input(fs);
201 jet_input.addVetoOnThisFinalState(elecs);
202 jet_input.addVetoOnThisFinalState(muons);
203 declare(FastJets(jet_input, JetAlg::ANTIKT, 0.4), "jets");
204
205 // // Book histograms
206 book(_h["pT_4l"], 2, 1, 1);
207 book(_h["pT_leading_dilepton"], 8, 1, 1);
208 book(_h["pT_subleading_dilepton"], 14, 1, 1);
209 book(_h["pT_lepton1"], 20, 1, 1);
210 book(_h["pT_lepton2"], 26, 1, 1);
211 book(_h["pT_lepton3"], 32, 1, 1);
212 book(_h["pT_lepton4"], 38, 1, 1);
213 book(_h["absy_4l"], 44, 1, 1);
214 book(_h["deltay_dileptons"], 50, 1, 1);
215 book(_h["deltaphi_dileptons"], 56, 1, 1);
216 book(_h["N_jets"], 62, 1, 1);
217 book(_h["N_central_jets"], 68, 1, 1);
218 book(_h["N_jets60"], 74, 1, 1);
219 book(_h["mass_dijet"], 80, 1, 1);
220 book(_h["deltay_dijet"], 86, 1, 1);
221 book(_h["scalarpTsum_jets"], 92, 1, 1);
222 book(_h["abseta_jet1"], 98, 1, 1);
223 book(_h["abseta_jet2"], 104, 1, 1);
224 book(_h["pT_jet1"], 110, 1, 1);
225 book(_h["pT_jet2"], 116, 1, 1);
226 }
227
228
229 /// Perform the per-event analysis
230 void analyze(Event const & event) {
231 const auto& baseline_electrons = apply<LeptonFinder>(event, "electrons").dressedLeptons();
232 const auto& baseline_muons = apply<LeptonFinder>(event, "muons").dressedLeptons();
233
234 // Form all possible quadruplets passing hierarchical lepton pT cuts
235 const auto quadruplets = getQuadruplets(baseline_electrons, baseline_muons);
236
237 if (quadruplets.empty()) vetoEvent;
238
239 // Select the best quadruplet, the one minimising the total distance from the Z pole mass
240 auto const quadruplet = selectQuadruplet(quadruplets);
241
242 // Event selection on the best quadruplet
243 if (!passesDileptonMinimumMassRequirement(quadruplet)) vetoEvent;
244 if (!passesLeptonDeltaRRequirements(quadruplet)) vetoEvent;
245 if (!inRange(quadruplet.leadingDilepton.momentum().mass(), 66*GeV, 116*GeV)) vetoEvent;
246 if (!inRange(quadruplet.subleadingDilepton.momentum().mass(), 66*GeV, 116*GeV)) vetoEvent;
247
248 // Select jets
249 Jets alljets = apply<JetFinder>(event, "jets").jetsByPt(Cuts::pT > 30*GeV);
250 for (const DressedLepton& lep : quadruplet.leptonsSortedByPt)
251 idiscard(alljets, deltaRLess(lep, 0.4));
252 const Jets jets = alljets;
253 const Jets centralJets = select(jets, Cuts::abseta < 2.4);
254 const Jets pt60Jets = select(jets, Cuts::pT > 60*GeV);
255
256 const auto& leadingDilepton = quadruplet.leadingDilepton.momentum();
257 const auto& subleadingDilepton = quadruplet.subleadingDilepton.momentum();
258
259 _h["pT_4l"]->fill((leadingDilepton + subleadingDilepton).pt()/GeV);
260 _h["pT_leading_dilepton"]->fill(leadingDilepton.pt()/GeV);
261 _h["pT_subleading_dilepton"]->fill(subleadingDilepton.pt()/GeV);
262 _h["pT_lepton1"]->fill(quadruplet.leptonsSortedByPt.at(0).pt()/GeV);
263 _h["pT_lepton2"]->fill(quadruplet.leptonsSortedByPt.at(1).pt()/GeV);
264 _h["pT_lepton3"]->fill(quadruplet.leptonsSortedByPt.at(2).pt()/GeV);
265 _h["pT_lepton4"]->fill(quadruplet.leptonsSortedByPt.at(3).pt()/GeV);
266 _h["absy_4l"]->fill((leadingDilepton + subleadingDilepton).absrapidity());
267 _h["deltay_dileptons"]->fill(fabs(leadingDilepton.rapidity() - subleadingDilepton.rapidity()));
268 _h["deltaphi_dileptons"]->fill(deltaPhi(leadingDilepton, subleadingDilepton)/pi);
269 _h["N_jets"]->fill(jets.size());
270 _h["N_central_jets"]->fill(centralJets.size());
271 _h["N_jets60"]->fill(pt60Jets.size());
272
273 // If at least one jet present
274 if (jets.empty()) vetoEvent;
275 _h["scalarpTsum_jets"]->fill(sum(jets, pT, 0.)/GeV);
276 _h["abseta_jet1"]->fill(jets.front().abseta());
277 _h["pT_jet1"]->fill(jets.front().pt()/GeV);
278
279 // If at least two jets present
280 if (jets.size() < 2) vetoEvent;
281 _h["mass_dijet"]->fill((jets.at(0).mom() + jets.at(1).mom()).mass()/GeV);
282 _h["deltay_dijet"]->fill(fabs(jets.at(0).rapidity() - jets.at(1).rapidity()));
283 _h["abseta_jet2"]->fill(jets.at(1).abseta());
284 _h["pT_jet2"]->fill(jets.at(1).pt()/GeV);
285 }
286
287
288 /// Normalise histograms etc., after the run
289 void finalize() {
290 // Normalise histograms to cross section
291 const double sf = crossSectionPerEvent() / femtobarn;
292 scale(_h, sf);
293 }
294 /// @}
295
296 private:
297 /// @name Histograms
298 /// @{
299 map<string, Histo1DPtr> _h;
300 static constexpr double Z_mass = 91.1876;
301 /// @}
302 };
303
304
305 RIVET_DECLARE_PLUGIN(ATLAS_2017_I1625109);
306}
|