1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/VetoedFinalState.hh"
  6#include "Rivet/Projections/InvMassFinalState.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief CDF Run II Z + b-jet cross-section measurement
 12  class CDF_2008_I806082 : public Analysis {
 13  public:
 14
 15    RIVET_DEFAULT_ANALYSIS_CTOR(CDF_2008_I806082);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    void init() {
 22      // Set up projections
 23      const FinalState fs((Cuts::etaIn(-3.2, 3.2)));
 24      declare(fs, "FS");
 25      // Create a final state with any e+e- or mu+mu- pair with
 26      // invariant mass 76 -> 106 GeV and ET > 18 (Z decay products)
 27      vector<pair<PdgId,PdgId> > vids;
 28      vids.push_back(make_pair(PID::ELECTRON, PID::POSITRON));
 29      vids.push_back(make_pair(PID::MUON, PID::ANTIMUON));
 30      FinalState fs2((Cuts::etaIn(-3.2, 3.2)));
 31      InvMassFinalState invfs(fs2, vids, 76*GeV, 106*GeV);
 32      declare(invfs, "INVFS");
 33      // Make a final state without the Z decay products for jet clustering
 34      VetoedFinalState vfs(fs);
 35      vfs.addVetoOnThisFinalState(invfs);
 36      declare(vfs, "VFS");
 37      declare(FastJets(vfs, JetAlg::CDFMIDPOINT, 0.7), "Jets");
 38
 39      // Book histograms
 40      book(_dStot    ,1, 1, 1);
 41      book(_dSdET    ,2, 1, 1);
 42      book(_dSdETA   ,3, 1, 1);
 43      book(_dSdZpT   ,4, 1, 1);
 44      book(_dSdNJet  ,5, 1, 1);
 45      book(_dSdNbJet ,6, 1, 1);
 46
 47      book(_sumWeightSelected,"sumWeightSelected");
 48    }
 49
 50
 51    // Do the analysis
 52    void analyze(const Event& event) {
 53      // Check we have an l+l- pair that passes the kinematic cuts
 54      // Get the Z decay products (mu+mu- or e+e- pair)
 55      const InvMassFinalState& invMassFinalState = apply<InvMassFinalState>(event, "INVFS");
 56      const Particles&  ZDecayProducts =  invMassFinalState.particles();
 57
 58      // make sure we have 2 Z decay products (mumu or ee)
 59      if (ZDecayProducts.size() < 2) vetoEvent;
 60      //new cuts
 61      double Lep1Pt = ZDecayProducts[0].perp();
 62      double Lep2Pt = ZDecayProducts[1].perp();
 63      double Lep1Eta = fabs(ZDecayProducts[0].rapidity());
 64      double Lep2Eta = fabs(ZDecayProducts[1].rapidity());
 65
 66      if (Lep1Eta > _LepEtaCut || Lep2Eta > _LepEtaCut) vetoEvent;
 67
 68      if (ZDecayProducts[0].abspid()==13 &&
 69          ((Lep1Eta > 1.5 || Lep2Eta > 1.5) || (Lep1Eta > 1.0 && Lep2Eta > 1.0))) {
 70        vetoEvent;
 71      }
 72
 73      if (Lep1Pt > Lep2Pt) {
 74        if (Lep1Pt < _Lep1PtCut || Lep2Pt < _Lep2PtCut) vetoEvent;
 75      }
 76      else {
 77        if (Lep1Pt < _Lep2PtCut || Lep2Pt < _Lep1PtCut) vetoEvent;
 78      }
 79
 80      _sumWeightSelected->fill();
 81      /// @todo: write out a warning if there are more than two decay products
 82      FourMomentum Zmom = ZDecayProducts[0].momentum() +  ZDecayProducts[1].momentum();
 83
 84      // Put all b-quarks in a vector
 85      /// @todo Use a b-hadron search rather than b-quarks for tagging
 86      Particles bquarks;
 87      for(ConstGenParticlePtr p: HepMCUtils::particles(event.genEvent())) {
 88        if (std::abs(p->pdg_id()) == PID::BQUARK) {
 89          bquarks += Particle(*p);
 90        }
 91      }
 92
 93      // Get jets
 94      const FastJets& jetpro = apply<FastJets>(event, "Jets");
 95      MSG_DEBUG("Jet multiplicity before any pT cut = " << jetpro.size());
 96
 97      const PseudoJets& jets = jetpro.pseudojetsByPt();
 98      MSG_DEBUG("jetlist size = " << jets.size());
 99
100      int numBJet = 0;
101      int numJet  = 0;
102      // for each b-jet plot the ET and the eta of the jet, normalise to the total cross section at the end
103      // for each event plot N jet and pT(Z), normalise to the total cross section at the end
104      for (PseudoJets::const_iterator jt = jets.begin(); jt != jets.end(); ++jt) {
105        // select jets that pass the kinematic cuts
106        if (jt->perp() > _JetPtCut && fabs(jt->rapidity()) <= _JetEtaCut) {
107          numJet++;
108          // does the jet contain a b-quark?
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            _dSdET->fill(jt->perp());
119            _dSdETA->fill(fabs(jt->rapidity()));
120          }
121        }
122      } // end loop around jets
123
124      // wasn't asking for b-jets before!!!!
125      if(numJet > 0 && numBJet > 0) _dSdNJet->fill(numJet);
126      if(numBJet > 0) {
127        _dStot->fill(1960);
128        _dSdNbJet->fill(numBJet);
129        _dSdZpT->fill(Zmom.pT());
130      }
131    }
132
133
134    // Finalize
135    void finalize() {
136      // normalise histograms
137      // scale by 1 / the sum-of-weights of events that pass the Z cuts
138      // since the cross sections are normalized to the inclusive
139      // Z cross sections.
140      double Scale = 1.0;
141      if (_sumWeightSelected->val() != 0.0) Scale = 1.0/dbl(*_sumWeightSelected);
142      scale(_dStot,Scale);
143      scale(_dSdET,Scale);
144      scale(_dSdETA,Scale);
145      scale(_dSdNJet,Scale);
146      scale(_dSdNbJet,Scale);
147      scale(_dSdZpT,Scale);
148    }
149
150    /// @}
151
152
153  private:
154
155    /// @name Cuts
156    /// @{
157    const double _Rjet = 0.7;
158    const double _JetPtCut = 20;
159    const double _JetEtaCut = 1.5;
160    const double _Lep1PtCut = 18;
161    const double _Lep2PtCut = 10;
162    const double _LepEtaCut = 3.2;
163    /// @}
164
165    /// Counter
166    CounterPtr _sumWeightSelected;
167
168    /// @name Histograms
169    /// @{
170    BinnedHistoPtr<int> _dStot, _dSdNbJet, _dSdNJet;
171    Histo1DPtr _dSdET, _dSdETA, _dSdZpT;
172    /// @}
173
174  };
175
176
177
178  RIVET_DECLARE_ALIASED_PLUGIN(CDF_2008_I806082, CDF_2008_S8095620);
179
180}