1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/IdentifiedFinalState.hh"
  5#include "Rivet/Projections/VetoedFinalState.hh"
  6#include "Rivet/Projections/FastJets.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief Measurement differential Z/\f$ \gamma^* \f$ + jet + \f$ X \f$ cross sections
 12  ///
 13  /// @author Frank Siegert
 14  class CDF_2008_I768451 : public Analysis {
 15  public:
 16
 17    RIVET_DEFAULT_ANALYSIS_CTOR(CDF_2008_I768451);
 18
 19
 20    /// @name Analysis methods
 21    /// @{
 22
 23    /// Book histograms
 24    void init() {
 25      // Full final state
 26      FinalState fs((Cuts::etaIn(-5.0, 5.0)));
 27      declare(fs, "FS");
 28
 29      // Leading electrons in tracking acceptance
 30      IdentifiedFinalState elfs(Cuts::abseta < 5 && Cuts::pT > 25*GeV);
 31      elfs.acceptIdPair(PID::ELECTRON);
 32      declare(elfs, "LeadingElectrons");
 33
 34      book(_h_jet_multiplicity ,1, 1, 1);
 35      book(_h_jet_pT_cross_section_incl_1jet ,2, 1, 1);
 36      book(_h_jet_pT_cross_section_incl_2jet ,3, 1, 1);
 37    }
 38
 39
 40    /// Do the analysis
 41    void analyze(const Event & event) {
 42      // Skip if the event is empty
 43      const FinalState& fs = apply<FinalState>(event, "FS");
 44      if (fs.empty()) {
 45        MSG_DEBUG("Skipping event " << numEvents() << " because no final state pair found");
 46        vetoEvent;
 47      }
 48
 49      // Find the Z candidates
 50      const FinalState & electronfs = apply<FinalState>(event, "LeadingElectrons");
 51      std::vector<std::pair<Particle, Particle> > Z_candidates;
 52      Particles all_els=electronfs.particles();
 53      for (size_t i=0; i<all_els.size(); ++i) {
 54        for (size_t j=i+1; j<all_els.size(); ++j) {
 55          bool candidate=true;
 56          double mZ = FourMomentum(all_els[i].momentum()+all_els[j].momentum()).mass()/GeV;
 57          if (mZ < 66.0 || mZ > 116.0) {
 58            candidate = false;
 59          }
 60          double abs_eta_0 = fabs(all_els[i].eta());
 61          double abs_eta_1 = fabs(all_els[j].eta());
 62          if (abs_eta_1 < abs_eta_0) {
 63            double tmp = abs_eta_0;
 64            abs_eta_0 = abs_eta_1;
 65            abs_eta_1 = tmp;
 66          }
 67          if (abs_eta_0 > 1.0) {
 68            candidate = false;
 69          }
 70          if (!(abs_eta_1 < 1.0 || (inRange(abs_eta_1, 1.2, 2.8)))) {
 71            candidate = false;
 72          }
 73          if (candidate) {
 74            Z_candidates.push_back(make_pair(all_els[i], all_els[j]));
 75          }
 76        }
 77      }
 78      if (Z_candidates.size() != 1) {
 79        MSG_DEBUG("Skipping event " << numEvents() << " because no unique electron pair found ");
 80        vetoEvent;
 81      }
 82
 83      // Now build the jets on a FS without the electrons from the Z (including QED radiation)
 84      Particles jetparts;
 85      for (const Particle& p : fs.particles()) {
 86        bool copy = true;
 87        if (p.pid() == PID::PHOTON) {
 88          FourMomentum p_e0 = Z_candidates[0].first.momentum();
 89          FourMomentum p_e1 = Z_candidates[0].second.momentum();
 90          FourMomentum p_P = p.momentum();
 91          if (deltaR(p_e0, p_P) < 0.2) copy = false;
 92          if (deltaR(p_e1, p_P) < 0.2) copy = false;
 93        } else {
 94          if (HepMCUtils::uniqueId(p.genParticle()) == HepMCUtils::uniqueId(Z_candidates[0].first.genParticle())) copy = false;
 95          if (HepMCUtils::uniqueId(p.genParticle()) == HepMCUtils::uniqueId(Z_candidates[0].second.genParticle())) copy = false;
 96        }
 97        if (copy) jetparts.push_back(p);
 98      }
 99
100      // Proceed to lepton dressing
101      const PseudoJets pjs = mkPseudoJets(jetparts);
102      const auto jplugin = make_shared<fastjet::CDFMidPointPlugin>(0.7, 0.5, 1.0);
103      const Jets jets_all = mkJets(fastjet::ClusterSequence(pjs, jplugin.get()).inclusive_jets());
104      const Jets jets_cut = sortByPt(select(jets_all, Cuts::pT > 30*GeV && Cuts::abseta < 2.1));
105      // FastJets jetpro(JetAlg::CDFMIDPOINT, 0.7);
106      // jetpro.calc(jetparts);
107      // // Take jets with pt > 30, |eta| < 2.1:
108      // const Jets& jets = jetpro.jets();
109      // Jets jets_cut;
110      // for (const Jet& j, jets) {
111      //   if (j.pT()/GeV > 30.0 && j.abseta() < 2.1) {
112      //     jets_cut.push_back(j);
113      //   }
114      // }
115      // // Sort by pT:
116      // sort(jets_cut.begin(), jets_cut.end(), cmpMomByPt);
117
118      // Return if there are no jets:
119      MSG_DEBUG("Num jets above 30 GeV = " << jets_cut.size());
120      if (jets_cut.empty()) {
121        MSG_DEBUG("No jets pass cuts ");
122        vetoEvent;
123      }
124
125      // Cut on Delta R between Z electrons and *all* jets
126      for (const Jet& j : jets_cut) {
127        if (deltaR(Z_candidates[0].first, j) < 0.7) vetoEvent;
128        if (deltaR(Z_candidates[0].second, j) < 0.7) vetoEvent;
129      }
130
131      // Fill histograms
132      for (size_t njet=1; njet<=jets_cut.size(); ++njet) {
133        _h_jet_multiplicity->fill(njet);
134      }
135      for (const Jet& j : jets_cut) {
136        if (jets_cut.size() > 0) {
137          _h_jet_pT_cross_section_incl_1jet->fill(j.pT());
138        }
139        if (jets_cut.size() > 1) {
140          _h_jet_pT_cross_section_incl_2jet->fill(j.pT());
141        }
142      }
143    }
144
145
146    /// Rescale histos
147    void finalize() {
148      const double invlumi = crossSection()/femtobarn/sumOfWeights();
149      scale(_h_jet_multiplicity, invlumi);
150      scale(_h_jet_pT_cross_section_incl_1jet, invlumi);
151      scale(_h_jet_pT_cross_section_incl_2jet, invlumi);
152    }
153
154    /// @}
155
156
157  private:
158
159    /// @name Histograms
160    /// @{
161    Histo1DPtr _h_jet_multiplicity;
162    Histo1DPtr _h_jet_pT_cross_section_incl_1jet;
163    Histo1DPtr _h_jet_pT_cross_section_incl_2jet;
164    /// @}
165
166  };
167
168
169
170  RIVET_DECLARE_ALIASED_PLUGIN(CDF_2008_I768451, CDF_2008_S7540469);
171
172}