1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/VetoedFinalState.hh"
  5#include "Rivet/Projections/ChargedFinalState.hh"
  6#include "Rivet/Projections/InvMassFinalState.hh"
  7#include "Rivet/Projections/FastJets.hh"
  8#include "Rivet/Projections/ChargedLeptons.hh"
  9
 10namespace Rivet {
 11
 12
 13  /// @brief CDF Run II analysis: jet \f$ p_T \f$ and \f$ \eta \f$
 14  ///   distributions in Z + (b) jet production
 15  ///
 16  /// @author Lars Sonnenschein
 17  ///
 18  /// This CDF analysis provides \f$ p_T \f$ and \f$ \eta \f$ distributions of
 19  /// jets in Z + (b) jet production, before and after tagging.
 20  class CDF_2006_I717572 : public Analysis {
 21  public:
 22
 23    /// Constructor
 24    RIVET_DEFAULT_ANALYSIS_CTOR(CDF_2006_I717572);
 25
 26
 27    /// @name Analysis methods
 28    /// @{
 29
 30    void init() {
 31      const FinalState fs(Cuts::abseta < 3.6);
 32      declare(fs, "FS");
 33
 34      // Create a final state with any e+e- or mu+mu- pair with
 35      // invariant mass 76 -> 106 GeV and ET > 20 (Z decay products)
 36      vector<pair<PdgId,PdgId> > vids;
 37      vids.push_back(make_pair(PID::ELECTRON, PID::POSITRON));
 38      vids.push_back(make_pair(PID::MUON, PID::ANTIMUON));
 39      FinalState fs2((Cuts::etaIn(-3.6, 3.6)));
 40      InvMassFinalState invfs(fs2, vids, 66*GeV, 116*GeV);
 41      declare(invfs, "INVFS");
 42
 43      // Make a final state without the Z decay products for jet clustering
 44      VetoedFinalState vfs(fs);
 45      vfs.addVetoOnThisFinalState(invfs);
 46      declare(vfs, "VFS");
 47      declare(FastJets(vfs, JetAlg::CDFMIDPOINT, 0.7), "Jets");
 48
 49      // Book histograms
 50      book(_sigmaBJet ,1, 1, 1);
 51      book(_ratioBJetToZ ,2, 1, 1);
 52      book(_ratioBJetToJet ,3, 1, 1);
 53      book(_sumWeightsWithZ, "sumWeightsWithZ");
 54      book(_sumWeightsWithZJet, "sumWeightsWithZJet");
 55    }
 56
 57
 58    /// Do the analysis
 59    void analyze(const Event& event) {
 60      // Check we have an l+l- pair that passes the kinematic cuts
 61      // Get the Z decay products (mu+mu- or e+e- pair)
 62      const InvMassFinalState& invMassFinalState = apply<InvMassFinalState>(event, "INVFS");
 63      const Particles&  ZDecayProducts =  invMassFinalState.particles();
 64
 65      // Make sure we have at least 2 Z decay products (mumu or ee)
 66      if (ZDecayProducts.size() < 2) vetoEvent;
 67      //
 68      double Lep1Pt = ZDecayProducts[0].pT();
 69      double Lep2Pt = ZDecayProducts[1].pT();
 70      double Lep1Eta = ZDecayProducts[0].absrap(); ///< @todo This is y... should be abseta()?
 71      double Lep2Eta = ZDecayProducts[1].absrap(); ///< @todo This is y... should be abseta()?
 72
 73      if (Lep1Eta > _LepEtaCut && Lep2Eta > _LepEtaCut) vetoEvent;
 74      if (ZDecayProducts[0].abspid()==13 && Lep1Eta > 1. && Lep2Eta > 1.) vetoEvent;
 75      if (Lep1Pt < _Lep1PtCut && Lep2Pt < _Lep2PtCut) vetoEvent;
 76
 77      _sumWeightsWithZ->fill();
 78
 79      /// @todo Write out a warning if there are more than two decay products
 80      FourMomentum Zmom = ZDecayProducts[0].momentum() +  ZDecayProducts[1].momentum();
 81
 82      // Put all b-quarks in a vector
 83      /// @todo Use jet contents rather than accessing quarks directly
 84      Particles bquarks;
 85      /// @todo is this HepMC wrangling necessary?
 86      for (ConstGenParticlePtr p: HepMCUtils::particles(event.genEvent())) {
 87        if ( std::abs(p->pdg_id()) == PID::BQUARK ) {
 88          bquarks.push_back(Particle(*p));
 89        }
 90      }
 91
 92      // Get jets
 93      const FastJets& jetpro = apply<FastJets>(event, "Jets");
 94      MSG_DEBUG("Jet multiplicity before any pT cut = " << jetpro.size());
 95
 96      const PseudoJets& jets = jetpro.pseudojetsByPt();
 97      MSG_DEBUG("jetlist size = " << jets.size());
 98
 99      int numBJet = 0;
100      int numJet  = 0;
101      // for each b-jet plot the ET and the eta of the jet, normalise to the total cross section at the end
102      // for each event plot N jet and pT(Z), normalise to the total cross section at the end
103      for (PseudoJets::const_iterator jt = jets.begin(); jt != jets.end(); ++jt) {
104        // select jets that pass the kinematic cuts
105        if (jt->perp() > _JetPtCut && fabs(jt->rapidity()) <= _JetEtaCut) {
106          ++numJet;
107          // Does the jet contain a b-quark?
108          /// @todo Use jet contents rather than accessing quarks directly
109          bool bjet = false;
110          for (const Particle& bquark : bquarks) {
111            if (deltaR(jt->rapidity(), jt->phi(), bquark.rapidity(), bquark.phi()) <= _Rjet) {
112              bjet = true;
113              break;
114            }
115          } // end loop around b-jets
116          if (bjet) {
117            numBJet++;
118          }
119        }
120      } // end loop around jets
121
122      if (numJet > 0) _sumWeightsWithZJet->fill();
123      if (numBJet > 0) {
124        _sigmaBJet->fill(1960);
125        _ratioBJetToZ->fill(1960);
126        _ratioBJetToJet->fill(1960);
127      }
128
129    }
130
131
132    /// Finalize
133    void finalize() {
134      MSG_DEBUG("Total sum of weights = " << sumOfWeights());
135      MSG_DEBUG("Sum of weights for Z production in mass range = " << dbl(*_sumWeightsWithZ));
136      MSG_DEBUG("Sum of weights for Z+jet production in mass range = " << dbl(*_sumWeightsWithZJet));
137
138      scale(_sigmaBJet, crossSection()/picobarn/sumOfWeights());
139      scale(_ratioBJetToZ, 1.0/ *_sumWeightsWithZ);
140      scale(_ratioBJetToJet, 1.0/ *_sumWeightsWithZJet);
141    }
142
143    /// @}
144
145
146  private:
147
148    /// @name Cuts and counters
149    /// @{
150    const double _Rjet = 0.7;
151    const double _JetPtCut = 20.;
152    const double _JetEtaCut = 1.5;
153    const double _Lep1PtCut = 18.;
154    const double _Lep2PtCut = 10.;
155    const double _LepEtaCut = 1.1;
156    /// @}
157
158    /// @name Cuts and counters
159    /// @{
160    CounterPtr _sumWeightsWithZ;
161    CounterPtr _sumWeightsWithZJet;
162    /// @}
163
164    /// @name Histograms
165    /// @{
166    BinnedHistoPtr<int> _sigmaBJet, _ratioBJetToZ, _ratioBJetToJet;
167    /// @}
168
169  };
170
171
172
173  RIVET_DECLARE_ALIASED_PLUGIN(CDF_2006_I717572, CDF_2006_S6653332);
174
175}