1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Lambda_b and B0 production at 7 and 8 TeV
  9  class LHCB_2016_I1391317 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(LHCB_2016_I1391317);
 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      // CMS energy
 23      int iloc=-1;
 24      if (isCompatibleWithSqrtS(7000)) {
 25	iloc = 0;
 26      }
 27      else if  (isCompatibleWithSqrtS(8000)) {
 28	iloc = 1;
 29      }
 30      else
 31	throw UserError("Centre-of-mass energy of the given input is neither 7 or 8 TeV.");
 32      // histograms
 33      book(_h_Lambda,{2.0,2.5,3.0,3.5,4.0,4.5});
 34      book(_h_B     ,{2.0,2.5,3.0,3.5,4.0,4.5});
 35      for(unsigned int iy=0;iy<5;++iy) {
 36        book(_h_Lambda->bin(iy+1),1+iloc,1,1+iy);
 37	book(_h_B     ->bin(iy+1),3+iloc,1,1+iy);
 38      }
 39      for(unsigned int ix=0;ix<2;++ix) {
 40	book(_h_pT[ix],"TMP/h_pT_"+toString(ix),refData(5,1,1+iloc));
 41	book(_h_y [ix],"TMP/h_y_ "+toString(ix),refData(6,1,1+iloc));
 42      }
 43    }
 44
 45
 46    /// Perform the per-event analysis
 47    void analyze(const Event& event) {
 48      // Final state of unstable particles to get particle spectra
 49      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 50      // loop over onium states
 51      for( const Particle & p : ufs.particles(Cuts::abspid==5122 || Cuts::abspid==511)) {
 52	// skip copies due mixing
 53	if(p.children().size()==1 && p.children()[0].abspid()==p.abspid()) continue;
 54	double rap=p.rapidity();
 55	if(rap<2. || rap>4.5) continue;
 56	double pT = p.perp();
 57	if(p.abspid()==5122) {
 58	  _h_Lambda->fill(rap,pT);
 59	  if(p.pid()>0) {
 60	    _h_pT[0]->fill(pT);
 61	    if(pT<20.) _h_y [0]->fill(rap);
 62	  }
 63	  else {
 64	    _h_pT[1]->fill(pT);
 65	    if(pT<20.) _h_y [1]->fill(rap);
 66	  }
 67	}
 68	else {
 69	  _h_B->fill(rap,pT);
 70	}
 71      }
 72    }
 73
 74
 75    /// Normalise histograms etc., after the run
 76    void finalize() {
 77      // CMS energy
 78      int iloc=-1;
 79      if (isCompatibleWithSqrtS(7000)) {
 80	iloc = 0;
 81      }
 82      else if  (isCompatibleWithSqrtS(8000)) {
 83	iloc = 1;
 84      }
 85      // branching ratios
 86      vector<double> br={3.2e-4,1.27e-3};
 87      // 0.5 particle/antiparticle
 88      double factor = 0.5*crossSection()/picobarn/sumOfWeights();
 89      scale(_h_Lambda,br[0]*factor);
 90      divByGroupWidth(_h_Lambda);
 91      scale(_h_B     ,br[1]*factor);
 92      divByGroupWidth(_h_B);
 93      Estimate1DPtr tmp;
 94      book(tmp,5,1,1+iloc);
 95      asymm(_h_pT[0],_h_pT[1],tmp);
 96      tmp->scale(100.);
 97      book(tmp,6,1,1+iloc);
 98      asymm(_h_y [0],_h_y [1],tmp);
 99      tmp->scale(100.);
100    }
101
102    /// @}
103
104
105    /// @name Histograms
106    /// @{
107    Histo1DGroupPtr _h_Lambda,_h_B;
108    Histo1DPtr _h_pT[2],_h_y[2];
109    /// @}
110
111
112  };
113
114
115  RIVET_DECLARE_PLUGIN(LHCB_2016_I1391317);
116
117}