1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief psi(2S) production at 7 and 13 TeV
  9  class LHCB_2019_I1748712 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(LHCB_2019_I1748712);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      // projections
 22      declare(UnstableParticles(), "UFS");
 23      // histograms
 24      _is=-1;
 25      if (isCompatibleWithSqrtS(7000)) {
 26      	_is=1;
 27      }
 28      else if (isCompatibleWithSqrtS(13000)) {
 29      	_is=0;
 30      }
 31      else  {
 32      	throw Error("Invalid CMS energy for LHCB_2019_I1748712");
 33      }
 34      for(unsigned int ix=0;ix<3;++ix) {
 35        book(_h_psi[ix],{2.0,2.5,3.0,3.5,4.0,4.5});
 36        for(unsigned int iy=0;iy<5;++iy) {
 37          if(ix<2)
 38            book(_h_psi[ix]->bin(iy+1),1+ix+2*_is,1,iy+1);
 39          else
 40            book(_h_psi[ix]->bin(iy+1),"TMP/psi_"+toString(iy),refData(1+2*_is,1,iy+1));
 41        }
 42      }
 43      for(unsigned int ix=0;ix<2;++ix) {
 44	book(_h_pT[ix],5+2*_is,1,ix+1);
 45	book(_h_y [ix],6+2*_is,1,ix+1);
 46	if(_is==0) {
 47	  book(_h_pT_J[ix],"TMP/Jpsi_pT_"+toString(ix),refData(11,1,ix+1));
 48	  book(_h_pT_2[ix],"TMP/psi_pT_"+toString(ix),refData(11,1,ix+1));
 49	  book(_h_y_J [ix] ,"TMP/Jpsi_y_" +toString(ix),refData(12,1,ix+1));
 50	}
 51      }
 52    }
 53
 54    /// Perform the per-event analysis
 55    void analyze(const Event& event) {
 56      // Final state of unstable particles to get particle spectra
 57      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 58      // J/psi
 59      for (const Particle& p : ufs.particles(Cuts::pid==443 || Cuts::pid==100443)) {
 60	// prompt/non-prompt
 61	bool nonPrompt = p.fromBottom();
 62        double absrap = p.absrap();
 63        double xp = p.perp();
 64	if(absrap<2. || absrap>4.5) continue;
 65	if     (_is==0 && (xp<2   || xp>20.)) continue;
 66	else if(_is==1 && (xp<3.5 || xp>14.)) continue;
 67	if (p.pid()==100443) {
 68	  _h_psi[nonPrompt]->fill(absrap,xp);
 69	  _h_psi[2        ]->fill(absrap,xp);
 70	  _h_pT[nonPrompt]->fill(xp);
 71	  _h_y [nonPrompt]->fill(absrap);
 72	  if(_h_pT_2[0]) _h_pT_2[nonPrompt]->fill(xp);
 73	}
 74	else if (_h_y_J[0]) {
 75	  _h_pT_J[nonPrompt]->fill(xp);
 76	  _h_y_J [nonPrompt]->fill(absrap);
 77	}
 78      }
 79    }
 80
 81
 82    /// Normalise histograms etc., after the run
 83    void finalize() {
 84      // 1/2 due rapidity folding +/-
 85      double factor = 0.5*crossSection()/nanobarn/sumOfWeights();
 86      for(unsigned int ix=0;ix<3;++ix) {
 87	scale(_h_psi[ix],factor);
 88        divByGroupWidth(_h_psi[ix]);
 89      }
 90      Estimate1DPtr tmp;
 91      for(unsigned int ix=0;ix<_h_psi[1]->numBins();++ix) {
 92	book(tmp,9+_is,1,1+ix);
 93	divide(_h_psi[1]->bin(ix+1),_h_psi[2]->bin(ix+1),tmp);
 94	tmp->scale(100.);
 95      }
 96      for(unsigned int ix=0;ix<2;++ix) {
 97	scale(_h_pT[ix],factor);
 98	scale(_h_y [ix],factor);
 99      }
100      // ratio psi(2s)/J/psi only at 13 TeV
101      if(_is==0) {
102	for(unsigned int ix=0;ix<2;++ix) {
103	  scale(_h_pT_J[ix],factor);
104	  scale(_h_pT_2[ix],factor);
105	  scale(_h_y_J [ix],factor);
106	  book(tmp,11,1,ix+1);
107	  divide(_h_pT_2[ix],_h_pT_J[ix],tmp);
108	  book(tmp,12,1,ix+1);
109	  divide(_h_y[ix],_h_y_J[ix],tmp);
110	}
111      }
112    }
113
114    /// @}
115
116
117    /// @name Histograms
118    /// @{
119    Histo1DGroupPtr _h_psi[3];
120    Histo1DPtr _h_pT[2],_h_pT_2[2],_h_y[2];
121    Histo1DPtr _h_pT_J[2],_h_y_J[2];
122    int _is;
123    /// @}
124
125
126  };
127
128
129  RIVET_DECLARE_PLUGIN(LHCB_2019_I1748712);
130
131}