Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/UnstableParticles.hh"

namespace Rivet {


  /// @brief Hadron Spectra in $e^+e^-$ collisions at 29 GeV
  class HRS_1987_I215848 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(HRS_1987_I215848);


    /// @name Analysis methods
    ///@{

    /// Book histograms and initialise projections before the run
    void init() {
      declare(UnstableParticles(), "UFS");
      // hists
      book(_d["pi"],  2, 1, 1);
      book(_d["Kp"],  3, 1, 1);
      book(_d["p"],   4, 1, 1);
      book(_h["K0"],  5, 1, 1);
      book(_h["lam"], 6, 1, 1);
      _axes["pi"]=YODA::Axis<double>({0.035800994999999995,0.049226100999999994,
          0.06281088,0.076470283,0.090170398,0.103895139,0.11763585,0.13138755});
      _axes["Kp"]=YODA::Axis<double>({0.048458616,0.059073984,0.070793603,0.083151849,
          0.095902214,0.108906933,0.122084929,0.135385471});
      _axes["p"]=YODA::Axis<double>({0.073322849,0.08073250800000001,0.08966457400000001,
          0.099710693,0.110567797,0.122019328,0.133913114,0.14614108599999998});
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      if (_edges.empty()) {
        for (const auto& item : _d) {
          _edges[item.first] = item.second->xEdges();
        }
      }
      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
      for (const Particle & p : ufs.particles()) {
        const double xE = 2.*p.E()/sqrtS();
        const double beta = p.p3().mod() / p.E();
        if (p.pid()==130 || p.pid()==310) {
          _h["K0"]->fill(xE,1./beta);
        }
        else if (p.abspid()==321) {
          discfill("Kp",xE, 1./beta);
        }
        else if (p.abspid()==211) {
          discfill("pi",xE, 1./beta);
        }
        else if (p.abspid()==2212) {
          discfill("p",xE, 1./beta);
        }
        else if (p.abspid()==3122) {
          _h["lam"]->fill(xE, 1./beta);
        }
      }
    }

    void discfill(const string& name, const double value, const double weight) {
      string edge = "OTHER";
      const size_t idx = _axes[name].index(value);
      if (idx && idx <= _edges[name].size())  edge = _edges[name][idx-1];
      _d[name]->fill(edge,weight);
    }

    /// Normalise histograms etc., after the run
    void finalize() {
      const double sf = crossSection()*sqr(sqrtS())/microbarn/sumOfWeights();
      scale(_h, sf);
      scale(_d, sf);
      for( auto & hist : _d) {
        for(auto & b: hist.second->bins()) {
          const size_t idx = b.index();
          b.scaleW(1./_axes[hist.first].width(idx));
        }
      }
    }

    ///@}


    /// @name Histograms
    ///@{
    map<string, Histo1DPtr> _h;
    map<string, BinnedHistoPtr<string> > _d;
    map<string, YODA::Axis<double> > _axes;
    map<string, vector<string> > _edges;
    ///@}


  };


  RIVET_DECLARE_PLUGIN(HRS_1987_I215848);

}