1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/IdentifiedFinalState.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/ChargedFinalState.hh"
  6
  7namespace Rivet {
  8
  9
 10
 11
 12  /// @brief MC validation analysis for photons
 13  /// @todo Rename to MC_DRESSEDPHOTONS, or add these plots to the generic particle analysis photons
 14  class MC_PHOTONS : public Analysis {
 15  public:
 16
 17    /// @name Constructors etc.
 18    /// @{
 19
 20    /// Constructor
 21    MC_PHOTONS()
 22      : Analysis("MC_PHOTONS")
 23    {    }
 24
 25    /// @}
 26
 27
 28    /// @name Analysis methods
 29    /// @{
 30
 31    /// Book histograms and initialise projections before the run
 32    void init() {
 33      // set FS cuts from input options
 34      const double etalcut = getOption<double>("ABSETALMAX", 5.);
 35      const double ptlcut = getOption<double>("PTLMIN", 10.);
 36
 37      IdentifiedFinalState leptons(Cuts::abseta < etalcut && Cuts::pT > ptlcut*GeV);
 38      leptons.acceptChLeptons();
 39      declare(leptons, "lFS");
 40
 41      // set photon cuts from input options
 42      const double etagamcut = getOption<double>("ABSETAGAMMAX", 5.0);
 43
 44      IdentifiedFinalState photons(Cuts::abseta < etagamcut);
 45      photons.acceptId(PID::PHOTON);
 46      declare(photons, "gammaFS");
 47
 48      book(_h_Ptgamma ,"Ptgamma", logspace(50, 0.01, 30));
 49      book(_h_Egamma ,"Egamma", logspace(50, 0.01, 200));
 50      book(_h_sumPtgamma ,"sumPtgamma", 50, 0, 100);
 51      book(_h_sumEgamma ,"sumEgamma", 50, 0, (sqrtS()>0.?sqrtS():14000.)/GeV/5.0);
 52      book(_h_DelR ,"DeltaR", 50, 0, 2);
 53      book(_h_DelR_weighted ,"DeltaR_ptweighted", 50, 0, 2);
 54      book(_h_DelR_R ,"DeltaR_R", 50, 0, 2);
 55      book(_h_DelR_R_weighted ,"DeltaR_R_ptweighted", 50, 0, 2);
 56      book(_p_DelR_vs_pTl ,"DeltaR_vs_pTlep", 50, 10, 120);
 57      book(_p_DelR_weighted_vs_pTl ,"DeltaR_ptweighted_vs_pTlep", 50, 10, 120);
 58      book(_p_DelR_R_vs_pTl ,"DeltaR_R_vs_pTlep", 50, 10, 120);
 59      book(_p_DelR_R_weighted_vs_pTl ,"DeltaR_R_ptweighted_vs_pTlep", 50, 10, 120);
 60      book(_p_sumPtgamma_vs_pTl ,"sumPtGamma_vs_pTlep", 50, 10, 120);
 61    }
 62
 63
 64    /// Perform the per-event analysis
 65    void analyze(const Event& event) {
 66
 67      /// Get photons and leptons
 68      const Particles& photons = apply<FinalState>(event, "gammaFS").particles();
 69      MSG_DEBUG("Photon multiplicity = " << photons.size());
 70      const Particles& leptons = apply<FinalState>(event, "lFS").particles();
 71      MSG_DEBUG("Photon multiplicity = " << leptons.size());
 72
 73      // Initialise a map of sumPtgamma for each lepton
 74      map<size_t, double> sumpT_per_lep;
 75      for (size_t il = 0; il < leptons.size(); ++il) sumpT_per_lep[il] = 0;
 76
 77      // Calculate photon energies and transverse momenta
 78      double sumPtgamma(0), sumEgamma(0);
 79      for (const Particle& p : photons) {
 80        // Individual and summed pTs and energies
 81        double pTgamma = p.pT()/GeV;
 82        double Egamma = p.E()/GeV;
 83        _h_Ptgamma->fill(pTgamma);
 84        _h_Egamma->fill(Egamma);
 85        sumPtgamma += pTgamma;
 86        sumEgamma += Egamma;
 87
 88        // Calculate delta R with respect to the nearest lepton
 89        int ilep = -1;
 90        double delR = 10000;
 91        for (size_t il = 0; il < leptons.size(); ++il) {
 92          const double tmpdelR = deltaR(leptons[il].momentum(), p.momentum());
 93          if (tmpdelR < delR) {
 94            ilep = il;
 95            delR = tmpdelR;
 96          }
 97        }
 98        if (ilep != -1) {
 99          _h_DelR->fill(delR);
100          _h_DelR_weighted->fill(delR, pTgamma/GeV);
101          _h_DelR_R->fill(delR, 1.0/(delR+1e-5));
102          _h_DelR_R_weighted->fill(delR, pTgamma/GeV/(delR+1e-5));
103          _p_DelR_vs_pTl->fill(leptons[ilep].pT()/GeV, delR);
104          _p_DelR_weighted_vs_pTl->fill(leptons[ilep].pT()/GeV, delR, pTgamma/GeV);
105          _p_DelR_R_vs_pTl->fill(leptons[ilep].pT()/GeV, delR, 1.0/(delR+1e-5));
106          _p_DelR_R_weighted_vs_pTl->fill(leptons[ilep].pT()/GeV, delR, pTgamma/GeV/(delR+1e-5));
107          sumpT_per_lep[ilep] += pTgamma;
108        }
109      }
110
111      // Histogram whole-event photon HT/energy
112      _h_sumPtgamma->fill(sumPtgamma/GeV);
113      _h_sumEgamma->fill(sumEgamma/GeV);
114
115      // Histogram per-lepton sum(pT)
116      for (size_t il = 0; il < leptons.size(); ++il) {
117        _p_sumPtgamma_vs_pTl->fill(leptons[il].pT()/GeV, sumpT_per_lep[il]/GeV);
118      }
119
120    }
121
122
123    /// Normalise histograms etc., after the run
124    void finalize() {
125      normalize(_h_Ptgamma);
126      normalize(_h_Egamma);
127      normalize(_h_sumPtgamma);
128      normalize(_h_sumEgamma);
129      normalize(_h_DelR);
130      normalize(_h_DelR_weighted);
131      normalize(_h_DelR_R);
132      normalize(_h_DelR_R_weighted);
133    }
134
135    /// @}
136
137
138  private:
139
140    /// @name Histograms
141    /// @{
142    Histo1DPtr _h_Ptgamma, _h_Egamma;
143    Histo1DPtr _h_sumPtgamma, _h_sumEgamma;
144    Histo1DPtr _h_DelR, _h_DelR_weighted;
145    Histo1DPtr _h_DelR_R, _h_DelR_R_weighted;
146    Profile1DPtr _p_DelR_vs_pTl, _p_DelR_weighted_vs_pTl;
147    Profile1DPtr _p_DelR_R_vs_pTl, _p_DelR_R_weighted_vs_pTl;
148    Profile1DPtr _p_sumPtgamma_vs_pTl;
149    /// @}
150
151  };
152
153
154  RIVET_DECLARE_PLUGIN(MC_PHOTONS);
155
156
157}