← Back to branch summary

~ubuntu-branches/ubuntu/saucy/rivet/saucy-proposed

« back to all changes in this revision

Viewing changes to src/Analyses/MC_PHOTONS.cc

  • Committer: Package Import Robot
  • Author(s): Lifeng Sun
  • Date: 2013-05-07 09:24:27 UTC
  • mfrom: (1.2.1) (2.1.2 experimental)
  • Revision ID: package-import@ubuntu.com-20130507092427-rpfijl4mn3y87ek7
Tags: 1.8.3-1
New upstream release.

Show diffs side-by-side

added added

removed removed

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