1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/VetoedFinalState.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6
  7#include "fastjet/contrib/SoftDrop.hh"
  8
  9namespace Rivet {
 10
 11
 12
 13  // Soft-drop and ungroomed jet mass measurement
 14  class CMS_2018_I1682495 : public Analysis {
 15  public:
 16
 17    /// @name Constructors etc.
 18    /// @{
 19
 20    /// Constructor
 21    CMS_2018_I1682495()
 22      : Analysis("CMS_2018_I1682495"),
 23        _softdrop(fjcontrib::SoftDrop(0, 0.1, 0.8) ) // parameters are beta, zcut, R0
 24    {    }
 25
 26    /// @}
 27
 28
 29    /// @name Analysis methods
 30    /// @{
 31
 32    /// Book histograms and initialise projections before the run
 33    void init() {
 34      // define a projection that keeps all the particles up to |eta|=5
 35      const FinalState fs(Cuts::abseta < 5.);
 36
 37      // use FastJet, anti-kt(R=0.8) to do the clustering
 38      declare(FastJets(fs, JetAlg::ANTIKT, 0.8), "JetsAK8");
 39
 40      // Histograms
 41      for (size_t i = 0; i < N_PT_BINS_dj; ++i ) {
 42        book(_h_ungroomedJetMass_dj[i][0], i+1+0*N_PT_BINS_dj, 1, 1); // Ungroomed mass, absolute
 43        book(_h_sdJetMass_dj[i][0],        i+1+1*N_PT_BINS_dj, 1, 1); // Groomed mass, absolute
 44        book(_h_ungroomedJetMass_dj[i][1], i+1+2*N_PT_BINS_dj, 1, 1); // Ungroomed mass, normalized
 45        book(_h_sdJetMass_dj[i][1],        i+1+3*N_PT_BINS_dj, 1, 1); // Groomed mass, normalized
 46      }
 47    }
 48
 49
 50    // Find the pT histogram bin index for value pt (in GeV), to hack a 2D histogram equivalent
 51    /// @todo Use a YODA axis/finder alg when available
 52    size_t findPtBin(double ptJ) {
 53      for (size_t ibin = 0; ibin < N_PT_BINS_dj; ++ibin) {
 54        if (inRange(ptJ, ptBins_dj[ibin], ptBins_dj[ibin+1])) return ibin;
 55      }
 56      return N_PT_BINS_dj;
 57    }
 58
 59
 60    /// Perform the per-event analysis
 61    void analyze(const Event& event) {
 62
 63      // Look at events with >= 2 jets
 64      auto jetsAK8 = apply<FastJets>(event, "JetsAK8").jetsByPt(Cuts::pT > 200*GeV and Cuts::abseta < 2.4);
 65      if (jetsAK8.size() < 2) vetoEvent;
 66
 67      // Get the leading two jets
 68      const fastjet::PseudoJet& j0 = jetsAK8[0].pseudojet();
 69      const fastjet::PseudoJet& j1 = jetsAK8[1].pseudojet();
 70
 71      // Calculate delta phi and the pt asymmetry
 72      double deltaPhi = Rivet::deltaPhi( j0.phi(), j1.phi() );
 73      double ptasym = (j0.pt() - j1.pt()) / (j0.pt() + j1.pt());
 74      if (deltaPhi < 2.0 ) vetoEvent;
 75      if (ptasym > 0.3) vetoEvent;
 76
 77      // Find the appropriate pT bins and fill the histogram
 78      const size_t njetBin0 = findPtBin(j0.pt()/GeV);
 79      const size_t njetBin1 = findPtBin(j1.pt()/GeV);
 80      if (njetBin0 < N_PT_BINS_dj && njetBin1 < N_PT_BINS_dj) {
 81        for ( size_t jbin = 0; jbin < N_CATEGORIES; jbin++ ){
 82          _h_ungroomedJetMass_dj[njetBin0][jbin]->fill(j0.m()/GeV);
 83          _h_ungroomedJetMass_dj[njetBin1][jbin]->fill(j1.m()/GeV);
 84        }
 85      }
 86
 87      // Now run the substructure algs...
 88      fastjet::PseudoJet sd0 = _softdrop(j0);
 89      fastjet::PseudoJet sd1 = _softdrop(j1);
 90      // ... and repeat
 91      if (njetBin0 < N_PT_BINS_dj && njetBin1 < N_PT_BINS_dj) {
 92        for ( size_t jbin = 0; jbin < N_CATEGORIES; jbin++ ){
 93          _h_sdJetMass_dj[njetBin0][jbin]->fill(sd0.m()/GeV);
 94          _h_sdJetMass_dj[njetBin1][jbin]->fill(sd1.m()/GeV);
 95        }
 96      }
 97    }
 98
 99
100    /// Normalise histograms etc., after the run
101    void finalize() {
102      // Normalize the normalized cross section histograms to unity,
103      for (size_t i = 0; i < N_PT_BINS_dj; ++i) {
104        normalize(_h_ungroomedJetMass_dj[i][1]);
105        normalize(_h_sdJetMass_dj[i][1]);
106      }
107      // Normalize the absolute cross section histograms to xs * lumi.
108      for (size_t i = 0; i < N_PT_BINS_dj; ++i) {
109        scale(_h_ungroomedJetMass_dj[i][0],   crossSection()/picobarn / sumOfWeights() / (ptBins_dj[i+1]-ptBins_dj[i]) );
110        scale(_h_sdJetMass_dj[i][0],          crossSection()/picobarn / sumOfWeights() / (ptBins_dj[i+1]-ptBins_dj[i]) );
111      }
112    }
113    /// @}
114
115
116  private:
117
118    /// @name FastJet grooming tools (configured in constructor init list)
119    /// @{
120    const fjcontrib::SoftDrop _softdrop;
121    /// @}
122
123
124    /// @name Histograms
125    /// @{
126    enum { PT_200_260_dj=0,
127           PT_260_350_dj,
128           PT_350_460_dj,
129           PT_460_550_dj,
130           PT_550_650_dj,
131           PT_650_760_dj,
132           PT_760_900_dj,
133           PT_900_1000_dj,
134           PT_1000_1100_dj,
135           PT_1100_1200_dj,
136           PT_1200_1300_dj,
137           PT_1300_Inf_dj,
138           N_PT_BINS_dj };
139    static const int N_CATEGORIES=2;
140    const double ptBins_dj[N_PT_BINS_dj+1]= { 200., 260., 350., 460., 550., 650., 760., 900., 1000., 1100., 1200., 1300., 13000.};
141
142    Histo1DPtr _h_ungroomedJet0pt, _h_ungroomedJet1pt;
143    Histo1DPtr _h_sdJet0pt, _h_sdJet1pt;
144    // Here, store both the absolute (index 0) and normalized (index 1) cross sections.
145    Histo1DPtr _h_ungroomedJetMass_dj[N_PT_BINS_dj][2];
146    Histo1DPtr _h_sdJetMass_dj[N_PT_BINS_dj][2];
147    /// @}
148
149
150  };
151
152
153  RIVET_DECLARE_PLUGIN(CMS_2018_I1682495);
154
155
156}