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      Scatter1D R = *_c_hadrons/ *_c_muons;
 65      double              rval = R.point(0).x();
 66      pair<double,double> rerr = R.point(0).xErrs();
 67      double fact = crossSection()/ sumOfWeights() /picobarn;
 68      double sig_h = _c_hadrons->val()*fact;
 69      double err_h = _c_hadrons->err()*fact;
 70      double sig_m = _c_muons  ->val()*fact;
 71      double err_m = _c_muons  ->err()*fact;
 72      Scatter2D temphisto(refData(1, 1, 1));
 73      Scatter2DPtr hadrons;
 74      book(hadrons, "sigma_hadrons");
 75      Scatter2DPtr muons;
 76      book(muons, "sigma_muons"  );
 77      Scatter2DPtr     mult;
 78      book(mult, 1, 1, 1);
 79      for (size_t b = 0; b < temphisto.numPoints(); b++) {
 80	const double x  = temphisto.point(b).x();
 81	pair<double,double> ex = temphisto.point(b).xErrs();
 82	pair<double,double> ex2 = ex;
 83	if(ex2.first ==0.) ex2. first=0.0001;
 84	if(ex2.second==0.) ex2.second=0.0001;
 85	if (inRange(sqrtS()/GeV, x-ex2.first, x+ex2.second)) {
 86	  mult   ->addPoint(x, rval, ex, rerr);
 87	  hadrons->addPoint(x, sig_h, ex, make_pair(err_h,err_h));
 88	  muons  ->addPoint(x, sig_m, ex, make_pair(err_m,err_m));
 89	}
 90	else {
 91	  mult   ->addPoint(x, 0., ex, make_pair(0.,.0));
 92	  hadrons->addPoint(x, 0., ex, make_pair(0.,.0));
 93	  muons  ->addPoint(x, 0., ex, make_pair(0.,.0));
 94	}
 95      }
 96    }
 97
 98    //@}
 99
100
101    /// @name Histograms
102    //@{
103    CounterPtr _c_hadrons, _c_muons;
104    //@}
105
106
107  };
108
109
110  // The hook for the plugin system
111  RIVET_DECLARE_PLUGIN(BABAR_2009_I797507);
112
113
114}