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- > J/psi X(3940)
 10  class BELLE_2007_I686580 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2007_I686580);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // set the PDG code
 23      _pid = getOption<int>("PID", 9010441);
 24      // projections
 25      declare("FS",FinalState());
 26      declare("UFS",UnstableParticles(Cuts::pid==443 || Cuts::pid==_pid));
 27      // histograms
 28      book(_h,1,1,1);
 29    }
 30
 31    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount,
 32		      unsigned int & nCharged) {
 33      for( const Particle &child : p.children()) {
 34        if (child.children().empty()) {
 35          --nRes[child.pid()];
 36          --ncount;
 37          if(PID::isCharged(p.pid())) ++nCharged;
 38        }
 39        else {
 40          findChildren(child,nRes,ncount,nCharged);
 41        }
 42      }
 43    }
 44
 45    /// Perform the per-event analysis
 46    void analyze(const Event& event) {
 47      // final state particles
 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      // loop over J/psi
 56      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 57      bool matched=false;
 58      for (const Particle& p : ufs.particles(Cuts::pid==443)) {
 59        if (p.children().empty()) continue;
 60        map<long,int> nRes = nCount;
 61        int ncount = ntotal;
 62        unsigned int nCharged=0;
 63        findChildren(p,nRes,ncount,nCharged);
 64        // X(3940)
 65      	for (const Particle& p2 : ufs.particles(Cuts::pid==_pid)) {
 66      	  map<long,int> nResB = nRes;
 67      	  int ncountB = ncount;
 68      	  unsigned int nChargedB=0;
 69      	  findChildren(p2,nResB,ncountB,nChargedB);
 70      	  if (ncountB!=0) continue;
 71      	  matched = true;
 72      	  for (const auto& val : nResB) {
 73      	    if (val.second!=0) {
 74      	      matched = false;
 75      	      break;
 76      	    }
 77      	  }
 78      	  if (matched) {
 79            if (nCharged>2) _h->fill("10.6"s);
 80            break;
 81          }
 82        }
 83        if (matched) break;
 84      }
 85    }
 86
 87
 88    /// Normalise histograms etc., after the run
 89    void finalize() {
 90      scale(_h, crossSection()/ sumOfWeights() /femtobarn);
 91    }
 92
 93    /// @}
 94
 95
 96    /// @name Histograms
 97    /// @{
 98    int _pid;
 99    BinnedHistoPtr<string> _h;
100    /// @}
101
102
103  };
104
105
106  RIVET_DECLARE_PLUGIN(BELLE_2007_I686580);
107
108}