1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Add a short analysis description here
  9  class GAMMAGAMMA_1979_I141722 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(GAMMAGAMMA_1979_I141722);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      // Initialise and register projections
 22      declare(FinalState(), "FS");
 23
 24      // Book histograms
 25      book(_c_hadrons, "/TMP/sigma_hadrons");
 26      book(_c_muons, "/TMP/sigma_muons");
 27      book(_c_charged, "/TMP/Ncharged");
 28      book(_c_neutral, "/TMP/Nneutral");
 29      book(_nHadrons, "/TMP/NHadrons");
 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),ncharged(0),nneutral(0);
 39      for (const Particle& p : fs.particles()) {
 40	nCount[p.pid()] += 1;
 41	++ntotal;
 42	if(PID::isCharged(p.pid()))
 43	  ncharged += 1;
 44	else
 45	  nneutral += 1;
 46      }
 47      // mu+mu- + photons
 48      if(nCount[-13]==1 and nCount[13]==1 &&
 49	 ntotal==2+nCount[22])
 50	_c_muons->fill();
 51      // everything else
 52      else {
 53	if(ntotal==2) vetoEvent;
 54	_c_hadrons->fill();
 55	_c_charged->fill(ncharged);
 56	_c_neutral->fill(nneutral);
 57	_nHadrons->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      scale(_c_charged, 1./_nHadrons->sumW());
 84      scale(_c_neutral, 1./_nHadrons->sumW());
 85      for (unsigned int iy=1; iy<3;++iy) {
 86        double aver(0.),error(0.);
 87        if(iy==1) {
 88          aver  = _c_charged->val();
 89          error = _c_charged->err();
 90        }
 91        else {
 92          aver  = _c_neutral->val();
 93          error = _c_neutral->err();
 94        }
 95        Estimate1DPtr mult;
 96        book(mult, 2, 1, iy);
 97        for (auto& b : mult->bins()) {
 98          if (inRange(sqrtS()/GeV, b.xMin(), b.xMax())) {
 99            b.set(aver, error);
100          }
101        }
102      }
103    }
104
105    /// @}
106
107
108    /// @name Histograms
109    /// @{
110    CounterPtr _c_hadrons, _c_muons,_c_neutral,_c_charged;
111    CounterPtr _nHadrons;
112    /// @}
113
114
115  };
116
117
118  RIVET_DECLARE_PLUGIN(GAMMAGAMMA_1979_I141722);
119
120
121}