1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5
  6namespace Rivet {
  7
  8
  9  /// @brief J/psi and psi(2S) at 1.8 TeV
 10  class CDF_1997_I440101 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(CDF_1997_I440101);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      declare(UnstableParticles(), "UFS");
 23      for (unsigned int ix=0; ix<2; ++ix) {
 24        book(_h_total[ix], 1, 1, 1+ix);
 25        for (unsigned int iy=0; iy<2; ++iy) {
 26          book(_h_psi[ix][iy], 2+2*ix+iy, 1, 1);
 27        }
 28      }
 29    }
 30
 31
 32    /// Perform the per-event analysis
 33    void analyze(const Event& event) {
 34
 35      if (_edges.empty()) {
 36        _edges.resize(2);
 37        for(unsigned int ix=0;ix<2;++ix) {
 38          _edges[ix].resize(2);
 39          for(unsigned int iy=0;iy<2;++iy)
 40            _edges[ix][iy] = _h_psi[ix][iy]->xEdges();
 41        }
 42      }
 43
 44      // Final state of unstable particles to get particle spectra
 45      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 46
 47      for (const Particle& p : ufs.particles(Cuts::pid==443 or Cuts::pid==100443)) {
 48        const double abseta = p.abseta();
 49        const double xp = p.perp();
 50          if (xp<5. || abseta>0.6) continue;
 51          bool prompt = !p.fromBottom();
 52          unsigned int ipsi=p.pid()/100000;
 53          _h_total[ipsi]->fill("< 0.6"s);
 54          _h_psi[prompt][ipsi]->fill(disc(xp, prompt, ipsi));
 55      }
 56    }
 57
 58    string disc(const double value, unsigned int prompt, unsigned int ipsi) const {
 59      size_t idx = _axis[ipsi].index(value);
 60      if (0 < idx && idx <= _axis[ipsi].numBins())
 61        return _edges[prompt][ipsi][idx-1];
 62      return "OTHER"s;
 63    }
 64
 65
 66    /// Normalise histograms etc., after the run
 67    void finalize() {
 68      // normalisation factor
 69      const double factor = crossSection()/nanobarn/sumOfWeights();
 70      // br to muons PDG 2021 (psi2s is e+e- due large errors on mu+mu-)
 71      const vector<double> br = {0.05961,0.00793};
 72      for (unsigned int ix=0; ix<2; ++ix) {
 73        scale(_h_total[ix], factor*br[ix]);
 74        for (unsigned int iy=0; iy<2; ++iy) {
 75          scale(_h_psi[ix][iy], factor*br[iy]);
 76          for(auto & b :_h_psi[ix][iy]->bins()) {
 77            b.scaleW(1./_axis[iy].width(b.index()));
 78          }
 79        }
 80      }
 81    }
 82
 83    /// @}
 84
 85
 86    /// @name Histograms
 87    /// @{
 88    BinnedHistoPtr<string> _h_total[2], _h_psi[2][2];
 89    YODA::Axis<double> _axis[2] = { YODA::Axis<double>({5.,5.5,6.,6.5,7.,8.,9.,10.,12.,14.,17.,20.}),
 90                                    YODA::Axis<double>({5.,6.,7.,9.,12.,16.}) };
 91    vector<vector<vector<string>>> _edges;
 92    /// @}
 93
 94
 95  };
 96
 97
 98  RIVET_DECLARE_PLUGIN(CDF_1997_I440101);
 99
100}