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 e+ e- > Omega- Omegabar+
 10  class BESIII_2022_I2611486 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2022_I2611486);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // projections
 23      declare(FinalState(), "FS");
 24      declare(UnstableParticles(Cuts::abspid==3334), "UFS");
 25      // histos
 26      book(_sigma,1,1,1);
 27      for (const string& en : _sigma.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
 49    /// Perform the per-event analysis
 50    void analyze(const Event& event) {
 51      const FinalState& fs = apply<FinalState>(event, "FS");
 52      // total hadronic and muonic cross sections
 53      map<long,int> nCount;
 54      int ntotal(0);
 55      for (const Particle& p : fs.particles()) {
 56      	nCount[p.pid()] += 1;
 57      	++ntotal;
 58      }
 59      // find the Omegas
 60      const FinalState& ufs = apply<UnstableParticles>(event, "UFS");
 61      for (const Particle& p1 : ufs.particles()) {
 62        bool matched = false;
 63        // find the children
 64        map<long,int> nRes = nCount;
 65        int ncount = ntotal;
 66        findChildren(p1,nRes,ncount);
 67        for(const Particle& p2 : ufs.particles()) {
 68          if(p2.pid() != -p1.pid()) continue;
 69          map<long,int> nRes2 = nRes;
 70          int ncount2 = ncount;
 71          findChildren(p2,nRes2,ncount2);
 72          if (ncount2!=0) continue;
 73          matched=true;
 74          for (const auto& val : nRes2) {
 75            if (val.second!=0) {
 76              matched = false;
 77              break;
 78            }
 79          }
 80          if (matched) {
 81            _sigma->fill(_ecms);
 82            break;
 83          }
 84        }
 85        if(matched) break;
 86      }
 87    }
 88
 89
 90    /// Normalise histograms etc., after the run
 91    void finalize() {
 92      scale(_sigma, crossSection()/ sumOfWeights() /femtobarn);
 93    }
 94
 95    /// @}
 96
 97
 98    /// @name Histograms
 99    /// @{
100    BinnedHistoPtr<string> _sigma;
101    string _ecms;
102    /// @}
103
104
105  };
106
107
108  RIVET_DECLARE_PLUGIN(BESIII_2022_I2611486);
109
110}