1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/PromptFinalState.hh"
  5#include "Rivet/Projections/VetoedFinalState.hh"
  6#include "Rivet/Projections/MissingMomentum.hh"
  7#include "Rivet/Projections/LeptonFinder.hh"
  8#include "Rivet/Projections/FastJets.hh"
  9#include "Rivet/Projections/HeavyHadrons.hh"
 10
 11namespace Rivet {
 12
 13
 14  /// @brief ATLAS W+b measurement
 15  class ATLAS_2013_I1219109: public Analysis {
 16  public:
 17
 18    /// Constructor
 19    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2013_I1219109);
 20
 21
 22    void init() {
 23
 24      // Get options from the option system
 25      _mode = 0;
 26      if ( getOption("LMODE") == "EL" ) _mode = 1;
 27      if ( getOption("LMODE") == "MU" ) _mode = 2;
 28
 29
 30      // Initialise and register projections
 31      declare("MET", MissingMomentum());
 32      Cut cuts = Cuts::abseta < 2.5 && Cuts::pT > 25*GeV;
 33      LeptonFinder ef(0.1, cuts && Cuts::abspid == PID::ELECTRON);
 34      declare(ef, "Elecs");
 35      LeptonFinder mf(0.1, cuts && Cuts::abspid == PID::MUON);
 36      declare(mf, "Muons");
 37
 38      // Jets
 39      VetoedFinalState jet_fs;
 40      jet_fs.addVetoOnThisFinalState(ef);
 41      jet_fs.addVetoOnThisFinalState(mf);
 42      FastJets fj(jet_fs, JetAlg::ANTIKT, 0.4);
 43      fj.useInvisibles();
 44      declare(fj, "Jets");
 45      declare(HeavyHadrons(Cuts::abseta < 2.5 && Cuts::pT > 5*GeV), "BHadrons");
 46
 47      // Book histograms
 48      book(_njet     ,1, 1, 1); // dSigma / dNjet
 49      book(_jet1_bPt ,3, 1, 1); // dSigma / dBjetPt for Njet = 1
 50      book(_jet2_bPt ,8, 1, 1); // dSigma / dBjetPt for Njet = 2
 51    }
 52
 53
 54    void analyze(const Event& event) {
 55
 56      // MET cut
 57      const P4& pmiss = apply<MissingMom>(event, "MET").missingMom();
 58      if (pmiss.pT() < 25*GeV) vetoEvent;
 59
 60      // Retrieve W boson candidate
 61      const Particles& es = apply<LeptonFinder>(event, "Elecs").particles();
 62      const int iefound = closestMatchIndex(es, pmiss, Kin::mass, 80.4*GeV);
 63      const Particles& mus = apply<LeptonFinder>(event, "Muons").particles();
 64      const int imfound = closestMatchIndex(mus, pmiss, Kin::mass, 80.4*GeV);
 65
 66      // Event selection and identify W lepton
 67      size_t nWel = (iefound >= 0);
 68      size_t nWmu = (imfound >= 0);
 69      if (_mode == 0 && !((nWmu == 1 && !nWel) || (!nWmu && nWel == 1)))  vetoEvent; // one W->munu OR W->elnu candidate, otherwise veto
 70      if (_mode == 1 && !(!nWmu && nWel == 1))  vetoEvent; // one W->elnu candidate, otherwise veto
 71      if (_mode == 2 && !(nWmu == 1 && !nWel))  vetoEvent; // one W->munu candidate, otherwise veto
 72      const Particle& lepton = (nWel > 0) ? es[iefound] : mus[imfound];
 73
 74      // mT cut
 75      /// @todo Shouldn't this go earlier, in the lepton filtering?
 76      if (mT(lepton, pmiss) < 60*GeV) vetoEvent;
 77
 78      // Count good jets, check if good jet contains B hadron
 79      const Particles& bHadrons = apply<HeavyHadrons>(event, "BHadrons").bHadrons();
 80      const Jets& jets = apply<JetFinder>(event, "Jets").jetsByPt(Cuts::pT > 25*GeV);
 81      int goodjets = 0, bjets = 0;
 82      double bPt = 0.;
 83      for (const Jet& j : jets) {
 84        if ( (j.abseta() < 2.1) && (deltaR(lepton, j) > 0.5) ) {
 85          // this jet passes the selection!
 86          ++goodjets;
 87          // Perform manual dR matching
 88          for (const Particle& b : bHadrons) {
 89            if ( deltaR(j, b) < 0.3 ) {
 90              // jet matched to B hadron!
 91              if (!bPt) bPt = j.pT()/GeV; // leading b-jet pT
 92              ++bjets; // count number of b-jets
 93              break;
 94            }
 95          }
 96        }
 97      }
 98      if (goodjets > 2) vetoEvent; // at most two jets
 99      if (!bjets) vetoEvent; // at least one of them b-tagged
100
101      const double ncomb = 3.0; //< ?
102      _njet->fill(goodjets);
103      _njet->fill(ncomb);
104
105      if (goodjets == 1) _jet1_bPt->fill(bPt);
106      else if (goodjets == 2) _jet2_bPt->fill(bPt);
107    }
108
109
110    void finalize() {
111      const double sf = _mode? 1.0 : 0.5;
112      const double xs_pb = sf * crossSection() / picobarn  / sumOfWeights();
113      const double xs_fb = sf * crossSection() / femtobarn / sumOfWeights();
114      scale(_njet,     xs_pb);
115      scale(_jet1_bPt, xs_fb);
116      scale(_jet2_bPt, xs_fb);
117    }
118
119
120  private:
121
122    size_t _mode;
123
124    Histo1DPtr _njet;
125    Histo1DPtr _jet1_bPt;
126    Histo1DPtr _jet2_bPt;
127
128  };
129
130
131  RIVET_DECLARE_PLUGIN(ATLAS_2013_I1219109);
132
133}