1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief R measurement
  9  class TASSO_1982_I176887 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(TASSO_1982_I176887);
 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      for(unsigned int ix=0;ix<2;++ix) {
 26        book(_c_hadrons[ix], "/TMP/sigma_hadrons_"+toString(2*ix), refData<YODA::BinnedEstimate<string>>(1+2*ix,1,1));
 27        book(_c_muons[ix],   "/TMP/sigma_muons_"  +toString(2*ix), refData<YODA::BinnedEstimate<string>>(1+2*ix,1,1));
 28        for (const string& en : _c_hadrons[ix].binning().edges<0>()) {
 29          const size_t idx = en.find("-");
 30          if(idx!=std::string::npos) {
 31            const double emin = std::stod(en.substr(0,idx));
 32            const double emax = std::stod(en.substr(idx+1,string::npos));
 33            if(inRange(sqrtS()/GeV, emin, emax)) {
 34              _ecms[ix] = en;
 35              break;
 36            }
 37          }
 38          else {
 39            const double end = std::stod(en)*GeV;
 40            if (isCompatibleWithSqrtS(end)) {
 41              _ecms[ix] = en;
 42              break;
 43            }
 44          }
 45        }
 46      }
 47      book(_c_hadronsI, "/TMP/sigma_hadrons_1", refData<YODA::BinnedEstimate<int>>(2,1,1));
 48      book(_c_muonsI,   "/TMP/sigma_muons_1"  , refData<YODA::BinnedEstimate<int>>(2,1,1));
 49    }
 50
 51
 52    /// Perform the per-event analysis
 53    void analyze(const Event& event) {
 54      const FinalState& fs = apply<FinalState>(event, "FS");
 55
 56      map<long,int> nCount;
 57      int ntotal(0);
 58      for (const Particle& p : fs.particles()) {
 59	nCount[p.pid()] += 1;
 60	++ntotal;
 61      }
 62      // mu+mu- + photons
 63      if(nCount[-13]==1 and nCount[13]==1 &&
 64	 ntotal==2+nCount[22]) {
 65	for(unsigned int ix=0;ix<2;++ix) _c_muons[ix]->fill(_ecms[ix]);
 66        _c_muonsI->fill(round(sqrtS()/GeV));
 67      }
 68      // everything else
 69      else {
 70	for(unsigned int ix=0;ix<2;++ix) _c_hadrons[ix]->fill(_ecms[ix]);
 71        _c_hadronsI->fill(round(sqrtS()/GeV));
 72      }
 73    }
 74
 75
 76    /// Normalise histograms etc., after the run
 77    void finalize() {
 78      for(unsigned int ix=0;ix<2;++ix) {
 79        BinnedEstimatePtr<string> mult;
 80        book(mult, 1+2*ix, 1, 1);
 81        divide(_c_hadrons[ix],_c_muons[ix],mult);
 82      }
 83      BinnedEstimatePtr<int> multI;
 84      book(multI, 2, 1, 1);
 85      divide(_c_hadronsI,_c_muonsI,multI);
 86    }
 87    /// @}
 88
 89
 90    /// @name Histograms
 91    /// @{
 92    BinnedHistoPtr<string> _c_hadrons[2], _c_muons[2];
 93    string _ecms[2];
 94    BinnedHistoPtr<int> _c_hadronsI, _c_muonsI;
 95    /// @}
 96
 97
 98  };
 99
100
101  RIVET_DECLARE_PLUGIN(TASSO_1982_I176887);
102
103
104}