Rivet analyses

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

namespace Rivet {


  /// @brief J/psi production at 13 TeV (central rapidity)
  class ALICE_2021_I1898832 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(ALICE_2021_I1898832);


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

    /// Book histograms and initialise projections before the run
    void init() {
      declare(UnstableParticles(Cuts::pid==443), "UFS");
      book(_h_pT,1,1,1);
      book(_h_y1,2,1,1);
      book(_h_y2,3,1,1);
      _axis = YODA::Axis<double>({ 0., 0.5, 1., 2., 3., 4., 5., 7., 10., 15., 20., 25., 30., 40.});
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      if (_edges.empty())  _edges = _h_pT->xEdges();
      // loop over J/Psi
      for (const Particle& p : apply<UnstableParticles>(event, "UFS").particles()) {
        const double absrap = p.absrap();
        if (absrap>0.9) continue;

        _h_pT->fill(disc(p.perp()));

        _h_y1->fill("0.0"s);

        string yedge("0.7");
        if (absrap < 0.2) yedge = "0.1";
        else if (absrap < 0.5) yedge = "0.3";

        _h_y2->fill(yedge);
      }
    }

    string disc(const double value) const {
      string edge("OTHER");
      const size_t idx = _axis.index(value);
      if (0 < idx && idx <= _axis.numBins())  return _edges[idx-1];
      return edge;
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      double fact = crossSection()/microbarn/sumOfWeights();
      // 1.8 for the rapidity -0.9 to 0.9 in the double differential
      scale(_h_pT,fact/1.8);
      for(auto & b : _h_pT->bins()) {
        b.scaleW(1./_axis.width(b.index()));
      }
      // divide by range -0.9 to 0.9 as quoted as dsigma/dy
      scale(_h_y1,fact/1.8);
      // 0.5 as fold + and - rapidity
      scale(_h_y2,0.5*fact);
      for(auto & b : _h_y2->bins()) {
        if     (b.index()==1) b.scaleW(1./0.2);
        else if(b.index()==2) b.scaleW(1./0.3);
        else                  b.scaleW(1./0.4);
      }
    }

    ///@}


    /// @name Histograms
    ///@{
    BinnedHistoPtr<string> _h_pT,_h_y1,_h_y2;
    vector<string> _edges;
    YODA::Axis<double> _axis;
    ///@}


  };


  RIVET_DECLARE_PLUGIN(ALICE_2021_I1898832);

}