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      IdentifiedFinalState leptons(Cuts::abseta < 5.0 && Cuts::pT > 10*GeV);
 34      leptons.acceptChLeptons();
 35      declare(leptons, "lFS");
 36
 37      IdentifiedFinalState photons(Cuts::abseta < 5.0);
 38      photons.acceptId(PID::PHOTON);
 39      declare(photons, "gammaFS");
 40
 41      book(_h_Ptgamma ,"Ptgamma", logspace(50, 0.01, 30));
 42      book(_h_Egamma ,"Egamma", logspace(50, 0.01, 200));
 43      book(_h_sumPtgamma ,"sumPtgamma", 50, 0, 100);
 44      book(_h_sumEgamma ,"sumEgamma", 50, 0, (sqrtS()>0.?sqrtS():14000.)/GeV/5.0);
 45      book(_h_DelR ,"DeltaR", 50, 0, 2);
 46      book(_h_DelR_weighted ,"DeltaR_ptweighted", 50, 0, 2);
 47      book(_h_DelR_R ,"DeltaR_R", 50, 0, 2);
 48      book(_h_DelR_R_weighted ,"DeltaR_R_ptweighted", 50, 0, 2);
 49      book(_p_DelR_vs_pTl ,"DeltaR_vs_pTlep", 50, 10, 120);
 50      book(_p_DelR_weighted_vs_pTl ,"DeltaR_ptweighted_vs_pTlep", 50, 10, 120);
 51      book(_p_DelR_R_vs_pTl ,"DeltaR_R_vs_pTlep", 50, 10, 120);
 52      book(_p_DelR_R_weighted_vs_pTl ,"DeltaR_R_ptweighted_vs_pTlep", 50, 10, 120);
 53      book(_p_sumPtgamma_vs_pTl ,"sumPtGamma_vs_pTlep", 50, 10, 120);
 54    }
 55
 56
 57    /// Perform the per-event analysis
 58    void analyze(const Event& event) {
 59      const double weight = 1.0;
 60
 61      /// Get photons and leptons
 62      const Particles& photons = apply<FinalState>(event, "gammaFS").particles();
 63      MSG_DEBUG("Photon multiplicity = " << photons.size());
 64      const Particles& leptons = apply<FinalState>(event, "lFS").particles();
 65      MSG_DEBUG("Photon multiplicity = " << leptons.size());
 66
 67      // Initialise a map of sumPtgamma for each lepton
 68      map<size_t, double> sumpT_per_lep;
 69      for (size_t il = 0; il < leptons.size(); ++il) sumpT_per_lep[il] = 0;
 70
 71      // Calculate photon energies and transverse momenta
 72      double sumPtgamma(0), sumEgamma(0);
 73      for (const Particle& p : photons) {
 74        // Individual and summed pTs and energies
 75        double pTgamma = p.pT()/GeV;
 76        double Egamma = p.E()/GeV;
 77        _h_Ptgamma->fill(pTgamma, weight);
 78        _h_Egamma->fill(Egamma, weight);
 79        sumPtgamma += pTgamma;
 80        sumEgamma += Egamma;
 81
 82        // Calculate delta R with respect to the nearest lepton
 83        int ilep = -1;
 84        double delR = 10000;
 85        for (size_t il = 0; il < leptons.size(); ++il) {
 86          const double tmpdelR = deltaR(leptons[il].momentum(), p.momentum());
 87          if (tmpdelR < delR) {
 88            ilep = il;
 89            delR = tmpdelR;
 90          }
 91        }
 92        if (ilep != -1) {
 93          _h_DelR->fill(delR, weight);
 94          _h_DelR_weighted->fill(delR, weight*pTgamma/GeV);
 95          _h_DelR_R->fill(delR, weight/(delR+1e-5));
 96          _h_DelR_R_weighted->fill(delR, weight*pTgamma/GeV/(delR+1e-5));
 97          _p_DelR_vs_pTl->fill(leptons[ilep].pT()/GeV, delR, weight);
 98          _p_DelR_weighted_vs_pTl->fill(leptons[ilep].pT()/GeV, delR, weight*pTgamma/GeV);
 99          _p_DelR_R_vs_pTl->fill(leptons[ilep].pT()/GeV, delR, weight/(delR+1e-5));
100          _p_DelR_R_weighted_vs_pTl->fill(leptons[ilep].pT()/GeV, delR, weight*pTgamma/GeV/(delR+1e-5));
101          sumpT_per_lep[ilep] += pTgamma;
102        }
103      }
104
105      // Histogram whole-event photon HT/energy
106      _h_sumPtgamma->fill(sumPtgamma/GeV, weight);
107      _h_sumEgamma->fill(sumEgamma/GeV, weight);
108
109      // Histogram per-lepton sum(pT)
110      for (size_t il = 0; il < leptons.size(); ++il) {
111        _p_sumPtgamma_vs_pTl->fill(leptons[il].pT()/GeV, sumpT_per_lep[il]/GeV, weight);
112      }
113
114    }
115
116
117    /// Normalise histograms etc., after the run
118    void finalize() {
119      normalize(_h_Ptgamma);
120      normalize(_h_Egamma);
121      normalize(_h_sumPtgamma);
122      normalize(_h_sumEgamma);
123      normalize(_h_DelR);
124      normalize(_h_DelR_weighted);
125      normalize(_h_DelR_R);
126      normalize(_h_DelR_R_weighted);
127    }
128
129    //@}
130
131
132  private:
133
134    /// @name Histograms
135    //@{
136    Histo1DPtr _h_Ptgamma, _h_Egamma;
137    Histo1DPtr _h_sumPtgamma, _h_sumEgamma;
138    Histo1DPtr _h_DelR, _h_DelR_weighted;
139    Histo1DPtr _h_DelR_R, _h_DelR_R_weighted;
140    Profile1DPtr _p_DelR_vs_pTl, _p_DelR_weighted_vs_pTl;
141    Profile1DPtr _p_DelR_R_vs_pTl, _p_DelR_R_weighted_vs_pTl;
142    Profile1DPtr _p_sumPtgamma_vs_pTl;
143    //@}
144
145  };
146
147
148  // The hook for the plugin system
149  RIVET_DECLARE_PLUGIN(MC_PHOTONS);
150
151
152}