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 measurement of R_b between 10.541 and 11.206 GeV
 10  class BABAR_2009_I797507 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2009_I797507);
 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_muons  , "/TMP/sigma_muons");
 30    }
 31
 32
 33    /// Perform the per-event analysis
 34    void analyze(const Event& event) {
 35      const FinalState& fs = apply<FinalState>(event, "FS");
 36      map<long,int> nCount;
 37      int ntotal(0);
 38      for (const Particle& p : fs.particles()) {
 39	nCount[p.pid()] += 1;
 40	++ntotal;
 41      }
 42      bool isBottom(false);
 43      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 44      for (const Particle& p : ufs.particles()) {
 45        if (PID::isBottomHadron(p.pid())) {
 46	  isBottom = true;
 47	  break;
 48	}
 49      }
 50
 51      // mu+mu- + photons
 52      if(nCount[-13]==1 and nCount[13]==1 &&
 53	 ntotal==2+nCount[22])
 54	_c_muons->fill();
 55      // everything else
 56      else if(isBottom) {
 57	_c_hadrons->fill();
 58      }
 59    }
 60
 61
 62    /// Normalise histograms etc., after the run
 63    void finalize() {
 64      Estimate0D R = *_c_hadrons/ *_c_muons;
 65      double fact = crossSection()/ sumOfWeights() /picobarn;
 66      double sig_h = _c_hadrons->val()*fact;
 67      double err_h = _c_hadrons->err()*fact;
 68      double sig_m = _c_muons  ->val()*fact;
 69      double err_m = _c_muons  ->err()*fact;
 70      Estimate1DPtr hadrons;
 71      book(hadrons, "sigma_hadrons");
 72      Estimate1DPtr muons;
 73      book(muons, "sigma_muons");
 74      Estimate1DPtr     mult;
 75      book(mult, 1, 1, 1);
 76      for (auto& b : mult->bins()) {
 77        if (inRange(sqrtS()/GeV, b.xMin(), b.xMax())) {
 78          b.set(R.val(), R.errPos());
 79          hadrons->bin(b.index()).set(sig_h, err_h);
 80          muons->bin(b.index()).set(sig_m, err_m);
 81        }
 82      }
 83    }
 84
 85    /// @}
 86
 87
 88    /// @name Histograms
 89    /// @{
 90    CounterPtr _c_hadrons, _c_muons;
 91    /// @}
 92
 93
 94  };
 95
 96
 97  RIVET_DECLARE_PLUGIN(BABAR_2009_I797507);
 98
 99
100}