Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/InvisibleFinalState.hh"
#include "Rivet/Projections/MissingMomentum.hh"
#include "Rivet/Projections/Smearing.hh"

namespace Rivet {


  /// @brief MC validation analysis for truth and reco MET and SET observables
  class MC_MET : public Analysis {
  public:

    RIVET_DEFAULT_ANALYSIS_CTOR(MC_MET);


    void init() {

      const FinalState inclfs;
      const FinalState calofs(Cuts::abseta < 5);
      const MissingMomentum inclmm(inclfs);
      const MissingMomentum calomm(calofs);

      declare(inclmm, "InclMET");
      declare(calomm, "CaloMET");
      declare(SmearedMET(calomm, MET_SMEARPARAMS_ATLAS_RUN2), "RecoMET");

      declare(InvisibleFinalState(), "InvisibleFS");
      declare(InvisibleFinalState(OnlyDirect::YES), "DirectInvisibleFS");

      const double METMAX = sqrtS()/5, SETMAX = sqrtS()/3;

      book(_h["met_incl"], "met_incl", logspace(50, 10, METMAX/GeV));
      book(_h["met_calo"], "met_calo", logspace(50, 10, METMAX/GeV));
      book(_h["met_reco"], "met_reco", logspace(50, 10, METMAX/GeV));

      book(_h["mpt_incl"], "mpt_incl", logspace(50, 10, METMAX/GeV));
      book(_h["mpt_calo"], "mpt_calo", logspace(50, 10, METMAX/GeV));
      book(_h["mpt_reco"], "mpt_reco", logspace(50, 10, METMAX/GeV));

      book(_h["set_incl"], "set_incl", logspace(50, 10, SETMAX/GeV));
      book(_h["set_calo"], "set_calo", logspace(50, 10, SETMAX/GeV));
      book(_h["set_reco"], "set_reco", logspace(50, 10, SETMAX/GeV));

      book(_h["pT_inv"],       "pT_inv",       logspace(50, 10, METMAX/GeV));
      book(_h["pT_directinv"], "pT_directinv", logspace(50, 10, METMAX/GeV));

      book(_h["mass_inv"],       "mass_inv",       logspace(100, 10, METMAX/GeV));
      book(_h["mass_directinv"], "mass_directinv", logspace(100, 10, METMAX/GeV));

      book(_h["rap_inv"],       "rap_inv",       50, -5., 5.);
      book(_h["rap_directinv"], "rap_directinv", 50, -5., 5.);
    }


    void analyze(const Event& event) {

      // Inclusive truth MET, etc.
      const METFinder& mmincl = apply<MissingMomentum>(event, "InclMET");
      _h["met_incl"]->fill(mmincl.met()/GeV);
      _h["mpt_incl"]->fill(mmincl.mpt()/GeV);
      _h["set_incl"]->fill(mmincl.set()/GeV);

      // Calo-acceptance truth MET, etc.
      const METFinder& mmcalo = apply<MissingMomentum>(event, "CaloMET");
      _h["met_calo"]->fill(mmcalo.met()/GeV);
      _h["mpt_calo"]->fill(mmcalo.mpt()/GeV);
      _h["set_calo"]->fill(mmcalo.set()/GeV);

      // (ATLAS) reco MET, etc.
      const METFinder& mmreco = apply<SmearedMET>(event, "RecoMET");
      _h["met_reco"]->fill(mmreco.met()/GeV);
      _h["mpt_reco"]->fill(mmreco.mpt()/GeV);
      _h["set_reco"]->fill(mmreco.set()/GeV);

      // Compute alternative MET from invisible particles, including from hadron decays
      const ParticleFinder& allinvisibles = apply<FinalState>(event, "InvisibleFS");
      if (!allinvisibles.empty()) {
        FourMomentum invsum = sumP4(allinvisibles.particles());
        _h["pT_inv"]->fill(invsum.pT()/GeV);
        _h["mass_inv"]->fill(invsum.mass()/GeV);
        _h["rap_inv"]->fill(invsum.rapidity());
      }

      // Compute alternative MET from direct, invisible particles only
      const ParticleFinder& directinvisibles = apply<FinalState>(event, "DirectInvisibleFS");
      if (!directinvisibles.empty()) {
        FourMomentum invsum = sumP4(directinvisibles.particles());
        _h["pT_directinv"]->fill(invsum.pT()/GeV);
        _h["mass_directinv"]->fill(invsum.mass()/GeV);
        _h["rap_directinv"]->fill(invsum.rapidity());
      }

    }


    void finalize() {
      const double sf = crossSectionPerEvent()/picobarn;
      scale(_h, sf);
    }


  private:

    map<string, Histo1DPtr> _h;

  };


  RIVET_DECLARE_PLUGIN(MC_MET);

}