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/IdentifiedFinalState.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.0*GeV && Cuts::abseta < 2.5 && !(Cuts::absetaIn(1.4, 1.6));
 22      IdentifiedFinalState photons(c_photons);
 23      photons.acceptId(PID::PHOTON);
 24      declare(photons, "PHOTFS");
 25
 26      Cut c_muons   = Cuts::pT > 9*GeV && Cuts::abseta < 2.4;
 27      IdentifiedFinalState muons(c_muons);
 28      muons.acceptIdPair(PID::MUON);
 29      declare(muons, "MUFS");
 30
 31
 32      book(_hist_pho_et           ,1, 1, 1);  // photon transverse energy
 33      book(_hist_pho_et_wide      ,1, 2, 1);  // photon transverse energy (0.5 < dr < 3.0)
 34      book(_hist_pho_et_close     ,1, 3, 1);  // photon transverse energy (0.05 < dr < 0.5)
 35      book(_hist_pho_et_lqt       ,1, 4, 1);  // photon transverse energy (q_T < 10)
 36      book(_hist_pho_et_hqt       ,1, 5, 1);  // photon transverse energy (q_T > 50)
 37      book(_hist_pho_dr           ,2, 1, 1);  // delta_R
 38      book(_hist_pho_dr_lqt       ,2, 2, 1);  // delta_R (q_T < 10)
 39      book(_hist_pho_dr_hqt       ,2, 3, 1);  // delta_R  (q_T > 50)
 40    }
 41
 42
 43    // Perform the per-event analysis
 44    void analyze(const Event& event) {
 45
 46      const Particles muons = apply<IdentifiedFinalState>(event, "MUFS").particlesByPt();
 47
 48      if (muons.size() < 2) vetoEvent;
 49      if (muons[0].pT()/GeV < 31) vetoEvent;
 50      if (muons[0].charge()*muons[1].charge() > 0) vetoEvent;
 51      const double mZ = (muons[0].momentum() + muons[1].momentum()).mass();
 52      if (!inRange(mZ, 30*GeV, 87*GeV)) vetoEvent;
 53
 54      const Particles photons = apply<IdentifiedFinalState>(event, "PHOTFS").particlesByPt();
 55      // We want the photon with the highest pT that does not come from a decay
 56      for (const Particle& p : photons) {
 57        if (!p.isDirect()) continue;
 58        if (!p.isStable()) continue;
 59
 60        const double dR = std::min(deltaR(p, muons[0]), deltaR(p, muons[1]) );
 61        if (!inRange(dR, 0.05, 3.0)) continue;
 62
 63        // Calculate the three-body (mu,mu,gamma) transverse momentum
 64        const double qT = (muons[0].mom() + muons[1].mom() + p.mom()).pT();
 65
 66        // Fill the analysis histograms
 67        _hist_pho_et->fill(p.pT()/GeV, 1.0);
 68        _hist_pho_dr->fill(dR, 1.0);
 69
 70        (dR <= 0.5 ? _hist_pho_et_close : _hist_pho_et_wide)->fill(p.pT()/GeV, 1.0);
 71
 72        if (qT / GeV < 10.) {
 73          _hist_pho_et_lqt->fill(p.pT()/GeV, 1.0);
 74          _hist_pho_dr_lqt->fill(dR, 1.0);
 75        }
 76
 77        if (qT / GeV > 50.) {
 78          _hist_pho_et_hqt->fill(p.pT()/GeV, 1.0);
 79          _hist_pho_dr_hqt->fill(dR, 1.0);
 80        }
 81
 82        break; // Exit the loop since we found the highest pT lepton already
 83      }
 84    }
 85
 86
 87    /// Normalise histograms etc., after the run
 88    void finalize() {
 89      scale(_hist_pho_et,       crossSection() / sumOfWeights());
 90      scale(_hist_pho_et_wide,  crossSection() / sumOfWeights());
 91      scale(_hist_pho_et_close, crossSection() / sumOfWeights());
 92      scale(_hist_pho_et_lqt,   crossSection() / sumOfWeights());
 93      scale(_hist_pho_et_hqt,   crossSection() / sumOfWeights());
 94      scale(_hist_pho_dr,       crossSection() / sumOfWeights());
 95      scale(_hist_pho_dr_lqt,   crossSection() / sumOfWeights());
 96      scale(_hist_pho_dr_hqt,   crossSection() / sumOfWeights());
 97    }
 98
 99
100  private:
101
102    Histo1DPtr _hist_pho_et;
103    Histo1DPtr _hist_pho_et_wide, _hist_pho_et_close;
104    Histo1DPtr _hist_pho_et_lqt,  _hist_pho_et_hqt;
105    Histo1DPtr _hist_pho_dr;
106    Histo1DPtr _hist_pho_dr_lqt, _hist_pho_dr_hqt;
107
108  };
109
110
111  RIVET_DECLARE_PLUGIN(CMS_2015_I1346843);
112
113}