1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/UnstableParticles.hh"
 4
 5namespace Rivet {
 6
 7
 8  /// @brief J/psi at 5.02 TeV
 9  class CMS_2017_I1512296 : public Analysis {
10  public:
11
12    /// Constructor
13    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2017_I1512296);
14
15
16    /// @name Analysis methods
17    /// @{
18
19    /// Book histograms and initialise projections before the run
20    void init() {
21      declare(UnstableParticles(), "UFS");
22      // histograms
23      for (int ix=0; ix<2; ++ix) {
24        book(_h_JPsi_pT[ix], {-2.4,-1.93,-1.5,-0.9,0.,0.9,1.5,1.93,2.4});
25        for (int iy=0; iy<8; ++iy) {
26          Histo1DPtr tmp;
27          if (iy<4) {
28            book(_h_JPsi_pT[ix]->bin(iy+1), 2+13*ix, 1, 4-iy);
29          }
30          else {
31            book(_h_JPsi_pT[ix]->bin(iy+1), 1+13*ix, 1, 8-iy);
32          }
33        }
34      }
35      for (unsigned int ix=0; ix<2; ++ix) {
36        book(_h_JPsi_y[ix], {6.5,10.,30.});
37        for (auto& b : _h_JPsi_y[ix]->bins()) {
38          book(b, 5+13*ix, 1, b.index());
39        }
40      }
41    }
42
43
44    /// Perform the per-event analysis
45    void analyze(const Event& event) {
46
47      // Final state of unstable particles to get particle spectra
48      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
49
50      for (const Particle& p : ufs.particles(Cuts::pid==443)) {
51        // prompt/non-prompt
52        bool nonPrompt = p.fromBottom();
53        const double rap = p.rap();
54        const double xp = p.perp();
55        _h_JPsi_pT[nonPrompt]->fill(rap,xp);
56        _h_JPsi_y [nonPrompt]->fill(xp,rap);
57      }
58    }
59
60
61    /// Normalise histograms etc., after the run
62    void finalize() {
63      // br to muons PDG 2021
64      const double br = 0.05961;
65      const double factor = br*crossSection()/microbarn/sumOfWeights();
66      const vector<double> pTbins={6.5,10.,30.};
67      for (unsigned int ix=0; ix<2; ++ix) {
68        scale(_h_JPsi_pT[ix], factor);
69        divByGroupWidth(_h_JPsi_pT[ix]);
70        // in y just single differential, don't divide by group width
71        scale(_h_JPsi_y[ix], factor);
72      }
73    }
74
75    /// @}
76
77
78    /// @name Histograms
79    /// @{
80    Histo1DGroupPtr _h_JPsi_pT[2],_h_JPsi_y[2];
81    /// @}
82
83
84  };
85
86
87  RIVET_DECLARE_PLUGIN(CMS_2017_I1512296);
88
89}