Rivet analyses

title: MC_VERTICES

Monte Carlo validation observables for interaction vertices

Experiment: ()

Status: VALIDATED

Authors: - Ynyr Harris - Christian Gutschow

References: none listed

Beams: * *

Beam energies: ANY

Run details: none listed

Vertices can be filtered by HepMC status code using the STATUS option. Signal and MPI vertices are always plotted.

Source code:MC_VERTICES.cc

```c++

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

} ```