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 Xi- Xibar+ cross section
 10  class BESIII_2020_I1823448 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2020_I1823448);
 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      book(_nXi, 1, 1, 4);
 26
 27      for (const string& en : _nXi.binning().edges<0>()) {
 28        const double end = std::stod(en)*GeV;
 29        if (isCompatibleWithSqrtS(end)) {
 30          _ecms = en;
 31          break;
 32        }
 33      }
 34      if(_ecms.empty()) MSG_ERROR("Beam energy incompatible with analysis.");
 35    }
 36
 37    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 38      for(const Particle &child : p.children()) {
 39	if(child.children().empty()) {
 40	  nRes[child.pid()]-=1;
 41	  --ncount;
 42	}
 43	else
 44	  findChildren(child,nRes,ncount);
 45      }
 46    }
 47
 48    /// Perform the per-event analysis
 49    void analyze(const Event& event) {
 50      const FinalState& fs = apply<FinalState>(event, "FS");
 51      // total hadronic and muonic cross sections
 52      map<long,int> nCount;
 53      int ntotal(0);
 54      for (const Particle& p : fs.particles()) {
 55	nCount[p.pid()] += 1;
 56	++ntotal;
 57      }
 58      // find the Xis
 59      const FinalState& ufs = apply<UnstableParticles>(event, "UFS");
 60      for(unsigned int ix=0;ix<ufs.particles().size();++ix) {
 61	const Particle& p1 = ufs.particles()[ix];
 62	if(abs(p1.pid())!=3312) continue;
 63	bool matched = false;
 64	// check fs
 65	bool fs = true;
 66	for(const Particle & child : p1.children()) {
 67	  if(child.pid()==p1.pid()) {
 68	    fs = false;
 69	    break;
 70	  }
 71	}
 72	if(!fs) continue;
 73	// find the children
 74	map<long,int> nRes = nCount;
 75	int ncount = ntotal;
 76	findChildren(p1,nRes,ncount);
 77	for(unsigned int iy=ix+1;iy<ufs.particles().size();++iy) {
 78	  const Particle& p2 = ufs.particles()[iy];
 79	  if(abs(p2.pid())!=3312) continue;
 80	  // check fs
 81	  bool fs = true;
 82	  for(const Particle & child : p2.children()) {
 83	    if(child.pid()==p2.pid()) {
 84	      fs = false;
 85	      break;
 86	    }
 87	  }
 88	  if(!fs) continue;
 89	  map<long,int> nRes2 = nRes;
 90	  int ncount2 = ncount;
 91	  findChildren(p2,nRes2,ncount2);
 92	  if(ncount2!=0) continue;
 93	  matched=true;
 94	  for(auto const & val : nRes2) {
 95	    if(val.second!=0) {
 96	      matched = false;
 97	      break;
 98	    }
 99	  }
100	  if(matched) {
101	    _nXi->fill(_ecms);
102	    break;
103	  }
104	}
105	if(matched) break;
106      }
107    }
108
109
110    /// Normalise histograms etc., after the run
111    void finalize() {
112      scale(_nXi, crossSection()/ sumOfWeights() /picobarn);
113    }
114
115    ///@}
116
117
118    /// @name Histograms
119    ///@{
120    BinnedHistoPtr<string> _nXi;
121    string _ecms;
122    ///@}
123
124
125  };
126
127
128  RIVET_DECLARE_PLUGIN(BESIII_2020_I1823448);
129
130}