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 Cross section for $e^+e^-\to K^0_SK^\pm\pi^\mp\pi^0$ and $K^0_SK^\pm\pi^\mp\eta$ between 3.90 to 4.60 GeV
 10  class BESIII_2018_I1699641 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2018_I1699641);
 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
 26      // Book histograms
 27      book(_cKKpipi, "TMP/2Kpipi" );
 28      book(_cKKpieta, "TMP/2Kpieta");
 29    }
 30
 31
 32    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 33      for (const Particle &child : p.children()) {
 34        if(child.children().empty()) {
 35          nRes[child.pid()]-=1;
 36          --ncount;
 37        }
 38        else
 39          findChildren(child,nRes,ncount);
 40      }
 41    }
 42
 43    /// Perform the per-event analysis
 44    void analyze(const Event& event) {
 45      const FinalState& fs = apply<FinalState>(event, "FS");
 46
 47      map<long,int> nCount;
 48      int ntotal(0);
 49      for (const Particle& p : fs.particles()) {
 50        nCount[p.pid()] += 1;
 51        ++ntotal;
 52      }
 53      // K K pi pi
 54      if(ntotal==4 && nCount[310]==1 && nCount[111]==1 &&
 55         ((nCount[ 321]==1 &&nCount[-211]==1) ||
 56          (nCount[-321]==1 &&nCount[ 211]==1) ))
 57        _cKKpipi->fill();
 58      // eta resonance
 59      const FinalState& ufs = apply<FinalState>(event, "UFS");
 60      for (const Particle& p : ufs.particles()) {
 61        if(p.children().empty()) continue;
 62        if(p.pid()!=221) continue;
 63        map<long,int> nRes=nCount;
 64        int ncount = ntotal;
 65        findChildren(p,nRes,ncount);
 66        if(ncount!=3) continue;
 67        bool matched=true;
 68        for(auto const & val : nRes) {
 69          if(abs(val.first)==321 || abs(val.first)==211) {
 70            continue;
 71          }
 72          else if(abs(val.first)==310) {
 73            if(val.second!=1) {
 74              matched = false;
 75              break;
 76            }
 77          }
 78          else if(val.second!=0) {
 79            matched = false;
 80            break;
 81          }
 82        }
 83        if(matched==false) continue;
 84        if((nCount[ 321] == 1 && nCount[-211] ==1) ||
 85           (nCount[-321] == 1 && nCount[ 211] ==1))
 86          _cKKpieta->fill();
 87      }
 88    }
 89
 90
 91    /// Normalise histograms etc., after the run
 92    void finalize() {
 93      for (unsigned int ix=1;ix<3;++ix) {
 94        double sigma = 0., error = 0.;
 95        if(ix==1) {
 96          sigma = _cKKpipi->val();
 97          error = _cKKpipi->err();
 98        }
 99        else if(ix==2) {
100          sigma = _cKKpieta->val();
101          error = _cKKpieta->err();
102        }
103        sigma *= crossSection()/ sumOfWeights() /picobarn;
104        error *= crossSection()/ sumOfWeights() /picobarn;
105        Estimate1DPtr  mult;
106        book(mult, ix, 1, 1);
107        for (auto& b : mult->bins()) {
108          if (inRange(sqrtS()/GeV, b.xMin(), b.xMax())) {
109            b.set(sigma, error);
110          }
111        }
112      }
113    }
114
115    /// @}
116
117
118    /// @name Histograms
119    /// @{
120    CounterPtr _cKKpipi,_cKKpieta;
121    /// @}
122
123  };
124
125
126  RIVET_DECLARE_PLUGIN(BESIII_2018_I1699641);
127
128
129}