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