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  class ATLAS_2017_I1514251 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2017_I1514251);
 16
 17    /// Book histograms and initialise projections before the run
 18    void init() {
 19
 20      // Get options from the new option system
 21      // default to combined.
 22      _mode = 2;
 23      if ( getOption("LMODE") == "EL" ) _mode = 0;
 24      if ( getOption("LMODE") == "MU" ) _mode = 1;
 25      if ( getOption("LMODE") == "EMU" ) _mode = 2;
 26
 27      Cut cuts = Cuts::pT > 25*GeV && Cuts::abseta < 2.5;
 28      DileptonFinder zeefinder(91.2*GeV, 0.1, cuts && Cuts::abspid == PID::ELECTRON, Cuts::massIn(71*GeV, 111*GeV));
 29      DileptonFinder zmumufinder(91.2*GeV, 0.1, cuts && Cuts::abspid == PID::MUON, Cuts::massIn(71*GeV, 111*GeV));
 30      declare(zeefinder, "zeefinder");
 31      declare(zmumufinder, "zmumufinder");
 32
 33      // Define veto FS in order to prevent Z-decay products entering the jet algorithm
 34      VetoedFinalState had_fs;
 35      had_fs.addVetoOnThisFinalState(zeefinder);
 36      had_fs.addVetoOnThisFinalState(zmumufinder);
 37      FastJets jets(had_fs, JetAlg::ANTIKT, 0.4, JetMuons::ALL, JetInvisibles::DECAY);
 38      declare(jets, "jets");
 39
 40      // individual channels
 41      book(_h_Njets_excl,  _mode + 1, 1, 1);
 42      book(_h_Njets,       _mode + 4, 1, 1);
 43      book(_h_Njets_Ratio, _mode + 7, 1, 1);
 44
 45      book(_h_leading_jet_pT_eq1jet, _mode + 10, 1, 1);
 46      book(_h_leading_jet_pT       , _mode + 13, 1, 1);
 47      book(_h_leading_jet_pT_2jet  , _mode + 16, 1, 1);
 48      book(_h_leading_jet_pT_3jet  , _mode + 19, 1, 1);
 49      book(_h_leading_jet_pT_4jet  , _mode + 22, 1, 1);
 50      book(_h_leading_jet_rap      , _mode + 25, 1, 1);
 51      book(_h_HT                   , _mode + 28, 1, 1);
 52      book(_h_jet_dphi             , _mode + 31, 1, 1);
 53      book(_h_jet_mass             , _mode + 34, 1, 1);
 54
 55    }
 56
 57
 58
 59    /// Perform the per-event analysis
 60    void analyze(const Event& event) {
 61
 62      const DileptonFinder& zeefinder = apply<DileptonFinder>(event, "zeefinder");
 63      const DileptonFinder& zmumufinder = apply<DileptonFinder>(event, "zmumufinder");
 64
 65      const Particles& zees = zeefinder.bosons();
 66      const Particles& zmumus = zmumufinder.bosons();
 67
 68      //Veto Z->mumu in electron mode, and vice versa:
 69      if (_mode==0 && (zees.size()!=1 || zmumus.size() ) )  vetoEvent;
 70
 71      if (_mode==1 && (zees.size() || zmumus.size()!=1 ) )  vetoEvent;
 72
 73      if (zees.size() + zmumus.size() != 1) {
 74        // Running in combined mode, we did not find exactly one Z. Not good.
 75        MSG_DEBUG("Did not find exactly one good Z candidate");
 76        vetoEvent;
 77      }
 78
 79      // Find the (dressed!) leptons
 80      const Particles& leptons = zees.size() ? zeefinder.constituents() : zmumufinder.constituents();
 81      if (leptons.size() != 2) vetoEvent;
 82
 83      Jets jets =  apply<JetFinder>(event, "jets").jetsByPt(Cuts::pT > 30*GeV && Cuts::absrap < 2.5);
 84
 85      bool veto = false;
 86      for(const Jet& j : jets)  {
 87        for(const Particle& l : leptons) { veto |= deltaR(j, l) < 0.4; }
 88      }
 89      if (veto) vetoEvent;
 90
 91      double HT=0;
 92      for(const Particle& l : leptons) { HT += l.pT(); }
 93
 94      const size_t Njets = jets.size();
 95      _h_Njets_excl->fill(Njets);
 96      for(size_t i = 0; i <= Njets; ++i) { _h_Njets->fill(i);	}
 97
 98      if (Njets < 1)  vetoEvent;
 99
100
101      for(size_t i = 0; i < Njets; ++i) { HT += jets[i].pT(); }
102      const double pT = jets[0].pT();
103      const double rap = jets[0].rapidity();
104
105      _h_HT->fill(HT);
106      _h_leading_jet_rap->fill(fabs(rap));
107      _h_leading_jet_pT->fill(pT);
108      if (Njets == 1)  _h_leading_jet_pT_eq1jet->fill(pT);
109      if (Njets > 1) {
110        _h_leading_jet_pT_2jet->fill(pT);
111        _h_jet_dphi->fill( deltaPhi(jets[0], jets[1]));
112        _h_jet_mass->fill( (jets[0].momentum()+jets[1].momentum()).mass() );
113      }
114
115      if (Njets > 2)  _h_leading_jet_pT_3jet->fill(pT);
116      if (Njets > 3)  _h_leading_jet_pT_4jet->fill(pT);
117
118    }
119
120    void finalize() {
121      for (size_t i = 1; i < _h_Njets->numBins(); ++i) {
122        double  n = _h_Njets->bin(i + 1).sumW();
123        double dN = _h_Njets->bin(i + 1).sumW2();
124        double  d = _h_Njets->bin(i).sumW();
125        double dD = _h_Njets->bin(i).sumW2();
126        double r = safediv(n, d);
127        double e = sqrt( safediv(r * (1 - r), d) );
128        if ( _h_Njets->effNumEntries() != _h_Njets->numEntries() ) {
129          // use F. James's approximation for weighted events:
130          e = sqrt( safediv((1 - 2 * r) * dN + r * r * dD, d * d) );
131        }
132        _h_Njets_Ratio->bin(i).set(r, e);
133      }
134
135      // when running in combined mode, need to average to get lepton xsec
136      double normfac = crossSectionPerEvent();
137      if (_mode == 2) normfac = 0.5*normfac;
138
139      scale(_h_Njets,                  normfac );
140      scale(_h_Njets_excl,             normfac );
141      scale(_h_HT,                     normfac );
142      scale(_h_leading_jet_rap,        normfac );
143      scale(_h_leading_jet_pT,         normfac );
144      scale(_h_leading_jet_pT_eq1jet,  normfac );
145      scale(_h_leading_jet_pT_2jet,    normfac );
146      scale(_h_leading_jet_pT_3jet,    normfac );
147      scale(_h_leading_jet_pT_4jet,    normfac );
148      scale(_h_jet_dphi,               normfac );
149      scale(_h_jet_mass,               normfac );
150
151    }
152
153
154  protected:
155
156    size_t _mode;
157
158
159  private:
160
161    Estimate1DPtr _h_Njets_Ratio;
162    Histo1DPtr   _h_Njets;
163    Estimate1DPtr _h_Njets_excl_Ratio;
164    Histo1DPtr   _h_Njets_excl;
165    Histo1DPtr   _h_HT;
166    Histo1DPtr   _h_leading_jet_rap;
167    Histo1DPtr   _h_leading_jet_pT;
168    Histo1DPtr   _h_leading_jet_pT_eq1jet;
169    Histo1DPtr   _h_leading_jet_pT_2jet;
170    Histo1DPtr   _h_leading_jet_pT_3jet;
171    Histo1DPtr   _h_leading_jet_pT_4jet;
172    Histo1DPtr   _h_jet_dphi;
173    Histo1DPtr   _h_jet_mass;
174
175  };
176
177
178  RIVET_DECLARE_PLUGIN(ATLAS_2017_I1514251);
179
180}