Rivet analyses

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

namespace Rivet {


  /// @brief  gamma, pi0 and eta spectra at 14, 22.5 and 34.4 GeV
  class JADE_1985_I213948 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(JADE_1985_I213948);


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

    /// Book histograms and initialise projections before the run
    void init() {
      declare(Beam(), "Beams");
      declare(FinalState(), "FS");
      declare(UnstableParticles(), "UFS");
      // book histos
      size_t ih = 2;
      for (double eVal : allowedEnergies()) {
        const string en = toString(round(eVal/MeV));
        if (isCompatibleWithSqrtS(eVal, 1e-3))  _sqs = en;
        book(_h[en+"gamma"], ih+1, 1, 1);
        book(_h[en+"pi0"],   ih+4, 1, 1);
        if (en == "34500"s) {
          book(_h[en+"eta"], 7, 1, 1);
          _axes[en+"gamma"] = vector<double>{0.00, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10,
                                             0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.20,
                                             0.22, 0.24, 0.26, 0.28, 0.30, 0.32, 0.34, 0.36, 0.38, 0.40,
                                             0.45, 0.50, 0.55, 0.60, 0.70, 0.80, 0.90}; // d01
          _axes[en+"pi0"] = vector<double>{0.011, 0.023, 0.035, 0.0465, 0.058, 0.0695, 0.081,
                                           0.093, 0.1045, 0.119, 0.1395, 0.1625, 0.1915, 0.2265}; // d04
          _axes[en+"eta"] = vector<double>{0.040, 0.056, 0.088, 0.274}; // d07
        }
        else if (en == "22500"s) {
          _axes[en+"gamma"] = vector<double>{0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10,
                                             0.12, 0.14, 0.16, 0.18, 0.20, 0.25, 0.40, 0.60, 1.00}; // d02
          _axes[en+"pi0"] = vector<double>{0.018, 0.036, 0.054, 0.072, 0.108, 0.144, 0.216, 0.360}; // d05
        }
        else {
          _axes[en+"gamma"] = vector<double>{0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10,
                                             0.12, 0.14, 0.16, 0.18, 0.20, 0.25, 0.40, 0.60, 1.00}; //d03
          _axes[en+"pi0"] = vector<double>{0.030, 0.058, 0.086, 0.114, 0.170, 0.226, 0.338, 0.562}; // d06
        }
        --ih;
      }
      raiseBeamErrorIf(_sqs.empty());
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      if (_edges.empty()) {
        _edges["gamma"] = _h[_sqs+"gamma"]->xEdges();
        _edges["pi0"] = _h[_sqs+"pi0"]->xEdges();
        if (_sqs == "34500"s)  _edges["eta"] = _h[_sqs+"eta"]->xEdges();
      }
      // Get beams and average beam momentum
      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
      const double meanBeamMom = 0.5*(beams.first.p3().mod()+beams.second.p3().mod());
      // gamma
      const FinalState& fs = apply<FinalState>(event, "FS");
      for (const Particle& p : fs.particles(Cuts::pid==22)) {
        const double xE = p.E()/meanBeamMom;
        _h[_sqs+"gamma"]->fill(map2string(xE, "gamma"));
      }
      // pi0, eta
      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
      for (const Particle& p : ufs.particles(Cuts::pid==111 || Cuts::pid==221)) {
        double xE = p.E()/meanBeamMom;
        if (p.pid()==111) {
          _h[_sqs+"pi0"]->fill(map2string(xE, "pi0"));
        }
        else if (_sqs == "34500"s) {
          _h[_sqs+"eta"]->fill(map2string(xE, "eta"));
        }
      }
    }

    string map2string(const double val, const string& obs) const {
      const size_t idx = _axes.at(_sqs+obs).index(val);
      if (idx && idx <= _axes.at(_sqs+obs).numBins())  return _edges.at(obs)[idx-1];
      return "OTHER"s;
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      for (auto& item : _h) {
        const double en = stod(item.first.substr(0,5));
        scale(item.second, crossSection()*sqr(en*MeV)/microbarn/sumOfWeights());
        for(auto& b: item.second->bins()) {
          b.scaleW(1./_axes[item.first].width(b.index()));
        }
      }
    }

    /// @}


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


  };


  RIVET_DECLARE_PLUGIN(JADE_1985_I213948);
}