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- -> p KS0 nbar K-
 10  class BESIII_2018_I1681638 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2018_I1681638);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      declare(FinalState(), "FS");
 23      declare(UnstableParticles(Cuts::pid==102134), "UFS");
 24      // counters
 25      for (unsigned int ix=0; ix<3; ++ix) {
 26        book(_sigma[ix], 1, 1, 1+ix);
 27      }
 28
 29      for (const string& en : _sigma[0].binning().edges<0>()) {
 30        const double end = std::stod(en)*GeV;
 31        if (isCompatibleWithSqrtS(end)) {
 32          _ecms = en;
 33          break;
 34        }
 35      }
 36      if(_ecms.empty()) MSG_ERROR("Beam energy incompatible with analysis.");
 37
 38    }
 39
 40    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 41      for (const Particle &child : p.children()) {
 42        if (child.children().empty()) {
 43          --nRes[child.pid()];
 44          --ncount;
 45        }
 46        else {
 47          findChildren(child,nRes,ncount);
 48        }
 49      }
 50    }
 51
 52    /// Perform the per-event analysis
 53    void analyze(const Event& event) {
 54      // find the final-state particles
 55      const FinalState& fs = apply<FinalState>(event, "FS");
 56      map<long,int> nCount;
 57      int ntotal(0);
 58      for (const Particle& p : fs.particles()) {
 59        nCount[p.pid()] += 1;
 60        ++ntotal;
 61      }
 62      // FS
 63      if (ntotal==4&&nCount[310]==1) {
 64        if ((nCount[ 2212]==1 && nCount[-2112]==1 && nCount[-321]==1) ||
 65           (nCount[-2212]==1 && nCount[ 2112]==1 && nCount[ 321]==1)) {
 66          _sigma[0]->fill(_ecms);
 67        }
 68      }
 69      for (const Particle& p : apply<FinalState>(event, "UFS").particles()) {
 70        if (p.children().empty()) continue;
 71        map<long,int> nRes = nCount;
 72        int ncount = ntotal;
 73        findChildren(p,nRes,ncount);
 74        if (ncount!=2) continue;
 75        bool matched = true;
 76        int sign = p.pid()/p.abspid();
 77        for (const auto& val : nRes) {
 78          if (val.first==-sign*2112 || val.first==310) {
 79            if (val.second !=1) {
 80              matched = false;
 81              break;
 82            }
 83          }
 84          else if (val.second!=0) {
 85            matched = false;
 86            break;
 87          }
 88        }
 89        if (matched) {
 90          _sigma[1]->fill(_ecms);
 91          break;
 92        }
 93        matched = true;
 94        for (const auto& val : nRes) {
 95          if(val.first==-sign*2212 ||
 96             val.first==sign*321) {
 97            if (val.second !=1) {
 98              matched = false;
 99              break;
100            }
101          }
102          else if (val.second!=0) {
103            matched = false;
104            break;
105          }
106        }
107        if (matched) {
108          _sigma[2]->fill(_ecms);
109          break;
110        }
111      }
112    }
113
114
115    /// Normalise histograms etc., after the run
116    void finalize() {
117      scale(_sigma, crossSection()/picobarn/sumOfWeights());
118    }
119
120    /// @}
121
122
123    /// @name Histograms
124    /// @{
125    BinnedHistoPtr<string> _sigma[3];
126    string _ecms;
127    /// @}
128
129
130  };
131
132
133  RIVET_DECLARE_PLUGIN(BESIII_2018_I1681638);
134
135}