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, FastJets::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      const double weight = 1.0;
 55
 56      // Retrieve clustered jets, sorted by pT, with a minimum pT cut
 57      Jets jets20 = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 20*GeV && Cuts::absrap < 4.7);
 58      Jets jets35 = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 35*GeV && Cuts::absrap < 4.7);
 59
 60      if (jets35.size() < 2) return;
 61
 62      // Loop over jet pairs
 63      double deltaY_MN = 0.0;
 64      for (size_t ij1 = 0; ij1 < jets35.size(); ++ij1) {
 65        for (size_t ij2 = ij1 + 1; ij2 < jets35.size(); ++ij2) {
 66          const double deltaY = fabs(jets35[ij1].rapidity() - jets35[ij2].rapidity());
 67          // Exclusive dijet case:
 68          if (jets35.size() == 2) {
 69            _h["exclusive"]->fill(deltaY, weight);
 70            //Exclusive with veto 20 GeV dijet case:
 71            if (jets20.size() == 2) {
 72              _h["exclusive_veto"]->fill(deltaY, weight);
 73            }
 74          }
 75          // Inclusive jets case:
 76          _h["inclusive"]->fill(deltaY, weight);
 77          // Mueller-Navelet:
 78          if (deltaY > deltaY_MN) deltaY_MN = deltaY;
 79        }
 80      }
 81      // Fill histogram with MN dijets Delta y
 82      _h["MN"]->fill(deltaY_MN, weight);
 83
 84    }
 85
 86
 87    /// Normalise histograms etc., after the run
 88    void finalize() {
 89
 90      // Calculate ratios
 91      efficiency(_h["exclusive"], _h["inclusive"], _s["R_incl"]);
 92      efficiency(_h["exclusive"], _h["MN"], _s["R_MN"]);
 93      efficiency(_h["exclusive_veto"], _h["inclusive"], _s["R_incl_veto"]);
 94      efficiency(_h["exclusive_veto"], _h["MN"], _s["R_MN_veto"]);
 95
 96      transformY(*_s["R_incl"], _invert);
 97      transformY(*_s["R_MN"], _invert);
 98      transformY(*_s["R_incl_veto"], _invert);
 99      transformY(*_s["R_MN_veto"], _invert);
100
101
102      scale(_h["inclusive"], crossSection()/picobarn/sumOfWeights()); // norm to generated cross-section in pb
103      scale(_h["MN"], crossSection()/picobarn/sumOfWeights());        // norm to generated cross-section in pb
104
105    }
106
107    ///@}
108
109
110    /// @name Histograms
111    ///@{
112    map<string, Histo1DPtr> _h;
113    map<string, Scatter2DPtr>_s;
114    ///@}
115    private:
116    
117    /// Reciprocal function with div-by-zero protection, for inverting the efficiency measure
118    static double _invert(double x) { return (x > 0) ? 1/x : 0; }
119
120
121  };
122
123
124  RIVET_DECLARE_PLUGIN(CMS_2021_I1963239);
125
126}