1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Upsilon production at 7,8 TeV
  9  class LHCB_2015_I1392456 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(LHCB_2015_I1392456);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      declare(UnstableParticles(), "UFS");
 22      int iloc=-1;
 23      if (isCompatibleWithSqrtS(7000)) {
 24	iloc = 0;
 25      }
 26      else if  (isCompatibleWithSqrtS(8000)) {
 27	iloc = 10;
 28      }
 29      else
 30	throw UserError("Centre-of-mass energy of the given input is neither 7 or 8 TeV.");
 31
 32      for(unsigned int iups=0;iups<3;++iups) {
 33        book(_h_Ups[iups],{2.0,2.5,3.0,3.5,4.0,4.5});
 34	for(unsigned int iy=0;iy<5;++iy) {
 35          book(_h_Ups[iups]->bin(iy+1),1+iups+iloc,1,1+iy);
 36	}
 37	book(_h_Ups_pT[iups],4+iloc,1,1+iups);
 38	book(_h_Ups_y [iups],5+iloc,1,1+iups);
 39      }
 40    }
 41
 42
 43    /// Perform the per-event analysis
 44    void analyze(const Event& event) {
 45      // Final state of unstable particles to get particle spectra
 46      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 47      for (const Particle& p : ufs.particles(Cuts::pid==553 or Cuts::pid==100553 or Cuts::pid==200553)) {
 48        double absrap = p.absrap();
 49        double xp = p.perp();
 50	if(absrap<2. || absrap>4.5) continue;
 51	unsigned int iups=p.pid()/100000;
 52	_h_Ups_pT[iups]  ->fill(xp);
 53	if(xp<30.) _h_Ups_y[iups]   ->fill(absrap);
 54	_h_Ups[iups]     ->fill(absrap,xp);	
 55      }
 56    }
 57
 58
 59    /// Normalise histograms etc., after the run
 60    void finalize() {
 61      // 1/2 due rapidity folding +/-
 62      double factor = 0.5*crossSection()/picobarn/sumOfWeights();
 63      // branching ratios for upsilon
 64      vector<double> br = {0.0248,0.0193,0.0218};
 65      for(unsigned int iups=0;iups<3;++iups) {
 66	scale(_h_Ups_pT  [iups],factor*br[iups]);
 67	scale(_h_Ups_y   [iups],factor/1000.*br[iups]);
 68	scale(_h_Ups[iups],factor*br[iups]);
 69        divByGroupWidth(_h_Ups[iups]);
 70      }
 71      unsigned int iloc = isCompatibleWithSqrtS(8000) ? 10 : 0;
 72      for(unsigned int ix=0;ix<_h_Ups[0]->numBins();++ix) {
 73	Estimate1DPtr tmp;
 74	for(unsigned int iups=1;iups<3;++iups) {
 75	  book(tmp,6+iups+iloc,1,1+ix);
 76	  divide(_h_Ups[iups]->bin(ix+1),_h_Ups[0]->bin(ix+1),tmp);
 77	}
 78      }
 79      for(unsigned int iups=1;iups<3;++iups) {
 80	Estimate1DPtr tmp;
 81	book(tmp,9+iloc,1,iups);
 82	divide(_h_Ups_pT[iups],_h_Ups_pT[0],tmp);
 83	book(tmp,10+iloc,1,iups);
 84	divide(_h_Ups_y[iups],_h_Ups_y[0],tmp);
 85      }
 86      Estimate1DPtr tmp;
 87      book(tmp,9+iloc,1,3);
 88      divide(_h_Ups_pT[2],_h_Ups_pT[1],tmp);
 89      book(tmp,10+iloc,1,3);
 90      divide(_h_Ups_y[2],_h_Ups_y[1],tmp);
 91    }
 92
 93    /// @}
 94
 95
 96    /// @name Histograms
 97    /// @{
 98    Histo1DPtr _h_Ups_pT[3], _h_Ups_y[3];
 99    Histo1DGroupPtr _h_Ups[3];
100    /// @}
101
102
103  };
104
105
106  RIVET_DECLARE_PLUGIN(LHCB_2015_I1392456);
107
108}