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