1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief chi_b production at 7,8 TeV
  9  class LHCB_2014_I1308738 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(LHCB_2014_I1308738);
 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      unsigned int iloc=0;
 23      if (isCompatibleWithSqrtS(7000)) {
 24	iloc = 1;
 25      }
 26      else if  (isCompatibleWithSqrtS(8000)) {
 27	iloc = 2;
 28      }
 29      else
 30	throw UserError("Centre-of-mass energy of the given input is neither 7 or 8 TeV.");
 31      // decays to 1S
 32      for(unsigned int ix=0;ix<3;++ix)
 33	book(_h_pT_1S[ix],"TMP/Ups_1_"+toString(ix),refData(3*iloc-2,1,1+ix));
 34      for(unsigned int ix=0;ix<2;++ix)
 35	book(_h_pT_1S[ix+3],"TMP/Ups_1_"+toString(ix+3),refData(3*iloc-2,1,2+ix));
 36      // decays to 2S
 37      for(unsigned int ix=0;ix<2;++ix)
 38	for(unsigned int iy=0;iy<2;++iy)
 39	  book(_h_pT_2S[ix+2*iy],"TMP/Ups_2_"+toString(ix+2*iy),refData(3*iloc-1,1,2+ix));
 40      // decays to 3S
 41      for(unsigned int ix=0;ix<2;++ix)
 42	book(_h_pT_3S[ix],"TMP/Ups_3_"+toString(ix),refData(3*iloc,1,3));
 43    }
 44
 45
 46    /// Perform the per-event analysis
 47    void analyze(const Event& event) {
 48
 49      // Final state of unstable particles to get particle spectra
 50      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 51
 52      // first the upsilon states for the denominator
 53      for (const Particle& p : ufs.particles(Cuts::pid==553 ||
 54					     Cuts::pid==100553 ||
 55					     Cuts::pid==200553)) {
 56	// rapidity cut
 57	double absrap=p.absrap();
 58	if(absrap<2. || absrap>4.5) continue;
 59	double xp = p.perp();
 60	if(p.pid()==553) {
 61	  for(unsigned int ix=0;ix<2;++ix) _h_pT_1S[3+ix]->fill(xp);
 62	}
 63	else if (p.pid()==100553) {
 64	  for(unsigned int ix=0;ix<2;++ix) _h_pT_2S[2+ix]->fill(xp);
 65	}
 66	else if (p.pid()==200553) {
 67	  _h_pT_3S[1]->fill(xp);
 68	}
 69      }
 70      // P states
 71      for (const Particle& p : ufs.particles(Cuts::pid==10551 || Cuts::pid==110551 || Cuts::pid==210551 ||
 72      					     Cuts::pid==20553 || Cuts::pid==120553 || Cuts::pid==220553 ||
 73      					     Cuts::pid==  555 || Cuts::pid==100555 || Cuts::pid==200555 )) {
 74	Particle Upsilon;
 75	if     (p.children()[0].pid()==22 && abs(p.children()[1].pid())%100000==553) {
 76	  Upsilon=p.children()[1];
 77	}
 78	else if(p.children()[1].pid()==22 && abs(p.children()[0].pid())%100000==553) {
 79	  Upsilon=p.children()[0];
 80	}
 81	else
 82	  continue;
 83	double absrap=Upsilon.absrap();
 84	if(absrap<2. || absrap>4.5) continue;
 85	unsigned int iups = p.children()[0].pid()/100000;
 86	unsigned int ichi = p.pid()/100000;
 87	double xp  = Upsilon.perp();
 88	if(iups==0)
 89	  _h_pT_1S[ichi]->fill(xp);
 90	else if(iups==1 && ichi>0)
 91	  _h_pT_2S[ichi-1]->fill(xp);
 92	else if(iups==2 && ichi==2)
 93	  _h_pT_3S[0]->fill(xp);
 94      }
 95    }
 96
 97
 98    /// Normalise histograms etc., after the run
 99    void finalize() {
100      unsigned int iloc=0;
101      if (isCompatibleWithSqrtS(7000)) {
102	iloc = 1;
103      }
104      else if  (isCompatibleWithSqrtS(8000)) {
105	iloc = 2;
106      }
107      // calculate the ratios
108      // for 1S
109      for(unsigned int iy=0;iy<3;++iy) {
110	Estimate1DPtr tmp;
111	book(tmp,3*iloc-2,1,iy+1);
112	if(iy<2)
113	  efficiency(_h_pT_1S[iy],_h_pT_1S[3],tmp);
114	else
115	  efficiency(_h_pT_1S[iy],_h_pT_1S[4],tmp);
116	tmp->scale(100.);
117      }
118      Estimate1DPtr tmp;
119      book(tmp,3*iloc-1,1,2);
120      efficiency(_h_pT_2S[0],_h_pT_2S[2],tmp);
121      tmp->scale(100.);
122      book(tmp,3*iloc-1,1,3);
123      efficiency(_h_pT_2S[1],_h_pT_2S[3],tmp);
124      tmp->scale(100.);
125      book(tmp,3*iloc,1,3);
126      efficiency(_h_pT_3S[0],_h_pT_3S[1],tmp);
127      tmp->scale(100.);
128    }
129
130    /// @}
131
132
133    /// @name Histograms
134    /// @{
135    Histo1DPtr _h_pT_1S[5],_h_pT_2S[4],_h_pT_3S[2];
136    /// @}
137
138
139  };
140
141
142  RIVET_DECLARE_PLUGIN(LHCB_2014_I1308738);
143
144}