Rivet analyses referenceCMS_2017_I1594909Search for SUSY in multijet events with missing transverse momentum in $pp$ collisions at 13 TeVExperiment: CMS (LHC) Inspire ID: 1594909 Status: VALIDATED Authors:
Beam energies: (6500.0, 6500.0) GeV Run details:
A search for supersymmetry based on multijet events with large missing transverse momentum produced in proton-proton collisions at a center-of-mass energy of $\sqrt{s} = 13\,\text{TeV}$. The data, corresponding to an integrated luminosity of 35.9/fb, were collected with the CMS detector at the CERN LHC in 2016. The analysis utilizes four-dimensional exclusive search regions defined in terms of the number of jets, the number of tagged bottom-quark jets, the scalar sum of jet transverse momenta, and the magnitude of the vector sum of jet transverse momenta. This coding presents the fully detailed 174 signal regions, as well as 12 aggregate signal regions. Source code: CMS_2017_I1594909.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/FastJets.hh"
5#include "Rivet/Projections/ChargedFinalState.hh"
6#include "Rivet/Projections/VisibleFinalState.hh"
7#include "Rivet/Projections/MissingMomentum.hh"
8#include "Rivet/Projections/SmearedParticles.hh"
9#include "Rivet/Projections/SmearedJets.hh"
10#include "Rivet/Projections/SmearedMET.hh"
11#include "Rivet/Tools/Cutflow.hh"
12#include <tuple>
13
14namespace Rivet {
15
16
17 /// CMS search for SUSY with multijet + MET signatures in 36/fb of 13 TeV pp data
18 class CMS_2017_I1594909 : public Analysis {
19 public:
20
21 /// Constructor
22 RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2017_I1594909);
23
24
25 /// @name Analysis methods
26 //@{
27
28 /// Book histograms and initialise projections before the run
29 void init() {
30
31 VisibleFinalState pfall;
32 declare(pfall, "PFAll");
33 ChargedFinalState pfchg(Cuts::abseta < 2.5);
34 declare(pfchg, "PFChg");
35
36 FastJets jets(FinalState(Cuts::abseta < 4.9), FastJets::ANTIKT, 0.4);
37 SmearedJets recojets(jets, JET_SMEAR_CMS_RUN2, [](const Jet& j){ return j.bTagged() ? 0.55 : j.cTagged() ? 0.12 : 0.016; });
38 declare(recojets, "Jets");
39
40 FinalState electrons(Cuts::abspid == PID::ELECTRON && Cuts::abseta < 2.5);
41 SmearedParticles recoelectrons(electrons, ELECTRON_EFF_CMS_RUN2);
42 declare(recoelectrons, "Electrons");
43
44 FinalState muons(Cuts::abspid == PID::MUON && Cuts::abseta < 2.4);
45 SmearedParticles recomuons(muons, MUON_EFF_CMS_RUN2);
46 declare(recomuons, "Muons");
47
48 VisibleFinalState calofs(Cuts::abseta < 4.9 && Cuts::abspid != PID::MUON);
49 MissingMomentum met(calofs);
50 SmearedMET recomet(met, MET_SMEAR_CMS_RUN2);
51 declare(recomet, "MET");
52
53
54 // Book counters, into a map of 3 indices since the global index is not obvious to calculate
55 size_t i = 0;
56 for (int j = 1; j <= 5; ++j) {
57 for (int b = 1; b <= 4; ++b) {
58 if (j == 1 && b == 4) continue;
59 for (int k = 1; k <= 10; ++k) {
60 if (j > 3 && (k == 1 || k == 4)) continue;
61 stringstream s; s << "count_" << (i+1); // << "_" << j << b << k;
62 book(_counts[std::make_tuple(j,b,k)], s.str());
63 i += 1;
64 }
65 }
66 }
67 MSG_DEBUG("Booked " << i << " signal regions (should be 174)");
68 // Aggregate SR counters
69 for (size_t i = 0; i < 12; ++i)
70 book(_counts_agg[i], "count_agg_" + toString(i+1));
71
72
73 // Book cut-flow
74 _flow = Cutflow("Presel", {"Njet>=2", "HT>300", "HTmiss>300",
75 "Nmuon=0", "Nmuisotrk=0", "Nelec=0", "Nelisotrk=0", "Nhadisotrk=0",
76 "dPhi_miss,j1>0.5", "dPhi_miss,j2>0.5", "dPhi_miss,j3>0.3", "dPhi_miss,j4>0.3"
77 });
78
79 }
80
81
82 /// Perform the per-event analysis
83 void analyze(const Event& event) {
84
85 _flow.fillinit();
86
87 // Find leptons and isolation particles
88 const Particles elecs = apply<ParticleFinder>(event, "Electrons").particlesByPt();
89 const Particles mus = apply<ParticleFinder>(event, "Muons").particlesByPt();
90 const Particles pfall = apply<ParticleFinder>(event, "PFAll").particlesByPt();
91 const Particles pfiso = filter_select(pfall, [](const Particle& p){ return p.isHadron() || p.pid() == PID::PHOTON; });
92
93 // Find isolated leptons
94 const Particles isoleps = filter_select(elecs+mus, [&](const Particle& l){
95 const double dR = l.pT() < 50*GeV ? 0.2 : l.pT() < 200*GeV ? 10*GeV/l.pT() : 0.05;
96 const double sumpt = sum(filter_select(pfiso, deltaRLess(l, dR)), Kin::pT, 0.0);
97 return sumpt/l.pT() < (l.abspid() == PID::ELECTRON ? 0.1 : 0.2); //< different I criteria for e and mu
98 });
99
100 // Find other isolated tracks
101 const Particles pfchg = apply<ParticleFinder>(event, "PFChg").particlesByPt();
102 const Particles isochgs = filter_select(pfchg, [&](const Particle& t){
103 if (t.abseta() > 2.4) return false;
104 if (any(isoleps, deltaRLess(t, 0.01))) return false; //< don't count isolated leptons here
105 const double sumpt = sum(filter_select(pfchg, deltaRLess(t, 0.3)), Kin::pT, -t.pT());
106 return sumpt/t.pT() < ((t.abspid() == PID::ELECTRON || t.abspid() == PID::MUON) ? 0.2 : 0.1);
107 });
108
109 // Find and isolate jets
110 const Jets jets = apply<JetAlg>(event, "Jets").jetsByPt(Cuts::pT > 30*GeV);
111 const Jets cjets = filter_select(jets, Cuts::abseta < 2.4);
112 const Jets isojets = cjets; //discardIfAnyDeltaRLess(cjets, elecs+mus, 0.4);
113 const int njets = isojets.size();
114 const Jets isobjets = filter_select(isojets, hasBTag());
115 const int nbjets = isobjets.size();
116 MSG_DEBUG("Njets = " << jets.size() << ", Nisojets = " << njets << ", Nbjets = " << nbjets);
117
118 // Calculate HT, HTmiss, and pTmiss quantities
119 const double ht = sum(jets, Kin::pT, 0.0);
120 const Vector3 vhtmiss = -sum(jets, pTvec, Vector3());
121 const double htmiss = vhtmiss.perp();
122 const Vector3& vptmiss = -apply<SmearedMET>(event, "MET").vectorEt();
123 const double ptmiss = vptmiss.perp();
124 MSG_DEBUG("HT = " << ht/GeV << " GeV, HTmiss = " << htmiss/GeV << " GeV");
125
126
127 /////////////////////////////////////
128 // Event selection
129
130 // Njet cut
131 if (njets < 2) vetoEvent;
132 _flow.fill(1);
133 // HT cut
134 if (ht < 300*GeV) vetoEvent;
135 _flow.fill(2);
136 // HTmiss cut
137 if (htmiss < 300*GeV) vetoEvent;
138 _flow.fill(3);
139
140 // Isolated leptons cut
141 if (!filter_select(isoleps, Cuts::pT > 10*GeV).empty()) vetoEvent;
142 // Isolated tracks cut
143 for (const Particle& t : isochgs) {
144 const double mT = sqrt(2*t.pT()*ptmiss * (1 - cos(deltaPhi(t, vptmiss))) );
145 if (mT < 100*GeV) continue;
146 const double pTmax = (t.abspid() == PID::ELECTRON || t.abspid() == PID::MUON) ? 5*GeV : 10*GeV;
147 if (t.pT() > pTmax) vetoEvent;
148 }
149 //
150 // // Inefficiently separated version of isolation cuts for detailed cutflow debugging
151 // // Muon cut
152 // if (!filter_select(isoleps, Cuts::pT > 10*GeV && Cuts::abspid == PID::MUON).empty()) vetoEvent;
153 // _flow.fill(4);
154 // // Muon isotrk cut
155 // for (const Particle& t : filter_select(isochgs, Cuts::abspid == PID::MUON)) {
156 // const double mT = sqrt(2*t.pT()*ptmiss * (1 - cos(deltaPhi(t, vptmiss))) );
157 // if (mT > 100*GeV && t.pT() > 5*GeV) vetoEvent;
158 // }
159 // _flow.fill(5);
160 // // Electron cut
161 // if (!filter_select(isoleps, Cuts::pT > 10*GeV && Cuts::abspid == PID::ELECTRON).empty()) vetoEvent;
162 // _flow.fill(6);
163 // // Electron isotrk cut
164 // for (const Particle& t : filter_select(isochgs, Cuts::abspid == PID::ELECTRON)) {
165 // const double mT = sqrt(2*t.pT()*ptmiss * (1 - cos(deltaPhi(t, vptmiss))) );
166 // if (mT > 100*GeV && t.pT() > 5*GeV) vetoEvent;
167 // }
168 // _flow.fill(7);
169 // // Hadron isotrk cut
170 // for (const Particle& t : filter_select(isochgs, Cuts::abspid != PID::ELECTRON && Cuts::abspid != PID::MUON)) {
171 // const double mT = sqrt(2*t.pT()*ptmiss * (1 - cos(deltaPhi(t, vptmiss))) );
172 // if (mT > 100*GeV && t.pT() > 10*GeV) vetoEvent;
173 // }
174 _flow.fill(8);
175
176
177 // dPhi(jet,HTmiss) cuts
178 if (deltaPhi(vhtmiss, isojets[0]) < 0.5) vetoEvent;
179 _flow.fill(9);
180 if (deltaPhi(vhtmiss, isojets[1]) < 0.5) vetoEvent;
181 _flow.fill(10);
182 if (njets >= 3 && deltaPhi(vhtmiss, isojets[2]) < 0.3) vetoEvent;
183 _flow.fill(11);
184 if (njets >= 4 && deltaPhi(vhtmiss, isojets[3]) < 0.3) vetoEvent;
185 _flow.fill(12);
186
187
188 /////////////////////////////////////
189 // Find SR index and fill counter
190
191 const double w = 1.0;
192
193 const int idx_j = binIndex(njets, vector<int>{2,3,5,7,9}, true);
194 const int idx_b = binIndex(nbjets, vector<int>{0,1,2,3}, true);
195 int idx_k = -1;
196 if (inRange(htmiss/GeV, 300, 350)) {
197 idx_k = ht < 500*GeV ? 1 : ht < 1000*GeV ? 2 : 3;
198 } else if (inRange(htmiss/GeV, 350, 500) && ht > 350*GeV) {
199 idx_k = ht < 500*GeV ? 4 : ht < 1000*GeV ? 5 : 6;
200 } else if (inRange(htmiss/GeV, 500, 750) && ht > 500*GeV) {
201 idx_k = ht < 1000*GeV ? 7 : 8;
202 } else if (htmiss/GeV > 750 && ht > 750*GeV) {
203 idx_k = ht < 1500*GeV ? 9 : 10;
204 }
205
206 // Fill via 3-tuple index
207 if (idx_j >= 0 && idx_b >= 0 && idx_k >= 0) {
208 const auto idx = std::make_tuple(idx_j+1,idx_b+1,idx_k);
209 if (has_key(_counts, idx)) _counts[idx]->fill(w);
210 }
211
212
213 /////////////////////////////////////
214 // Aggregate SRs
215
216 // Region Njet Nb-jet HT [GeV] HTmiss [GeV] Parton multiplicity Heavy flavor ? ∆m
217 if (njets >= 2 && nbjets == 0 && ht >= 500*GeV && htmiss >= 500*GeV) _counts_agg[0]->fill(w);
218 if (njets >= 3 && nbjets == 0 && ht >= 1500*GeV && htmiss >= 750*GeV) _counts_agg[1]->fill(w);
219 if (njets >= 5 && nbjets == 0 && ht >= 500*GeV && htmiss >= 500*GeV) _counts_agg[2]->fill(w);
220 if (njets >= 5 && nbjets == 0 && ht >= 1500*GeV && htmiss >= 750*GeV) _counts_agg[3]->fill(w);
221 if (njets >= 9 && nbjets == 0 && ht >= 1500*GeV && htmiss >= 750*GeV) _counts_agg[4]->fill(w);
222 if (njets >= 2 && nbjets >= 2 && ht >= 500*GeV && htmiss >= 500*GeV) _counts_agg[5]->fill(w);
223 if (njets >= 3 && nbjets >= 1 && ht >= 750*GeV && htmiss >= 750*GeV) _counts_agg[6]->fill(w);
224 if (njets >= 5 && nbjets >= 3 && ht >= 500*GeV && htmiss >= 500*GeV) _counts_agg[7]->fill(w);
225 if (njets >= 5 && nbjets >= 2 && ht >= 1500*GeV && htmiss >= 750*GeV) _counts_agg[8]->fill(w);
226 if (njets >= 9 && nbjets >= 3 && ht >= 750*GeV && htmiss >= 750*GeV) _counts_agg[9]->fill(w);
227 if (njets >= 7 && nbjets >= 1 && ht >= 300*GeV && htmiss >= 300*GeV) _counts_agg[10]->fill(w);
228 if (njets >= 5 && nbjets >= 1 && ht >= 750*GeV && htmiss >= 750*GeV) _counts_agg[11]->fill(w);
229
230 }
231
232
233 /// Normalise histograms etc., after the run
234 void finalize() {
235
236 const double norm = 35.9*crossSection()/femtobarn;
237 const double sf = norm/sumOfWeights();
238
239 for (auto& idx_cptr : _counts)
240 scale(idx_cptr.second, sf);
241 for (CounterPtr& cptr : _counts_agg)
242 scale(cptr, sf);
243
244 _flow.scale(sf);
245 MSG_INFO("CUTFLOWS:\n\n" << _flow);
246
247 }
248
249 //@}
250
251
252 private:
253
254 Cutflow _flow;
255
256 map<tuple<int,int,int>, CounterPtr> _counts;
257 CounterPtr _counts_agg[12];
258
259 };
260
261
262 // The hook for the plugin system
263 RIVET_DECLARE_PLUGIN(CMS_2017_I1594909);
264
265
266}
|