1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief e+e- > bbar and Upsilon(1,2,3) pi+pi-
 10  class BELLE_2015_I1336624 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2015_I1336624);
 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
 27      // Book histograms
 28      book(_c_hadrons, "TMP/c_hadrons", refData<YODA::BinnedEstimate<string>>(2,1,1));
 29      for(unsigned int ix=0;ix<3;++ix)
 30        book(_c_Ups[ix], "/TMP/c_Ups_"+toString(ix+1),refData<YODA::BinnedEstimate<string>>(1,1,1+ix));
 31      book(_c_muons[0], "/TMP/sigma_muons_1",refData<YODA::BinnedEstimate<string>>(2,1,1));
 32      book(_c_muons[1], "/TMP/sigma_muons_2",refData<YODA::BinnedEstimate<string>>(1,1,1));
 33      for(unsigned int ix=0;ix<2;++ix) {
 34        for (const string& en : _c_muons[ix].binning().edges<0>()) {
 35          const double end = std::stod(en)*MeV;
 36          if (isCompatibleWithSqrtS(end)) {
 37            _ecms[ix] = en;
 38            break;
 39          }
 40        }
 41      }
 42      if (_ecms[0].empty() && _ecms[1].empty())
 43        MSG_ERROR("Beam energy incompatible with analysis.");
 44    }
 45
 46    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 47      for (const Particle &child : p.children()) {
 48	if(child.children().empty()) {
 49	  --nRes[child.pid()];
 50	  --ncount;
 51	}
 52	else
 53	  findChildren(child,nRes,ncount);
 54      }
 55    }
 56
 57    /// Perform the per-event analysis
 58    void analyze(const Event& event) {
 59      // analyse the final state
 60      const FinalState& fs = apply<FinalState>(event, "FS");
 61      map<long,int> nCount;
 62      int ntotal(0);
 63      for (const Particle& p : fs.particles()) {
 64	nCount[p.pid()] += 1;
 65	++ntotal;
 66      }
 67      // intermediates
 68      bool isBottom(false);
 69      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 70      for (const Particle& p : ufs.particles()) {
 71	// check for bottom hadrons
 72        if (PID::isBottomHadron(p.pid()) &&
 73            (p.abspid()%1000)/10!=55) {
 74	  isBottom = true;
 75	  break;
 76	}
 77	// upsilon + pi+pi-
 78	if(p.children().empty()) continue;
 79	if(p.pid() !=   553  &&
 80	   p.pid() != 100553 &&
 81	   p.pid() != 200553 ) continue;
 82	map<long,int> nRes = nCount;
 83	int ncount = ntotal;
 84	findChildren(p,nRes,ncount);
 85	if(ncount!=2) continue;
 86	bool matched = true;
 87	for(auto const & val : nRes) {
 88	  if(abs(val.first)==211) {
 89	    continue;
 90	  }
 91	  else if(val.second!=0) {
 92	    matched = false;
 93	    break;
 94	  }
 95	}
 96	if(matched) {
 97	  if(nRes[211]==1 && nRes[-211]==1 ) {
 98	    if(p.pid()==553)
 99	      _c_Ups[0]->fill(_ecms[1]);
100	    if(p.pid()==100553)
101	      _c_Ups[1]->fill(_ecms[1]);
102	    if(p.pid()==200553)
103	      _c_Ups[2]->fill(_ecms[1]);
104	  }
105	}
106      }
107      // mu+mu- + photons
108      if(nCount[-13]==1 and nCount[13]==1 &&
109	 ntotal==2+nCount[22]) {
110	_c_muons[0]->fill(_ecms[0]);
111	_c_muons[1]->fill(_ecms[1]);
112      }
113      // open bottom
114      else if(isBottom) {
115	_c_hadrons->fill(_ecms[0]);
116      }
117    }
118
119
120    /// Normalise histograms etc., after the run
121    void finalize() {
122      BinnedEstimatePtr<string> tmp;
123      book(tmp,2,1,1);
124      divide(_c_hadrons,_c_muons[0],tmp);
125      for(unsigned int ix=0;ix<3;++ix) {
126        book(tmp,1,1,1+ix);
127        divide(_c_Ups[ix],_c_muons[1],tmp);
128      }
129    }
130
131    /// @}
132
133
134    /// @name Histograms
135    /// @{
136    BinnedHistoPtr<string> _c_hadrons, _c_muons[2], _c_Ups[3];
137    string _ecms[2];
138    /// @}
139
140
141  };
142
143
144  RIVET_DECLARE_PLUGIN(BELLE_2015_I1336624);
145
146}