Rivet analyses referenceATLAS_2016_I14582700-lepton SUSY search with 3.2/fb of 13 TeV $pp$ dataExperiment: ATLAS (LHC) Inspire ID: 1458270 Status: VALIDATED, UNOFFICIAL Authors:
Beam energies: (6500.0, 6500.0) GeV Run details:
ATLAS 0-lepton SUSY search using 3.2/fb of LHC $pp$ data at 13 TeV, recorded in 2015. The event selection is via signal regions requiring from 2-6 high-energy jets, and significant missing transverse energy. Detailed info: http://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PAPERS/SUSY-2015-06/ Source code: ATLAS_2016_I1458270.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/PromptFinalState.hh"
5#include "Rivet/Projections/FastJets.hh"
6#include "Rivet/Projections/Sphericity.hh"
7#include "Rivet/Projections/SmearedParticles.hh"
8#include "Rivet/Projections/SmearedJets.hh"
9#include "Rivet/Projections/SmearedMET.hh"
10#include "Rivet/Tools/Cutflow.hh"
11
12namespace Rivet {
13
14
15 /// @brief ATLAS 0-lepton SUSY search with 3.2/fb of 13 TeV pp data
16 class ATLAS_2016_I1458270 : public Analysis {
17 public:
18
19 /// Constructor
20 RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2016_I1458270);
21
22
23 /// @name Analysis methods
24 /// @{
25
26 /// Book histograms and initialise projections before the run
27 void init() {
28
29 // Initialise and register projections
30 FinalState calofs(Cuts::abseta < 4.8);
31 FastJets fj(calofs, JetAlg::ANTIKT, 0.4);
32 declare(fj, "TruthJets");
33 declare(SmearedJets(fj, JET_SMEAR_ATLAS_RUN2, JET_BTAG_ATLAS_RUN2_MV2C20), "RecoJets");
34
35 MissingMomentum mm(calofs);
36 declare(mm, "TruthMET");
37 declare(SmearedMET(mm, MET_SMEAR_ATLAS_RUN2), "RecoMET");
38
39 PromptFinalState es(Cuts::abseta < 2.47 && Cuts::abspid == PID::ELECTRON, TauDecaysAs::PROMPT, MuDecaysAs::PROMPT);
40 declare(es, "TruthElectrons");
41 declare(SmearedParticles(es, ELECTRON_RECOEFF_ATLAS_RUN2, ELECTRON_SMEAR_ATLAS_RUN2), "RecoElectrons");
42
43 PromptFinalState mus(Cuts::abseta < 2.7 && Cuts::abspid == PID::MUON, TauDecaysAs::PROMPT);
44 declare(mus, "TruthMuons");
45 declare(SmearedParticles(mus, MUON_EFF_ATLAS_RUN2, MUON_SMEAR_ATLAS_RUN2), "RecoMuons");
46
47
48 // Book histograms/counters
49 book(_h_2jl, "2jl");
50 book(_h_2jm, "2jm");
51 book(_h_2jt, "2jt");
52 book(_h_4jt, "4jt");
53 book(_h_5j , "5j");
54 book(_h_6jm, "6jm");
55 book(_h_6jt, "6jt");
56
57 book(_hMeff_2jl, 4,1,1);
58 book(_hMeff_2jm, 5,1,1);
59 book(_hMeff_2jt, 6,1,1);
60 book(_hMeff_4jt, 7,1,1);
61 book(_hMeff_5j , 8,1,1);
62 book(_hMeff_6jm, 9,1,1);
63 book(_hMeff_6jt, 10,1,1);
64
65 book(_h_temp_Meff_2jl, "_temp_Meff_2jl",refData( 4,1,1));
66 book(_h_temp_Meff_2jm, "_temp_Meff_2jm",refData( 5,1,1));
67 book(_h_temp_Meff_2jt, "_temp_Meff_2jt",refData( 6,1,1));
68 book(_h_temp_Meff_4jt, "_temp_Meff_4jt",refData( 7,1,1));
69 book(_h_temp_Meff_5j , "_temp_Meff_5j" ,refData( 8,1,1));
70 book(_h_temp_Meff_6jm, "_temp_Meff_6jm",refData( 9,1,1));
71 book(_h_temp_Meff_6jt, "_temp_Meff_6jt",refData(10,1,1));
72
73
74
75
76 // Book cut-flows
77 const vector<string> cuts2j = {"Pre-sel+MET+pT1", "Njet", "Dphi_min(j,MET)",
78 "pT2", "MET/sqrtHT", "m_eff(incl)"};
79 const vector<string> cutsXj = {"Pre-sel+MET+pT1", "Njet", "Dphi_min(j,MET)",
80 "pT2", "pT4", "Aplanarity", "MET/m_eff(Nj)", "m_eff(incl)"};
81 book(_flows, {"CF-2jl", "CF-2jm", "CF-2jt", "CF-4jt", "CF-5j", "CF-6jm", "CF-6jt"});
82 for (auto& b : _flows->bins()) {
83 if (b.index() < 4) book(b, b.xEdge(), cuts2j);
84 else book(b, b.xEdge(), cutsXj);
85 }
86 }
87
88 /// Perform the per-event analysis
89 void analyze(const Event& event) {
90
91 _flows->groupfillinit();
92
93 // Same MET cut for all signal regions
94 //const Vector3 vmet = -apply<MissingMomentum>(event, "TruthMET").vectorEt();
95 const Vector3 vmet = -apply<SmearedMET>(event, "RecoMET").vectorEt();
96 const double met = vmet.mod();
97 if (met < 200*GeV) vetoEvent;
98
99 // Get baseline electrons, muons, and jets
100 Particles elecs = apply<ParticleFinder>(event, "RecoElectrons").particles(Cuts::pT > 10*GeV);
101 Particles muons = apply<ParticleFinder>(event, "RecoMuons").particles(Cuts::pT > 10*GeV);
102 Jets jets = apply<JetFinder>(event, "RecoJets").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 2.8); ///< @todo Pile-up subtraction
103
104 // Jet/electron/muons overlap removal and selection
105 // Remove any |eta| < 2.8 jet within dR = 0.2 of a baseline electron
106 for (const Particle& e : elecs)
107 idiscard(jets, deltaRLess(e, 0.2, RAPIDITY));
108 // Remove any electron or muon with dR < 0.4 of a remaining (Nch > 3) jet
109 for (const Jet& j : jets) {
110 /// @todo Add track efficiency random ing
111 idiscard(elecs, deltaRLess(j, 0.4, RAPIDITY));
112 if (j.particles(Cuts::abscharge > 0 && Cuts::pT > 500*MeV).size() >= 3)
113 idiscard(muons, deltaRLess(j, 0.4, RAPIDITY));
114 }
115 // Discard the softer of any electrons within dR < 0.05
116 for (size_t i = 0; i < elecs.size(); ++i) {
117 const Particle& e1 = elecs[i];
118 /// @todo Would be nice to pass a "tail view" for the filtering, but awkward without range API / iterator guts
119 idiscard(elecs, [&](const Particle& e2){ return e2.pT() < e1.pT() && deltaR(e1,e2) < 0.05; });
120 }
121
122 // Loose electron selection
123 iselect(elecs, ParticleEffFilter(ELECTRON_EFF_ATLAS_RUN2_LOOSE));
124
125 // Veto the event if there are any remaining baseline leptons
126 if (!elecs.empty()) vetoEvent;
127 if (!muons.empty()) vetoEvent;
128
129 // Signal jets have pT > 50 GeV
130 const Jets jets50 = select(jets, Cuts::pT > 50*GeV);
131 if (jets50.size() < 2) vetoEvent;
132 vector<double> jetpts; transform(jets, jetpts, pT);
133 vector<double> jetpts50; transform(jets50, jetpts50, pT);
134 const double j1pt = jetpts50[0];
135 const double j2pt = jetpts50[1];
136 if (j1pt < 200*GeV) vetoEvent;
137
138 // Construct multi-jet observables
139 const double ht = sum(jetpts, 0.0);
140 const double met_sqrt_ht = met / sqrt(ht);
141 const double meff_incl = sum(jetpts50, met);
142
143 // Get dphis between MET and jets
144 vector<double> dphimets50; transform(jets50, dphimets50, deltaPhiWRT(vmet));
145 const double min_dphi_met_3 = min(head(dphimets50, 3));
146 MSG_DEBUG(dphimets50 << ", " << min_dphi_met_3);
147
148 // Jet aplanarity
149 Sphericity sph; sph.calc(jets);
150 const double aplanarity = sph.aplanarity();
151
152
153 // Fill SR counters
154 // 2-jet SRs
155 if (_flows->fillnext("CF-2jl", {true, true, min_dphi_met_3 > 0.8, j2pt > 200*GeV,
156 met_sqrt_ht > 15*sqrt(GeV), meff_incl > 1200*GeV})) _h_2jl->fill();
157 if (_flows->fillnext("CF-2jm", {j1pt > 300*GeV, true, min_dphi_met_3 > 0.4, j2pt > 50*GeV,
158 met_sqrt_ht > 15*sqrt(GeV), meff_incl > 1600*GeV})) _h_2jm->fill();
159 if (_flows->fillnext("CF-2jt", {true, true, min_dphi_met_3 > 0.8, j2pt > 200*GeV,
160 met_sqrt_ht > 20*sqrt(GeV), meff_incl > 2000*GeV})) _h_2jt->fill();
161
162 // Fill SR Meff Histo1Ds
163 // 2-jet SRs
164 if ((min_dphi_met_3 > 0.8) && (j2pt > 200*GeV) && (met_sqrt_ht > 15*sqrt(GeV))) _h_temp_Meff_2jl->fill(meff_incl);
165 if ((j1pt > 300*GeV) && (min_dphi_met_3 > 0.4) && (j2pt > 50*GeV) && (met_sqrt_ht > 15*sqrt(GeV))) _h_temp_Meff_2jm->fill(meff_incl);
166 if ((min_dphi_met_3 > 0.8) && (j2pt > 200*GeV ) && (met_sqrt_ht > 20*sqrt(GeV))) _h_temp_Meff_2jt->fill(meff_incl);
167
168 // Upper multiplicity SRs
169 const double j4pt = jets50.size() > 3 ? jetpts50[3] : -1;
170 const double j5pt = jets50.size() > 4 ? jetpts50[4] : -1;
171 const double j6pt = jets50.size() > 5 ? jetpts50[5] : -1;
172 const double meff_4 = jets50.size() > 3 ? sum(head(jetpts50, 4), met) : -1;
173 const double meff_5 = jets50.size() > 4 ? meff_4 + jetpts50[4] : -1;
174 const double meff_6 = jets50.size() > 5 ? meff_5 + jetpts50[5] : -1;
175 const double met_meff_4 = met / meff_4;
176 const double met_meff_5 = met / meff_5;
177 const double met_meff_6 = met / meff_6;
178 const double min_dphi_met_more = jets50.size() > 3 ? min(tail(dphimets50, -3)) : -1;
179
180
181 if (_flows->fillnext("CF-4jt", {true, jets50.size() >= 4, min_dphi_met_3 > 0.4 && min_dphi_met_more > 0.2,
182 jetpts[1] > 100*GeV, j4pt > 100*GeV, aplanarity > 0.04, met_meff_4 > 0.20, meff_incl > 2200*GeV}))
183 _h_4jt->fill();
184 if (_flows->fillnext("CF-5j", {true, jets50.size() >= 5, min_dphi_met_3 > 0.4 && min_dphi_met_more > 0.2,
185 jetpts[1] > 100*GeV, j4pt > 100*GeV && j5pt > 50*GeV, aplanarity > 0.04, met_meff_5 > 0.25, meff_incl > 1600*GeV}))
186 _h_5j->fill();
187 if (_flows->fillnext("CF-6jm", {true, jets50.size() >= 6, min_dphi_met_3 > 0.4 && min_dphi_met_more > 0.2,
188 jetpts[1] > 100*GeV, j4pt > 100*GeV && j6pt > 50*GeV, aplanarity > 0.04, met_meff_6 > 0.25, meff_incl > 1600*GeV}))
189 _h_6jm->fill();
190 if (_flows->fillnext("CF-6jt", {true, jets50.size() >= 6, min_dphi_met_3 > 0.4 && min_dphi_met_more > 0.2,
191 jetpts[1] > 100*GeV, j4pt > 100*GeV && j6pt > 50*GeV, aplanarity > 0.04, met_meff_6 > 0.20, meff_incl > 2000*GeV}))
192 _h_6jt->fill();
193
194 // Fill SR Meff Histo1Ds
195 // Upper multiplicity SRs
196 if (((jets50.size() >= 4) && (min_dphi_met_3 > 0.4) && (min_dphi_met_more > 0.2) && (jetpts[1] > 100*GeV) && (j4pt > 100*GeV) && (aplanarity > 0.04) && (met_meff_4 > 0.20))) _h_temp_Meff_4jt->fill(meff_incl);
197 if (((jets50.size() >= 5) && (min_dphi_met_3 > 0.4) && (min_dphi_met_more > 0.2 ) &&
198 (jetpts[1] > 100*GeV) && (j4pt > 100*GeV) && (j5pt > 50*GeV) && (aplanarity > 0.04) && (met_meff_5 > 0.25))) _h_temp_Meff_5j->fill(meff_incl);
199 if (((jets50.size() >= 6) && (min_dphi_met_3 > 0.4) && (min_dphi_met_more > 0.2) &&
200 (jetpts[1] > 100*GeV) && (j4pt > 100*GeV) && (j6pt > 50*GeV) && (aplanarity > 0.04) && (met_meff_6 > 0.25))) _h_temp_Meff_6jm->fill(meff_incl);
201 if (((jets50.size() >= 6) && (min_dphi_met_3 > 0.4) && (min_dphi_met_more > 0.2) &&
202 (jetpts[1] > 100*GeV) && (j4pt > 100*GeV) && (j6pt > 50*GeV) && (aplanarity > 0.04) && (met_meff_6 > 0.20))) _h_temp_Meff_6jt->fill(meff_incl);
203
204 }
205
206
207 /// Normalise histograms etc., after the run
208 void finalize() {
209
210 const double sf = 3.2*crossSection()/femtobarn/sumOfWeights();
211 scale(_h_2jl, sf); scale(_h_2jm, sf); scale(_h_2jt, sf);
212 scale(_h_4jt, sf); scale(_h_5j, sf);
213 scale(_h_6jm, sf); scale(_h_6jt, sf);
214
215 scale(_h_temp_Meff_2jl, sf); scale(_h_temp_Meff_2jm, sf); scale(_h_temp_Meff_2jt, sf);
216 scale(_h_temp_Meff_4jt, sf); scale(_h_temp_Meff_5j, sf);
217 scale(_h_temp_Meff_6jm, sf); scale(_h_temp_Meff_6jt, sf);
218
219
220 // the HEPData entry corresponding to this does not divide their distributions
221 // by bin width... so to avoid this we need to convert to Estimate1D which is not divided by bw
222 barchart(_h_temp_Meff_2jl,_hMeff_2jl);
223 barchart(_h_temp_Meff_2jm,_hMeff_2jm);
224 barchart(_h_temp_Meff_2jt,_hMeff_2jt);
225 barchart(_h_temp_Meff_4jt,_hMeff_4jt);
226 barchart(_h_temp_Meff_5j ,_hMeff_5j );
227 barchart(_h_temp_Meff_6jm,_hMeff_6jm);
228 barchart(_h_temp_Meff_6jt,_hMeff_6jt);
229
230 }
231
232 /// @}
233
234
235 private:
236
237 /// @name Histograms
238 /// @{
239 CounterPtr _h_2jl, _h_2jm, _h_2jt;
240 CounterPtr _h_4jt, _h_5j;
241 CounterPtr _h_6jm, _h_6jt;
242
243 Estimate1DPtr _hMeff_2jl, _hMeff_2jm, _hMeff_2jt;
244 Estimate1DPtr _hMeff_4jt, _hMeff_5j;
245 Estimate1DPtr _hMeff_6jm, _hMeff_6jt;
246
247 Histo1DPtr _h_temp_Meff_2jl, _h_temp_Meff_2jm, _h_temp_Meff_2jt;
248 Histo1DPtr _h_temp_Meff_4jt, _h_temp_Meff_5j;
249 Histo1DPtr _h_temp_Meff_6jm, _h_temp_Meff_6jt;
250 /// @}
251
252 /// Cut-flows
253 CutflowsPtr _flows;
254
255 };
256
257
258
259 RIVET_DECLARE_PLUGIN(ATLAS_2016_I1458270);
260
261
262}
|