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