1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/DileptonFinder.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/VetoedFinalState.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// Z + jets in pp at 13 TeV
 11  ///
 12  /// @note This base class contains a "mode" variable for combined, e, and mu channel derived classes
 13  class ATLAS_2015_CONF_2015_041 : public Analysis {
 14  public:
 15
 16    /// @name Constructors etc.
 17    /// @{
 18
 19    /// Constructor
 20    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2015_CONF_2015_041);
 21
 22    /// @}
 23
 24
 25    /// Book histograms and initialise projections before the run
 26    void init() {
 27
 28
 29      _mode = 0;
 30      if ( getOption("LMODE") == "EL" ) _mode = 1;
 31      else if ( getOption("LMODE") == "MU" ) _mode = 2;
 32
 33      Cut cuts = Cuts::pT > 25*GeV && Cuts::abseta < 2.5;
 34      DileptonFinder eefinder(91.2*GeV, 0.1, cuts && Cuts::abspid == PID::ELECTRON,
 35                             Cuts::massIn(66*GeV, 116*GeV));
 36      DileptonFinder mmfinder(91.2*GeV, 0.1, cuts && Cuts::abspid == PID::MUON,
 37                             Cuts::massIn(66*GeV, 116*GeV));
 38      declare(eefinder, "eefinder");
 39      declare(mmfinder, "mmfinder");
 40
 41      // Define veto FS in order to prevent Z-decay products entering the jet algorithm
 42      VetoedFinalState had_fs;
 43      had_fs.addVetoOnThisFinalState(eefinder);
 44      had_fs.addVetoOnThisFinalState(mmfinder);
 45      FastJets jets(had_fs, JetAlg::ANTIKT, 0.4, JetMuons::ALL, JetInvisibles::DECAY);
 46      declare(jets, "jets");
 47
 48      // combination
 49      book(_h["Njets_comb"], 1, 1, 1);
 50      book(_e["ratio_comb"], 2, 1, 1);
 51      // individual channels
 52      if (_mode == 0 || _mode == 1) {
 53        book(_h["Njets_el"], 1, 2, 1);
 54        book(_e["ratio_el"], 2, 2, 1);
 55      }
 56      if (_mode == 0 || _mode == 2) {
 57        book(_h["Njets_mu"], 1, 2, 2);
 58        book(_e["ratio_mu"], 2, 2, 2);
 59      }
 60
 61    }
 62
 63
 64    /// Perform the per-event analysis
 65    void analyze(const Event& event) {
 66
 67      const Particles& elecs = apply<DileptonFinder>(event, "eefinder").constituents();
 68      const Particles& muons = apply<DileptonFinder>(event, "mmfinder").constituents();
 69
 70      if (elecs.size() + muons.size() != 2)  vetoEvent;
 71
 72      Jets jets = apply<JetFinder>(event, "jets").jetsByPt(Cuts::pT > 30*GeV && Cuts::absrap < 2.5);
 73      idiscardIfAnyDeltaRLess(jets, elecs, 0.4);
 74      idiscardIfAnyDeltaRLess(jets, muons, 0.4);
 75
 76      fillHistos(elecs, jets);
 77      fillHistos(muons, jets);
 78    }
 79
 80    void fillHistos(const Particles& leptons, const Jets& jets) {
 81      if (leptons.size() != 2)  return;
 82
 83      const size_t njets = jets.size();
 84      for (size_t i = 0; i <= njets; ++i) {
 85        if (_mode == 0 || _mode == 1) {
 86          _h["Njets_comb"]->fill(i + 0.5);
 87          _h["Njets_el"]->fill(i + 0.5);
 88        }
 89        if (_mode == 0 || _mode == 2) {
 90          _h["Njets_comb"]->fill(i + 0.5);
 91          _h["Njets_mu"]->fill(i + 0.5);
 92        }
 93      }
 94
 95    }
 96
 97    void finalize() {
 98      scale(_h, crossSectionPerEvent());
 99      if (_mode == 0)  scale(_h["Njets_comb"], 0.5); // average of el + mu
100
101      vector<string> channels = { "comb" };
102      if (_mode == 0 || _mode == 1)  channels.push_back("el");
103      if (_mode == 0 || _mode == 2)  channels.push_back("mu");
104      YODA::Estimate num, den;
105      for (size_t i = 1; i < 5; ++i) {
106        for (const string& ch : channels) {
107          num.set(_h["Njets_"+ch]->bin(i+1).sumW(), _h["Njets_"+ch]->bin(i+1).errW());
108          den.set(_h["Njets_"+ch]->bin(i).sumW(),   _h["Njets_"+ch]->bin(i).errW());
109          _e["ratio_"+ch]->bin(i) = YODA::divide(num, den);
110        }
111      }
112
113    }
114
115    /// @}
116
117
118  private:
119
120    size_t _mode;
121
122    map<string, Estimate1DPtr> _e;
123    map<string, Histo1DPtr> _h;
124  };
125
126
127  RIVET_DECLARE_PLUGIN(ATLAS_2015_CONF_2015_041);
128}