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