1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/DileptonFinder.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6#include "Rivet/Projections/HeavyHadrons.hh"
  7#include "Rivet/Projections/VetoedFinalState.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// Electroweak Wjj production at 8 TeV
 13  class ATLAS_2014_I1306294 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2014_I1306294);
 18
 19
 20    /// @name Analysis methods
 21    /// @{
 22
 23    /// Book histograms and initialise projections before the run
 24    void init() {
 25
 26      // Get options from the new option system
 27      _mode = 1;
 28      if ( getOption("LMODE") == "EL" ) _mode = 1;
 29      if ( getOption("LMODE") == "MU" ) _mode = 2;
 30
 31      Cut cuts = Cuts::abseta < 2.5 && Cuts::pT > 20*GeV;
 32      DileptonFinder zfinder(91.2*GeV, 0.1, cuts && Cuts::abspid == (_mode == 1 ? PID::ELECTRON : PID::MUON),
 33                             Cuts::massIn(76.0*GeV, 106.0*GeV),
 34                             LeptonOrigin::ALL, PhotonOrigin::NODECAY);
 35      declare(zfinder, "DileptonFinder");
 36
 37      VetoedFinalState jet_fs;
 38      jet_fs.addVetoOnThisFinalState(getProjection<DileptonFinder>("DileptonFinder"));
 39      FastJets jetpro1(jet_fs, JetAlg::ANTIKT, 0.4, JetMuons::ALL, JetInvisibles::ALL);
 40      declare(jetpro1, "AntiKtJets04");
 41      declare(HeavyHadrons(), "BHadrons");
 42
 43      // Histograms with data binning
 44      book(_h_bjet_Pt      , 3, 1, 1);
 45      book(_h_bjet_Y       , 5, 1, 1);
 46      book(_h_bjet_Yboost  , 7, 1, 1);
 47      book(_h_bjet_DY20    , 9, 1, 1);
 48      book(_h_bjet_ZdPhi20 ,11, 1, 1);
 49      book(_h_bjet_ZdR20   ,13, 1, 1);
 50      book(_h_bjet_ZPt     ,15, 1, 1);
 51      book(_h_bjet_ZY      ,17, 1, 1);
 52      book(_h_2bjet_dR     ,21, 1, 1);
 53      book(_h_2bjet_Mbb    ,23, 1, 1);
 54      book(_h_2bjet_ZPt    ,25, 1, 1);
 55      book(_h_2bjet_ZY     ,27, 1, 1);
 56    }
 57
 58
 59    /// Perform the per-event analysis
 60    void analyze(const Event& e) {
 61
 62      // Check we have a Z:
 63      const DileptonFinder& zfinder = apply<DileptonFinder>(e, "DileptonFinder");
 64      if (zfinder.bosons().size() != 1) vetoEvent;
 65
 66      const Particles boson_s =  zfinder.bosons();
 67      const Particle boson_f =  boson_s[0];
 68      const Particles zleps   =  zfinder.constituents();
 69
 70      // Stop processing the event if no true b-partons or hadrons are found
 71      const Particles allBs = apply<HeavyHadrons>(e, "BHadrons").bHadrons(Cuts::pT > 5.0*GeV);
 72      Particles stableBs = select(allBs, Cuts::abseta < 2.5);
 73      if (stableBs.empty()) vetoEvent;
 74
 75      // Get the b-jets
 76      const Jets& jets = apply<JetFinder>(e, "AntiKtJets04").jetsByPt(Cuts::pT >20.0*GeV && Cuts::abseta <2.4);
 77      Jets b_jets;
 78      for (const Jet& jet : jets) {
 79        //veto overlaps with Z leptons:
 80        bool veto = false;
 81        for (const Particle& zlep : zleps) {
 82          if (deltaR(jet, zlep) < 0.5) veto = true;
 83        }
 84        if (veto) continue;
 85
 86        for (const Particle& bhadron : stableBs) {
 87          if (deltaR(jet, bhadron) <= 0.3) {
 88            b_jets.push_back(jet);
 89            break; // match
 90          }
 91	}
 92      }
 93
 94      // Make sure we have at least 1
 95      if (b_jets.empty()) vetoEvent;
 96
 97      // Fill the plots
 98      const double ZpT = boson_f.pT()/GeV;
 99      const double ZY  = boson_f.absrap();
100
101      _h_bjet_ZPt->fill(ZpT);
102      _h_bjet_ZY ->fill(ZY);
103
104      for (const Jet& jet : b_jets) {
105        _h_bjet_Pt->fill(jet.pT()/GeV);
106        _h_bjet_Y ->fill(jet.absrap());
107
108        const double Yboost = 0.5 * fabs(boson_f.rapidity() + jet.rapidity());
109
110        _h_bjet_Yboost->fill(Yboost);
111
112        if(ZpT > 20.) {
113
114          const double ZBDY   = fabs( boson_f.rapidity() - jet.rapidity() );
115          const double ZBDPHI = fabs( deltaPhi(jet.phi(), boson_f.phi()) );
116          const double ZBDR   = deltaR(jet, boson_f, RAPIDITY);
117          _h_bjet_DY20->fill(   ZBDY);
118          _h_bjet_ZdPhi20->fill(ZBDPHI);
119          _h_bjet_ZdR20->fill(  ZBDR);
120        }
121
122      } //loop over b-jets
123
124      if (b_jets.size() < 2) return;
125
126      _h_2bjet_ZPt->fill(ZpT);
127      _h_2bjet_ZY ->fill(ZY);
128
129      const double BBDR = deltaR(b_jets[0], b_jets[1], RAPIDITY);
130      const double Mbb  = (b_jets[0].momentum() + b_jets[1].momentum()).mass();
131
132      _h_2bjet_dR ->fill(BBDR);
133      _h_2bjet_Mbb->fill(Mbb);
134
135    } // end of analysis loop
136
137
138    /// Normalise histograms etc., after the run
139    void finalize() {
140
141      const double normfac = crossSection()/picobarn / sumOfWeights();
142
143      scale( _h_bjet_Pt,      normfac);
144      scale( _h_bjet_Y,       normfac);
145      scale( _h_bjet_Yboost,  normfac);
146      scale( _h_bjet_DY20,    normfac);
147      scale( _h_bjet_ZdPhi20, normfac);
148      scale( _h_bjet_ZdR20,   normfac);
149      scale( _h_bjet_ZPt,     normfac);
150      scale( _h_bjet_ZY,      normfac);
151      scale( _h_2bjet_dR,     normfac);
152      scale( _h_2bjet_Mbb,    normfac);
153      scale( _h_2bjet_ZPt,    normfac);
154      scale( _h_2bjet_ZY,     normfac);
155    }
156
157    /// @}
158
159
160  protected:
161
162    // Data members like post-cuts event weight counters go here
163    size_t _mode;
164
165
166  private:
167
168    Histo1DPtr _h_bjet_Pt;
169    Histo1DPtr _h_bjet_Y;
170    Histo1DPtr _h_bjet_Yboost;
171    Histo1DPtr _h_bjet_DY20;
172    Histo1DPtr _h_bjet_ZdPhi20;
173    Histo1DPtr _h_bjet_ZdR20;
174    Histo1DPtr _h_bjet_ZPt;
175    Histo1DPtr _h_bjet_ZY;
176    Histo1DPtr _h_2bjet_dR;
177    Histo1DPtr _h_2bjet_Mbb;
178    Histo1DPtr _h_2bjet_ZPt;
179    Histo1DPtr _h_2bjet_ZY;
180
181  };
182
183
184  RIVET_DECLARE_PLUGIN(ATLAS_2014_I1306294);
185
186}