1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/LeptonFinder.hh"
  6#include "Rivet/Projections/MissingMomentum.hh"
  7#include "Rivet/Projections/PromptFinalState.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief Production cross-sections of WZ and same-sign WW with two jets in pp collisions at 13 TeV
 13  class CMS_2020_I1794169 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2020_I1794169);
 18
 19
 20    /// @name Analysis methods
 21    /// @{
 22
 23    /// Book histograms and initialise projections before the run
 24    void init() {
 25
 26      // Initialise and register projections
 27      _mode = 0;
 28      if ( getOption("LMODE") == "WZ" ) _mode = 1;
 29
 30      // The basic final-state projection:
 31      // all final-state particles within
 32      // the given eta acceptance
 33      const FinalState fs(Cuts::abseta < 4.9);
 34      const FinalState fsjet4p7(Cuts::abseta < 4.7);
 35
 36      // The final-state particles declared above are clustered using FastJet with
 37      // the anti-kT algorithm and a jet-radius parameter 0.4
 38      // muons and neutrinos are excluded from the clustering
 39      FastJets jet4p7fs(fsjet4p7, JetAlg::ANTIKT, 0.4);
 40      declare(jet4p7fs, "jets4p7");
 41
 42      // FinalState of prompt photons and bare muons and electrons in the event
 43      PromptFinalState photons(Cuts::abspid == PID::PHOTON);
 44      PromptFinalState bare_leps(Cuts::abspid == PID::MUON || Cuts::abspid == PID::ELECTRON);
 45      bare_leps.acceptTauDecays(false);
 46
 47      // Dress the prompt bare leptons with prompt photons within dR < 0.1,
 48      // and apply some fiducial cuts on the dressed leptons
 49      Cut lepton_cuts = Cuts::abseta < 2.5 && Cuts::pT > 20*GeV;
 50      LeptonFinder dressed_leps(bare_leps, photons, 0.1, lepton_cuts);
 51      declare(dressed_leps, "leptons");
 52
 53      // Missing momentum
 54      declare(MissingMomentum(fs), "MET");
 55
 56      // Book histograms
 57      book(_h_WW_mjj   ,  9, 1, 1);
 58      book(_h_WW_mll   , 11, 1, 1);
 59      book(_h_WW_ptlmax, 13, 1, 1);
 60      book(_h_WZ_mjj   , 15, 1, 1);
 61
 62    }
 63
 64    /// Perform the per-event analysis
 65    void analyze(const Event& event) {
 66
 67      // Retrieve dressed leptons, sorted by pT
 68      Particles leptons = apply<LeptonFinder>(event, "leptons").particles();
 69
 70      // Apply a #leptons requirement
 71      if (leptons.size() <= 1 || leptons.size() >= 4) return;
 72
 73      // Retrieve clustered jets, sorted by pT, with a minimum pT cut
 74      Jets jets50 = apply<FastJets>(event, "jets4p7").jetsByPt(Cuts::pT > 50*GeV);
 75
 76      // Remove all jets within dR < 0.4 of a dressed lepton
 77      idiscardIfAnyDeltaRLess(jets50, leptons, 0.4);
 78
 79      // Apply a njets >= 2 cut
 80      if (jets50.size() < 2) return;
 81
 82      FourMomentum dijetCand = jets50[0].momentum() + jets50[1].momentum();
 83      double deltaEtaJJ = std::abs(jets50[0].eta() - jets50[1].eta());
 84
 85      // Apply a mjj > 500 and detajj > 2.5 cuts
 86      if (dijetCand.mass() <= 500*GeV || deltaEtaJJ <= 2.5) return;
 87
 88      // W+W+ selection
 89      if (leptons.size() == 2 && leptons[0].pid() * leptons[1].pid() > 0 && _mode == 0) {
 90        FourMomentum dilCand = leptons[0].momentum() + leptons[1].momentum();
 91        if (dilCand.mass() > 20*GeV) {
 92          double ptlmax = leptons[0].pt(); double ptlmin = leptons[1].pt();
 93          if (ptlmax < ptlmin) {
 94            ptlmax = leptons[1].pt(); ptlmin = leptons[0].pt();
 95          }
 96
 97          _h_WW_mjj   ->fill(min(dijetCand.mass()/GeV, 2999.999));
 98          _h_WW_mll   ->fill(min(dilCand.mass()/GeV, 499.999));
 99          _h_WW_ptlmax->fill(min(ptlmax/GeV, 299.999));
100
101        }
102      }
103
104      // WZ selection
105      else if (leptons.size() == 3 && _mode == 1) {
106        double mllZ = 10000; int iW = -1;
107        if (leptons[0].pid() * leptons[1].pid() < 0 && leptons[0].abspid() == leptons[1].abspid() &&
108            fabs((leptons[0].momentum() + leptons[1].momentum()).mass() - 91.1876*GeV) < fabs(mllZ - 91.1876*GeV)) {
109          mllZ = (leptons[0].momentum() + leptons[1].momentum()).mass(); iW = 2;
110        }
111
112        if (leptons[0].pid() * leptons[2].pid() < 0 && leptons[0].abspid() == leptons[2].abspid() &&
113            fabs((leptons[0].momentum() + leptons[2].momentum()).mass() - 91.1876*GeV) < fabs(mllZ - 91.1876*GeV)) {
114          mllZ = (leptons[0].momentum() + leptons[2].momentum()).mass(); iW = 1;
115        }
116
117        if (leptons[1].pid() * leptons[2].pid() < 0 && leptons[1].abspid() == leptons[2].abspid() &&
118            fabs((leptons[1].momentum() + leptons[2].momentum()).mass() - 91.1876*GeV) < fabs(mllZ - 91.1876*GeV)) {
119          mllZ = (leptons[1].momentum() + leptons[2].momentum()).mass(); iW = 0;
120        }
121
122        // Plot
123        if (iW >= 0 && fabs(mllZ - 91.1876*GeV) < 15*GeV) {
124          _h_WZ_mjj->fill(min(dijetCand.mass()/GeV, 2999.999));
125        }
126      }
127
128    }
129
130
131    /// @todo Replace with barchart()
132    void normalizeToSum(Histo1DPtr hist) {
133      for (size_t i = 0; i < hist->numBins(); ++i) {
134        float width = hist->bin(i).xWidth();
135        hist->bin(i).scaleW(width != 0 ? width : 1.);
136      }
137      if (hist->integral() > 0) scale(hist, 1./hist->integral());
138    }
139
140
141    /// Normalise histograms etc., after the run
142    void finalize() {
143      double norm = (sumOfWeights() != 0) ? crossSection()/femtobarn/sumOfWeights() : 1.0;
144      scale(_h_WW_mjj   , norm);
145      scale(_h_WW_mll   , norm);
146      scale(_h_WW_ptlmax, norm);
147      scale(_h_WZ_mjj   , norm);
148    }
149
150
151  private:
152
153    /// Lepton-mode flag
154    size_t _mode;
155
156    /// @name Histograms
157    /// @{
158    Histo1DPtr _h_WW_mjj, _h_WW_mll, _h_WW_ptlmax, _h_WZ_mjj;
159    /// @}
160
161  };
162
163
164  RIVET_DECLARE_PLUGIN(CMS_2020_I1794169);
165
166}