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 pi0, K* Kbar K2 Kbar
 10  class BELLE_2013_I1252555 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2013_I1252555);
 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      // counters
 26      for (unsigned int ix=0; ix<5; ++ix) {
 27        book(_sigma[ix],1+ix,1,1);
 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    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 40      for (const Particle &child : p.children()) {
 41        if (child.children().empty()) {
 42          nRes[child.pid()]-=1;
 43          --ncount;
 44        }
 45        else {
 46          findChildren(child,nRes,ncount);
 47        }
 48      }
 49    }
 50
 51    /// Perform the per-event analysis
 52    void analyze(const Event& event) {
 53      // final state particles
 54      const FinalState& fs = apply<FinalState>(event, "FS");
 55      map<long,int> nCount;
 56      int ntotal(0);
 57      for (const Particle& p : fs.particles()) {
 58        nCount[p.pid()] += 1;
 59        ++ntotal;
 60      }
 61      // resonances
 62      for (const Particle & p : apply<UnstableParticles>(event, "UFS").particles(Cuts::abspid==223 ||
 63										Cuts::abspid==313 || Cuts::abspid==323 ||	Cuts::abspid==315 || Cuts::abspid==325)) {
 64        map<long,int> nRes = nCount;
 65        int ncount = ntotal;
 66        findChildren(p,nRes,ncount);
 67        if (ncount!=1) continue;
 68        unsigned int imode=0;
 69        long pid = 111;
 70        if (p.abspid()==313) {
 71          imode=1;
 72          pid=310;
 73        }
 74        else if (p.abspid()==323) {
 75          imode=2;
 76          pid = p.pid()>0 ? -321 : 321;
 77        }
 78        else if (p.abspid()==315) {
 79          imode=3;
 80          pid=310;
 81        }
 82        else if (p.abspid()==325) {
 83          imode=4;
 84          pid = p.pid()>0 ? -321 : 321;
 85        }
 86        bool matched=true;
 87        for (const auto& val : nRes) {
 88          if ((pid==310 && (val.first==310 || val.second==130)) || val.first==pid) {
 89            if (val.second!=1) {
 90              matched = false;
 91              break;
 92            }
 93          }
 94          else if (val.second!=0) {
 95            matched = false;
 96            break;
 97          }
 98        }
 99        if (matched) {
100          _sigma[imode]->fill(_ecms);
101          break;
102        }
103      }
104    }
105
106
107    /// Normalise histograms etc., after the run
108    void finalize() {
109      scale(_sigma, crossSection()/ sumOfWeights() /femtobarn);
110    }
111
112    /// @}
113
114
115    /// @name Histograms
116    /// @{
117    BinnedHistoPtr<string> _sigma[5];
118    string _ecms;
119    /// @}
120
121
122  };
123
124
125  RIVET_DECLARE_PLUGIN(BELLE_2013_I1252555);
126
127}