1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/WFinder.hh"
  5#include "Rivet/Projections/VetoedFinalState.hh"
  6#include "Rivet/Projections/FastJets.hh"
  7#include "Rivet/Projections/HeavyHadrons.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief ATLAS W+b measurement
 13  class ATLAS_2013_I1219109: public Analysis {
 14  public:
 15
 16    ///@brief: Electroweak Wjj production at 8 TeV
 17    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2013_I1219109);
 18    //@}
 19
 20    void init() {
 21
 22      // Get options from the new option system
 23      _mode = 0;
 24      if ( getOption("LMODE") == "EL" ) _mode = 1;
 25      if ( getOption("LMODE") == "MU" ) _mode = 2;
 26
 27      const FinalState fs;
 28
 29      Cut cuts = Cuts::abseta < 2.5 && Cuts::pT >= 25*GeV;
 30
 31      // W finder for electrons and muons
 32      WFinder wf_mu(fs, cuts, PID::MUON, 0.0*GeV, DBL_MAX, 0.0, 0.1,
 33                 WFinder::ChargedLeptons::PROMPT, WFinder::ClusterPhotons::NODECAY, WFinder::AddPhotons::NO, WFinder::MassWindow::MT);
 34      WFinder wf_el(fs, cuts, PID::ELECTRON, 0.0*GeV, DBL_MAX, 0.0, 0.1,
 35                 WFinder::ChargedLeptons::PROMPT, WFinder::ClusterPhotons::NODECAY, WFinder::AddPhotons::NO, WFinder::MassWindow::MT);
 36      declare(wf_mu, "WFmu");
 37      declare(wf_el, "WFel");
 38
 39      // jets
 40      VetoedFinalState jet_fs(fs);
 41      jet_fs.addVetoOnThisFinalState(wf_el);
 42      jet_fs.addVetoOnThisFinalState(wf_mu);
 43      FastJets fj(jet_fs, FastJets::ANTIKT, 0.4);
 44      fj.useInvisibles();
 45      declare(fj, "Jets");
 46      declare(HeavyHadrons(Cuts::abseta < 2.5 && Cuts::pT > 5*GeV), "BHadrons");
 47
 48
 49      // book histograms
 50      book(_njet     ,1, 1, 1); // dSigma / dNjet
 51      book(_jet1_bPt ,3, 1, 1); // dSigma / dBjetPt for Njet = 1
 52      book(_jet2_bPt ,8, 1, 1); // dSigma / dBjetPt for Njet = 2
 53
 54    }
 55
 56
 57    void analyze(const Event& event) {
 58
 59      //  retrieve W boson candidate
 60      const WFinder& wf_mu = apply<WFinder>(event, "WFmu");
 61      const WFinder& wf_el = apply<WFinder>(event, "WFel");
 62
 63      size_t nWmu = wf_mu.size();
 64      size_t nWel = wf_el.size();
 65
 66      if (_mode == 0 && !((nWmu == 1 && !nWel) || (!nWmu && nWel == 1)))  vetoEvent; // one W->munu OR W->elnu candidate, otherwise veto
 67      if (_mode == 1 && !(!nWmu && nWel == 1))  vetoEvent; // one W->elnu candidate, otherwise veto
 68      if (_mode == 2 && !(nWmu == 1 && !nWel))  vetoEvent; // one W->munu candidate, otherwise veto
 69
 70
 71      if (   (nWmu? wf_mu : wf_el).bosons().size() != 1 )  vetoEvent; // only one W boson candidate
 72      if ( !((nWmu? wf_mu : wf_el).mT() > 60.0*GeV) )      vetoEvent;
 73      //const Particle& Wboson  = wf.boson();
 74
 75
 76      // retrieve constituent neutrino
 77      const Particle& neutrino = (nWmu? wf_mu : wf_el).constituentNeutrino();
 78      if( !(neutrino.pT() > 25*GeV) )  vetoEvent;
 79
 80      // retrieve constituent lepton
 81      const Particle& lepton = (nWmu? wf_mu : wf_el).constituentLepton();
 82
 83      // count good jets, check if good jet contains B hadron
 84      const Particles& bHadrons = apply<HeavyHadrons>(event, "BHadrons").bHadrons();
 85      const Jets& jets = apply<JetAlg>(event, "Jets").jetsByPt(25*GeV);
 86      int goodjets = 0, bjets = 0;
 87      double bPt = 0.;
 88      for(const Jet& j : jets) {
 89        if( (j.abseta() < 2.1) && (deltaR(lepton, j) > 0.5) ) {
 90          // this jet passes the selection!
 91          ++goodjets;
 92          // j.bTagged() uses ghost association which is
 93          // more elegant, but not what has been used in
 94          // this analysis originally, will match B had-
 95          // rons in eta-phi space instead
 96          for(const Particle& b : bHadrons) {
 97            if( deltaR(j, b) < 0.3 ) {
 98              // jet matched to B hadron!
 99              if(!bPt)  bPt = j.pT() * GeV; // leading b-jet pT
100              ++bjets; // count number of b-jets
101              break;
102            }
103          }
104        }
105      }
106      if( goodjets > 2 )  vetoEvent; // at most two jets
107      if( !bjets )  vetoEvent; // at least one of them b-tagged
108
109      double njets = double(goodjets);
110      double ncomb = 3.0;
111      _njet->fill(njets);
112      _njet->fill(ncomb);
113
114      if(     goodjets == 1)  _jet1_bPt->fill(bPt);
115      else if(goodjets == 2)  _jet2_bPt->fill(bPt);
116    }
117
118
119    void finalize() {
120      const double sf = _mode? 1.0 : 0.5;
121      const double xs_pb = sf * crossSection() / picobarn  / sumOfWeights();
122      const double xs_fb = sf * crossSection() / femtobarn / sumOfWeights();
123      scale(_njet,     xs_pb);
124      scale(_jet1_bPt, xs_fb);
125      scale(_jet2_bPt, xs_fb);
126    }
127
128  protected:
129
130    size_t _mode;
131
132  private:
133
134    Histo1DPtr _njet;
135    Histo1DPtr _jet1_bPt;
136    Histo1DPtr _jet2_bPt;
137
138    //bool _isMuon;
139
140  };
141
142
143  // The hook for the plugin system
144  RIVET_DECLARE_PLUGIN(ATLAS_2013_I1219109);
145
146}