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);

}