1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5
  6namespace Rivet {
  7
  8
  9  /// @brief xi_cc++ production
 10  class LHCB_2020_I1760788 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(LHCB_2020_I1760788);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // projections
 23      declare(UnstableParticles(), "UFS");
 24      // histograms
 25      book(_n_lambda_c,"TMP/nLambdac");
 26      for(unsigned int ix=0;ix<2;++ix)
 27	book(_n_Xicc[ix],"TMP/nXicc_"+toString(ix));
 28      book(_h_xi,1,1,1);
 29    }
 30    
 31    void findDecayProducts(Particle mother, double sign, Particles & lambdac, Particles & Km, Particles & pip, unsigned int & nstable) {
 32      for(const Particle & p: mother.children()) {
 33	if(p.pid()==4122*sign)
 34	  lambdac.push_back(p);
 35	else if(p.pid()==-321*sign)
 36	  Km.push_back(p);
 37	else if(p.pid()==211*sign)
 38	  pip.push_back(p);
 39	else if(p.pid()==111 || p.children().empty())
 40	  ++nstable;
 41	else 
 42	  findDecayProducts(p,sign,lambdac,Km,pip,nstable);
 43      }
 44    }
 45
 46
 47    /// Perform the per-event analysis
 48    void analyze(const Event& event) {
 49      if(_edges.empty()) _edges = _h_xi->xEdges();
 50      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 51      for (const Particle& p : ufs.particles(Cuts::abspid==4122 || Cuts::abspid==4422)) {
 52	if(p.abspid()==4422) {
 53	  Particles lambdac,Km,pip;
 54	  unsigned int nstable=0;
 55	  double sign = p.pid()>0 ? 1. : -1.;
 56	  findDecayProducts(p,sign,lambdac,Km,pip,nstable);
 57	  if(lambdac.size()==1 && Km.size()==1 && pip.size()==2 && nstable ==0)
 58	    _n_Xicc[0]->fill();
 59	  _n_Xicc[1]->fill();
 60	}
 61	// pT and rapidity cuts
 62	double pT=p.perp();
 63	if(pT<4. || pT>15.) continue;
 64	double y=p.absrap();
 65	if(y<2. || y> 4.5) continue;
 66	if (p.abspid()==4422) {
 67	  for(const string & tau : _edges) _h_xi->fill(tau);
 68	}
 69	else
 70	  _n_lambda_c->fill();
 71      }
 72    }
 73
 74
 75    /// Normalise histograms etc., after the run
 76    void finalize() {
 77      // scale by Xi_cc br
 78      if (_n_Xicc[1]->effNumEntries()>0.)
 79      	scale(_h_xi, *_n_Xicc[0]/ *_n_Xicc[1]);
 80      // and by lambda_c rate for normalisation (10^6 as rate in units 10^-4)
 81      if (_n_lambda_c->effNumEntries()>0.)
 82	scale(_h_xi,1e4/ *_n_lambda_c);
 83    }
 84
 85    /// @}
 86
 87
 88    /// @name Histograms
 89    /// @{
 90    BinnedHistoPtr<string> _h_xi;
 91    CounterPtr _n_lambda_c, _n_Xicc[2];
 92    vector<string> _edges;
 93    /// @}
 94
 95
 96  };
 97
 98
 99  RIVET_DECLARE_PLUGIN(LHCB_2020_I1760788);
100
101}