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 Add a short analysis description here
 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/sigma_hadrons");
 29      book(_c_1S     , "/TMP/1S");
 30      book(_c_2S     , "/TMP/2S");
 31      book(_c_3S     , "/TMP/3S");
 32      book(_c_muons, "/TMP/sigma_muons");
 33    }
 34
 35    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 36      for (const Particle &child : p.children()) {
 37	if(child.children().empty()) {
 38	  --nRes[child.pid()];
 39	  --ncount;
 40	}
 41	else
 42	  findChildren(child,nRes,ncount);
 43      }
 44    }
 45
 46    /// Perform the per-event analysis
 47    void analyze(const Event& event) {
 48      // analyse the final state
 49      const FinalState& fs = apply<FinalState>(event, "FS");
 50      map<long,int> nCount;
 51      int ntotal(0);
 52      for (const Particle& p : fs.particles()) {
 53	nCount[p.pid()] += 1;
 54	++ntotal;
 55      }
 56      // intermediates
 57      bool isBottom(false);
 58      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 59      for (const Particle& p : ufs.particles()) {
 60	// check for bottom hadrons
 61        if (PID::isBottomHadron(p.pid())) {
 62	  isBottom = true;
 63	  break;
 64	}
 65	// upsilon + pi+pi-
 66	if(p.children().empty()) continue;
 67	if(p.pid() !=   553  &&
 68	   p.pid() != 100553 &&
 69	   p.pid() != 200553 ) continue;
 70	map<long,int> nRes = nCount;
 71	int ncount = ntotal;
 72	findChildren(p,nRes,ncount);
 73	if(ncount!=2) continue;
 74	bool matched = true;
 75	for(auto const & val : nRes) {
 76	  if(abs(val.first)==211) {
 77	    continue;
 78	  }
 79	  else if(val.second!=0) {
 80	    matched = false;
 81	    break;
 82	  }
 83	}
 84	if(matched) {
 85	  if(nRes[211]==1 && nRes[-211]==1 ) {
 86	    if(p.pid()==553)
 87	      _c_1S->fill();
 88	    if(p.pid()==100553)
 89	      _c_2S->fill();
 90	    if(p.pid()==200553)
 91	      _c_3S->fill();
 92	  }
 93	}
 94      }
 95      // mu+mu- + photons
 96      if(nCount[-13]==1 and nCount[13]==1 &&
 97	 ntotal==2+nCount[22])
 98	_c_muons->fill();
 99      // open bottom
100      else if(isBottom) {
101	_c_hadrons->fill();
102      }
103    }
104
105
106    /// Normalise histograms etc., after the run
107    void finalize() {
108      // loop over histos to be filled
109      for(unsigned int ix=1;ix<3;++ix) {
110        Estimate0D R;
111        for (unsigned int iy=1;iy<4;++iy) {
112          if (ix==2 && iy!=1) continue;
113          if (ix==1) {
114            if (iy==1) {
115              R = *_c_1S/ *_c_muons;
116            }
117            else if (iy==2) {
118              R = *_c_2S/ *_c_muons;
119            }
120            else {
121              R = *_c_3S/ *_c_muons;
122            }
123          }
124          else if(ix==2) {
125            R = *_c_hadrons/ *_c_muons;
126          }
127          Estimate1DPtr     mult;
128          book(mult, ix, 1, iy);
129          for (auto& b : mult->bins()) {
130            if (inRange(sqrtS()/MeV, b.xMin(), b.xMax())) {
131              b.set(R.val(), R.errPos());
132            }
133          }
134        }
135      }
136    }
137
138    /// @}
139
140
141    /// @name Histograms
142    /// @{
143    CounterPtr _c_hadrons, _c_muons, _c_1S, _c_2S, _c_3S;
144    /// @}
145
146
147  };
148
149
150  RIVET_DECLARE_PLUGIN(BELLE_2015_I1336624);
151
152}