1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/FinalState.hh"
 4#include "Rivet/Projections/LeptonFinder.hh"
 5#include "Rivet/Projections/MissingMomentum.hh"
 6#include "Rivet/Projections/PromptFinalState.hh"
 7
 8namespace Rivet {
 9
10  /// @brief WWgamma production in emu final state at 13 TeV
11  class CMS_2023_I2709669 : public Analysis {
12  public:
13    /// Constructor
14    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2023_I2709669);
15
16    /// @name Analysis methods
17    /// @{
18    void init() {
19      // Initialise and register projections
20
21      // The basic final-state projection:
22      // all final-state particles within
23      // the given eta acceptance
24
25      const FinalState fs(Cuts::abseta < 4.9);
26
27      // FinalState of prompt photons and bare muons and electrons in the event
28      PromptFinalState photons(Cuts::abspid == PID::PHOTON);
29      PromptFinalState bare_electrons(Cuts::abspid == PID::ELECTRON, TauDecaysAs::PROMPT);
30      PromptFinalState bare_muons(Cuts::abspid == PID::MUON, TauDecaysAs::PROMPT);
31      // Dress the prompt bare leptons with prompt photons within dR < 0.1,
32      // and apply some fiducial cuts on the dressed leptons
33
34      Cut muon_cuts = Cuts::abseta < 2.4 && Cuts::pT > 20 * GeV;
35      Cut electron_cuts = Cuts::abseta < 2.5 && Cuts::pT > 25 * GeV;
36      LeptonFinder dressed_electrons(bare_electrons, photons, 0.1, electron_cuts);
37      LeptonFinder dressed_muons(bare_muons, photons, 0.1, muon_cuts);
38      declare(dressed_electrons, "electrons");
39      declare(dressed_muons, "muons");
40      declare(photons, "photons");
41
42      // Missing momentum
43      declare(MissingMomentum(fs), "MET");
44
45      // Book histograms
46      book(_h_xsec_exp, 1, 1, 1);
47    }
48
49    /// Perform the per-event analysis
50    void analyze(const Event& event) {
51      // Retrieve dressed leptons, sorted by pT
52      Particles muons = apply<LeptonFinder>(event, "muons").particlesByPt();
53      Particles electrons = apply<LeptonFinder>(event, "electrons").particlesByPt();
54
55      Particles photons =
56          apply<PromptFinalState>(event, "photons").particles(Cuts::abseta < 2.5 && Cuts::pT > 20 * GeV);
57      idiscardIfAnyDeltaRLess(photons, muons, 0.5);
58      idiscardIfAnyDeltaRLess(photons, electrons, 0.5);
59      idiscardIfAnyDeltaRLess(electrons, muons, 0.5);
60
61      // Retrieve Etmiss variable
62      const double Etmiss = apply<MissingMomentum>(event, "MET").missingEt();
63
64      if ((muons.size() < 1) || (electrons.size() < 1))  vetoEvent;
65      if (photons.size() != 1)  vetoEvent;
66      Particle muon = muons.at(0);
67      Particle electron = electrons.at(0);
68      if (muon.pid() * electron.pid() > 0)  vetoEvent;
69      Particle photon = photons.at(0);
70      FourMomentum LL = muon.momentum() + electron.momentum();
71      if (LL.mass() < 10 * GeV)  vetoEvent;
72      if (LL.pt() < 15 * GeV)    vetoEvent;
73
74      FourMomentum EtMiss = apply<MissingMomentum>(event, "MET").missingMomentum();
75      const double dphi = deltaPhi(LL, EtMiss);
76      const double mT = sqrt(2 * LL.pT() * EtMiss.pT() * (1 - cos(dphi)));
77
78      if (mT < 10*GeV)      vetoEvent;
79      if (Etmiss < 20*GeV)  vetoEvent;
80      _h_xsec_exp->fill("Results"s);
81    }
82
83    /// Normalise histograms etc., after the run
84    void finalize() {
85      const double norm = crossSection() / femtobarn / sumOfWeights();
86
87      scale(_h_xsec_exp, norm);
88    }
89
90    BinnedHistoPtr<string> _h_xsec_exp;
91  };
92
93  RIVET_DECLARE_PLUGIN(CMS_2023_I2709669);
94
95}  // namespace Rivet