1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/DressedLeptons.hh"
  6#include "Rivet/Projections/MissingMomentum.hh"
  7#include "Rivet/Projections/PromptFinalState.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief Measurements of W+W- boson pair production in pp collisions at 13 TeV
 13  class CMS_2020_I1814328 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2020_I1814328);
 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
 28      // The basic final-state projection:
 29      // all final-state particles within
 30      // the given eta acceptance
 31      const FinalState fs(Cuts::abseta < 4.9);
 32      const FinalState fsjet4p5(Cuts::abseta < 4.5);
 33      const FinalState fsjet2p5(Cuts::abseta < 2.5);
 34
 35      // The final-state particles declared above are clustered using FastJet with
 36      // the anti-kT algorithm and a jet-radius parameter 0.4
 37      // muons and neutrinos are excluded from the clustering
 38      FastJets jet4p5fs(fsjet4p5, FastJets::ANTIKT, 0.4);
 39      declare(jet4p5fs, "jets4p5");
 40      FastJets jet2p5fs(fsjet2p5, FastJets::ANTIKT, 0.4);
 41      declare(jet2p5fs, "jets2p5");
 42
 43      // FinalState of prompt photons and bare muons and electrons in the event
 44      PromptFinalState photons(Cuts::abspid == PID::PHOTON);
 45      PromptFinalState bare_leps(Cuts::abspid == PID::MUON || Cuts::abspid == PID::ELECTRON);
 46      bare_leps.acceptTauDecays(false);
 47
 48      // Dress the prompt bare leptons with prompt photons within dR < 0.1,
 49      // and apply some fiducial cuts on the dressed leptons
 50      Cut lepton_cuts = Cuts::abseta < 2.5 && Cuts::pT > 25*GeV;
 51      DressedLeptons dressed_leps(photons, bare_leps, 0.1, lepton_cuts);
 52      declare(dressed_leps, "leptons");
 53
 54      // Missing momentum
 55      declare(MissingMomentum(fs), "MET");
 56
 57      // Book histograms
 58      book(_h_WW_njets_norm , 2, 1, 1);
 59      book(_h_WW_mll_norm   , 4, 1, 1);
 60      book(_h_WW_ptlmax_norm, 5, 1, 1);
 61      book(_h_WW_ptlmin_norm, 6, 1, 1);
 62      book(_h_WW_dphill_norm, 7, 1, 1);
 63      book(_h_WW_njet0      , 8, 1, 1);
 64
 65    }
 66
 67    /// Perform the per-event analysis
 68    void analyze(const Event& event) {
 69
 70      // Apply a missing-momentum cut
 71      if (apply<MissingMomentum>(event, "MET").missingPt() < 20*GeV) return;
 72
 73      // Retrieve dressed leptons, sorted by pT
 74      vector<DressedLepton> leptons = apply<DressedLeptons>(event, "leptons").dressedLeptons();
 75
 76      // Retrieve clustered jets, sorted by pT, with a minimum pT cut
 77      Jets jets25 = apply<FastJets>(event, "jets4p5").jetsByPt(Cuts::pT > 25*GeV);
 78      Jets jets30 = apply<FastJets>(event, "jets4p5").jetsByPt(Cuts::pT > 30*GeV);
 79      Jets jets35 = apply<FastJets>(event, "jets4p5").jetsByPt(Cuts::pT > 35*GeV);
 80      Jets jets45 = apply<FastJets>(event, "jets4p5").jetsByPt(Cuts::pT > 45*GeV);
 81      Jets jets60 = apply<FastJets>(event, "jets4p5").jetsByPt(Cuts::pT > 60*GeV);
 82      Jets jetsNj = apply<FastJets>(event, "jets2p5").jetsByPt(Cuts::pT > 30*GeV);
 83
 84      // Remove all jets within dR < 0.4 of a dressed lepton
 85      idiscardIfAnyDeltaRLess(jets25, leptons, 0.4);
 86      idiscardIfAnyDeltaRLess(jets30, leptons, 0.4);
 87      idiscardIfAnyDeltaRLess(jets35, leptons, 0.4);
 88      idiscardIfAnyDeltaRLess(jets45, leptons, 0.4);
 89      idiscardIfAnyDeltaRLess(jets60, leptons, 0.4);
 90      idiscardIfAnyDeltaRLess(jetsNj, leptons, 0.4);
 91
 92      if (leptons.size() == 2 && leptons[0].pid() * leptons[1].pid() < 0) {
 93        FourMomentum dilCand = leptons[0].momentum() + leptons[1].momentum();
 94        if (dilCand.mass() > 20*GeV && dilCand.pT() > 30*GeV){
 95          double ptlmax = leptons[0].pT(); double ptlmin = leptons[1].pT();
 96          if (ptlmax < ptlmin) {
 97            ptlmax = leptons[1].pT(); ptlmin = leptons[0].pT();
 98          }
 99
100          if (leptons[0].abspid() != leptons[1].abspid()) {
101            _h_WW_njets_norm ->fill(min((double)jetsNj.size()+1, 2.999));
102            _h_WW_mll_norm   ->fill(min(dilCand.mass()/GeV, 1499.999));
103            _h_WW_ptlmax_norm->fill(min(ptlmax/GeV, 399.999));
104            _h_WW_ptlmin_norm->fill(min(ptlmin/GeV, 149.999));
105            _h_WW_dphill_norm->fill(deltaPhi(leptons[0], leptons[1]));
106          }
107
108          if (jets25.size() == 0) _h_WW_njet0->fill(1.0);
109          if (jets30.size() == 0) _h_WW_njet0->fill(2.0);
110          if (jets35.size() == 0) _h_WW_njet0->fill(3.0);
111          if (jets45.size() == 0) _h_WW_njet0->fill(4.0);
112          if (jets60.size() == 0) _h_WW_njet0->fill(5.0);
113
114        }
115      }
116
117    }
118
119
120    /// @todo Replace with barchart()
121    void normalizeToSum(Histo1DPtr hist) {
122      double sum = 0.;
123      for (size_t i = 0; i < hist->numBins(); ++i) {
124        sum += hist->bin(i).height();
125      }
126      scale(hist, 1./sum);
127    }
128
129
130    /// Normalise histograms etc., after the run
131    void finalize() {
132
133      double norm = (sumOfWeights() != 0) ? crossSection()/picobarn/sumOfWeights() : 1.0;
134
135      normalizeToSum(_h_WW_njets_norm );
136      normalizeToSum(_h_WW_mll_norm   );
137      normalizeToSum(_h_WW_ptlmax_norm);
138      normalizeToSum(_h_WW_ptlmin_norm);
139      normalizeToSum(_h_WW_dphill_norm);
140
141      scale(_h_WW_njet0, norm);
142
143    }
144
145    /// @}
146
147
148    /// @name Histograms
149    /// @{
150    Histo1DPtr _h_WW_njets_norm;
151    Histo1DPtr _h_WW_mll_norm, _h_WW_ptlmax_norm, _h_WW_ptlmin_norm, _h_WW_dphill_norm;
152    Histo1DPtr _h_WW_njet0;
153    /// @}
154
155  };
156
157
158
159  RIVET_DECLARE_PLUGIN(CMS_2020_I1814328);
160
161}