Rivet analyses referenceATLAS_2017_I1517194Electroweak $Wjj$ production at 8 TeVExperiment: ATLAS (LHC) Inspire ID: 1517194 Status: VALIDATED Authors:
Beam energies: (4000.0, 4000.0) GeV Run details:
Measurements of the electroweak production of a $W$ boson in association with two jets at high dijet invariant mass are performed using $\sqrt{s} = 7$ and 8 TeV proton-proton collision data produced by the Large Hadron Collider, corresponding respectively to 4.7 and 20.2 fb$^{-1}$ of integrated luminosity collected by the ATLAS detector. The measurements are sensitive to the production of a $W$ boson via a triple-gauge-boson vertex and include both the fiducial and differential cross sections of the electroweak process. Source code: ATLAS_2017_I1517194.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/PromptFinalState.hh"
4#include "Rivet/Projections/MissingMomentum.hh"
5#include "Rivet/Projections/FastJets.hh"
6#include "Rivet/Projections/LeptonFinder.hh"
7
8namespace Rivet {
9
10
11 namespace FOO {
12 double mass(const Particle&, const P4&) { return 0.0; }
13 // double mass(const P4&, double) { return 0.0; }
14 }
15
16 /// @brief Electroweak Wjj production at 8 TeV
17 class ATLAS_2017_I1517194 : public Analysis {
18 public:
19
20 /// Constructor
21 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2017_I1517194);
22
23
24 /// @name Analysis methods
25 /// @{
26
27 /// Book histograms and initialise projections before the run
28 void init() {
29
30 // Get options from the option system
31 _mode = 0;
32 if ( getOption("LMODE") == "EL" ) _mode = 0;
33 if ( getOption("LMODE") == "MU" ) _mode = 1;
34
35 // W Selection
36 MissingMomentum missmom(FinalState(Cuts::eta < 5.0));
37 declare(missmom, "MET");
38 LeptonFinder lf(0.1, Cuts::rap < 2.5 && Cuts::pT > 25*GeV &&
39 Cuts::abspid == (_mode? PID::MUON : PID::ELECTRON));
40 declare(lf, "Leptons");
41
42 VetoedFinalState vfs;
43 vfs.vetoFinalState(lf);
44 FastJets jets(vfs, JetAlg::ANTIKT, 0.4, JetMuons::DECAY, JetInvisibles::ALL); //< !!
45 declare(jets, "Jets_w");
46
47 const vector<string> phase_spaces = { "highmass15", "antiLC", "signal10",
48 "highmass10", "inclusive", "highmass20",
49 "antiLCantiJC", "antiJC", "signal", };
50 const vector<string> variables = { "dijetmass", "dijetpt", "dphi12", "dy12", "j1pt", "JC", "LC", "ngapjets" };
51 size_t hepdataid = 10;
52 for (size_t group = 0; group < 4; ++group) {
53 for ( size_t ps=0; ps < phase_spaces.size(); ++ps ) {
54 for ( size_t var=0; var < variables.size(); ++var ) {
55 if (group < 2) {
56 if ((ps == 0 || ps == 2 || ps == 3 || ps == 5) && var == 0) continue;
57 if ((ps == 1 || ps == 2 || ps > 5) && var > 4) continue;
58 if (group == 1 && ps == 7 && var == 3) continue;
59 }
60 else {
61 if (ps == 1 || ps == 4 || ps > 5) continue;
62 if ((ps == 0 || ps == 5) && var < 2) continue;
63 if (ps == 2 && var > 4) continue;
64 if ((ps == 0 || ps == 3 || ps == 5) && var == 5) continue;
65 if (group == 2) {
66 if ((ps == 0 || ps == 5) && var == 3) continue;
67 if (ps == 3 && var == 1) continue;
68 }
69 else {
70 if ((ps == 0 || ps == 3 || ps == 5) && var == 6) continue;
71 if ((ps == 2 || ps == 3) && var < 2) continue;
72 }
73 }
74
75 ++hepdataid;
76
77 string label = variables[var]+"_"+phase_spaces[ps];
78 //string pre = group > 1? "ew_" : "";
79 //std::cout << "rivet -- " << pre << label << suff << " " << hepdataid << std::endl;
80 string suff = group % 2? "" : "_norm";
81 if (group > 1) book(_hists["ew_" + label + suff], hepdataid, 1, 1);
82 else book(_hists[label + suff], hepdataid, 1, 1);
83 }
84 }
85 }
86 }
87
88
89 /// Perform the per-event analysis
90 void analyze(const Event& event) {
91
92 const P4& pmiss = apply<MissingMomentum>(event, "MET").missingMom();
93 const Particles& ls = apply<LeptonFinder>(event, "Leptons").particles();
94 const int ifound = closestMatchIndex(ls, pmiss, FOO::mass, 80.4*GeV);
95
96 if (ifound < 0) vetoEvent;
97 const FourMomentum lepton(ls[ifound]), neutrino(pmiss);
98
99 const FastJets& jetpro = apply<FastJets>(event, "Jets_w");
100 vector<FourMomentum> jets;
101 for (const Jet& jet : jetpro.jetsByPt(Cuts::pT > 30*GeV && Cuts::absrap < 4.4)) {
102 if ( fabs(deltaR(jet, lepton)) < 0.3 ) continue;
103 jets.push_back(jet.momentum());
104 }
105
106 if (jets.size() < 2) vetoEvent;
107
108 double mT = sqrt( 2.0*lepton.pT()*neutrino.Et()*( 1.0-cos( lepton.phi()-neutrino.phi() ) ) );
109 double DeltaRap = fabs( jets[0].rapidity()-jets[1].rapidity() );
110 double dijet_mass = FourMomentum(jets[0]+jets[1]).mass();
111 size_t nojets = jets.size();
112
113 if (pmiss.pT() < 25*GeV) vetoEvent;
114 if (jets[0].pT() < 80*GeV) vetoEvent;
115 if (jets[1].pT() < 60*GeV) vetoEvent;
116 if (dijet_mass < 500*GeV) vetoEvent;
117 if (DeltaRap < 2.0) vetoEvent;
118 if (mT < 40*GeV) vetoEvent;
119
120 // By now, the event has passed all VBF cuts
121 double DijetPt = FourMomentum(jets[0]+jets[1]).pT();
122 double dphi = fabs(jets[0].phi() - jets[1].phi());
123 double DeltaPhi = ( dphi<=pi ) ? dphi/pi : (2.*pi-dphi)/pi;
124 double jet_1_rap = jets[0].rapidity();
125 double jet_2_rap = jets[1].rapidity();
126
127 // Njets in Gap Control Regions info
128 int njetsingap = 0;
129 bool firstIsForward = jet_1_rap > jet_2_rap;
130 int jF = (firstIsForward) ? 0 : 1; // sets most positive jet to be forward index
131 int jB = (firstIsForward) ? 1 : 0; // sets most negative jet to be backward index
132 for (size_t j = 2; j < nojets; ++j) {
133 if ( (jets[j].rapidity()<jets[jF].rapidity()) && (jets[j].rapidity()>jets[jB].rapidity()) ) njetsingap++;
134 }
135
136 // Third+ Jet Centrality Cut (Vetos any event that has a jet between raps of the two leading jets)
137 bool passJC = false;
138 std::vector<double> JCvals;
139 JCvals.clear();
140 if ( nojets < 3) passJC = true;
141 else if ( nojets > 2 ) {
142 passJC = true;
143 for (size_t j = 2; j < nojets; ++j) {
144 double jet_3_rap = jets[j].rapidity();
145 double jet3gap = fabs(( jet_3_rap - ((jet_1_rap + jet_2_rap)/2.0))/(jet_1_rap - jet_2_rap));
146 JCvals.push_back(jet3gap);
147 if ( jet3gap < 0.4 ) { passJC = false; }
148 }
149 }
150
151 double lepton_rap = lepton.rapidity();
152 double lep_cent = fabs((lepton_rap - ((jet_1_rap + jet_2_rap)/2.0) )/(jet_1_rap - jet_2_rap));
153 bool passLC = (lep_cent < 0.4);
154
155 map<string, bool> phaseSpaces;
156 phaseSpaces["inclusive"] = true;
157 phaseSpaces["highmass10"] = (dijet_mass>1000.0*GeV);
158 phaseSpaces["highmass15"] = (dijet_mass>1500.0*GeV);
159 phaseSpaces["highmass20"] = (dijet_mass>2000.0*GeV);
160 phaseSpaces["antiLC"] = ( !passLC && passJC );
161 phaseSpaces["antiJC"] = ( passLC && !passJC );
162 phaseSpaces["antiLCantiJC"] = ( !passLC && !passJC );
163 phaseSpaces["signal"] = ( passLC && passJC );
164 phaseSpaces["signal10"] = ( (dijet_mass>1000.0*GeV) && passLC && passJC );
165
166 for (const auto& ps : phaseSpaces) {
167 if (!ps.second) continue;
168 fillHisto("dijetmass_"+ps.first, dijet_mass);
169 fillHisto("dijetpt_"+ps.first, DijetPt);
170 fillHisto("dy12_"+ps.first, fabs(jet_1_rap-jet_2_rap));
171 fillHisto("dphi12_"+ps.first, DeltaPhi);
172 fillHisto("j1pt_"+ps.first, jets[0].pT());
173 if (ps.first == "inclusive" || ps.first.find("highmass") != string::npos) {
174 fillHisto("LC_"+ps.first, lep_cent);
175 fillHisto("ngapjets_"+ps.first, njetsingap);
176 for (auto& jc : JCvals) {
177 fillHisto("JC_"+ps.first, jc);
178 }
179 }
180 }
181
182 }
183
184
185 /// Normalise histograms etc., after the run
186 void finalize() {
187 double factor = crossSection()/sumOfWeights()/femtobarn; // refData is in fb
188 for (auto& key_hist : _hists) {
189 scale(key_hist.second, factor);
190 if (key_hist.first.find("_norm") != string::npos) normalize(key_hist.second);
191 }
192 }
193
194 /// @}
195
196
197 void fillHisto(const string& label, const double value) {
198 if (_hists.find(label) != _hists.end()) { // QCD+EW, absolute
199 _hists[label]->fill(value);
200 }
201 if (_hists.find(label + "_norm") != _hists.end()) { // QCD+EW, normalised
202 _hists[label + "_norm"]->fill(value);
203 }
204 if (_hists.find("ew_" + label) != _hists.end()) { // EW-only, absolute
205 _hists["ew_" + label]->fill(value);
206 }
207 if (_hists.find("ew_" + label + "_norm") != _hists.end()) { // EW-only, normalised
208 _hists["ew_" + label + "_norm"]->fill(value);
209 }
210 }
211
212
213 protected:
214
215 size_t _mode;
216
217 private:
218
219 map<string, Histo1DPtr> _hists;
220
221 };
222
223
224 RIVET_DECLARE_PLUGIN(ATLAS_2017_I1517194);
225
226}
|