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- > hb pi+pi-
 10  class BELLE_2016_I1389855 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2016_I1389855);
 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==10553 ||
 24                                Cuts::pid==110553), "UFS");
 25
 26      book(_nhb1, 1, 1, 3);
 27      book(_nhb2, 1, 1, 4);
 28
 29      // Check which indices match the Ecm energy
 30      YODA::BinnedEstimate<int> tmp = refData<YODA::BinnedEstimate<int>>(1, 1, 1); 
 31      for (const auto& est : tmp.bins()) {
 32        if (abs(sqrtS()/MeV-est.val())<est.totalErrAvg()) {
 33          edges.push_back(est.xEdge());
 34        }
 35      }
 36      if (edges.empty())
 37        MSG_ERROR("Beam energy incompatible with analysis.");
 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
 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      const FinalState& ufs = apply<FinalState>(event, "UFS");
 63      for (const Particle& p : ufs.particles()) {
 64        if (p.children().empty()) continue;
 65        map<long,int> nRes = nCount;
 66        int ncount = ntotal;
 67        findChildren(p,nRes,ncount);
 68        // omega pi+pi-
 69        if(ncount!=2) continue;
 70        bool matched = true;
 71        for (const auto& val : nRes) {
 72          if (abs(val.first)==211) {
 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) {
 84          for (const int Ecm : edges) {
 85            if (p.pid() == 10553) {
 86              _nhb1->fill(Ecm);
 87            }
 88            else {
 89              _nhb2->fill(Ecm);
 90            }
 91          }
 92          break;
 93        }
 94      }
 95    }
 96
 97
 98    /// Normalise histograms etc., after the run
 99    void finalize() {
100      const double sf = crossSection()/sumOfWeights()/picobarn;
101      scale({_nhb1, _nhb2}, sf);
102    }
103    /// @}
104
105
106    /// @name Histograms
107    /// @{
108    BinnedHistoPtr<int> _nhb1, _nhb2;
109    vector<int> edges;
110    /// @}
111
112  };
113
114  RIVET_DECLARE_PLUGIN(BELLE_2016_I1389855);
115}