1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ZFinder.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// D0 measurement of the ratio \f$ \sigma(Z/\gamma^* + n \text{ jets})/\sigma(Z/\gamma^*) \f$
 10  class D0_2008_S6879055 : public Analysis {
 11  public:
 12
 13    RIVET_DEFAULT_ANALYSIS_CTOR(D0_2008_S6879055);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    // Book histograms
 20    void init() {
 21      FinalState fs;
 22      ZFinder zfinder(fs, Cuts::open(), PID::ELECTRON,
 23                      40*GeV, 200*GeV, 0.2, ZFinder::ClusterPhotons::NODECAY, ZFinder::AddPhotons::YES);
 24      declare(zfinder, "ZFinder");
 25
 26      FastJets conefinder(zfinder.remainingFinalState(), FastJets::D0ILCONE, 0.5);
 27      declare(conefinder, "ConeFinder");
 28
 29      book(_crossSectionRatio ,1, 1, 1);
 30      book(_pTjet1 ,2, 1, 1);
 31      book(_pTjet2 ,3, 1, 1);
 32      book(_pTjet3 ,4, 1, 1);
 33    }
 34
 35
 36    /// Do the analysis
 37    void analyze(const Event& event) {
 38      const ZFinder& zfinder = apply<ZFinder>(event, "ZFinder");
 39      if (zfinder.bosons().size()!=1) {
 40        vetoEvent;
 41      }
 42
 43      FourMomentum e0 = zfinder.constituents()[0].momentum();
 44      FourMomentum e1 = zfinder.constituents()[1].momentum();
 45      const double e0eta = e0.eta();
 46      const double e0phi = e0.phi();
 47      const double e1eta = e1.eta();
 48      const double e1phi = e1.phi();
 49
 50      vector<FourMomentum> finaljet_list;
 51      for (const Jet& j : apply<JetAlg>(event, "ConeFinder").jetsByPt(20*GeV)) {
 52        const double jeta = j.eta();
 53        const double jphi = j.phi();
 54        if (fabs(jeta) < 2.5) {
 55          if (deltaR(e0eta, e0phi, jeta, jphi) > 0.4 &&
 56              deltaR(e1eta, e1phi, jeta, jphi) > 0.4) {
 57            finaljet_list.push_back(j.momentum());
 58          }
 59        }
 60      }
 61
 62      // For normalisation of crossSection data (includes events with no jets passing cuts)
 63      _crossSectionRatio->fill(0);
 64
 65      // Fill jet pT and multiplicities
 66      if (finaljet_list.size() >= 1) {
 67        _crossSectionRatio->fill(1);
 68        _pTjet1->fill(finaljet_list[0].pT());
 69      }
 70      if (finaljet_list.size() >= 2) {
 71        _crossSectionRatio->fill(2);
 72        _pTjet2->fill(finaljet_list[1].pT());
 73      }
 74      if (finaljet_list.size() >= 3) {
 75        _crossSectionRatio->fill(3);
 76        _pTjet3->fill(finaljet_list[2].pT());
 77      }
 78      if (finaljet_list.size() >= 4) {
 79        _crossSectionRatio->fill(4);
 80      }
 81    }
 82
 83
 84    /// Finalize
 85    void finalize() {
 86      // Now divide by the inclusive result
 87      scale(_crossSectionRatio,1/_crossSectionRatio->bin(0).area());
 88
 89      // Normalise jet pTs to integrals of data
 90      // @note There is no other way to do this, because these quantities are not detector-corrected
 91      /// @todo Use integrals of refData()?
 92      normalize(_pTjet1, 10439); // fixed norm OK
 93      normalize(_pTjet2, 1461.5); // fixed norm OK
 94      normalize(_pTjet3, 217); // fixed norm OK
 95    }
 96
 97    /// @}
 98
 99
100  private:
101
102    /// @name Histograms
103    /// @{
104    Histo1DPtr _crossSectionRatio;
105    Histo1DPtr _pTjet1;
106    Histo1DPtr _pTjet2;
107    Histo1DPtr _pTjet3;
108    /// @}
109
110  };
111
112
113
114  RIVET_DECLARE_ALIASED_PLUGIN(D0_2008_S6879055, D0_2008_I724239);
115
116}