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 J/\psi \eta$ at energies between 3.82 and 4.6 GeV
 10  class BESIII_2015_I1355215 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2015_I1355215);
 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(), "UFS");
 24      book(_nJPsi, 1, 1, 10);
 25      for (const string& en : _nJPsi.binning().edges<0>()) {
 26        const double end = std::stod(en)*GeV;
 27        if (isCompatibleWithSqrtS(end)) {
 28          _ecms = en;
 29          break;
 30        }
 31      }
 32      if (_ecms.empty()) MSG_ERROR("Beam energy incompatible with analysis.");
 33    }
 34
 35    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 36      for (const Particle &child : p.children()) {
 37        if(child.children().empty()) {
 38          --nRes[child.pid()];
 39          --ncount;
 40        }
 41        else
 42          findChildren(child,nRes,ncount);
 43      }
 44    }
 45
 46    /// Perform the per-event analysis
 47    void analyze(const Event& event) {
 48      const FinalState& fs = apply<FinalState>(event, "FS");
 49      map<long,int> nCount;
 50      int ntotal(0);
 51      for (const Particle& p : fs.particles()) {
 52        nCount[p.pid()] += 1;
 53        ++ntotal;
 54      }
 55      const FinalState& ufs = apply<FinalState>(event, "UFS");
 56      for (const Particle& p1 : ufs.particles()) {
 57        if(p1.children().empty()) continue;
 58        bool matched=false;
 59        // find the j/psi
 60        if(p1.pid()!=443) continue;
 61        map<long,int> nRes = nCount;
 62        int ncount = ntotal;
 63        findChildren(p1,nRes,ncount);
 64
 65        for (const Particle& p2 : ufs.particles()) {
 66          if(p2.children().empty()) continue;
 67          // find the j/psi
 68          if(p2.pid()!=331) 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(auto const & val : nRes2) {
 75            if(val.second!=0) {
 76              matched = false;
 77              break;
 78            }
 79          }
 80          if(matched) {
 81            _nJPsi->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(_nJPsi,crossSection()/ sumOfWeights() /picobarn);
 93    }
 94
 95    /// @}
 96
 97
 98    /// @name Histograms
 99    /// @{
100    BinnedHistoPtr<string> _nJPsi;
101    string _ecms;
102    /// @}
103
104  };
105
106
107  RIVET_DECLARE_PLUGIN(BESIII_2015_I1355215);
108
109
110}