Rivet analyses

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

namespace Rivet {


  /// @brief onium production at 5.02 TeV
  class ATLAS_2018_I1622737 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2018_I1622737);


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

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

      // binning in y
      const vector<double> ybins={0.,0.75,1.5,2.0};
      // book histos
      for (unsigned int ix=0; ix<7; ++ix) {
        book(_h_Onium[ix], ybins);
        for (auto& b : _h_Onium[ix]->bins()) {
          book(b, ix+1, 1, b.index());
        }
        if (ix == 0 || ix == 2) {
          _axes[ix] = YODA::Axis<double>({ 8.,9.,10.,11.,12.,13.,14.,16.,18.,22.,30.,40.});
        }
        else if (ix == 1 || ix == 3) {
          _axes[ix] = YODA::Axis<double>({ 8.,10.,12.,16.,22.,40.});
        }
        else if (ix == 4) {
          _axes[ix] = YODA::Axis<double>({ 0.,2.,4.,6.,8.,10.,14.,20.,40.});
        }
        else {
          _axes[ix] = YODA::Axis<double>({ 0.,6.,10.,15.,40.});
        }
      }
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {

      if (_edges[0].empty()) {
        for (unsigned int ix=0; ix<7; ++ix) {
          _edges[ix] = _h_Onium[ix]->bin(1)->xEdges();
        }
      }

      // Final state of unstable particles to get particle spectra
      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");

      for (const Particle& p : ufs.particles()) {
        const double absrap = p.absrap();
        const double xp = p.perp();
        if (p.pid()==443 || p.pid()==100443) {
          // prompt/non-prompt
          const bool prompt = !p.fromBottom();
          if (p.pid()==443) {
            _h_Onium[2*prompt]->fill(absrap, disc(xp, 2*prompt));
          }
          else {
            _h_Onium[1+2*prompt]->fill(absrap, disc(xp, 1+2*prompt));
          }
        }
        else if(p.pid()==553)     _h_Onium[4]->fill(absrap,disc(xp,4));
        else if(p.pid()==100553)  _h_Onium[5]->fill(absrap,disc(xp,5));
        else if(p.pid()==200553)  _h_Onium[6]->fill(absrap,disc(xp,6));
      }
    }

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


    /// Normalise histograms etc., after the run
    void finalize() {
      // 0.5 due rapidity folding +/-
      const double factor = 0.5*crossSection()/nanobarn/sumOfWeights();
      // br to muons PDG 2021 (psi2s is e+e- due large errors on mu+mu-)
      const vector<double> br = {0.05961,0.00793,0.05961,0.00793,0.0248,0.0193,0.0218};
      for (unsigned int ix=0; ix<7; ++ix) {
        scale(_h_Onium[ix], factor*br[ix]);
        divByGroupWidth(_h_Onium[ix]);
      }
    }

    /// @}


    /// @name Histograms
    /// @{
    HistoGroupPtr<double,string> _h_Onium[7];
    YODA::Axis<double> _axes[7];
    vector<string> _edges[7];
    /// @}


  };


  RIVET_DECLARE_PLUGIN(ATLAS_2018_I1622737);

}