1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/FinalState.hh"
 4#include "Rivet/Projections/UnstableParticles.hh"
 5
 6namespace Rivet {
 7
 8
 9  /// @brief J/psi production at 13 TeV (central rapidity)
10  class ALICE_2021_I1898832 : public Analysis {
11  public:
12
13    /// Constructor
14    RIVET_DEFAULT_ANALYSIS_CTOR(ALICE_2021_I1898832);
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), "UFS");
23      book(_h_pT,1,1,1);
24      book(_h_y1,2,1,1);
25      book(_h_y2,3,1,1);
26      _axis = YODA::Axis<double>({ 0., 0.5, 1., 2., 3., 4., 5., 7., 10., 15., 20., 25., 30., 40.});
27    }
28
29
30    /// Perform the per-event analysis
31    void analyze(const Event& event) {
32      if (_edges.empty())  _edges = _h_pT->xEdges();
33      // loop over J/Psi
34      for (const Particle& p : apply<UnstableParticles>(event, "UFS").particles()) {
35        const double absrap = p.absrap();
36        if (absrap>0.9) continue;
37
38        _h_pT->fill(disc(p.perp()));
39
40        _h_y1->fill("0.0"s);
41
42        string yedge("0.7");
43        if (absrap < 0.2) yedge = "0.1";
44        else if (absrap < 0.5) yedge = "0.3";
45
46        _h_y2->fill(yedge);
47      }
48    }
49
50    string disc(const double value) const {
51      string edge("OTHER");
52      const size_t idx = _axis.index(value);
53      if (0 < idx && idx <= _axis.numBins())  return _edges[idx-1];
54      return edge;
55    }
56
57
58    /// Normalise histograms etc., after the run
59    void finalize() {
60      double fact = crossSection()/microbarn/sumOfWeights();
61      // 1.8 for the rapidity -0.9 to 0.9 in the double differential
62      scale(_h_pT,fact/1.8);
63      for(auto & b : _h_pT->bins()) {
64        b.scaleW(1./_axis.width(b.index()));
65      }
66      // divide by range -0.9 to 0.9 as quoted as dsigma/dy
67      scale(_h_y1,fact/1.8);
68      // 0.5 as fold + and - rapidity
69      scale(_h_y2,0.5*fact);
70      for(auto & b : _h_y2->bins()) {
71        if     (b.index()==1) b.scaleW(1./0.2);
72        else if(b.index()==2) b.scaleW(1./0.3);
73        else                  b.scaleW(1./0.4);
74      }
75    }
76
77    ///@}
78
79
80    /// @name Histograms
81    ///@{
82    BinnedHistoPtr<string> _h_pT,_h_y1,_h_y2;
83    vector<string> _edges;
84    YODA::Axis<double> _axis;
85    ///@}
86
87
88  };
89
90
91  RIVET_DECLARE_PLUGIN(ALICE_2021_I1898832);
92
93}