Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
#include "Rivet/Projections/Thrust.hh"

namespace Rivet {


  /// @brief p pbar correlations
  class DELPHI_2000_I531568 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(DELPHI_2000_I531568);


    /// @name Analysis methods
    /// @{

    /// Book histograms and initialise projections before the run
    void init() {

      // Initialise and register projections
      const ChargedFinalState cfs;
      declare(cfs, "CFS");
      declare(Thrust(cfs), "Thrust");

      // book histos
      book(_h_pMp, "_n_pMp", 8, 0., 1.);
      book(_h_sum, "_n_sum", 8, 0., 1.);
    }

    double findPP(const size_t mode, const Vector3& axis, const Particles& part, size_t& pp) const {
      pp = 0;
      const size_t n = part.size();
      double dy = numeric_limits<double>::infinity();
      // compute rapidity with respect to axis
      vector<pair<double, size_t>> keys;
      keys.reserve(n);
      for (size_t i = 0; i < n; ++i) {
        const double mom = dot(axis, part[i].p3());
        const double energy = part[i].E();
        double rap = 0.5 * std::log((energy + mom) / (energy - mom));
        if (!isfinite(rap)) continue; //< skip nans
        if (mode == 0 && rap > 0.)
          keys.emplace_back(rap, i);
        else if (mode == 1 && fuzzyLessEquals(rap, 0.0))
          keys.emplace_back(rap, i);
      }
      // sort by rapidity, break ties by original index for determinism
      std::stable_sort(keys.begin(), keys.end(), [](const auto& a, const auto& b) {
        if (fuzzyLessEquals(a.first, b.first)) return true;
        if (fuzzyLessEquals(b.first, a.first)) return false;
        return a.second < b.second;
      });
      // rebuild particles and rapidities in sorted order
      Particles particles;
      particles.reserve(n);
      vector<double> rapvals;
      rapvals.reserve(n);
      for (const auto& [rap, idx] : keys) {
        particles += part[idx];
        rapvals.push_back(rap);
      }
      // number of particles between the proton and antiproton
      auto itp1 = std::find_if(particles.begin(), particles.end(),
                               [](const Particle& p) { return p.abspid() == PID::PROTON; });
      auto itp2 = std::find_if(itp1 + 1, particles.end(),
                               [](const Particle& p) { return p.abspid() == PID::PROTON; });
      size_t rank = std::distance(itp1, itp2);
      // number of particles between the proton and antiproton
      if (rank > 4) return dy;
      size_t idx1 = itp1 - particles.begin();
      size_t idx2 = itp2 - particles.begin();
      // proton/antiproton next to each other, distance to nearest meson near them
      if (rank == 1) {
        if (idx1 > 0) {
          pp = 1;
          dy = min(2. / 3. * std::abs(rapvals[idx1] - rapvals[idx1 - 1]), dy);
        }
        if (idx2 + 1 < particles.size()) {
          pp = 1;
          dy = min(2. / 3. * std::abs(rapvals[idx2] - rapvals[idx2 + 1]), dy);
        }
      }
      else {
        double ycent = 0.5 * (rapvals[idx1] + rapvals[idx2]);
        double ymin = numeric_limits<double>::infinity();
        double ymes = numeric_limits<double>::infinity();
        for (size_t im = idx1 + 1; im != idx2; ++im) {
          double rap = rapvals[im];
          double test = std::abs(ycent - rap);
          if (test < ymin) {
            ymes = rap;
            ymin = test;
          }
        }
        pp = 2;
        dy = std::min(std::abs(rapvals[idx1] - ymes), std::abs(rapvals[idx2] - ymes));
      }
      return dy;
    }

    /// Perform the per-event analysis
    void analyze(const Event& event) {
      const Vector3 axis = apply<Thrust>(event, "Thrust").thrustAxis();
      const Particles chargedParticles = apply<ChargedFinalState>(event, "CFS").particlesByPt();
      unsigned int np[2] = {0, 0}, npbar[2] = {0, 0};
      for (const Particle& p : chargedParticles) {
        if (p.abspid() != PID::PROTON) continue;
        const double mom = dot(axis, p.p3());
        const double energy = p.E();
        const double rap = 0.5 * std::log((energy + mom) / (energy - mom));
        unsigned int irap = rap > 0 ? 0 : 1;
        if (p.pid() > 0)
          ++np[irap];
        else
          ++npbar[irap];
      }
      MSG_DEBUG("np[0]=" << np[0] << " np[1]=" << np[1]);
      MSG_DEBUG("npbar[0]=" << npbar[0] << " npbar[1]=" << npbar[1]);
      if (np[0] == 1 && npbar[0] == 1) {
        size_t pp = 0;
        double dy = findPP(0, axis, chargedParticles, pp);
        if (pp == 1)
          _h_sum->fill(dy);
        else if (pp == 2) {
          _h_sum->fill(dy);
          _h_pMp->fill(dy);
        }
      }
      if (np[1] == 1 && npbar[1] == 1) {
        size_t pp = 0;
        double dy = findPP(1, axis, chargedParticles, pp);
        if (pp == 1)
          _h_sum->fill(dy);
        else if (pp == 2) {
          _h_sum->fill(dy);
          _h_pMp->fill(dy);
        }
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      Estimate1DPtr h_r;
      book(h_r, 1, 1, 1);
      divide(_h_pMp, _h_sum, h_r);
    }

    /// @}


    /// @name Histograms
    /// @{
    Histo1DPtr _h_sum, _h_pMp;
    /// @}
  };

  RIVET_DECLARE_PLUGIN(DELPHI_2000_I531568);
}