rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

MC_PHOTONS

Monte Carlo validation observables for general photons
Experiment: ()
Status: VALIDATED
Authors:
  • Steve Lloyd
  • Andy Buckley
No references listed
Beams: * *
Beam energies: ANY
Run details:
  • Any event type, but there are many observables for photons associated to (semi-)hard leptons.

Observables for testing general unisolated photon properties, especially those associated with charged leptons (e and mu).

Source code: MC_PHOTONS.cc
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/IdentifiedFinalState.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/ChargedFinalState.hh"

namespace Rivet {

  


  /// @brief MC validation analysis for photons
  /// @todo Rename to MC_DRESSEDPHOTONS, or add these plots to the generic particle analysis photons
  class MC_PHOTONS : public Analysis {
  public:

    /// @name Constructors etc.
    //@{

    /// Constructor
    MC_PHOTONS()
      : Analysis("MC_PHOTONS")
    {    }

    //@}


    /// @name Analysis methods
    //@{

    /// Book histograms and initialise projections before the run
    void init() {
      IdentifiedFinalState leptons(Cuts::abseta < 5.0 && Cuts::pT > 10*GeV);
      leptons.acceptChLeptons();
      declare(leptons, "lFS");

      IdentifiedFinalState photons(Cuts::abseta < 5.0);
      photons.acceptId(PID::PHOTON);
      declare(photons, "gammaFS");

      _h_Ptgamma = bookHisto1D("Ptgamma", logspace(50, 0.01, 30));
      _h_Egamma = bookHisto1D("Egamma", logspace(50, 0.01, 200));
      _h_sumPtgamma = bookHisto1D("sumPtgamma", 50, 0, 100);
      _h_sumEgamma = bookHisto1D("sumEgamma", 50, 0, (sqrtS()>0.?sqrtS():14000.)/GeV/5.0);
      _h_DelR = bookHisto1D("DeltaR", 50, 0, 2);
      _h_DelR_weighted = bookHisto1D("DeltaR_ptweighted", 50, 0, 2);
      _h_DelR_R = bookHisto1D("DeltaR_R", 50, 0, 2);
      _h_DelR_R_weighted = bookHisto1D("DeltaR_R_ptweighted", 50, 0, 2);
      _p_DelR_vs_pTl = bookProfile1D("DeltaR_vs_pTlep", 50, 10, 120);
      _p_DelR_weighted_vs_pTl = bookProfile1D("DeltaR_ptweighted_vs_pTlep", 50, 10, 120);
      _p_DelR_R_vs_pTl = bookProfile1D("DeltaR_R_vs_pTlep", 50, 10, 120);
      _p_DelR_R_weighted_vs_pTl = bookProfile1D("DeltaR_R_ptweighted_vs_pTlep", 50, 10, 120);
      _p_sumPtgamma_vs_pTl = bookProfile1D("sumPtGamma_vs_pTlep", 50, 10, 120);
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      const double weight = event.weight();

      /// Get photons and leptons
      const Particles& photons = apply<FinalState>(event, "gammaFS").particles();
      MSG_DEBUG("Photon multiplicity = " << photons.size());
      const Particles& leptons = apply<FinalState>(event, "lFS").particles();
      MSG_DEBUG("Photon multiplicity = " << leptons.size());

      // Initialise a map of sumPtgamma for each lepton
      map<size_t, double> sumpT_per_lep;
      for (size_t il = 0; il < leptons.size(); ++il) sumpT_per_lep[il] = 0;

      // Calculate photon energies and transverse momenta
      double sumPtgamma(0), sumEgamma(0);
      foreach (const Particle& p, photons) {
        // Individual and summed pTs and energies
        double pTgamma = p.pT()/GeV;
        double Egamma = p.E()/GeV;
        _h_Ptgamma->fill(pTgamma, weight);
        _h_Egamma->fill(Egamma, weight);
        sumPtgamma += pTgamma;
        sumEgamma += Egamma;

        // Calculate delta R with respect to the nearest lepton
        int ilep = -1;
        double delR = 10000;
        for (size_t il = 0; il < leptons.size(); ++il) {
          const double tmpdelR = deltaR(leptons[il].momentum(), p.momentum());
          if (tmpdelR < delR) {
            ilep = il;
            delR = tmpdelR;
          }
        }
        if (ilep != -1) {
          _h_DelR->fill(delR, weight);
          _h_DelR_weighted->fill(delR, weight*pTgamma/GeV);
          _h_DelR_R->fill(delR, weight/(delR+1e-5));
          _h_DelR_R_weighted->fill(delR, weight*pTgamma/GeV/(delR+1e-5));
          _p_DelR_vs_pTl->fill(leptons[ilep].pT()/GeV, delR, weight);
          _p_DelR_weighted_vs_pTl->fill(leptons[ilep].pT()/GeV, delR, weight*pTgamma/GeV);
          _p_DelR_R_vs_pTl->fill(leptons[ilep].pT()/GeV, delR, weight/(delR+1e-5));
          _p_DelR_R_weighted_vs_pTl->fill(leptons[ilep].pT()/GeV, delR, weight*pTgamma/GeV/(delR+1e-5));
          sumpT_per_lep[ilep] += pTgamma;
        }
      }

      // Histogram whole-event photon HT/energy
      _h_sumPtgamma->fill(sumPtgamma/GeV, weight);
      _h_sumEgamma->fill(sumEgamma/GeV, weight);

      // Histogram per-lepton sum(pT)
      for (size_t il = 0; il < leptons.size(); ++il) {
        _p_sumPtgamma_vs_pTl->fill(leptons[il].pT()/GeV, sumpT_per_lep[il]/GeV, weight);
      }

    }


    /// Normalise histograms etc., after the run
    void finalize() {
      normalize(_h_Ptgamma);
      normalize(_h_Egamma);
      normalize(_h_sumPtgamma);
      normalize(_h_sumEgamma);
      normalize(_h_DelR);
      normalize(_h_DelR_weighted);
      normalize(_h_DelR_R);
      normalize(_h_DelR_R_weighted);
    }

    //@}


  private:

    /// @name Histograms
    //@{
    Histo1DPtr _h_Ptgamma, _h_Egamma;
    Histo1DPtr _h_sumPtgamma, _h_sumEgamma;
    Histo1DPtr _h_DelR, _h_DelR_weighted;
    Histo1DPtr _h_DelR_R, _h_DelR_R_weighted;
    Profile1DPtr _p_DelR_vs_pTl, _p_DelR_weighted_vs_pTl;
    Profile1DPtr _p_DelR_R_vs_pTl, _p_DelR_R_weighted_vs_pTl;
    Profile1DPtr _p_sumPtgamma_vs_pTl;
    //@}

  };


  // The hook for the plugin system
  DECLARE_RIVET_PLUGIN(MC_PHOTONS);


}