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      if (isCompatibleWithSqrtS(13000)) {
 23        book(_h_JPsi_pT,1,1,1);
 24        book(_h_JPsi_y,2,1,1);
 25        book(_h_Psi2S_pT,3,1,1);
 26        book(_h_Psi2S_y ,4,1,1);
 27        book(_h_JPsi_pT2,"TMP/JPsi_pY",refData(5,1,1));
 28        book(_h_JPsi_y2 ,"TMP/JPsi_y", refData(6,1,1));
 29      }
 30      else if(isCompatibleWithSqrtS(5020) ) {
 31      	book(_h_JPsi_pT,7,1,1);
 32      	book(_h_JPsi_y,8,1,1);
 33      }
 34      else
 35        throw UserError("Centre-of-mass energy of the given input is neither 5020 nor 13000 GeV.");
 36    }
 37
 38
 39    /// Perform the per-event analysis
 40    void analyze(const Event& event) {
 41      // loop over J/Psi
 42      for (const Particle& p : apply<UnstableParticles>(event, "UFS").particles()) {
 43        // rapidity cut
 44        double absrap = p.absrap();
 45        if (absrap<2.5 || absrap>4) continue;
 46        double xp = p.perp();
 47        // J/Psi
 48        if (p.pid()==443) {
 49          if(_h_JPsi_pT2) {
 50            if (xp>30.) continue;
 51            _h_JPsi_pT->fill(xp);
 52            _h_JPsi_y->fill(absrap);
 53            if (xp<=16.) {
 54            _h_JPsi_pT2->fill(xp);
 55            _h_JPsi_y2->fill(absrap);
 56            }
 57          }
 58          else {
 59            if (xp>12.) continue;
 60            _h_JPsi_pT->fill(xp);
 61            _h_JPsi_y->fill(absrap);
 62          }
 63        }
 64        // psi(2S)
 65        else if(_h_Psi2S_pT) {
 66          if (xp>16.) continue;
 67          _h_Psi2S_pT->fill(xp);
 68          _h_Psi2S_y->fill(absrap);
 69        }
 70      }
 71    }
 72
 73
 74    /// Normalise histograms etc., after the run
 75    void finalize() {
 76      // factor 1/2 due folding +/- rap
 77      double fact = 0.5*crossSection()/nanobarn/sumOfWeights();
 78      // factor 1.5 for rapidity range 2.5-4
 79      scale(_h_JPsi_pT,fact/1.5);
 80      scale(_h_JPsi_y ,fact);
 81      if (_h_Psi2S_pT) {
 82        scale(_h_Psi2S_pT,fact/1.5);
 83        scale(_h_Psi2S_y ,fact);
 84      }
 85      if (_h_JPsi_pT2) {
 86        scale(_h_JPsi_pT2,fact/1.5);
 87        Estimate1DPtr tmp;
 88        book(tmp,5,1,1);
 89        divide(_h_Psi2S_pT,_h_JPsi_pT2,tmp);
 90      }
 91      if (_h_JPsi_y2) {
 92        scale(_h_JPsi_y2 ,fact);
 93        Estimate1DPtr tmp;
 94        book(tmp,6,1,1);
 95        divide(_h_Psi2S_y,_h_JPsi_y2,tmp);
 96      }
 97    }
 98
 99    ///@}
100
101
102    /// @name Histograms
103    ///@{
104    Histo1DPtr _h_JPsi_pT,_h_JPsi_y,_h_Psi2S_pT,_h_Psi2S_y;
105    Histo1DPtr _h_JPsi_pT2,_h_JPsi_y2;
106    ///@}
107
108
109  };
110
111
112  RIVET_DECLARE_PLUGIN(ALICE_2017_I1511865);
113
114}