Rivet analyses referenceCMS_2013_I1223519Searches for SUSY using $\alpha_T$ and $b$-quark multiplicity at 8 \text{TeV}Experiment: CMS (LHC) Inspire ID: 1223519 Status: UNVALIDATED Authors:
Beam energies: (4000.0, 4000.0) GeV Run details:
An inclusive search for supersymmetric processes that produce final states with jets and missing transverse energy in $pp$ collisions at 8 \text{TeV}. The data sample corresponds to an integrated luminosity of 11.7/fb collected by the CMS experiment. In this search, a dimensionless kinematic variable, $\alpha_T$, is used to discriminate between events with genuine and misreconstructed missing transverse energy. The search was based on an examination of the number of reconstructed jets per event, the scalar sum of transverse energies of these jets, and the number of these jets identified as originating from bottom quarks. No significant excess of events over the Standard Model expectation was found. Source code: CMS_2013_I1223519.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalStates.hh"
4#include "Rivet/Projections/FastJets.hh"
5#include "Rivet/Projections/MissingMomentum.hh"
6#include "Rivet/Projections/Smearing.hh"
7#include <bitset>
8
9namespace Rivet {
10
11
12 /// @brief Searches for SUSY using $\alpha_T$ and $b$-quark multiplicity at 8~\TeV
13 class CMS_2013_I1223519 : public Analysis {
14 public:
15
16 /// Constructor
17 RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2013_I1223519);
18
19
20 /// @name Analysis methods
21 /// @{
22
23 /// Book histograms and initialise projections before the run
24 void init() {
25
26 // Initialise and register projections
27 FinalState calofs(Cuts::abseta < 5.0);
28 declare(calofs, "Clusters");
29
30 MissingMomentum mm(calofs);
31 declare(mm, "TruthMET");
32 declare(SmearedMET(mm, MET_SMEAR_CMS_RUN2), "MET");
33
34 FastJets fj(calofs, JetAlg::ANTIKT, 0.5);
35 declare(fj, "TruthJets");
36 declare(SmearedJets(fj, JET_SMEAR_CMS_RUN2, [](const Jet& j) {
37 if (j.abseta() > 2.4) return 0.;
38 return j.bTagged() ? 0.65 : 0.01; }), "Jets"); ///< @note Charm mistag and exact b-tag eff not given
39
40 FinalState ys(Cuts::abspid == PID::PHOTON && Cuts::abseta < 5.0);
41 declare(ys, "TruthPhotons");
42 declare(SmearedParticles(ys, PHOTON_EFF_CMS_RUN2 /*, PHOTON_SMEAR_CMS_RUN2 */), "Photons");
43
44 FinalState es(Cuts::abspid == PID::ELECTRON && Cuts::abseta < 2.5);
45 declare(es, "TruthElectrons");
46 declare(SmearedParticles(es, ELECTRON_EFF_CMS_RUN2, ELECTRON_SMEAR_CMS_RUN2), "Electrons");
47
48 FinalState mus(Cuts::abspid == PID::MUON && Cuts::abseta < 2.4);
49 declare(mus, "TruthMuons");
50 declare(SmearedParticles(mus, MUON_EFF_CMS_RUN2, MUON_SMEAR_CMS_RUN2), "Muons");
51
52 ChargedFinalState cfs(Cuts::abseta < 2.5);
53 declare(cfs, "TruthTracks");
54 declare(SmearedParticles(cfs, TRK_EFF_CMS_RUN2), "Tracks");
55
56
57 // Book histograms
58 book(_h_alphaT23, "alphaT23", 15, 0, 3);
59 book(_h_alphaT4 , "alphaT4", 15, 0, 3);
60 /// @todo Add HT histograms
61
62 // Book counters
63 _h_srcounters.resize(8*7 + 3);
64 for (size_t inj = 0; inj < 2; ++inj) {
65 const size_t njmax = inj + 3;
66 for (size_t nb = 0; nb < njmax; ++nb) {
67 for (size_t iht = 0; iht < 8; ++iht) {
68 const size_t i = 8 * ((inj == 0 ? 0 : 3) + nb) + iht;
69 book(_h_srcounters[i], "srcount_j" + toString(njmax) + "_b" + toString(nb) + "_ht" + toString(iht+1));
70 }
71 }
72 }
73 // Special nj >= 4, nb >= 4 bins
74 for (size_t iht = 0; iht < 3; ++iht) {
75 book(_h_srcounters[8*7 + iht], "srcount_j4_b4_ht" + toString(iht+1));
76 }
77
78 }
79
80
81 /// Perform the per-event analysis
82 void analyze(const Event& event) {
83
84 // Get baseline photons, electrons & muons
85 Particles photons = apply<ParticleFinder>(event, "Photons").particles(Cuts::pT > 25*GeV);
86 Particles elecs = apply<ParticleFinder>(event, "Electrons").particles(Cuts::pT > 10*GeV);
87 Particles muons = apply<ParticleFinder>(event, "Muons").particles(Cuts::pT > 10*GeV);
88
89 // Electron/muon isolation (guesswork/copied from other CMS analysis -- paper is unspecific)
90 const Particles calofs = apply<ParticleFinder>(event, "Clusters").particles();
91 idiscard(photons, [&](const Particle& y) {
92 double ptsum = -y.pT();
93 for (const Particle& p : calofs)
94 if (deltaR(p,y) < 0.3) ptsum += p.pT();
95 return ptsum / y.pT() > 0.1;
96 });
97 idiscard(elecs, [&](const Particle& e) {
98 double ptsum = -e.pT();
99 for (const Particle& p : calofs)
100 if (deltaR(p,e) < 0.3) ptsum += p.pT();
101 return ptsum / e.pT() > 0.1;
102 });
103 idiscard(muons, [&](const Particle& m) {
104 double ptsum = -m.pT();
105 for (const Particle& p : calofs)
106 if (deltaR(p,m) < 0.3) ptsum += p.pT();
107 return ptsum / m.pT() > 0.2;
108 });
109
110 // Veto the event if there are any remaining baseline photons or leptons
111 if (!photons.empty()) vetoEvent;
112 if (!elecs.empty()) vetoEvent;
113 if (!muons.empty()) vetoEvent;
114
115
116 // Get jets and apply jet-based event-selection cuts
117 const JetFinder& jetproj = apply<JetFinder>(event, "Jets");
118 const Jets alljets = jetproj.jetsByPt(Cuts::abseta < 3.0 && Cuts::Et > 37*GeV); //< most inclusive jets requirement
119 if (select(alljets, Cuts::Et > 73*GeV).size() < 2) vetoEvent; //< most inclusive lead jets requirement
120
121 // Filter jets into different Et requirements & compute corresponding HTs
122 /// @note It's not clear if different HTs are used to choose the HT bins
123 const Jets jets37 = select(alljets, Cuts::Et > 37*GeV);
124 const Jets jets43 = select(jets37, Cuts::Et > 43*GeV);
125 const Jets jets50 = select(jets43, Cuts::Et > 50*GeV);
126 const double ht37 = sum(jets37, Kin::Et, 0.0);
127 const double ht43 = sum(jets43, Kin::Et, 0.0);
128 const double ht50 = sum(jets50, Kin::Et, 0.0);
129
130 // Find the relevant HT bin and apply leading jet event-selection cuts
131 static const vector<double> htcuts = { /* 275., 325., */ 375., 475., 575., 675., 775., 875.}; //< comment to avoid jets50 "fall-down"
132 const int iht = inRange(ht37, 275*GeV, 325*GeV) ? 0 : inRange(ht43, 325*GeV, 375*GeV) ? 1 : (2+binIndex(ht50, htcuts, true));
133 MSG_TRACE("HT = {" << ht37 << ", " << ht43 << ", " << ht50 << "} => IHT = " << iht);
134 if (iht < 0) vetoEvent;
135 if (iht == 1 && select(jets43, Cuts::Et > 78*GeV).size() < 2) vetoEvent;
136 if (iht >= 2 && select(jets50, Cuts::Et > 100*GeV).size() < 2) vetoEvent;
137
138 // Create references for uniform access to relevant set of jets & HT
139 const double etcut = iht == 0 ? 37. : iht == 1 ? 43. : 50.;
140 const double& ht = iht == 0 ? ht37 : iht == 1 ? ht43 : ht50;
141 const Jets& jets = iht == 0 ? jets37 : iht == 1 ? jets43 : jets50;
142 if (!jetproj.jets(Cuts::abseta > 3 && Cuts::Et > etcut*GeV).empty()) vetoEvent;
143 const size_t nj = jets.size();
144 const size_t nb = count_if(jets.begin(), jets.end(), [](const Jet& j) { return j.bTagged(Cuts::pT > 5*GeV); });
145
146 // Compute HTmiss = pT of 4-vector sum of jet momenta
147 const FourMomentum jsum = sum(jets, mom, FourMomentum());
148 const double htmiss = jsum.pT();
149
150 // Require HTmiss / ETmiss < 1.25
151 const double etmiss = apply<SmearedMET>(event, "MET").met();
152 if (htmiss/etmiss > 1.25) vetoEvent;
153
154 // Compute DeltaHT = minimum difference of "dijet" ETs, i.e. max(|1+2-3|, |1+3-2|, |2+3-1|)
155 double deltaht = -1;
156 vector<double> jetets; transform(jets, jetets, Kin::Et);
157 for (int i = 1; i < (1 << (jetets.size()-1)); ++i) { // count from 1 to 2**N-1, i.e. through all heterogeneous bitmasks with MSB(2**N)==0
158 const std::bitset<10> bits(i); /// @warning There'd better not be more than 10 jets...
159 const double htdiff = partition_diff(bits, jetets);
160 // MSG_INFO(bits.to_string() << " => " << htdiff);
161 if (deltaht < 0 || htdiff < deltaht) deltaht = htdiff;
162 }
163 MSG_DEBUG("dHT_bitmask = " << deltaht);
164
165 // Cross-check calculation in 2- and 3-jet cases
166 // if (jets.size() == 2) {
167 // MSG_INFO("dHT2 = " << fabs(jets[0].Et() - jets[1].Et()));
168 // } else if (jets.size() == 3) {
169 // double deltaht_01_2 = fabs(jets[0].Et()+jets[1].Et()-jets[2].Et());
170 // double deltaht_02_1 = fabs(jets[0].Et()+jets[2].Et()-jets[1].Et());
171 // double deltaht_12_0 = fabs(jets[1].Et()+jets[2].Et()-jets[0].Et());
172 // MSG_INFO("dHT3 = " << min({deltaht_01_2, deltaht_02_1, deltaht_12_0}));
173 // }
174
175 // Compute alphaT from the above
176 double alphaT = fabs(0.5*((ht-deltaht)/(sqrt((ht*ht)-(htmiss*htmiss)))));
177 if (alphaT < 0.55) vetoEvent;
178
179 /// @todo Need to include trigger efficiency sampling or weighting?
180
181 // Fill histograms
182 const size_t inj = nj < 4 ? 0 : 1;
183 const size_t inb = nb < 4 ? nb : 4;
184 if (iht >= 2)
185 (inj == 0 ? _h_alphaT23 : _h_alphaT4)->fill(alphaT);
186
187 // Fill the appropriate counter -- after working out the irregular SR bin index! *sigh*
188 size_t i = 8 * ((inj == 0 ? 0 : 3) + inb) + iht;
189 if (inj == 1 && inb == 4) i = 8*7 + (iht < 3 ? iht : 2);
190 MSG_INFO("inj = " << inj << ", inb = " << inb << ", i = " << i);
191 _h_srcounters[i]->fill();
192
193 }
194
195
196 /// Normalise histograms etc., after the run
197 void finalize() {
198
199 const double sf = crossSection()/femtobarn*11.7/sumOfWeights();
200 scale(_h_alphaT23, sf);
201 scale(_h_alphaT4, sf);
202 for (size_t i = 0; i < 8*7+3; ++i)
203 scale(_h_srcounters[i], sf);
204
205 }
206
207 /// @}
208
209
210 /// @name Utility functions for partitioning jet pTs into two groups and summing/diffing them
211 /// @{
212
213 /// Sum the given values into two subsets according to the provided bitmask
214 template <size_t N>
215 pair<double, double> partition_sum(const std::bitset<N>& mask,
216 const vector<double>& vals) const {
217 pair<double, double> rtn(0., 0.);
218 for (size_t i = 0; i < vals.size(); ++i) {
219 (!mask[vals.size()-1-i] ? rtn.first : rtn.second) += vals[i];
220 }
221 return rtn;
222 }
223
224 /// Return the difference between summed subsets according to the provided bitmask
225 template <size_t N>
226 double partition_diff(const std::bitset<N>& mask, const vector<double>& vals) const {
227 const pair<double, double> sums = partition_sum(mask, vals);
228 const double diff = fabs(sums.first - sums.second);
229 MSG_TRACE(mask.to_string() << ": " << sums.first << "/" << sums.second << " => " << diff);
230 return diff;
231 }
232
233 /// @}
234
235
236 /// @name Histograms
237 /// @{
238 Histo1DPtr _h_alphaT23, _h_alphaT4;
239 vector<CounterPtr> _h_srcounters;
240 /// @}
241
242
243 };
244
245
246 RIVET_DECLARE_PLUGIN(CMS_2013_I1223519);
247
248
249}
|