1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/UnstableParticles.hh"
 4
 5namespace Rivet {
 6
 7
 8  /// @brief B_s0 at 5.02 TeV
 9  class CMS_2019_I1697571 : public Analysis {
10  public:
11
12    /// Constructor
13    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2019_I1697571);
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      // histograms
24      book(_h_pT[0],1,1,1);
25      book(_h_pT[1],2,1,1);
26
27      _axes[0] = YODA::Axis<double>({7.,15.,20.,50.});
28      _axes[1] = YODA::Axis<double>({7.,15.,50.});
29    }
30
31
32    /// Perform the per-event analysis
33    void analyze(const Event& event) {
34      if (_edges[0].empty()) {
35        _edges[0] = _h_pT[0]->xEdges();
36        _edges[1] = _h_pT[1]->xEdges();
37      }
38      // Final state of unstable particles to get particle spectra
39      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
40      // loop over onium states
41      for (const Particle& p : ufs.particles(Cuts::abspid==531)) {
42        // skip copies due mixing
43        if (p.children().size()==1 && p.children()[0].abspid()==531) continue;
44        // rapidity and pT
45        const double y = p.absrap();
46        if (y>2.4) continue;
47        const double pT = p.perp();
48        if (pT<7. || pT>50.) continue;
49        for (unsigned int ix=0; ix<2; ++ix) {
50          const size_t idx = _axes[ix].index(pT);
51          string edge("OTHER");
52          if (idx && idx <= _edges[ix].size()) edge = _edges[ix][idx-1];
53          _h_pT[ix]->fill(edge);
54        }
55      }
56    }
57
58
59    /// Normalise histograms etc., after the run
60    void finalize() {
61      // 0.5 from particle/antiparticle
62      double fact = 0.5*crossSection()/picobarn/sumOfWeights();
63      for (unsigned int ix=0; ix<2; ++ix) {
64        scale(_h_pT[ix], fact);
65        for (auto& b : _h_pT[ix]->bins()) {
66          b.scaleW( 1.0/_axes[ix].width(b.index()) );
67        }
68      }
69    }
70
71    /// @}
72
73
74    /// @name Histograms
75    /// @{
76    BinnedHistoPtr<string> _h_pT[2];
77    YODA::Axis<double> _axes[2];
78    vector<string> _edges[2];
79    /// @}
80
81
82  };
83
84
85  RIVET_DECLARE_PLUGIN(CMS_2019_I1697571);
86
87}