rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

CMS_2013_I1223519

Searches for SUSY using $\alpha_T$ and $b$-quark multiplicity at 8 \text{TeV}
Experiment: CMS (LHC)
Inspire ID: 1223519
Status: UNVALIDATED
Authors:
  • Andy Buckley
References: Beams: p+ p+
Beam energies: (4000.0, 4000.0) GeV
Run details:
  • SM background or BSM physics model, depending on interpretation usage

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, FastJets::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      ifilter_discard(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      ifilter_discard(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      ifilter_discard(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 JetAlg& jetproj = apply<JetAlg>(event, "Jets");
118      const Jets alljets = jetproj.jetsByPt(Cuts::abseta < 3.0 && Cuts::Et > 37*GeV); //< most inclusive jets requirement
119      if (filter_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 = filter_select(alljets, Cuts::Et > 37*GeV);
124      const Jets jets43 = filter_select(jets37, Cuts::Et > 43*GeV);
125      const Jets jets50 = filter_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 && filter_select(jets43, Cuts::Et > 78*GeV).size() < 2) vetoEvent;
136      if (iht >= 2 && filter_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 double weight = 1.0;
183      const size_t inj = nj < 4 ? 0 : 1;
184      const size_t inb = nb < 4 ? nb : 4;
185      if (iht >= 2)
186        (inj == 0 ? _h_alphaT23 : _h_alphaT4)->fill(alphaT, weight);
187
188      // Fill the appropriate counter -- after working out the irregular SR bin index! *sigh*
189      size_t i = 8 * ((inj == 0 ? 0 : 3) + inb) + iht;
190      if (inj == 1 && inb == 4) i = 8*7 + (iht < 3 ? iht : 2);
191      MSG_INFO("inj = " << inj << ", inb = " << inb << ", i = " << i);
192      _h_srcounters[i]->fill(weight);
193
194    }
195
196
197    /// Normalise histograms etc., after the run
198    void finalize() {
199
200      const double sf = crossSection()/femtobarn*11.7/sumOfWeights();
201      scale(_h_alphaT23, sf);
202      scale(_h_alphaT4, sf);
203      for (size_t i = 0; i < 8*7+3; ++i)
204        scale(_h_srcounters[i], sf);
205
206    }
207
208    //@}
209
210
211    /// @name Utility functions for partitioning jet pTs into two groups and summing/diffing them
212    //@{
213
214    /// Sum the given values into two subsets according to the provided bitmask
215    template <size_t N>
216    pair<double, double> partition_sum(const std::bitset<N>& mask,
217                                       const vector<double>& vals) const {
218      pair<double, double> rtn(0., 0.);
219      for (size_t i = 0; i < vals.size(); ++i) {
220        (!mask[vals.size()-1-i] ? rtn.first : rtn.second) += vals[i];
221      }
222      return rtn;
223    }
224
225    /// Return the difference between summed subsets according to the provided bitmask
226    template <size_t N>
227    double partition_diff(const std::bitset<N>& mask, const vector<double>& vals) const {
228      const pair<double, double> sums = partition_sum(mask, vals);
229      const double diff = fabs(sums.first - sums.second);
230      MSG_TRACE(mask.to_string() << ": " << sums.first << "/" << sums.second << " => " << diff);
231      return diff;
232    }
233
234    //@}
235
236
237    /// @name Histograms
238    //@{
239    Histo1DPtr _h_alphaT23, _h_alphaT4;
240    vector<CounterPtr> _h_srcounters;
241    //@}
242
243
244  };
245
246
247  // The hook for the plugin system
248  RIVET_DECLARE_PLUGIN(CMS_2013_I1223519);
249
250
251}