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
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalStates.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/MissingMomentum.hh"
#include "Rivet/Projections/Smearing.hh"
#include <bitset>

namespace Rivet {


  /// @brief Add a short analysis description here
  class CMS_2013_I1223519 : public Analysis {
  public:

    /// Constructor
    DEFAULT_RIVET_ANALYSIS_CTOR(CMS_2013_I1223519);


    /// @name Analysis methods
    //@{

    /// Book histograms and initialise projections before the run
    void init() {

      // Initialise and register projections
      FinalState calofs(Cuts::abseta < 5.0);
      declare(calofs, "Clusters");

      MissingMomentum mm(calofs);
      declare(mm, "TruthMET");
      declare(SmearedMET(mm, MET_SMEAR_CMS_RUN2), "MET");

      FastJets fj(calofs, FastJets::ANTIKT, 0.5);
      declare(fj, "TruthJets");
      declare(SmearedJets(fj, JET_SMEAR_CMS_RUN2, [](const Jet& j) {
            if (j.abseta() > 2.4) return 0.;
            return j.bTagged() ? 0.65 : 0.01; }), "Jets"); ///< @note Charm mistag and exact b-tag eff not given

      FinalState ys(Cuts::abspid == PID::PHOTON && Cuts::abseta < 5.0);
      declare(ys, "TruthPhotons");
      declare(SmearedParticles(ys, PHOTON_EFF_CMS_RUN2 /*, PHOTON_SMEAR_CMS_RUN2 */), "Photons");

      FinalState es(Cuts::abspid == PID::ELECTRON && Cuts::abseta < 2.5);
      declare(es, "TruthElectrons");
      declare(SmearedParticles(es, ELECTRON_EFF_CMS_RUN2, ELECTRON_SMEAR_CMS_RUN2), "Electrons");

      FinalState mus(Cuts::abspid == PID::MUON && Cuts::abseta < 2.4);
      declare(mus, "TruthMuons");
      declare(SmearedParticles(mus, MUON_EFF_CMS_RUN2, MUON_SMEAR_CMS_RUN2), "Muons");

      ChargedFinalState cfs(Cuts::abseta < 2.5);
      declare(cfs, "TruthTracks");
      declare(SmearedParticles(cfs, TRK_EFF_CMS_RUN2), "Tracks");


      // Book histograms
      _h_alphaT23 = bookHisto1D("alphaT23", 15, 0, 3);
      _h_alphaT4 = bookHisto1D("alphaT4", 15, 0, 3);
      /// @todo Add HT histograms

      // Book counters
      _h_srcounters.resize(8*7 + 3);
      for (size_t inj = 0; inj < 2; ++inj) {
        const size_t njmax = inj + 3;
        for (size_t nb = 0; nb < njmax; ++nb) {
          for (size_t iht = 0; iht < 8; ++iht) {
            const size_t i = 8 * ((inj == 0 ? 0 : 3) + nb) + iht;
            _h_srcounters[i] = bookCounter("srcount_j" + toString(njmax) + "_b" + toString(nb) + "_ht" + toString(iht+1));
          }
        }
      }
      // Special nj >= 4, nb >= 4 bins
      for (size_t iht = 0; iht < 3; ++iht) {
        _h_srcounters[8*7 + iht] = bookCounter("srcount_j4_b4_ht" + toString(iht+1));
      }

    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {

      // Get baseline photons, electrons & muons
      Particles photons = apply<ParticleFinder>(event, "Photons").particles(Cuts::pT > 25*GeV);
      Particles elecs = apply<ParticleFinder>(event, "Electrons").particles(Cuts::pT > 10*GeV);
      Particles muons = apply<ParticleFinder>(event, "Muons").particles(Cuts::pT > 10*GeV);

      // Electron/muon isolation (guesswork/copied from other CMS analysis -- paper is unspecific)
      const Particles calofs = apply<ParticleFinder>(event, "Clusters").particles();
      ifilter_discard(photons, [&](const Particle& y) {
          double ptsum = -y.pT();
          for (const Particle& p : calofs)
            if (deltaR(p,y) < 0.3) ptsum += p.pT();
          return ptsum / y.pT() > 0.1;
        });
      ifilter_discard(elecs, [&](const Particle& e) {
          double ptsum = -e.pT();
          for (const Particle& p : calofs)
            if (deltaR(p,e) < 0.3) ptsum += p.pT();
          return ptsum / e.pT() > 0.1;
        });
      ifilter_discard(muons, [&](const Particle& m) {
          double ptsum = -m.pT();
          for (const Particle& p : calofs)
            if (deltaR(p,m) < 0.3) ptsum += p.pT();
          return ptsum / m.pT() > 0.2;
        });

      // Veto the event if there are any remaining baseline photons or leptons
      if (!photons.empty()) vetoEvent;
      if (!elecs.empty()) vetoEvent;
      if (!muons.empty()) vetoEvent;


      // Get jets and apply jet-based event-selection cuts
      const JetAlg& jetproj = apply<JetAlg>(event, "Jets");
      const Jets alljets = jetproj.jetsByPt(Cuts::abseta < 3.0 && Cuts::Et > 37*GeV); //< most inclusive jets requirement
      if (filter_select(alljets, Cuts::Et > 73*GeV).size() < 2) vetoEvent; //< most inclusive lead jets requirement

      // Filter jets into different Et requirements & compute corresponding HTs
      /// @note It's not clear if different HTs are used to choose the HT bins
      const Jets jets37 = filter_select(alljets, Cuts::Et > 37*GeV);
      const Jets jets43 = filter_select(jets37, Cuts::Et > 43*GeV);
      const Jets jets50 = filter_select(jets43, Cuts::Et > 50*GeV);
      const double ht37 = sum(jets37, Et, 0.0);
      const double ht43 = sum(jets43, Et, 0.0);
      const double ht50 = sum(jets50, Et, 0.0);

      // Find the relevant HT bin and apply leading jet event-selection cuts
      static const vector<double> htcuts = { /* 275., 325., */ 375., 475., 575., 675., 775., 875.}; //< comment to avoid jets50 "fall-down"
      const int iht = inRange(ht37, 275*GeV, 325*GeV) ? 0 : inRange(ht43, 325*GeV, 375*GeV) ? 1 : (2+binIndex(ht50, htcuts, true));
      MSG_TRACE("HT = {" << ht37 << ", " << ht43 << ", " << ht50 << "} => IHT = " << iht);
      if (iht < 0) vetoEvent;
      if (iht == 1 && filter_select(jets43, Cuts::Et > 78*GeV).size() < 2) vetoEvent;
      if (iht >= 2 && filter_select(jets50, Cuts::Et > 100*GeV).size() < 2) vetoEvent;

      // Create references for uniform access to relevant set of jets & HT
      const double etcut = iht == 0 ? 37. : iht == 1 ? 43. : 50.;
      const double& ht = iht == 0 ? ht37 : iht == 1 ? ht43 : ht50;
      const Jets& jets = iht == 0 ? jets37 : iht == 1 ? jets43 : jets50;
      if (!jetproj.jets(Cuts::abseta > 3 && Cuts::Et > etcut*GeV).empty()) vetoEvent;
      const size_t nj = jets.size();
      const size_t nb = count_if(jets.begin(), jets.end(), [](const Jet& j) { return j.bTagged(Cuts::pT > 5*GeV); });

      // Compute HTmiss = pT of 4-vector sum of jet momenta
      const FourMomentum jsum = sum(jets, mom, FourMomentum());
      const double htmiss = jsum.pT();

      // Require HTmiss / ETmiss < 1.25
      const double etmiss = apply<SmearedMET>(event, "MET").met();
      if (htmiss/etmiss > 1.25) vetoEvent;

      // Compute DeltaHT = minimum difference of "dijet" ETs, i.e. max(|1+2-3|, |1+3-2|, |2+3-1|)
      double deltaht = -1;
      vector<double> jetets; transform(jets, jetets, Et);
      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
        const bitset<10> bits(i); /// @warning There'd better not be more than 10 jets...
        const double htdiff = partition_diff(bits, jetets);
        // MSG_INFO(bits.to_string() << " => " << htdiff);
        if (deltaht < 0 || htdiff < deltaht) deltaht = htdiff;
      }
      MSG_DEBUG("dHT_bitmask = " << deltaht);

      // Cross-check calculation in 2- and 3-jet cases
      // if (jets.size() == 2) {
      //   MSG_INFO("dHT2 = " << fabs(jets[0].Et() - jets[1].Et()));
      // } else if (jets.size() == 3) {
      //   double deltaht_01_2 = fabs(jets[0].Et()+jets[1].Et()-jets[2].Et());
      //   double deltaht_02_1 = fabs(jets[0].Et()+jets[2].Et()-jets[1].Et());
      //   double deltaht_12_0 = fabs(jets[1].Et()+jets[2].Et()-jets[0].Et());
      //   MSG_INFO("dHT3 = " << min({deltaht_01_2, deltaht_02_1, deltaht_12_0}));
      // }

      // Compute alphaT from the above
      double alphaT = fabs(0.5*((ht-deltaht)/(sqrt((ht*ht)-(htmiss*htmiss)))));
      if (alphaT < 0.55) vetoEvent;

      /// @todo Need to include trigger efficiency sampling or weighting?

      // Fill histograms
      const double weight = event.weight();
      const size_t inj = nj < 4 ? 0 : 1;
      const size_t inb = nb < 4 ? nb : 4;
      if (iht >= 2)
        (inj == 0 ? _h_alphaT23 : _h_alphaT4)->fill(alphaT, weight);

      // Fill the appropriate counter -- after working out the irregular SR bin index! *sigh*
      size_t i = 8 * ((inj == 0 ? 0 : 3) + inb) + iht;
      if (inj == 1 && inb == 4) i = 8*7 + (iht < 3 ? iht : 2);
      MSG_INFO("inj = " << inj << ", inb = " << inb << ", i = " << i);
      _h_srcounters[i]->fill(weight);

    }


