Rivet analyses

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

namespace Rivet {


  /// @brief psi(2S) production at 7 and 13 TeV
  class LHCB_2019_I1748712 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(LHCB_2019_I1748712);


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

    /// Book histograms and initialise projections before the run
    void init() {
      // projections
      declare(UnstableParticles(), "UFS");

      // histograms
      for (double eVal : allowedEnergies()) {
        const int en = round(eVal);
        if (isCompatibleWithSqrtS(eVal)) _sqs = en;
        int ih = en == 7000;
        for (size_t ix = 0; ix < 3; ++ix) {
          book(_h_psi[en + ix], {2.0, 2.5, 3.0, 3.5, 4.0, 4.5});
          for (auto& b : _h_psi[en + ix]->bins()) {
            if (ix < 2)
              book(b, 1 + ix + 2 * ih, 1, b.index());
            else
              book(b, "TMP/psi_" + toString(en + b.index() - 1), refData(1 + 2 * ih, 1, b.index()));
          }
        }

        for (size_t ix = 0; ix < 2; ++ix) {
          book(_h_pT[en + ix], 5 + 2 * ih, 1, ix + 1);
          book(_h_y[en + ix], 6 + 2 * ih, 1, ix + 1);
          if (ih == 0) {
            book(_h_pT_J[en + ix], "TMP/Jpsi_pT_" + toString(en + ix), refData(11, 1, ix + 1));
            book(_h_pT_2[en + ix], "TMP/psi_pT_" + toString(en + ix), refData(11, 1, ix + 1));
            book(_h_y_J[en + ix], "TMP/Jpsi_y_" + toString(en + ix), refData(12, 1, ix + 1));
          }
        }
      }
      raiseBeamErrorIf(_sqs == 0);
    }

    /// Perform the per-event analysis
    void analyze(const Event& event) {
      // Final state of unstable particles to get particle spectra
      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
      // J/psi
      for (const Particle& p : ufs.particles(Cuts::pid == 443 || Cuts::pid == 100443)) {
        // prompt/non-prompt
        bool nonPrompt = p.fromBottom();
        double absrap = p.absrap();
        double xp = p.perp() / GeV;
        if (absrap < 2. || absrap > 4.5) continue;
        if (_sqs == 13000 && (xp < 2 || xp > 20.))
          continue;
        else if (_sqs == 7000 && (xp < 3.5 || xp > 14.))
          continue;
        if (p.pid() == 100443) {
          _h_psi[_sqs + nonPrompt]->fill(absrap, xp);
          _h_psi[_sqs + 2]->fill(absrap, xp);
          _h_pT[_sqs + nonPrompt]->fill(xp);
          _h_y[_sqs + nonPrompt]->fill(absrap);
          if (_sqs == 13000) _h_pT_2[_sqs + nonPrompt]->fill(xp);
        }
        else if (_sqs == 13000) {
          _h_pT_J[_sqs + nonPrompt]->fill(xp);
          _h_y_J[_sqs + nonPrompt]->fill(absrap);
        }
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      // 1/2 due rapidity folding +/-
      const double factor = 0.5 * crossSection() / nanobarn / sumOfWeights();
      scale(_h_psi, factor);
      divByGroupWidth(_h_psi);
      Estimate1DPtr tmp;
      for (double eVal : allowedEnergies()) {
        const int en(eVal);
        int ih = en == 7000;
        for (const auto& b : _h_psi[en + 1]->bins()) {
          book(tmp, 9 + ih, 1, b.index());
          divide(b, _h_psi[en + 2]->bin(b.index()), tmp);
          tmp->scale(100.);
        }
      }
      scale(_h_pT, factor);
      scale(_h_y, factor);
      scale(_h_pT_J, factor);
      scale(_h_pT_2, factor);
      scale(_h_y_J, factor);
      // ratio psi(2s)/J/psi only at 13 TeV
      size_t ih = 0;
      for (const auto& item : _h_pT_2) {
        book(tmp, 11, 1, ih + 1);
        divide(item.second, _h_pT_J[item.first], tmp);
        book(tmp, 12, 1, ih + 1);
        divide(_h_y[item.first], _h_y_J[item.first], tmp);
        ++ih;
      }
    }

    /// @}


    /// @name Histograms
    /// @{
    map<int, Histo1DGroupPtr> _h_psi;
    map<int, Histo1DPtr> _h_pT, _h_pT_2, _h_y, _h_pT_J, _h_y_J;
    int _sqs = 0;
    /// @}
  };


  RIVET_DECLARE_PLUGIN(LHCB_2019_I1748712);

}