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 Add a short analysis description here
 10  class CRYSTAL_BALL_1986_I238081 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(CRYSTAL_BALL_1986_I238081);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // Initialise and register projections
 23      declare(FinalState(), "FS");
 24      declare(UnstableParticles(), "UFS");
 25
 26      // Book histograms
 27      book(_c_hadrons, "/TMP/sigma_hadrons");
 28      book(_c_muons, "/TMP/sigma_muons");
 29      book(_c_D_star, "/TMP/sigma_D_star");
 30    }
 31
 32
 33    /// Perform the per-event analysis
 34    void analyze(const Event& event) {
 35      const FinalState& fs = apply<FinalState>(event, "FS");
 36
 37      map<long,int> nCount;
 38      int ntotal(0);
 39      for (const Particle& p : fs.particles()) {
 40	nCount[p.pid()] += 1;
 41	++ntotal;
 42      }
 43      // mu+mu- + photons
 44      if(nCount[-13]==1 and nCount[13]==1 &&
 45	 ntotal==2+nCount[22])
 46	_c_muons->fill();
 47      // everything else
 48      else
 49	_c_hadrons->fill();
 50
 51      const FinalState& ufs = apply<UnstableParticles>(event, "UFS");
 52      bool found = false;
 53      for (const Particle & p : ufs.particles()) {
 54	if(abs(p.pid())!=413 && abs(p.pid())!=423) continue;
 55	bool fs = true;
 56	for (const Particle & child : p.children()) {
 57	  if(child.pid()==p.pid()) {
 58	    fs = false;
 59	    break;
 60	  }
 61	}
 62	if(fs) {
 63	  found = true;
 64	  break;
 65	}
 66      }
 67      if(found)
 68	_c_D_star->fill();
 69    }
 70
 71
 72    /// Normalise histograms etc., after the run
 73    void finalize() {
 74      // R
 75      Estimate0D R = *_c_hadrons/ *_c_muons;
 76      double fact = crossSection()/ sumOfWeights() /picobarn;
 77      double sig_h = _c_hadrons->val()*fact;
 78      double err_h = _c_hadrons->err()*fact;
 79      double sig_m = _c_muons  ->val()*fact;
 80      double err_m = _c_muons  ->err()*fact;
 81      Estimate1DPtr hadrons;
 82      book(hadrons, "sigma_hadrons");
 83      Estimate1DPtr muons;
 84      book(muons, "sigma_muons");
 85      Estimate1DPtr mult;
 86      book(mult, 1, 1, 1);
 87      for (auto& b : mult->bins()) {
 88        if (inRange(sqrtS()/GeV, b.xMin(), b.xMax())) {
 89          b.set(R.val(), R.errPos());
 90          hadrons->bin(b.index()).set(sig_h, err_h);
 91          muons  ->bin(b.index()).set(sig_m, err_m);
 92        }
 93      }
 94      // D*
 95      fact = crossSection()/ sumOfWeights() /nanobarn;
 96      double sigma = _c_D_star->val()*fact;
 97      double error = _c_D_star->err()*fact;
 98      Estimate1DPtr mult2;
 99      book(mult2, 2, 1, 1);
100      for (auto& b : mult2->bins()) {
101      	if (inRange(sqrtS()/GeV, b.xMin(), b.xMax())) {
102      	  b.set(sigma, error);
103      	}
104      }
105    }
106
107    /// @}
108
109
110    /// @name Histograms
111    /// @{
112    CounterPtr _c_hadrons, _c_muons,_c_D_star;
113    /// @}
114
115
116  };
117
118
119  RIVET_DECLARE_PLUGIN(CRYSTAL_BALL_1986_I238081);
120
121
122}