Rivet analyses

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

namespace Rivet {


  /// @brief J/psi, and psi(2s) production at 5.02 and 13 TeV
  class ALICE_2017_I1511865 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(ALICE_2017_I1511865);


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

    /// Book histograms and initialise projections before the run
    void init() {
      declare(UnstableParticles(Cuts::pid == 443 || Cuts::pid == 100443), "UFS");

      for (double eVal : allowedEnergies()) {
        const string en = toString(round(eVal));
        if (isCompatibleWithSqrtS(eVal)) _sqs = en;

        size_t offset = (en == "5020"s) ? 6 : 0;
        book(_h[en + "JPsi_pT"], 1 + offset, 1, 1);
        book(_h[en + "JPsi_y"], 2 + offset, 1, 1);
      }
      raiseBeamErrorIf(_sqs.empty());

      book(_h["Psi2S_pT"], 3, 1, 1);
      book(_h["Psi2S_y"], 4, 1, 1);
      book(_h["JPsi_pT2"], "TMP/JPsi_pY", refData(5, 1, 1));
      book(_h["JPsi_y2"], "TMP/JPsi_y", refData(6, 1, 1));
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      // loop over J/Psi
      for (const Particle& p : apply<UnstableParticles>(event, "UFS").particles()) {
        // rapidity cut
        const double absrap = p.absrap();
        if (absrap < 2.5 || absrap > 4) continue;
        const double xp = p.perp();
        // J/Psi
        if (p.pid() == 443) {
          if (_sqs != "5020"s) {
            if (xp > 30.) continue;
            _h[_sqs + "JPsi_pT"]->fill(xp);
            _h[_sqs + "JPsi_y"]->fill(absrap);
            if (xp <= 16.) {
              _h["JPsi_pT2"]->fill(xp);
              _h["JPsi_y2"]->fill(absrap);
            }
          }
          else {
            if (xp > 12.) continue;
            _h[_sqs + "JPsi_pT"]->fill(xp);
            _h[_sqs + "JPsi_y"]->fill(absrap);
          }
        }
        // psi(2S)
        else if (_sqs != "5020"s) {
          if (xp > 16.) continue;
          _h["Psi2S_pT"]->fill(xp);
          _h["Psi2S_y"]->fill(absrap);
        }
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      // factor 1/2 due folding +/- rap
      const double fact = 0.5 * crossSection() / nanobarn / sumOfWeights();
      // factor 1.5 for rapidity range 2.5-4
      for (auto& item : _h) {
        if (item.first.find("_pT") != string::npos) {
          scale(item.second, fact / 1.5);
        }
        else {
          scale(item.second, fact);
        }
      }

      Estimate1DPtr tmp;
      book(tmp, 5, 1, 1);
      divide(_h["Psi2S_pT"], _h["JPsi_pT2"], tmp);
      book(tmp, 6, 1, 1);
      divide(_h["Psi2S_y"], _h["JPsi_y2"], tmp);
    }

    ///@}


    /// @name Histograms
    ///@{
    map<string, Histo1DPtr> _h;
    Histo1DPtr _h_JPsi_pT[2], _h_JPsi_y[2];
    Histo1DPtr _h_Psi2S_pT, _h_Psi2S_y;
    Histo1DPtr _h_JPsi_pT2, _h_JPsi_y2;

    string _sqs = "";
    ///@}
  };


  RIVET_DECLARE_PLUGIN(ALICE_2017_I1511865);

}