1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Hadron Spectra in $e^+e^-$ collisions at 29 GeV
  9  class HRS_1987_I215848 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(HRS_1987_I215848);
 14
 15
 16    /// @name Analysis methods
 17    ///@{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      declare(UnstableParticles(), "UFS");
 22      // hists
 23      book(_d["pi"],  2, 1, 1);
 24      book(_d["Kp"],  3, 1, 1);
 25      book(_d["p"],   4, 1, 1);
 26      book(_h["K0"],  5, 1, 1);
 27      book(_h["lam"], 6, 1, 1);
 28      _axes["pi"]=YODA::Axis<double>({0.035800994999999995,0.049226100999999994,
 29          0.06281088,0.076470283,0.090170398,0.103895139,0.11763585,0.13138755});
 30      _axes["Kp"]=YODA::Axis<double>({0.048458616,0.059073984,0.070793603,0.083151849,
 31          0.095902214,0.108906933,0.122084929,0.135385471});
 32      _axes["p"]=YODA::Axis<double>({0.073322849,0.08073250800000001,0.08966457400000001,
 33          0.099710693,0.110567797,0.122019328,0.133913114,0.14614108599999998});
 34    }
 35
 36
 37    /// Perform the per-event analysis
 38    void analyze(const Event& event) {
 39      if (_edges.empty()) {
 40        for (const auto& item : _d) {
 41          _edges[item.first] = item.second->xEdges();
 42        }
 43      }
 44      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 45      for (const Particle & p : ufs.particles()) {
 46        const double xE = 2.*p.E()/sqrtS();
 47        const double beta = p.p3().mod() / p.E();
 48        if (p.pid()==130 || p.pid()==310) {
 49          _h["K0"]->fill(xE,1./beta);
 50        }
 51        else if (p.abspid()==321) {
 52          discfill("Kp",xE, 1./beta);
 53        }
 54        else if (p.abspid()==211) {
 55          discfill("pi",xE, 1./beta);
 56        }
 57        else if (p.abspid()==2212) {
 58          discfill("p",xE, 1./beta);
 59        }
 60        else if (p.abspid()==3122) {
 61          _h["lam"]->fill(xE, 1./beta);
 62        }
 63      }
 64    }
 65
 66    void discfill(const string& name, const double value, const double weight) {
 67      string edge = "OTHER";
 68      const size_t idx = _axes[name].index(value);
 69      if (idx && idx <= _edges[name].size())  edge = _edges[name][idx-1];
 70      _d[name]->fill(edge,weight);
 71    }
 72
 73    /// Normalise histograms etc., after the run
 74    void finalize() {
 75      const double sf = crossSection()*sqr(sqrtS())/microbarn/sumOfWeights();
 76      scale(_h, sf);
 77      scale(_d, sf);
 78      for( auto & hist : _d) {
 79        for(auto & b: hist.second->bins()) {
 80          const size_t idx = b.index();
 81          b.scaleW(1./_axes[hist.first].width(idx));
 82        }
 83      }
 84    }
 85
 86    ///@}
 87
 88
 89    /// @name Histograms
 90    ///@{
 91    map<string, Histo1DPtr> _h;
 92    map<string, BinnedHistoPtr<string> > _d;
 93    map<string, YODA::Axis<double> > _axes;
 94    map<string, vector<string> > _edges;
 95    ///@}
 96
 97
 98  };
 99
100
101  RIVET_DECLARE_PLUGIN(HRS_1987_I215848);
102
103}