1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Quarkonium prouction at 13 TeV
  9  class CMS_2018_I1633431 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2018_I1633431);
 14
 15
 16    /// @name Analysis methods
 17    ///@{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      // projection
 22      declare(UnstableParticles(), "UFS");
 23      // J/psi/psi(2s) histos
 24      for (unsigned int ix=1; ix<3; ++ix) {
 25        book(_h_onium[ix-1], {0.,0.3,0.6,0.9,1.2});
 26        for (unsigned int iy=1; iy<5; ++iy) {
 27          Histo1DPtr tmp;
 28          book(_h_onium[ix-1]->bin(iy), ix, 1, iy);
 29        }
 30        book(_h_total[ix-1], ix, 1, 5);
 31      }
 32      // Upsilon
 33      for (unsigned int ix=3; ix<6; ++ix) {
 34        book(_h_onium[ix-1], {0.,0.6,1.2});
 35        for (unsigned int iy=1; iy<3; ++iy) {
 36          book(_h_onium[ix-1]->bin(iy), ix, 1, iy);
 37        }
 38        book(_h_total[ix-1], ix, 1, 3);
 39      }
 40      // histos for ratios
 41      for (unsigned int ix=0; ix<5; ++ix) {
 42        if (ix<2) {
 43          book(_h_ratio[ix],"TMP/h_ratio_"+to_str(ix), refData(11,1,1));
 44        }
 45        else {
 46          book(_h_ratio[ix],"TMP/h_ratio_"+to_str(ix), refData(12,1,1));
 47        }
 48      }
 49    }
 50
 51
 52    /// Perform the per-event analysis
 53    void analyze(const Event& event) {
 54      // Final state of unstable particles to get particle spectra
 55      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 56      // loop over onium states
 57      for (const Particle& p : ufs.particles(Cuts::pid==553 or Cuts::pid==100553 or
 58                                              Cuts::pid==200553 or
 59                                              Cuts::pid==443 or Cuts::pid==100443)) {
 60        // cuts on pT and rapidity
 61        const double y = p.absrap();
 62        const double pT = p.perp();
 63        if (y>1.2 || pT<20.) continue;
 64        unsigned int itype(0);
 65        if (p.pid()==443)         itype=0;
 66        else if (p.pid()==100443) itype=1;
 67        else if (p.pid()==553)    itype=2;
 68        else if (p.pid()==100553) itype=3;
 69        else if (p.pid()==200553) itype=4;
 70        // for J/pis and psi(2s) check if prompt
 71        if (itype<=1 && p.fromBottom()) continue;
 72        _h_onium[itype]->fill(y,pT);
 73        _h_total[itype]->fill(pT);
 74        _h_ratio[itype]->fill(pT);
 75      }
 76    }
 77
 78
 79    /// Normalise histograms etc., after the run
 80    void finalize() {
 81      // 0.5 due folded rapidity
 82      const double factor = 0.5*crossSection() / picobarn/ sumOfWeights();
 83      // branching ratios to muons
 84      const vector<double> brs={0.05961,0.00793,0.0248,0.0191,0.0218};
 85      for (unsigned int ix=0; ix<5; ++ix) {
 86        scale(_h_onium[ix], factor * brs[ix]);
 87        divByGroupWidth(_h_onium[ix]);
 88        scale(_h_total[ix],brs[ix]*factor/1.2);
 89        scale(_h_ratio[ix],brs[ix]*factor/1.2);
 90      }
 91      // ratio plots
 92      Estimate1DPtr tmp;
 93      book(tmp, 11, 1, 1);
 94      divide(_h_ratio[1], _h_ratio[0], tmp);
 95      for (unsigned int ix=0; ix<2; ++ix) {
 96        book(tmp,12,1,1+ix);
 97        divide(_h_ratio[3+ix], _h_ratio[2], tmp);
 98      }
 99    }
100
101    ///@}
102
103
104    /// @name Histograms
105    ///@{
106    Histo1DGroupPtr _h_onium[5];
107    Histo1DPtr _h_total[5];
108    Histo1DPtr _h_ratio[5];
109    ///@}
110
111
112  };
113
114
115  RIVET_DECLARE_PLUGIN(CMS_2018_I1633431);
116
117}