1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/FinalState.hh"
 4#include "Rivet/Projections/ChargedLeptons.hh"
 5#include "Rivet/Projections/NeutralFinalState.hh"
 6#include "Rivet/Projections/PromptFinalState.hh"
 7
 8namespace Rivet {
 9
10
11  /// Differential cross-section of FSR photons in Z decays
12  class CMS_2015_I1346843 : public Analysis {
13  public:
14
15    /// Constructor
16    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2015_I1346843);
17
18    /// Book histograms and initialise projections before the run
19    void init() {
20
21      Cut c_photons = Cuts::pT >= 5*GeV && Cuts::abseta < 2.5 && !(Cuts::absetaIn(1.4, 1.6));
22      PromptFinalState photons(Cuts::abspid == PID::PHOTON && c_photons, TauDecaysAs::PROMPT, MuDecaysAs::PROMPT);
23      declare(photons, "PHOTONFS");
24
25      Cut c_muons = Cuts::pT > 9*GeV && Cuts::abseta < 2.4;
26      PromptFinalState muons(Cuts::abspid == PID::MUON && c_muons);
27      declare(muons, "MUFS");
28
29      book(_h["pho_et"]           ,1, 1, 1);  // photon transverse energy
30      book(_h["pho_et_wide"]      ,1, 2, 1);  // photon transverse energy (0.5 < dr < 3.0)
31      book(_h["pho_et_close"]     ,1, 3, 1);  // photon transverse energy (0.05 < dr < 0.5)
32      book(_h["pho_et_lqt"]       ,1, 4, 1);  // photon transverse energy (q_T < 10)
33      book(_h["pho_et_hqt"]       ,1, 5, 1);  // photon transverse energy (q_T > 50)
34      book(_h["pho_dr"]           ,2, 1, 1);  // delta_R
35      book(_h["pho_dr_lqt"]       ,2, 2, 1);  // delta_R (q_T < 10)
36      book(_h["pho_dr_hqt"]       ,2, 3, 1);  // delta_R  (q_T > 50)
37    }
38
39
40    // Perform the per-event analysis
41    void analyze(const Event& event) {
42
43      const Particles muons = apply<PromptFinalState>(event, "MUFS").particlesByPt();
44
45      if (muons.size() < 2) vetoEvent;
46      if (muons[0].pT() < 31*GeV) vetoEvent;
47      if (muons[0].charge()*muons[1].charge() > 0) vetoEvent;
48      const double mZ = (muons[0].momentum() + muons[1].momentum()).mass();
49      if (!inRange(mZ, 30*GeV, 87*GeV)) vetoEvent;
50
51      const Particles photons = apply<PromptFinalState>(event, "PHOTONFS").particlesByPt();
52      // We want the photon with the highest pT that does not come from a decay
53      for (const Particle& p : photons) {
54
55        const double dR = std::min(deltaR(p, muons[0]), deltaR(p, muons[1]) );
56        if (!inRange(dR, 0.05, 3.0)) continue;
57
58        // Calculate the three-body (mu,mu,gamma) transverse momentum
59        const double qT = (muons[0].mom() + muons[1].mom() + p.mom()).pT();
60
61        // Fill the analysis histograms
62        _h["pho_et"]->fill(p.pT()/GeV, 1.0);
63        _h["pho_dr"]->fill(dR, 1.0);
64
65        _h[(dR <= 0.5 ? "pho_et_close" : "pho_et_wide")]->fill(p.pT()/GeV, 1.0);
66
67        if (qT / GeV < 10.) {
68          _h["pho_et_lqt"]->fill(p.pT()/GeV, 1.0);
69          _h["pho_dr_lqt"]->fill(dR, 1.0);
70        }
71
72        if (qT / GeV > 50.) {
73          _h["pho_et_hqt"]->fill(p.pT()/GeV, 1.0);
74          _h["pho_dr_hqt"]->fill(dR, 1.0);
75        }
76
77        break; // Exit the loop since we found the highest pT lepton already
78      }
79    }
80
81
82    /// Normalise histograms etc., after the run
83    void finalize() {
84      scale(_h, crossSection()/picobarn / sumOfWeights());
85    }
86
87
88  private:
89
90    map<std::string,Histo1DPtr> _h;
91
92  };
93
94
95  RIVET_DECLARE_PLUGIN(CMS_2015_I1346843);
96
97}