1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Xi_b- 7,8 and 13 TeV
  9  class LHCB_2019_I1716259 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(LHCB_2019_I1716259);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      // projections
 22      declare(UnstableParticles(), "UFS");
 23      // CMS energy
 24      if (isCompatibleWithSqrtS(7000) ||
 25	  isCompatibleWithSqrtS(8000)) {
 26	_rootS="7-8"s;
 27      }
 28      else if (isCompatibleWithSqrtS(13000)) {
 29	_rootS="13"s;
 30      }
 31      else
 32	throw UserError("Centre-of-mass energy of the given input is neither 7, 8, 13 TeV.");
 33      // histograms
 34      book(_h_xi ,"TMP/h_xi" ,refData<YODA::BinnedEstimate<string>>(1,1,1));
 35      book(_h_lam,"TMP/h_lam",refData<YODA::BinnedEstimate<string>>(1,1,1));
 36      for(unsigned int ix=0;ix<2;++ix) {
 37	book(_c_xi[ix] ,"TMP/c_xi_" +toString(ix+1));
 38	book(_c_lam[ix],"TMP/c_lam_"+toString(ix+1));
 39      }
 40    }
 41    
 42    void findDecayProducts(Particle mother, double sign, Particles & Lambda, Particles & Xi, Particles & Jpsi, unsigned int & nstable) {
 43      for(const Particle & p: mother.children()) {
 44	if(p.pid() == 3122*sign)
 45	  Lambda.push_back(p);
 46	else if(p.pid() == 3312*sign)
 47	  Xi.push_back(p);
 48	else if(p.pid()==443)
 49	  Jpsi.push_back(p);
 50	else if(p.pid()==111 || p.children().empty())
 51	  ++nstable;
 52	else 
 53	  findDecayProducts(p,sign,Lambda,Xi,Jpsi,nstable);
 54      }
 55    }
 56
 57    /// Perform the per-event analysis
 58    void analyze(const Event& event) {
 59      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 60      for (const Particle& p : ufs.particles(Cuts::abspid==5122 || Cuts::abspid==5132)) {
 61	// decay modes
 62	Particles Lambda,Xi,Jpsi;
 63	unsigned int nstable=0;
 64	double sign = p.pid()>0 ? 1. : -1.;
 65	findDecayProducts(p,sign,Lambda,Xi,Jpsi,nstable);
 66	if (p.abspid()==5112) {
 67	  _c_lam[1]->fill();
 68	  if(Lambda.size()==1 && Jpsi.size()==1 && nstable==0)
 69	    _c_lam[0]->fill();
 70	}
 71	else {
 72	  _c_xi[1]->fill();
 73	  if(Xi.size()==1 && Jpsi.size()==1 && nstable==0)
 74	    _c_xi[0]->fill();
 75	}
 76	// pT and rapidity cuts
 77	if(p.perp()>20.) continue;
 78	double eta=p.abseta();
 79	if(eta<2. || eta> 6.) continue;
 80	if (p.abspid()==5122) {
 81	  _h_lam->fill(_rootS);
 82	}
 83	else {
 84	  _h_xi->fill(_rootS);
 85	}
 86      }
 87    }
 88
 89
 90    /// Normalise histograms etc., after the run
 91    void finalize() {
 92      // first the simple ratio
 93      BinnedEstimatePtr<string> tmp;
 94      book(tmp,2,1,1);
 95      divide(_h_xi,_h_lam,tmp);
 96      // and the one with brs included
 97      book(tmp,1,1,1);
 98      if (_c_xi[1]->numEntries()>0.)
 99      	scale(_h_xi, *_c_xi[0]/ *_c_xi[1]);
100      if (_c_lam[1]->numEntries()>0.)
101      	scale(_h_lam, *_c_lam[0]/ *_c_lam[1]);
102      divide(_h_xi,_h_lam,tmp);
103    }
104
105    /// @}
106
107
108    /// @name Histograms
109    /// @{
110    BinnedHistoPtr<string> _h_xi,_h_lam;
111    CounterPtr _c_xi[2],_c_lam[2];
112    string _rootS;
113    /// @}
114
115
116  };
117
118
119  RIVET_DECLARE_PLUGIN(LHCB_2019_I1716259);
120
121}