1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/FinalState.hh"
 4#include "Rivet/Projections/InvisibleFinalState.hh"
 5#include "Rivet/Projections/MissingMomentum.hh"
 6
 7namespace Rivet {
 8
 9
10
11  /// @brief MC validation analysis for truth-MET measurement
12  /// @todo Add plots for MET based on prompt invisibles
13  class MC_MET : public Analysis {
14  public:
15
16    RIVET_DEFAULT_ANALYSIS_CTOR(MC_MET);
17
18
19    void init() {
20      const FinalState inclfs;
21      FinalState calofs(Cuts::abseta < 5);
22      declare(MissingMomentum(inclfs), "InclMET");
23      declare(MissingMomentum(calofs), "CaloMET");
24
25      declare(InvisibleFinalState(), "InvisibleFS");
26      declare(InvisibleFinalState(true), "PromptInvisibleFS");
27
28      book(_h["met_incl"], "met_incl", logspace(50, 10, sqrtS()/GeV/5));
29      book(_h["met_calo"], "met_calo", logspace(50, 10, sqrtS()/GeV/5));
30      book(_h["set_incl"], "set_incl", logspace(50, 10, sqrtS()/GeV/3));
31      book(_h["set_calo"], "set_calo", logspace(50, 10, sqrtS()/GeV/3));
32
33      book(_h["pT_inv"],   "pT_inv",   logspace(50, 10, sqrtS()/GeV/5));
34      book(_h["mass_inv"], "mass_inv", logspace(100, 10, sqrtS()/GeV/5));
35      book(_h["rap_inv"],  "rap_inv",   50, -5., 5.);
36
37      book(_h["pT_promptinv"],   "pT_promptinv",   logspace(50, 10, sqrtS()/GeV/5));
38      book(_h["mass_promptinv"], "mass_promptinv", logspace(100, 10, sqrtS()/GeV/5));
39      book(_h["rap_promptinv"],  "rap_promptinv",  50, -5., 5.);
40    }
41
42
43    void analyze(const Event& event) {
44
45      const MissingMomentum& mmincl = apply<MissingMomentum>(event, "InclMET");
46      _h["met_incl"]->fill(mmincl.met()/GeV);
47      _h["set_incl"]->fill(mmincl.set()/GeV);
48
49      const MissingMomentum& mmcalo = apply<MissingMomentum>(event, "CaloMET");
50      _h["met_calo"]->fill(mmcalo.met()/GeV);
51      _h["set_calo"]->fill(mmcalo.set()/GeV);
52
53      // Get the invisible final state particles
54      const Particles& invisibles = apply<InvisibleFinalState>(event, "InvisibleFS").particlesByPt();
55      const Particles& promptinvisibles = apply<InvisibleFinalState>(event, "PromptInvisibleFS").particlesByPt();
56
57      if (!invisibles.empty()) {
58        FourMomentum invsum;
59        for (const Particle& p : invisibles) {
60          invsum += p.momentum();
61        }
62        _h["pT_inv"]->fill(invsum.pT()/GeV);
63        _h["mass_inv"]->fill(invsum.mass()/GeV);
64        _h["rap_inv"]->fill(invsum.rapidity());
65      }
66
67      if (!promptinvisibles.empty()) {
68        FourMomentum promptinvsum;
69        for (const Particle& p : promptinvisibles) {
70          promptinvsum += p.momentum();
71        }
72        _h["pT_promptinv"]->fill(promptinvsum.pT()/GeV);
73        _h["mass_promptinv"]->fill(promptinvsum.mass()/GeV);
74        _h["rap_promptinv"]->fill(promptinvsum.rapidity());
75      }
76
77    }
78
79
80    void finalize() {
81      const double sf = crossSectionPerEvent()/picobarn;
82      scale(_h, sf);
83    }
84
85
86  private:
87
88    map<string, Histo1DPtr> _h;
89
90  };
91
92
93
94  // The hook for the plugin system
95  RIVET_DECLARE_PLUGIN(MC_MET);
96}