1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4#include "Rivet/Tools/HistoGroup.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief onium production at 5.02 TeV
 10  class ATLAS_2018_I1622737 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2018_I1622737);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      declare(UnstableParticles(Cuts::pid==443 || Cuts::pid==100443 ||
 23                                Cuts::pid==553 || Cuts::pid==100553 ||
 24                                Cuts::pid==200553), "UFS");
 25
 26      // binning in y
 27      const vector<double> ybins={0.,0.75,1.5,2.0};
 28      // book histos
 29      for (unsigned int ix=0; ix<7; ++ix) {
 30        book(_h_Onium[ix], ybins);
 31        for (auto& b : _h_Onium[ix]->bins()) {
 32          book(b, ix+1, 1, b.index());
 33        }
 34        if (ix == 0 || ix == 2) {
 35          _axes[ix] = YODA::Axis<double>({ 8.,9.,10.,11.,12.,13.,14.,16.,18.,22.,30.,40.});
 36        }
 37        else if (ix == 1 || ix == 3) {
 38          _axes[ix] = YODA::Axis<double>({ 8.,10.,12.,16.,22.,40.});
 39        }
 40        else if (ix == 4) {
 41          _axes[ix] = YODA::Axis<double>({ 0.,2.,4.,6.,8.,10.,14.,20.,40.});
 42        }
 43        else {
 44          _axes[ix] = YODA::Axis<double>({ 0.,6.,10.,15.,40.});
 45        }
 46      }
 47    }
 48
 49
 50    /// Perform the per-event analysis
 51    void analyze(const Event& event) {
 52
 53      if (_edges[0].empty()) {
 54        for (unsigned int ix=0; ix<7; ++ix) {
 55          _edges[ix] = _h_Onium[ix]->bin(1)->xEdges();
 56        }
 57      }
 58
 59      // Final state of unstable particles to get particle spectra
 60      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 61
 62      for (const Particle& p : ufs.particles()) {
 63        const double absrap = p.absrap();
 64        const double xp = p.perp();
 65        if (p.pid()==443 || p.pid()==100443) {
 66          // prompt/non-prompt
 67          const bool prompt = !p.fromBottom();
 68          if (p.pid()==443) {
 69            _h_Onium[2*prompt]->fill(absrap, disc(xp, 2*prompt));
 70          }
 71          else {
 72            _h_Onium[1+2*prompt]->fill(absrap, disc(xp, 1+2*prompt));
 73          }
 74        }
 75        else if(p.pid()==553)     _h_Onium[4]->fill(absrap,disc(xp,4));
 76        else if(p.pid()==100553)  _h_Onium[5]->fill(absrap,disc(xp,5));
 77        else if(p.pid()==200553)  _h_Onium[6]->fill(absrap,disc(xp,6));
 78      }
 79    }
 80
 81    string disc(const double value, size_t ix) const {
 82      size_t idx = _axes[ix].index(value);
 83      if (0 < idx && idx <= _axes[ix].numBins()) return _edges[ix][idx-1];
 84      return "OTHER"s;
 85    }
 86
 87
 88    /// Normalise histograms etc., after the run
 89    void finalize() {
 90      // 0.5 due rapidity folding +/-
 91      const double factor = 0.5*crossSection()/nanobarn/sumOfWeights();
 92      // br to muons PDG 2021 (psi2s is e+e- due large errors on mu+mu-)
 93      const vector<double> br = {0.05961,0.00793,0.05961,0.00793,0.0248,0.0193,0.0218};
 94      for (unsigned int ix=0; ix<7; ++ix) {
 95        scale(_h_Onium[ix], factor*br[ix]);
 96        for(auto & hist : _h_Onium[ix]->bins()) {
 97          for(auto & b : hist->bins()) {
 98            b.scaleW(1./_axes[ix].width(b.index()));
 99          }
100        }
101        divByGroupWidth(_h_Onium[ix]);
102      }
103    }
104
105    /// @}
106
107
108    /// @name Histograms
109    /// @{
110    HistoGroupPtr<double,string> _h_Onium[7];
111    YODA::Axis<double> _axes[7];
112    vector<string> _edges[7];
113    /// @}
114
115
116  };
117
118
119  RIVET_DECLARE_PLUGIN(ATLAS_2018_I1622737);
120
121}