    /// Normalise histograms etc., after the run
    void finalize() {

      const double sf = crossSection()/femtobarn*11.7/sumOfWeights();
      scale({_h_alphaT23,_h_alphaT4}, sf);
      for (size_t i = 0; i < 8*7+3; ++i)
        scale(_h_srcounters[i], sf);

    }

    //@}


    /// @name Utility functions for partitioning jet pTs into two groups and summing/diffing them
    //@{

    /// Sum the given values into two subsets according to the provided bitmask
    template <size_t N>
    pair<double, double> partition_sum(const bitset<N>& mask, const vector<double>& vals) const {
      pair<double, double> rtn(0., 0.);
      for (size_t i = 0; i < vals.size(); ++i) {
        (!mask[vals.size()-1-i] ? rtn.first : rtn.second) += vals[i];
      }
      return rtn;
    }

    /// Return the difference between summed subsets according to the provided bitmask
    template <size_t N>
    double partition_diff(const bitset<N>& mask, const vector<double>& vals) const {
      const pair<double, double> sums = partition_sum(mask, vals);
      const double diff = fabs(sums.first - sums.second);
      MSG_TRACE(mask.to_string() << ": " << sums.first << "/" << sums.second << " => " << diff);
      return diff;
    }

    //@}


    /// @name Histograms
    //@{
    Histo1DPtr _h_alphaT23, _h_alphaT4;
    vector<CounterPtr> _h_srcounters;
    //@}


  };


  // The hook for the plugin system
  DECLARE_RIVET_PLUGIN(CMS_2013_I1223519);


}