1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief B(*) bar{B}(*) exclusive cross section
 10  class BELLE_2021_I1859137 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2021_I1859137);
 15
 16
 17    /// @name Analysis methods
 18    ///@{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // Initialise and register projections
 23      declare(FinalState(), "FS");
 24      declare(UnstableParticles(), "UFS");
 25      // histograms
 26      for(unsigned int ix=0;ix<3;++ix)
 27        book(_sigma[ix],1,1,1+ix);
 28      for (const string& en : _sigma[0].binning().edges<0>()) {
 29        const double end = std::stod(en)*GeV;
 30        if (isCompatibleWithSqrtS(end)) {
 31          _ecms = en;
 32          break;
 33        }
 34      }
 35      if (_ecms.empty())
 36        MSG_ERROR("Beam energy incompatible with analysis.");
 37    }
 38
 39    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 40      for(const Particle &child : p.children()) {
 41	if(child.children().empty()) {
 42	  nRes[child.pid()]-=1;
 43	  --ncount;
 44	}
 45	else
 46	  findChildren(child,nRes,ncount);
 47      }
 48    }
 49
 50    /// Perform the per-event analysis
 51    void analyze(const Event& event) {
 52      const FinalState& fs = apply<FinalState>(event, "FS");
 53
 54      map<long,int> nCount;
 55      int ntotal(0);
 56      for (const Particle& p : fs.particles()) {
 57	nCount[p.pid()] += 1;
 58	++ntotal;
 59      }
 60      // extract botton hadrons
 61      Particles bHadrons=apply<FinalState>(event, "UFS").particles(Cuts::abspid==511 or Cuts::abspid==513 or
 62								   Cuts::abspid==521 or Cuts::abspid==523);
 63      for(unsigned int ix=0;ix<bHadrons.size();++ix) {
 64	long pix = bHadrons[ix].parents()[0].abspid();
 65	if(pix==511 || pix==413 || pix==521 || pix==523) continue;
 66	map<long,int> nRes = nCount;
 67	int ncount = ntotal;
 68	findChildren(bHadrons[ix],nRes,ncount);
 69	bool matched=false;
 70	for(unsigned int iy=ix+1;iy<bHadrons.size();++iy) {
 71	  long piy = bHadrons[ix].parents()[0].abspid();
 72	  if(piy==511 || piy==413 || piy==521 || piy==523) continue;
 73	  map<long,int> nRes2 = nRes;
 74	  int ncount2 = ncount;
 75	  findChildren(bHadrons[iy],nRes2,ncount2);
 76	  if(ncount2!=0) continue;
 77	  matched=true;
 78	  for(auto const & val : nRes2) {
 79	    if(val.second!=0) {
 80	      matched = false;
 81	      break;
 82	    }
 83	  }
 84	  if(matched) {
 85	    if(bHadrons[ix].abspid()==511 ||
 86	       bHadrons[ix].abspid()==521) {
 87	      if(bHadrons[iy].pid()==-bHadrons[ix].pid())
 88		_sigma[0]->fill(_ecms);
 89	      else
 90		_sigma[1]->fill(_ecms);
 91	    }
 92	    else if(bHadrons[ix].abspid()==513 ||
 93		    bHadrons[ix].abspid()==523) {
 94	      if(bHadrons[iy].pid()==-bHadrons[ix].pid())
 95		_sigma[2]->fill(_ecms);
 96	      else
 97		_sigma[1]->fill(_ecms);
 98	    }
 99	    break;
100	  }
101	}
102	if(matched) break;
103      }
104    }
105
106    /// Normalise histograms etc., after the run
107    void finalize() {
108      double fact = crossSection()/ sumOfWeights() /picobarn;
109      for(unsigned int ix=0;ix<3;++ix)
110        scale(_sigma[ix],fact);
111    }
112
113    ///@}
114
115
116    /// @name Histograms
117    ///@{
118    BinnedHistoPtr<string> _sigma[3];
119    string _ecms;
120    ///@}
121
122
123  };
124
125
126  RIVET_DECLARE_PLUGIN(BELLE_2021_I1859137);
127
128}