rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

CMS_2016_PAS_SUS_16_14

Search for supersymmetry in events with jets and missing transverse momentum at 13 \text{TeV}
Experiment: CMS (LHC)
Status: VALIDATED
Authors:
  • Andy Buckley
No references listed
Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • BSM physics signal events

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