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 Add a short analysis description here
 10  class BELLE_2007_I723333 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2007_I723333);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22
 23      // Initialise and register projections
 24      declare(FinalState(), "FS");
 25      declare(UnstableParticles(), "UFS");
 26      book(_nDSS, "/TMP/nDSS");
 27      book(_nDS, "/TMP/nDS");
 28    }
 29
 30    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 31      for (const Particle &child : p.children()) {
 32        if (child.children().empty()) {
 33          nRes[child.pid()]-=1;
 34          --ncount;
 35        }
 36        else {
 37          findChildren(child,nRes,ncount);
 38        }
 39      }
 40    }
 41
 42    /// Perform the per-event analysis
 43    void analyze(const Event& event) {
 44      const FinalState& fs = apply<FinalState>(event, "FS");
 45
 46      map<long,int> nCount;
 47      int ntotal(0);
 48      for (const Particle& p : fs.particles()) {
 49        nCount[p.pid()] += 1;
 50        ++ntotal;
 51      }
 52      const FinalState& ufs = apply<FinalState>(event, "UFS");
 53
 54      for (unsigned int ix=0; ix<ufs.particles().size(); ++ix) {
 55        const Particle& p1 = ufs.particles()[ix];
 56        if (abs(p1.pid())!=413) continue;
 57        map<long,int> nRes = nCount;
 58        int ncount = ntotal;
 59        findChildren(p1,nRes,ncount);
 60        bool matched=false;
 61        int sign = -p1.pid()/abs(p1.pid());
 62        for (unsigned int iy=0; iy<ufs.particles().size(); ++iy) {
 63          if (ix==iy) continue;
 64          const Particle& p2 = ufs.particles()[iy];
 65          if (!p2.parents().empty() && p2.parents()[0].pid()==p1.pid()) {
 66            continue;
 67          }
 68          if (p2.pid()!=sign*413 && p2.pid()!=sign*411) 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 (const auto& val : nRes2) {
 75            if (val.second!=0) {
 76              matched = false;
 77              break;
 78            }
 79          }
 80          if (matched) {
 81            sign = abs(p2.pid());
 82            break;
 83          }
 84        }
 85        if (matched) {
 86          if (sign==411) {
 87            _nDS->fill();
 88          }
 89          else if (sign==413) {
 90            _nDSS->fill();
 91          }
 92        }
 93      }
 94    }
 95
 96
 97    /// Normalise histograms etc., after the run
 98    void finalize() {
 99      scale({_nDSS,_nDS}, crossSection()/ sumOfWeights() /nanobarn);
100      for (unsigned int ix=1; ix<3; ++ix) {
101        Estimate1DPtr mult;
102        book(mult, ix, 1, 1);
103        for (auto& b : mult->bins()) {
104          if (inRange(sqrtS()/GeV, b.xMin(), b.xMax())) {
105            b.set( (ix==1? _nDSS : _nDS)->val(),
106                   (ix==1? _nDSS : _nDS)->err() );
107          }
108        }
109      }
110    }
111
112    /// @}
113
114
115    /// @name Histograms
116    /// @{
117    CounterPtr _nDSS,_nDS;
118    /// @}
119
120
121  };
122
123
124  RIVET_DECLARE_PLUGIN(BELLE_2007_I723333);
125
126
127}