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

}