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}