1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief  Measurement of inclusive and Mueller-Navelet dijet cross sections and their ratios at 2.76 TeV
 10  class CMS_2021_I1963239 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2021_I1963239);
 15
 16
 17    /// @name Analysis methods
 18    ///@{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22
 23      // Initialise and register projections
 24
 25      // The basic final-state projection:
 26      // all final-state particles within
 27      // the given eta acceptance
 28      const FinalState fs(Cuts::abseta < 5.2);
 29
 30      // The final-state particles declared above are clustered using FastJet with
 31      // the anti-kT algorithm and a jet-radius parameter 0.5
 32      FastJets jetfs(fs, JetAlg::ANTIKT, 0.5);
 33      declare(jetfs, "jets");
 34
 35      // Book histograms
 36      // specify custom binning
 37      // take binning from reference data using HEPData ID (digits in "d01-x01-y01" etc.)
 38
 39      book(_h["inclusive"], 7, 1, 1);
 40      book(_h["MN"], 8, 1, 1);
 41      book(_s["R_incl"], 9, 1, 1);
 42      book(_s["R_incl_veto"], 10, 1, 1);
 43      book(_s["R_MN"], 11, 1, 1);
 44      book(_s["R_MN_veto"], 12, 1, 1);
 45
 46      // Temporary histograms (directly instantiated)
 47      book(_h["exclusive"], "_exclusive", refData(7, 1, 1));
 48      book(_h["exclusive_veto"], "_exclusive_veto", refData(7, 1, 1));
 49    }
 50
 51
 52    /// Perform the per-event analysis
 53    void analyze(const Event& event) {
 54
 55      // Retrieve clustered jets, sorted by pT, with a minimum pT cut
 56      Jets jets20 = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 20*GeV && Cuts::absrap < 4.7);
 57      Jets jets35 = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 35*GeV && Cuts::absrap < 4.7);
 58
 59      if (jets35.size() < 2) return;
 60
 61      // Loop over jet pairs
 62      double deltaY_MN = 0.0;
 63      for (size_t ij1 = 0; ij1 < jets35.size(); ++ij1) {
 64        for (size_t ij2 = ij1 + 1; ij2 < jets35.size(); ++ij2) {
 65          const double deltaY = fabs(jets35[ij1].rapidity() - jets35[ij2].rapidity());
 66          // Exclusive dijet case:
 67          if (jets35.size() == 2) {
 68            _h["exclusive"]->fill(deltaY);
 69            //Exclusive with veto 20 GeV dijet case:
 70            if (jets20.size() == 2) {
 71              _h["exclusive_veto"]->fill(deltaY);
 72            }
 73          }
 74          // Inclusive jets case:
 75          _h["inclusive"]->fill(deltaY);
 76          // Mueller-Navelet:
 77          if (deltaY > deltaY_MN) deltaY_MN = deltaY;
 78        }
 79      }
 80      // Fill histogram with MN dijets Delta y
 81      _h["MN"]->fill(deltaY_MN);
 82
 83    }
 84
 85
 86    /// Normalise histograms etc., after the run
 87    void finalize() {
 88
 89      // Calculate ratios
 90      efficiency(_h["exclusive"], _h["inclusive"], _s["R_incl"]);
 91      efficiency(_h["exclusive"], _h["MN"], _s["R_MN"]);
 92      efficiency(_h["exclusive_veto"], _h["inclusive"], _s["R_incl_veto"]);
 93      efficiency(_h["exclusive_veto"], _h["MN"], _s["R_MN_veto"]);
 94
 95      transform(*_s["R_incl"], _invert);
 96      transform(*_s["R_MN"], _invert);
 97      transform(*_s["R_incl_veto"], _invert);
 98      transform(*_s["R_MN_veto"], _invert);
 99
100
101      scale(_h["inclusive"], crossSection()/picobarn/sumOfWeights()); // norm to generated cross-section in pb
102      scale(_h["MN"], crossSection()/picobarn/sumOfWeights());        // norm to generated cross-section in pb
103
104    }
105
106    ///@}
107
108
109    /// @name Histograms
110    ///@{
111    map<string, Histo1DPtr> _h;
112    map<string, Estimate1DPtr>_s;
113    ///@}
114    private:
115
116    /// Reciprocal function with div-by-zero protection, for inverting the efficiency measure
117    static double _invert(double x) { return (x > 0) ? 1/x : 0; }
118
119
120  };
121
122
123  RIVET_DECLARE_PLUGIN(CMS_2021_I1963239);
124
125}