1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief J/psi, and psi(2s) production at 5.02 and 13 TeV
  9  class ALICE_2017_I1511865 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(ALICE_2017_I1511865);
 14
 15
 16    /// @name Analysis methods
 17    ///@{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      declare(UnstableParticles(Cuts::pid==443 || Cuts::pid==100443), "UFS");
 22
 23      for (double eVal : allowedEnergies()) {
 24        const string en = toString(int(eVal/MeV));
 25        if (isCompatibleWithSqrtS(eVal))  _sqs = en;
 26
 27        size_t offset = (en == "5020"s)? 6 : 0;
 28        book(_h[en+"JPsi_pT"], 1+offset,1,1);
 29        book(_h[en+"JPsi_y"],  2+offset,1,1);
 30      }
 31      if (_sqs == "" && !merging()) {
 32        throw BeamError("Invalid beam energy for " + name() + "\n");
 33      }
 34      book(_h["Psi2S_pT"], 3,1,1);
 35      book(_h["Psi2S_y"],  4,1,1);
 36      book(_h["JPsi_pT2"], "TMP/JPsi_pY", refData(5,1,1));
 37      book(_h["JPsi_y2"],  "TMP/JPsi_y",  refData(6,1,1));
 38    }
 39
 40
 41    /// Perform the per-event analysis
 42    void analyze(const Event& event) {
 43      // loop over J/Psi
 44      for (const Particle& p : apply<UnstableParticles>(event, "UFS").particles()) {
 45        // rapidity cut
 46        const double absrap = p.absrap();
 47        if (absrap<2.5 || absrap>4) continue;
 48        const double xp = p.perp();
 49        // J/Psi
 50        if (p.pid()==443) {
 51          if (_sqs != "5020"s) {
 52            if (xp>30.) continue;
 53            _h[_sqs+"JPsi_pT"]->fill(xp);
 54            _h[_sqs+"JPsi_y"]->fill(absrap);
 55            if (xp<=16.) {
 56              _h["JPsi_pT2"]->fill(xp);
 57              _h["JPsi_y2"]->fill(absrap);
 58            }
 59          }
 60          else {
 61            if (xp>12.) continue;
 62            _h[_sqs+"JPsi_pT"]->fill(xp);
 63            _h[_sqs+"JPsi_y"]->fill(absrap);
 64          }
 65        }
 66        // psi(2S)
 67        else if (_sqs != "5020"s) {
 68          if (xp>16.) continue;
 69          _h["Psi2S_pT"]->fill(xp);
 70          _h["Psi2S_y"]->fill(absrap);
 71        }
 72      }
 73    }
 74
 75
 76    /// Normalise histograms etc., after the run
 77    void finalize() {
 78      // factor 1/2 due folding +/- rap
 79      const double fact = 0.5*crossSection()/nanobarn/sumOfWeights();
 80      // factor 1.5 for rapidity range 2.5-4
 81      for (auto& item : _h) {
 82        if (item.first.find("_pT") != string::npos) {
 83          scale(item.second, fact/1.5);
 84        }
 85        else {
 86          scale(item.second, fact);
 87        }
 88      }
 89
 90      Estimate1DPtr tmp;
 91      book(tmp,5,1,1);
 92      divide(_h["Psi2S_pT"], _h["JPsi_pT2"], tmp);
 93      book(tmp,6,1,1);
 94      divide(_h["Psi2S_y"], _h["JPsi_y2"], tmp);
 95    }
 96
 97    ///@}
 98
 99
100    /// @name Histograms
101    ///@{
102    map<string,Histo1DPtr> _h;
103    Histo1DPtr _h_JPsi_pT[2], _h_JPsi_y[2];
104    Histo1DPtr _h_Psi2S_pT,_h_Psi2S_y;
105    Histo1DPtr _h_JPsi_pT2,_h_JPsi_y2;
106
107    string _sqs = "";
108    ///@}
109
110
111  };
112
113
114  RIVET_DECLARE_PLUGIN(ALICE_2017_I1511865);
115
116}