Rivet analyses referenceCMS_2016_PAS_SUS_16_14Search for supersymmetry in events with jets and missing transverse momentum at 13 \text{TeV}Experiment: CMS (LHC) Status: VALIDATED Authors:
Beams: p+ p+ Beam energies: (6500.0, 6500.0) GeV Run details:
A search for supersymmetry in all-hadronic events with large missing transverse momentum, produced in proton--proton collisions at $\sqrt{s} = 13 \text{TeV}$. The data sample, corresponding to an integrated luminosity of 12.9/fb, was collected with the CMS detector at the CERN LHC in 2016. The data are examined in search regions of jet multiplicity, tagged bottom quark jet multiplicity, missing transverse momentum, and the scalar sum of jet transverse momenta. The observed numbers of events in all search regions are found to be consistent with the expectations from Standard Model processes. Source code: CMS_2016_PAS_SUS_16_14.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
11namespace Rivet {
12
13
14 /// CMS 2016 0-lepton SUSY search, from 13/fb PAS note
15 class CMS_2016_PAS_SUS_16_14 : public Analysis {
16 public:
17
18 /// Constructor
19 RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2016_PAS_SUS_16_14);
20
21
22 /// @name Analysis methods
23 /// @{
24
25 /// Book histograms and initialise projections before the run
26 void init() {
27
28 // Initialise and register projections
29 FinalState calofs(Cuts::abseta < 5.0);
30 FastJets fj(calofs, JetAlg::ANTIKT, 0.4);
31 declare(fj, "TruthJets");
32 declare(SmearedJets(fj, JET_SMEAR_CMS_RUN2, [](const Jet& j) {
33 if (j.abseta() > 2.5) return 0.;
34 return j.bTagged() ? 0.55 : j.cTagged() ? 0.12 : 0.016; }), "Jets");
35
36 FinalState es(Cuts::abspid == PID::ELECTRON && Cuts::abseta < 2.5);
37 declare(es, "TruthElectrons");
38 declare(SmearedParticles(es, ELECTRON_EFF_CMS_RUN2, ELECTRON_SMEAR_CMS_RUN2), "Electrons");
39
40 FinalState mus(Cuts::abspid == PID::MUON && Cuts::abseta < 2.4);
41 declare(mus, "TruthMuons");
42 declare(SmearedParticles(mus, MUON_EFF_CMS_RUN2, MUON_SMEAR_CMS_RUN2), "Muons");
43
44 FinalState isofs(Cuts::abseta < 3.0 && Cuts::abspid != PID::ELECTRON && Cuts::abspid != PID::MUON);
45 declare(isofs, "IsoFS");
46 FinalState cfs(Cuts::abseta < 2.5 && Cuts::abscharge != 0);
47 declare(cfs, "TruthTracks");
48 declare(SmearedParticles(cfs, TRK_EFF_CMS_RUN2), "Tracks");
49
50 // Book histograms/counters
51 _h_srcounts.resize(160);
52 for (size_t ij = 0; ij < 4; ++ij) {
53 for (size_t ib = 0; ib < 4; ++ib) {
54 for (size_t ih = 0; ih < 10; ++ih) {
55 const size_t i = 40*ij + 10*ib + ih;
56 book(_h_srcounts[i], toString(2*ij+3) + "j-" + toString(ib) + "b-" + toString(ih));
57 }
58 }
59 }
60 _h_srcountsagg.resize(12);
61 for (size_t ia = 0; ia < 12; ++ia) {
62 book(_h_srcountsagg[ia], "agg-" + toString(ia));
63 }
64
65 }
66
67
68 /// Perform the per-event analysis
69 void analyze(const Event& event) {
70
71 // Get jets and require Nj >= 3
72 const Jets jets24 = apply<JetFinder>(event, "Jets").jetsByPt(Cuts::pT > 30*GeV && Cuts::abseta < 2.4);
73 if (jets24.size() < 3) vetoEvent;
74
75 // HT cut
76 vector<double> jetpts24; transform(jets24, jetpts24, Kin::pT);
77 const double ht = sum(jetpts24, 0.0);
78 if (ht < 300*GeV) vetoEvent;
79
80 // HTmiss cut
81 const Jets jets50 = apply<JetFinder>(event, "Jets").jetsByPt(Cuts::pT > 30*GeV && Cuts::abseta < 5.0);
82 const FourMomentum htmissvec = -sum(jets24, mom, FourMomentum());
83 const double htmiss = htmissvec.pT();
84 if (htmissvec.pT() < 300*GeV) vetoEvent;
85
86
87 // Get baseline electrons & muons
88 Particles elecs = apply<ParticleFinder>(event, "Electrons").particles(Cuts::pT > 10*GeV);
89 Particles muons = apply<ParticleFinder>(event, "Muons").particles(Cuts::pT > 10*GeV);
90
91 // Electron/muon isolation
92 const Particles calofs = apply<ParticleFinder>(event, "IsoFS").particles();
93 idiscard(elecs, [&](const Particle& e) {
94 const double R = max(0.05, min(0.2, 10*GeV/e.pT()));
95 double ptsum = -e.pT();
96 for (const Particle& p : calofs)
97 if (deltaR(p,e) < R) ptsum += p.pT();
98 return ptsum / e.pT() > 0.1;
99 });
100 idiscard(muons, [&](const Particle& m) {
101 const double R = max(0.05, min(0.2, 10*GeV/m.pT()));
102 double ptsum = -m.pT();
103 for (const Particle& p : calofs)
104 if (deltaR(p,m) < R) ptsum += p.pT();
105 return ptsum / m.pT() > 0.2;
106 });
107
108 // Veto the event if there are any remaining baseline leptons
109 if (!elecs.empty()) vetoEvent;
110 if (!muons.empty()) vetoEvent;
111
112
113 // Get isolated tracks
114 Particles trks25 = apply<ParticleFinder>(event, "Tracks").particles();
115 idiscard(trks25, [&](const Particle& t) {
116 double ptsum = -t.pT();
117 for (const Particle& p : trks25)
118 if (deltaR(p,t) < 0.3) ptsum += p.pT();
119 return ptsum/t.pT() > ((t.abspid() == PID::ELECTRON || t.abspid() == PID::MUON) ? 0.2 : 0.1);
120 });
121 const Particles trks = select(trks25, Cuts::abseta < 2.4);
122
123 // Isolated track pT, pTmiss and mT cut
124 // mT^2 = m1^2 + m2^2 + 2(ET1 ET2 - pT1 . pT2))
125 // => mT0^2 = 2(ET1 |pT2| - pT1 . pT2)) for m1, m2 -> 0
126 FourMomentum ptmissvec = htmissvec; ///< @todo Can we do better? No e,mu left...
127 const double ptmiss = ptmissvec.pT();
128 for (const Particle& t : trks) {
129 const double ptcut = (t.abspid() == PID::ELECTRON || t.abspid() == PID::MUON) ? 5*GeV : 10*GeV;
130 const double mT = sqrt( t.mass2() + 2*(t.Et()*ptmiss - t.pT()*ptmiss*cos(deltaPhi(t,ptmissvec))) );
131 if (mT < 100*GeV && t.pT() < ptcut) vetoEvent;
132 }
133
134 // Lead jets isolation from Htmiss
135 if (deltaPhi(htmissvec, jets24[0]) < 0.5) vetoEvent;
136 if (deltaPhi(htmissvec, jets24[1]) < 0.5) vetoEvent;
137 if (deltaPhi(htmissvec, jets24[2]) < 0.3) vetoEvent;
138 if (jets24.size() >= 4 && deltaPhi(htmissvec, jets24[3]) < 0.3) vetoEvent;
139
140 // Count jet and b-jet multiplicities
141 const size_t nj = jets24.size();
142 size_t nbj = 0;
143 for (const Jet& j : jets24)
144 if (j.bTagged()) nbj += 1;
145
146
147 ////////
148
149
150 // Fill the aggregate signal regions first
151 if (nj >= 3 && nbj == 0 && ht > 500*GeV && htmiss > 500*GeV) _h_srcountsagg[ 0]->fill(1.0);
152 if (nj >= 3 && nbj == 0 && ht > 1500*GeV && htmiss > 750*GeV) _h_srcountsagg[ 1]->fill(1.0);
153 if (nj >= 5 && nbj == 0 && ht > 500*GeV && htmiss > 500*GeV) _h_srcountsagg[ 2]->fill(1.0);
154 if (nj >= 5 && nbj == 0 && ht > 1500*GeV && htmiss > 750*GeV) _h_srcountsagg[ 3]->fill(1.0);
155 if (nj >= 9 && nbj == 0 && ht > 1500*GeV && htmiss > 750*GeV) _h_srcountsagg[ 4]->fill(1.0);
156 if (nj >= 3 && nbj >= 2 && ht > 500*GeV && htmiss > 500*GeV) _h_srcountsagg[ 5]->fill(1.0);
157 if (nj >= 3 && nbj >= 1 && ht > 750*GeV && htmiss > 750*GeV) _h_srcountsagg[ 6]->fill(1.0);
158 if (nj >= 5 && nbj >= 3 && ht > 500*GeV && htmiss > 500*GeV) _h_srcountsagg[ 7]->fill(1.0);
159 if (nj >= 5 && nbj >= 2 && ht > 1500*GeV && htmiss > 750*GeV) _h_srcountsagg[ 8]->fill(1.0);
160 if (nj >= 9 && nbj >= 3 && ht > 750*GeV && htmiss > 750*GeV) _h_srcountsagg[ 9]->fill(1.0);
161 if (nj >= 7 && nbj >= 1 && ht > 300*GeV && htmiss > 300*GeV) _h_srcountsagg[10]->fill(1.0);
162 if (nj >= 5 && nbj >= 1 && ht > 750*GeV && htmiss > 750*GeV) _h_srcountsagg[11]->fill(1.0);
163
164
165 // Nj bin and Nbj bins
166 static const vector<double> njedges = {3., 5., 7., 9.};
167 const size_t inj = binIndex(nj, njedges, true);
168 static const vector<double> njbedges = {0., 1., 2., 3.};
169 const size_t inbj = binIndex(nbj, njbedges, true);
170 // HTmiss vs HT 2D bin
171 int iht = 0;
172 if (htmiss < 350*GeV) {
173 iht = ht < 500 ? 1 : ht < 1000 ? 2 : 3;
174 } else if (htmiss < 500*GeV && ht > 350*GeV) {
175 iht = ht < 500 ? 4 : ht < 1000 ? 5 : 6;
176 } else if (htmiss < 750*GeV && ht > 500*GeV) {
177 iht = ht < 1000 ? 7 : 8;
178 } else if (ht > 750*GeV) {
179 iht = ht < 1500 ? 9 : 10;
180 }
181 if (iht == 0) vetoEvent;
182 iht -= 1; //< change from the paper's indexing scheme to C++ zero-indexed
183 // Total bin number
184 const size_t ibin = 40*inj + 10*inbj + (size_t)iht;
185
186 // Fill SR counter
187 _h_srcounts[ibin]->fill(1.0);
188
189 }
190
191
192 /// Normalise counters after the run
193 void finalize() {
194
195 const double sf = 12.9*crossSection()/femtobarn/sumOfWeights();
196 scale(_h_srcounts, sf);
197 scale(_h_srcountsagg, sf);
198
199 }
200
201 /// @}
202
203
204 private:
205
206 /// @name Histograms
207 /// @{
208 vector<CounterPtr> _h_srcounts, _h_srcountsagg;
209 /// @}
210
211 };
212
213
214
215 RIVET_DECLARE_PLUGIN(CMS_2016_PAS_SUS_16_14);
216
217}
|