Rivet analyses

title: DELPHI_2000_I531568

Rapidity structure of $p\bar{p}$ pairs

Experiment: DELPHI (LEP)

Inspire ID: 531568

Status: VALIDATED

Authors: - Peter Richardson

References: - Phys.Lett.B 490 (2000) 61-70

Beams: e+ e-

Beam energies: (45.6, 45.6)GeV

Run details: - e+e- -> hadrons

Rapidity-rank structure of $p\bar{p}$ pairs, the relative occurrence of the rapidity-ordered configuration $pM\bar{p}$, where $M$ is a meson, and that of adjacent $p\bar{p}$ pairs is compared as a function of the rapidity w.r.t the thrust axis.

Source code:DELPHI_2000_I531568.cc

```c++ // -- 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); } ```