Rivet analyses referenceATLAS_2018_I1711223Observation of electroweak WZjj productionExperiment: ATLAS (LHC) Inspire ID: 1711223 Status: VALIDATED Authors:
Beams: p+ p+ Beam energies: (6500.0, 6500.0) GeV Run details:
An observation of electroweak $W^\pm Z$ production in association with two jets in proton-proton collisions is presented. The data collected by the ATLAS detector at the Large Hadron Collider in 2015 and 2016 at a centre-of-mass energy of $\sqrt{s}=13$ TeV are used, corresponding to an integrated luminosity of 36.1fb$^{-1}$. Events containing three identified leptons, either electrons or muons, and two jets are selected. The electroweak production of $W^\pm Z$ bosons in association with two jets is measured with an observed significance of 5.3 standard deviations. A fiducial cross-section for electroweak production including interference effects and for a single leptonic decay mode is measured to be $\sigma_{WZjj-EW}=0.57-0.13+0.14$(stat.)$-0.06+0.07$(syst.)fb. Total and differential fiducial cross-sections of the sum of $W^\pm Zjj$ electroweak and strong productions for several kinematic observables are also measured. Uses SM neutrino-lepton flavour matching and a resonant shape algorithm assuming the Standard Model, to match the MC-based correction to the fiducial region applied in the paper. This routine is therefore only valid under the assumption of the Standard Model and cannot be used for BSM reinterpretation Source code: ATLAS_2018_I1711223.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/FastJets.hh"
5#include "Rivet/Projections/IdentifiedFinalState.hh"
6#include "Rivet/Projections/PromptFinalState.hh"
7#include "Rivet/Projections/DressedLeptons.hh"
8#include "Rivet/Projections/VetoedFinalState.hh"
9
10
11namespace Rivet {
12
13
14 /// @brief Electroweak WZjj production cross section at 13 TeV
15 class ATLAS_2018_I1711223 : public Analysis {
16 public:
17
18 /// Constructor
19 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2018_I1711223);
20
21 /// @name Analysis methods
22 //@{
23
24 /// Book histograms and initialise projections before the run
25 void init() {
26
27 // Get photons to dress leptons
28 PromptFinalState photons(Cuts::abspid == PID::PHOTON);
29
30 // Electrons and muons in Fiducial PS
31 PromptFinalState leptons(Cuts::abspid == PID::ELECTRON || Cuts::abspid == PID::MUON);
32 leptons.acceptTauDecays(false);
33 DressedLeptons dressedleptons(photons, leptons, 0.1, Cuts::open(), true); // useDecayPhotons=true -- useJetClustering? auto-set to false?
34 declare(dressedleptons, "DressedLeptons");
35
36 // Prompt neutrinos (yikes!)
37 IdentifiedFinalState nu_id;
38 nu_id.acceptNeutrinos();
39 PromptFinalState neutrinos(nu_id);
40 neutrinos.acceptTauDecays(false);
41 declare(neutrinos, "Neutrinos");
42 MSG_WARNING("\033[91;1mLIMITED VALIDITY - check info file for details!\033[m");
43
44 //Jets
45
46 // Muons
47 PromptFinalState bare_mu(Cuts::abspid == PID::MUON, true); // true = use muons from prompt tau decays
48 DressedLeptons all_dressed_mu(photons, bare_mu, 0.1, Cuts::abseta < 5.0, true);
49
50 // Electrons
51 PromptFinalState bare_el(Cuts::abspid == PID::ELECTRON, true); // true = use electrons from prompt tau decays
52 DressedLeptons all_dressed_el(photons, bare_el, 0.1, Cuts::abseta < 5.0, true);
53
54 //Jet forming
55 VetoedFinalState vfs(FinalState(Cuts::abseta < 5));
56 vfs.addVetoOnThisFinalState(all_dressed_el);
57 vfs.addVetoOnThisFinalState(all_dressed_mu);
58
59 FastJets jets(vfs, FastJets::ANTIKT, 0.4, JetAlg::Muons::ALL, JetAlg::Invisibles::DECAY);
60 declare(jets, "Jets");
61
62 // Book auxiliary histograms
63 book(_h["MTWZ"], "_mTWZ", refData( 6, 1, 1));
64 book(_h["sumpt"], "_sumpT", refData( 8, 1, 1));
65 book(_h["dphiWZ"], "_dphiWZ", refData(10, 1, 1));
66 book(_h["Njets_VBS"], "_njets", refData(12, 1, 1));
67 book(_h["mjj"], "_mjj", refData(14, 1, 1));
68 book(_h["dyjj"], "_dRapjj", refData(16, 1, 1));
69 book(_h["dphijj"], "_dPhijj", refData(18, 1, 1));
70 book(_h["Njets_gap"], "_gapJets", refData(20, 1, 1));
71
72 // book output bar charts
73 book(_s["MTWZ"], 6, 1, 1);
74 book(_s["sumpt"], 8, 1, 1);
75 book(_s["dphiWZ"], 10, 1, 1);
76 book(_s["Njets_VBS"], 12, 1, 1);
77 book(_s["mjj"], 14, 1, 1);
78 book(_s["dyjj"], 16, 1, 1);
79 book(_s["dphijj"], 18, 1, 1);
80 book(_s["Njets_gap"], 20, 1, 1);
81
82 }
83
84
85 void analyze(const Event& event) {
86
87 const Particles& dressedleptons = apply<DressedLeptons>(event, "DressedLeptons").particlesByPt();
88 const Particles& neutrinos = apply<PromptFinalState>(event, "Neutrinos").particlesByPt();
89 Jets jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 4.5);
90
91 int i, j, k;
92 double MassZ01 = 0., MassZ02 = 0., MassZ12 = 0.;
93 double MassW0 = 0., MassW1 = 0., MassW2 = 0.;
94 double WeightZ1, WeightZ2, WeightZ3;
95 double WeightW1, WeightW2, WeightW3;
96 double M1, M2, M3;
97 double WeightTotal1, WeightTotal2, WeightTotal3;
98
99 //---Fiducial PS: assign leptons to W and Z bosons using Resonant shape algorithm
100 if (dressedleptons.size() < 3 || neutrinos.size() < 1) vetoEvent;
101
102 //--- count num of electrons and muons
103 int Nel = 0, Nmu = 0;
104 for (const Particle& l : dressedleptons) {
105 if (l.abspid() == 11) ++Nel;
106 if (l.abspid() == 13) ++Nmu;
107 }
108
109 int icomb=0;
110 // try Z pair of leptons 01
111 if ( (dressedleptons[0].pid() ==-(dressedleptons[1].pid())) && (dressedleptons[2].pid()*neutrinos[0].pid()< 0) && (dressedleptons[2].abspid()==neutrinos[0].abspid()-1)) {
112 MassZ01 = (dressedleptons[0].momentum() + dressedleptons[1].momentum()).mass();
113 MassW2 = (dressedleptons[2].momentum() + neutrinos[0].momentum()).mass();
114 icomb = 1;
115 }
116
117 // try Z pair of leptons 02
118 if ( (dressedleptons[0].pid()==-(dressedleptons[2].pid())) && (dressedleptons[1].pid()*neutrinos[0].pid()< 0) && (dressedleptons[1].abspid()==neutrinos[0].abspid()-1)) {
119 MassZ02 = (dressedleptons[0].momentum() + dressedleptons[2].momentum()).mass();
120 MassW1 = (dressedleptons[1].momentum() + neutrinos[0].momentum()).mass();
121 icomb = 2;
122 }
123 // try Z pair of leptons 12
124 if ( (dressedleptons[1].pid()==-(dressedleptons[2].pid())) && (dressedleptons[0].pid()*neutrinos[0].pid()< 0) && (dressedleptons[0].abspid()==neutrinos[0].abspid()-1)) {
125 MassZ12 = (dressedleptons[1].momentum() + dressedleptons[2].momentum()).mass();
126 MassW0 = (dressedleptons[0].momentum() + neutrinos[0].momentum()).mass();
127 icomb = 3;
128 }
129
130 if (icomb<=0) vetoEvent;
131
132
133 WeightZ1 = 1/(pow(MassZ01*MassZ01 - MZ_PDG*MZ_PDG,2) + pow(MZ_PDG*GammaZ_PDG,2));
134 WeightW1 = 1/(pow(MassW2*MassW2 - MW_PDG*MW_PDG,2) + pow(MW_PDG*GammaW_PDG,2));
135 WeightTotal1 = WeightZ1*WeightW1;
136 M1 = -1*WeightTotal1;
137
138 WeightZ2 = 1/(pow(MassZ02*MassZ02- MZ_PDG*MZ_PDG,2) + pow(MZ_PDG*GammaZ_PDG,2));
139 WeightW2 = 1/(pow(MassW1*MassW1- MW_PDG*MW_PDG,2) + pow(MW_PDG*GammaW_PDG,2));
140 WeightTotal2 = WeightZ2*WeightW2;
141 M2 = -1*WeightTotal2;
142
143 WeightZ3 = 1/(pow(MassZ12*MassZ12 - MZ_PDG*MZ_PDG,2) + pow(MZ_PDG*GammaZ_PDG,2));
144 WeightW3 = 1/(pow(MassW0*MassW0 - MW_PDG*MW_PDG,2) + pow(MW_PDG*GammaW_PDG,2));
145 WeightTotal3 = WeightZ3*WeightW3;
146 M3 = -1*WeightTotal3;
147
148 if( (M1 < M2 && M1 < M3) || (MassZ01 != 0 && MassW2 != 0 && MassZ02 == 0 && MassZ12 == 0) ) {
149 i = 0; j = 1; k = 2;
150 }
151 if((M2 < M1 && M2 < M3) || (MassZ02 != 0 && MassW1 != 0 && MassZ01 == 0 && MassZ12 == 0) ) {
152 i = 0; j = 2; k = 1;
153 }
154 if((M3 < M1 && M3 < M2) || (MassZ12 != 0 && MassW0 != 0 && MassZ01 == 0 && MassZ02 == 0) ) {
155 i = 1; j = 2; k = 0;
156 }
157
158 FourMomentum Zlepton1 = dressedleptons[i].momentum();
159 FourMomentum Zlepton2 = dressedleptons[j].momentum();
160 FourMomentum Wlepton = dressedleptons[k].momentum();
161 FourMomentum Zboson = dressedleptons[i].momentum()+dressedleptons[j].momentum();
162 FourMomentum Wboson = dressedleptons[k].momentum()+neutrinos[0].momentum();
163
164 double cosLepNeut;
165 double Wboson_mT = 0.;
166 double norm = Wlepton.pT() * neutrinos[0].pt();
167 if(norm != 0 ) {
168 cosLepNeut = ( Wlepton.px()*neutrinos[0].px() + Wlepton.py()*neutrinos[0].py() )/norm ;
169 if (1-cosLepNeut >= 0. ) Wboson_mT = sqrt( 2 * Wlepton.pT() * neutrinos[0].pt() * (1-cosLepNeut ) );
170 }
171
172 //---- CUTS (based on Table 1 WZ: 36.1 fb-1)----//
173 if (Wlepton.pT() <= 20*GeV || Zlepton1.pT() <= 15*GeV || Zlepton2.pT() <= 15*GeV) vetoEvent;
174 if (Wlepton.abseta() >= 2.5 || Zlepton1.abseta() >= 2.5 || Zlepton2.abseta() >= 2.5) vetoEvent;
175 if (fabs(Zboson.mass()/GeV - MZ_PDG) >= 10.) vetoEvent;
176 if (Wboson_mT <= 30*GeV) vetoEvent;
177 if (deltaR(Zlepton1, Zlepton2) <= 0.2) vetoEvent;
178 if (deltaR(Zlepton1, Wlepton) <= 0.3) vetoEvent;
179 if (deltaR(Zlepton2, Wlepton) <= 0.3) vetoEvent;
180
181 double WZ_pt = (Zlepton1.pt() + Zlepton2.pt() + Wlepton.pt() + neutrinos[0].pt())/GeV;
182 double WZ_px = (Zlepton1.px() + Zlepton2.px() + Wlepton.px() + neutrinos[0].px())/GeV;
183 double WZ_py = (Zlepton1.py() + Zlepton2.py() + Wlepton.py() + neutrinos[0].py())/GeV;
184 double mTWZ = sqrt( pow(WZ_pt, 2) - ( pow(WZ_px, 2) + pow(WZ_py,2) ) );
185 double sumptleptons = (Zlepton1.pt() + Zlepton2.pt() + Wlepton.pt())/GeV;
186 double dPhiWZTruth = acos(cos(Zboson.phi()-Wboson.phi()));
187
188
189
190 //---- Jet CUTS----//
191 ifilter_discard(jets, [&](const Jet& j) {
192 return deltaR(j, Zlepton1) < 0.3 || deltaR(j, Zlepton2) < 0.3 || deltaR(j, Wlepton) < 0.3;
193 });
194 if (jets.size() < 2) vetoEvent;
195 if (jets[0].pT() < 40*GeV) vetoEvent;
196
197 // Selection of the second jet as the second highest pT jet and in opposite hemisphere with the fisrt jet
198 FourMomentum jet_lead = jets[0].mom();
199 FourMomentum jet_sublead;
200 bool foundVBSJetPair = false;
201 for (const Jet& jet : jets) {
202 if(jet.pT() > 40*GeV && jet.eta()*jets[0].eta() < 0.) {
203 jet_sublead = jet.mom();
204 foundVBSJetPair = true;
205 break;
206 }
207 }
208 if (!foundVBSJetPair) vetoEvent;
209
210 const double mJJ = (jet_lead + jet_sublead).mass()/GeV;
211 const double dphi_jj = acos(cos(jet_lead.phi() - jet_sublead.phi()));
212 const double dyjj = fabs(jet_lead.rap() - jet_sublead.rap());
213
214 //Plots in the SR
215 if (mJJ < 500*GeV) vetoEvent;
216
217 const size_t njets40 = filter_select(jets, Cuts::pT > 40*GeV).size();
218 fillWithOverflow("Njets_VBS", njets40, 5.1);
219
220 const double y_min = std::min(jet_lead.rap(), jet_sublead.rap());
221 const double y_max = std::max(jet_lead.rap(), jet_sublead.rap());
222 const size_t njetsGap = count(jets, [&](const Jet& j) {
223 return (j.rap() > y_min && j.rap() < y_max);
224 });
225 fillWithOverflow("Njets_gap", njetsGap, 3.1);
226
227 fillWithOverflow("MTWZ", mTWZ, 551);
228 fillWithOverflow("sumpt", sumptleptons, 501);
229 fillWithOverflow("mjj", mJJ, 2001);
230
231 _h["dphiWZ"]->fill(dPhiWZTruth);
232 _h["dyjj"]->fill(dyjj);
233 _h["dphijj"]->fill(dphi_jj);
234
235 }
236
237
238 void fillWithOverflow(const string& tag, const double value, const double overflow) {
239 if (value < overflow) _h[tag]->fill(value);
240 else _h[tag]->fill(overflow);
241 }
242
243
244 void finalize() {
245
246 scale(_h, crossSectionPerEvent() / femtobarn);
247 // unfortunately, no differential cross-sections were measured in this analysis
248 for (auto &item : _h) barchart(item.second, _s[item.first]);
249
250 }
251
252
253 //@}
254
255 private:
256
257
258 /// @name Histograms
259 //@{
260
261 map<string, Histo1DPtr> _h;
262 map<string, Scatter2DPtr> _s;
263
264 //@}
265
266 double MZ_PDG = 91.1876;
267 double MW_PDG = 80.385;
268 double GammaZ_PDG = 2.4952;
269 double GammaW_PDG = 2.085;
270
271 };
272
273 // The hook for the plugin system
274 RIVET_DECLARE_PLUGIN(ATLAS_2018_I1711223);
275}
|