Rivet analyses

#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"

namespace Rivet {

  /// @brief Monte Carlo validation observables for interaction vertices
  class MC_VERTICES : public Analysis {
  public:

    RIVET_DEFAULT_ANALYSIS_CTOR(MC_VERTICES);

    void init() {

      _status = getOption<int>("STATUS", -1);

      declare(FinalState(), "fs");

      // Book histograms
      book(_h1["ip_mm"], "h_impact_parameter_wide", 50, -2.0, 2.0);
      book(_h1["ip_fm"], "h_impact_parameter_narrow", 50, -2.0, 2.0);

      book(_h2["hs"], "h_signal_vertex_xy", 50, -2.0, 2.0, 50, -2.0, 2.0);
      book(_h2["mpi"], "h_mpi_vertex_xy", 50, -2.0, 2.0, 50, -2.0, 2.0);

      const vector<double> logedges = logspace(50, 1e-2, 2.0);
      book(_h2["hs_log"], "h_signal_vertex_xy_log", logedges, logedges);
      book(_h2["mpi_log"], "h_mpi_vertex_xy_log", logedges, logedges);

      if (_status >= 0) {
        const string hname("h_status_" + toString(_status) + "_vertex_xy");
        book(_h2["status"], hname, 50, -2.0, 2.0, 50, -2.0, 2.0);
        book(_h2["status_log"], hname + "_vertex_xy_log", logedges, logedges);
      }
    }

    void analyze(const Event& event) {

      const HepMC3::GenEvent* genEvt = event.genEvent();

      for (const Particle& p : apply<FinalState>(event, "fs").particles()) {
        const double signedIP = p.impactParam(true);
        _h1["ip_mm"]->fill(signedIP / millimetre);
        _h1["ip_fm"]->fill(signedIP / fermi);
      }

      // If requested, filter by status and full corresponding subset
      if (_status >= 0) {
        for (ConstGenVertexPtr vtx : genEvt->vertices()) {
          if (_status != vtx->status()) continue;
          const HepMC3::FourVector pos = vtx->position();
          const double x = pos.x();
          const double y = pos.y();
          MSG_DEBUG("New vertex with status " << vtx->status() << " at position x = " << x << " y = " << y);
          _h2["status"]->fill(x / fermi, y / fermi);
          _h2["status_log"]->fill(std::abs(x) / fermi, std::abs(y) / fermi);
        }
      }

      // Access signal-process and MPI vertices
      auto spv_idx = genEvt->attribute<HepMC3::IntAttribute>("signal_process_vertex");
      auto mpi_idx = genEvt->attribute<HepMC3::IntAttribute>("mpi");

      if (spv_idx) {
        auto spv = genEvt->vertices().at(-spv_idx->value() - 1);
        const HepMC3::FourVector pos = spv->position();
        const double x = pos.x();
        const double y = pos.y();
        _h2["hs"]->fill(x / fermi, y / fermi);
        _h2["hs_log"]->fill(std::abs(x) / fermi, std::abs(y) / fermi);
      }

      if (mpi_idx) {
        auto mpi = genEvt->vertices().at(-mpi_idx->value() - 1);
        const HepMC3::FourVector pos = mpi->position();
        const double x = pos.x();
        const double y = pos.y();
        _h2["mpi"]->fill(x / fermi, y / fermi);
        _h2["mpi_log"]->fill(std::abs(x) / fermi, std::abs(y) / fermi);
      }
    }

    void finalize() {
      normalize(_h1);
      normalize(_h2);
    }

  private:

    map<string, Histo2DPtr> _h2;
    map<string, Histo1DPtr> _h1;

    int _status;
  };

  RIVET_DECLARE_PLUGIN(MC_VERTICES);

}