rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2016_I1458270

0-lepton SUSY search with 3.2/fb of 13 TeV $pp$ data
Experiment: ATLAS (LHC)
Inspire ID: 1458270
Status: VALIDATED, UNOFFICIAL
Authors:
  • Andy Buckley
  • Louie Corpe
References:
  • Eur.Phys.J. C76 (2016) no.7, 392
  • DOI:10.1140/epjc/s10052-016-4184-8
  • arXiv: 1605.03814
Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • BSM signal events

ATLAS 0-lepton SUSY search using 3.2/fb of LHC $pp$ data at 13 TeV, recorded in 2015. The event selection is via signal regions requiring from 2-6 high-energy jets, and significant missing transverse energy. Detailed info: http://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PAPERS/SUSY-2015-06/

Source code: ATLAS_2016_I1458270.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  /// @brief ATLAS 0-lepton SUSY search with 3.2/fb of 13 TeV pp data
 16  class ATLAS_2016_I1458270 : public Analysis {
 17  public:
 18
 19    /// Constructor
 20    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2016_I1458270);
 21
 22
 23    /// @name Analysis methods
 24    //@{
 25
 26    // method to turn Hist1D into Scatter... so we can write this out witout dividing by bin width
 27    // since the HEPData entry corresponding to this does not divide the refData by bin width!
 28    // Have requested they update their HEPData entry but until they do so, we use this workaround.
 29    Scatter2DPtr convertToScatterWithoutBinWidthDivision(Histo1DPtr input , Scatter2DPtr output ){
 30      for (size_t b = 0; b < input->numBins(); ++b) {
 31            const double x   = input->bin(b).xMid();
 32            const double ex  = input->bin(b).xWidth()/2.;
 33            const double val = input->bin(b).sumW();
 34            const double err = input->bin(b).relErr() * val;
 35            output->addPoint(x, val, ex, err);
 36       }
 37       return output;
 38    }
 39
 40    /// Book histograms and initialise projections before the run
 41    void init() {
 42
 43      // Initialise and register projections
 44      FinalState calofs(Cuts::abseta < 4.8);
 45      FastJets fj(calofs, FastJets::ANTIKT, 0.4);
 46      declare(fj, "TruthJets");
 47      declare(SmearedJets(fj, JET_SMEAR_ATLAS_RUN2, JET_BTAG_ATLAS_RUN2_MV2C20), "RecoJets");
 48
 49      MissingMomentum mm(calofs);
 50      declare(mm, "TruthMET");
 51      declare(SmearedMET(mm, MET_SMEAR_ATLAS_RUN2), "RecoMET");
 52
 53      PromptFinalState es(Cuts::abseta < 2.47 && Cuts::abspid == PID::ELECTRON, true, true);
 54      declare(es, "TruthElectrons");
 55      declare(SmearedParticles(es, ELECTRON_RECOEFF_ATLAS_RUN2, ELECTRON_SMEAR_ATLAS_RUN2), "RecoElectrons");
 56
 57      PromptFinalState mus(Cuts::abseta < 2.7 && Cuts::abspid == PID::MUON, true);
 58      declare(mus, "TruthMuons");
 59      declare(SmearedParticles(mus, MUON_EFF_ATLAS_RUN2, MUON_SMEAR_ATLAS_RUN2), "RecoMuons");
 60
 61
 62      // Book histograms/counters
 63      book(_h_2jl, "2jl");
 64      book(_h_2jm, "2jm");
 65      book(_h_2jt, "2jt");
 66      book(_h_4jt, "4jt");
 67      book(_h_5j , "5j");
 68      book(_h_6jm, "6jm");
 69      book(_h_6jt, "6jt");
 70
 71      book(_hMeff_2jl, 4,1,1);
 72      book(_hMeff_2jm, 5,1,1);
 73      book(_hMeff_2jt, 6,1,1);
 74      book(_hMeff_4jt, 7,1,1);
 75      book(_hMeff_5j , 8,1,1);
 76      book(_hMeff_6jm, 9,1,1);
 77      book(_hMeff_6jt, 10,1,1);
 78
 79      book(_h_temp_Meff_2jl, "_temp_Meff_2jl",refData( 4,1,1));
 80      book(_h_temp_Meff_2jm, "_temp_Meff_2jm",refData( 5,1,1));
 81      book(_h_temp_Meff_2jt, "_temp_Meff_2jt",refData( 6,1,1));
 82      book(_h_temp_Meff_4jt, "_temp_Meff_4jt",refData( 7,1,1));
 83      book(_h_temp_Meff_5j , "_temp_Meff_5j" ,refData( 8,1,1));
 84      book(_h_temp_Meff_6jm, "_temp_Meff_6jm",refData( 9,1,1));
 85      book(_h_temp_Meff_6jt, "_temp_Meff_6jt",refData(10,1,1));
 86
 87
 88
 89
 90      // Book cut-flows
 91      const vector<string> cuts2j = {"Pre-sel+MET+pT1", "Njet", "Dphi_min(j,MET)", "pT2", "MET/sqrtHT", "m_eff(incl)"};
 92      _flows.addCutflow("2jl", cuts2j);
 93      _flows.addCutflow("2jm", cuts2j);
 94      _flows.addCutflow("2jt", cuts2j);
 95      const vector<string> cutsXj = {"Pre-sel+MET+pT1", "Njet", "Dphi_min(j,MET)", "pT2", "pT4", "Aplanarity", "MET/m_eff(Nj)", "m_eff(incl)"};
 96      _flows.addCutflow("4jt", cutsXj);
 97      _flows.addCutflow("5j",  cutsXj);
 98      _flows.addCutflow("6jm", cutsXj);
 99      _flows.addCutflow("6jt", cutsXj);
100
101    }
102
103    /// Perform the per-event analysis
104    void analyze(const Event& event) {
105
106      _flows.fillinit();
107
108      // Same MET cut for all signal regions
109      //const Vector3 vmet = -apply<MissingMomentum>(event, "TruthMET").vectorEt();
110      const Vector3 vmet = -apply<SmearedMET>(event, "RecoMET").vectorEt();
111      const double met = vmet.mod();
112      if (met < 200*GeV) vetoEvent;
113
114      // Get baseline electrons, muons, and jets
115      Particles elecs = apply<ParticleFinder>(event, "RecoElectrons").particles(Cuts::pT > 10*GeV);
116      Particles muons = apply<ParticleFinder>(event, "RecoMuons").particles(Cuts::pT > 10*GeV);
117      Jets jets = apply<JetAlg>(event, "RecoJets").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 2.8); ///< @todo Pile-up subtraction
118
119      // Jet/electron/muons overlap removal and selection
120      // Remove any |eta| < 2.8 jet within dR = 0.2 of a baseline electron
121      for (const Particle& e : elecs)
122        ifilter_discard(jets, deltaRLess(e, 0.2, RAPIDITY));
123      // Remove any electron or muon with dR < 0.4 of a remaining (Nch > 3) jet
124      for (const Jet& j : jets) {
125        /// @todo Add track efficiency random filtering
126        ifilter_discard(elecs, deltaRLess(j, 0.4, RAPIDITY));
127        if (j.particles(Cuts::abscharge > 0 && Cuts::pT > 500*MeV).size() >= 3)
128          ifilter_discard(muons, deltaRLess(j, 0.4, RAPIDITY));
129      }
130      // Discard the softer of any electrons within dR < 0.05
131      for (size_t i = 0; i < elecs.size(); ++i) {
132        const Particle& e1 = elecs[i];
133        /// @todo Would be nice to pass a "tail view" for the filtering, but awkward without range API / iterator guts
134        ifilter_discard(elecs, [&](const Particle& e2){ return e2.pT() < e1.pT() && deltaR(e1,e2) < 0.05; });
135      }
136
137      // Loose electron selection
138      ifilter_select(elecs, ParticleEffFilter(ELECTRON_EFF_ATLAS_RUN2_LOOSE));
139
140      // Veto the event if there are any remaining baseline leptons
141      if (!elecs.empty()) vetoEvent;
142      if (!muons.empty()) vetoEvent;
143
144      // Signal jets have pT > 50 GeV
145      const Jets jets50 = filter_select(jets, Cuts::pT > 50*GeV);
146      if (jets50.size() < 2) vetoEvent;
147      vector<double> jetpts; transform(jets, jetpts, pT);
148      vector<double> jetpts50; transform(jets50, jetpts50, pT);
149      const double j1pt = jetpts50[0];
150      const double j2pt = jetpts50[1];
151      if (j1pt < 200*GeV) vetoEvent;
152
153      // Construct multi-jet observables
154      const double ht = sum(jetpts, 0.0);
155      const double met_sqrt_ht = met / sqrt(ht);
156      const double meff_incl = sum(jetpts50, met);
157
158      // Get dphis between MET and jets
159      vector<double> dphimets50; transform(jets50, dphimets50, deltaPhiWRT(vmet));
160      const double min_dphi_met_3 = min(head(dphimets50, 3));
161      MSG_DEBUG(dphimets50 << ", " << min_dphi_met_3);
162
163      // Jet aplanarity
164      Sphericity sph; sph.calc(jets);
165      const double aplanarity = sph.aplanarity();
166
167
168      // Fill SR counters
169      // 2-jet SRs
170      if (_flows["2jl"].filltail({true, true, min_dphi_met_3 > 0.8, j2pt > 200*GeV,
171            met_sqrt_ht > 15*sqrt(GeV), meff_incl > 1200*GeV})) _h_2jl->fill();
172      if (_flows["2jm"].filltail({j1pt > 300*GeV, true, min_dphi_met_3 > 0.4, j2pt > 50*GeV,
173            met_sqrt_ht > 15*sqrt(GeV), meff_incl > 1600*GeV})) _h_2jm->fill();
174      if (_flows["2jt"].filltail({true, true, min_dphi_met_3 > 0.8, j2pt > 200*GeV,
175            met_sqrt_ht > 20*sqrt(GeV), meff_incl > 2000*GeV})) _h_2jt->fill();
176
177      // Fill SR Meff Histo1Ds
178      // 2-jet SRs
179      if ((min_dphi_met_3 > 0.8) && (j2pt > 200*GeV) && (met_sqrt_ht > 15*sqrt(GeV))) _h_temp_Meff_2jl->fill(meff_incl);
180      if ((j1pt > 300*GeV)  && (min_dphi_met_3 > 0.4) && (j2pt > 50*GeV) && (met_sqrt_ht > 15*sqrt(GeV))) _h_temp_Meff_2jm->fill(meff_incl);
181      if ((min_dphi_met_3 > 0.8) && (j2pt > 200*GeV ) && (met_sqrt_ht > 20*sqrt(GeV))) _h_temp_Meff_2jt->fill(meff_incl);
182
183      // Upper multiplicity SRs
184      const double j4pt = jets50.size() > 3 ? jetpts50[3] : -1;
185      const double j5pt = jets50.size() > 4 ? jetpts50[4] : -1;
186      const double j6pt = jets50.size() > 5 ? jetpts50[5] : -1;
187      const double meff_4 = jets50.size() > 3 ? sum(head(jetpts50, 4), met) : -1;
188      const double meff_5 = jets50.size() > 4 ? meff_4 + jetpts50[4] : -1;
189      const double meff_6 = jets50.size() > 5 ? meff_5 + jetpts50[5] : -1;
190      const double met_meff_4 = met / meff_4;
191      const double met_meff_5 = met / meff_5;
192      const double met_meff_6 = met / meff_6;
193      const double min_dphi_met_more = jets50.size() > 3 ? min(tail(dphimets50, -3)) : -1;
194
195
196      if (_flows["4jt"].filltail({true, jets50.size() >= 4, min_dphi_met_3 > 0.4 && min_dphi_met_more > 0.2,
197            jetpts[1] > 100*GeV, j4pt > 100*GeV, aplanarity > 0.04, met_meff_4 > 0.20, meff_incl > 2200*GeV}))
198        _h_4jt->fill();
199      if (_flows["5j"].filltail({true, jets50.size() >= 5, min_dphi_met_3 > 0.4 && min_dphi_met_more > 0.2,
200            jetpts[1] > 100*GeV, j4pt > 100*GeV && j5pt > 50*GeV, aplanarity > 0.04, met_meff_5 > 0.25, meff_incl > 1600*GeV}))
201        _h_5j->fill();
202      if (_flows["6jm"].filltail({true, jets50.size() >= 6, min_dphi_met_3 > 0.4 && min_dphi_met_more > 0.2,
203            jetpts[1] > 100*GeV, j4pt > 100*GeV && j6pt > 50*GeV, aplanarity > 0.04, met_meff_6 > 0.25, meff_incl > 1600*GeV}))
204        _h_6jm->fill();
205      if (_flows["6jt"].filltail({true, jets50.size() >= 6, min_dphi_met_3 > 0.4 && min_dphi_met_more > 0.2,
206            jetpts[1] > 100*GeV, j4pt > 100*GeV && j6pt > 50*GeV, aplanarity > 0.04, met_meff_6 > 0.20, meff_incl > 2000*GeV}))
207        _h_6jt->fill();
208
209      // Fill SR Meff Histo1Ds
210      // Upper multiplicity SRs
211      if (((jets50.size() >= 4) && (min_dphi_met_3 > 0.4) && (min_dphi_met_more > 0.2) &&  (jetpts[1] > 100*GeV) && (j4pt > 100*GeV) && (aplanarity > 0.04) && (met_meff_4 > 0.20))) _h_temp_Meff_4jt->fill(meff_incl);
212      if (((jets50.size() >= 5) && (min_dphi_met_3 > 0.4) && (min_dphi_met_more > 0.2 ) &&
213            (jetpts[1] > 100*GeV) && (j4pt > 100*GeV) && (j5pt > 50*GeV) && (aplanarity > 0.04) && (met_meff_5 > 0.25))) _h_temp_Meff_5j->fill(meff_incl);
214      if (((jets50.size() >= 6) && (min_dphi_met_3 > 0.4) && (min_dphi_met_more > 0.2) &&
215            (jetpts[1] > 100*GeV) && (j4pt > 100*GeV) && (j6pt > 50*GeV) && (aplanarity > 0.04) && (met_meff_6 > 0.25))) _h_temp_Meff_6jm->fill(meff_incl);
216      if (((jets50.size() >= 6) && (min_dphi_met_3 > 0.4) && (min_dphi_met_more > 0.2) &&
217            (jetpts[1] > 100*GeV) && (j4pt > 100*GeV) && (j6pt > 50*GeV) && (aplanarity > 0.04) && (met_meff_6 > 0.20))) _h_temp_Meff_6jt->fill(meff_incl);
218
219    }
220
221
222    /// Normalise histograms etc., after the run
223    void finalize() {
224
225      const double sf = 3.2*crossSection()/femtobarn/sumOfWeights();
226      scale(_h_2jl, sf); scale(_h_2jm, sf); scale(_h_2jt, sf);
227      scale(_h_4jt, sf); scale(_h_5j, sf);
228      scale(_h_6jm, sf); scale(_h_6jt, sf);
229
230      scale(_h_temp_Meff_2jl, sf); scale(_h_temp_Meff_2jm, sf); scale(_h_temp_Meff_2jt, sf);
231      scale(_h_temp_Meff_4jt, sf); scale(_h_temp_Meff_5j, sf);
232      scale(_h_temp_Meff_6jm, sf); scale(_h_temp_Meff_6jt, sf);
233
234
235      // the HEPData entry corresponding to this does not divide their distributions
236      // by bin width... so to avoid this we need to convert to Scatter2D which is not divided by bw
237      _hMeff_2jl = convertToScatterWithoutBinWidthDivision(_h_temp_Meff_2jl,_hMeff_2jl);
238      _hMeff_2jm = convertToScatterWithoutBinWidthDivision(_h_temp_Meff_2jm,_hMeff_2jm);
239      _hMeff_2jt = convertToScatterWithoutBinWidthDivision(_h_temp_Meff_2jt,_hMeff_2jt);
240      _hMeff_4jt = convertToScatterWithoutBinWidthDivision(_h_temp_Meff_4jt,_hMeff_4jt);
241      _hMeff_5j  = convertToScatterWithoutBinWidthDivision(_h_temp_Meff_5j ,_hMeff_5j ) ;
242      _hMeff_6jm = convertToScatterWithoutBinWidthDivision(_h_temp_Meff_6jm,_hMeff_6jm);
243      _hMeff_6jt = convertToScatterWithoutBinWidthDivision(_h_temp_Meff_6jt,_hMeff_6jt);
244      MSG_INFO("CUTFLOWS:\n\n" << _flows);
245
246    }
247
248    //@}
249
250
251    private:
252
253    /// @name Histograms
254    //@{
255    CounterPtr _h_2jl, _h_2jm, _h_2jt;
256    CounterPtr _h_4jt, _h_5j;
257    CounterPtr _h_6jm, _h_6jt;
258
259    Scatter2DPtr  _hMeff_2jl, _hMeff_2jm, _hMeff_2jt;
260    Scatter2DPtr  _hMeff_4jt, _hMeff_5j;
261    Scatter2DPtr  _hMeff_6jm, _hMeff_6jt;
262
263    Histo1DPtr _h_temp_Meff_2jl, _h_temp_Meff_2jm, _h_temp_Meff_2jt;
264    Histo1DPtr _h_temp_Meff_4jt, _h_temp_Meff_5j;
265    Histo1DPtr _h_temp_Meff_6jm, _h_temp_Meff_6jt;
266    //@}
267
268    /// Cut-flows
269    Cutflows _flows;
270
271  };
272
273
274
275  // The hook for the plugin system
276  RIVET_DECLARE_PLUGIN(ATLAS_2016_I1458270);
277
278
279}