1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief chi_c at 7 TeV
  9  class LHCB_2013_I1242869 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(LHCB_2013_I1242869);
 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      for(unsigned int ichi=0;ichi<3;++ichi) {
 23	book(_h_chi[ichi],"TMP/h_CHI_"+toString(ichi),refData(1,1,1));
 24	book(_c_chi[ichi],"TMP/c_CHI_"+toString(ichi));
 25      }
 26      book(_r,2,1,1);
 27    }
 28
 29
 30    /// Perform the per-event analysis
 31    void analyze(const Event& event) {
 32      if(_edges.empty()) _edges = _r->xEdges();
 33      // Final state of unstable particles to get particle spectra
 34      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 35
 36      for (const Particle& p : ufs.particles(Cuts::pid==10441 ||
 37					     Cuts::pid==20443 ||
 38					     Cuts::pid==445)) {
 39	// prompt
 40	if(p.fromBottom()) continue;
 41	// J/psi /gamma mode
 42	if(p.children().size()!=2) continue;
 43	Particle Jpsi;
 44	if(p.children()[0].pid()==22 && p.children()[1].pid()==443) {
 45	  Jpsi=p.children()[1];
 46	}
 47	else if(p.children()[1].pid()==22 && p.children()[0].pid()==443) {
 48	  Jpsi=p.children()[0];
 49	}
 50	else
 51	  continue;
 52	double absrap=Jpsi.absrap();
 53	if(absrap<2. || absrap>4.5) continue;
 54	unsigned int ichi = 0;
 55	if(p.pid()==20443) ichi=1;
 56	else if(p.pid()==445) ichi=2;
 57	double xp=Jpsi.perp();
 58	_h_chi[ichi]->fill(xp);
 59	if(xp>4. && xp<20.) {
 60          _c_chi[ichi]->fill();
 61          if     (ichi==0) _r->fill(_edges[0]);
 62          else if(ichi==2) _r->fill(_edges[1]);
 63        }
 64      }
 65    }
 66
 67
 68    /// Normalise histograms etc., after the run
 69    void finalize() {
 70      // chi_c to gamma J/psi branching ratios from PDG 2021
 71      vector<double> br = {0.014,0.343,0.190};
 72      // divide out the branching ratio for mode used to get total rate
 73      for(unsigned int ichi=0;ichi<3;++ichi) {
 74	scale(_h_chi[ichi],1./br[ichi]);
 75	scale(_c_chi[ichi],1./br[ichi]);
 76      }
 77      // ratio chi_c2/chi_c2
 78      Estimate1DPtr tmp;
 79      book(tmp,1,1,1);
 80      divide(_h_chi[2],_h_chi[1],tmp);
 81      _r->bin(1).scaleW(double(1./br[0]/ _c_chi[2]->val()));
 82      _r->bin(2).scaleW(1./br[2]/ _c_chi[1]->val());
 83    }
 84
 85    /// @}
 86
 87
 88    /// @name Histograms
 89    /// @{
 90    Histo1DPtr _h_chi[3];
 91    BinnedHistoPtr<string> _r;
 92    CounterPtr _c_chi[3];
 93    vector<string> _edges;
 94    /// @}
 95
 96
 97  };
 98
 99
100  RIVET_DECLARE_PLUGIN(LHCB_2013_I1242869);
101
102}