Rivet analyses

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

namespace Rivet {


  /// @brief rho0, K*=, K*0 spectra at 29 GeV
  class HRS_1989_I276948 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(HRS_1989_I276948);


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

    /// Book histograms and initialise projections before the run
    void init() {
      declare(UnstableParticles(), "UFS");
      book(_h["rho"], 1, 1, 1);
      book(_h["Kstar0"], 2, 1, 1);
      book(_h["Kstarp"], 2, 1, 2);
      book(_h_sig_rho    ,3,1,1);
      book(_h_sig_Kstar0 ,4,1,1);
      book(_h_mult_rho   ,3,1,3);
      book(_h_mult_Kstar0,4,1,3);
      _axes["rho"] = YODA::Axis<double>({0.0725, 0.0925, 0.114, 0.136, 0.159, 0.183, 0.207, 0.231,
                                         0.2555, 0.28, 0.3045, 0.3295, 0.354, 0.3785, 0.4035, 0.4285,
                                         0.4595, 0.5155, 0.6025, 0.702, 0.802, 0.902});
      _axes["Kstar0"] = YODA::Axis<double>({0.077, 0.119, 0.163, 0.2095, 0.2575, 0.3065,
                                            0.3555, 0.405, 0.455, 0.5, 0.6, 0.7, 1.0});
      _axes["Kstarp"] = _axes["Kstar0"];
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      if (_edges.empty()) {
        for (const auto& item : _h) {
          _edges[item.first] = item.second->xEdges();
        }
      }
      UnstableParticles ufs = apply<UnstableParticles>(event,"UFS");
      for (const Particle& p : ufs.particles(Cuts::abspid==113 || Cuts::abspid==323 || Cuts::abspid==313)) {
        Vector3 mom3 = p.p3();
        const double energy = p.E();
        double modp = mom3.mod();
        double beta = modp/energy;
        double xE = 2.*modp/sqrtS();
        if (p.pid()==113) {
          discfill("rho", xE, 1./beta);
          _h_sig_rho->fill(29);
          _h_mult_rho->fill(29);
        }
        else if (p.pid()==313) {
          discfill("Kstar0", xE, 1./beta);
          _h_sig_Kstar0->fill(29);
          _h_mult_Kstar0->fill(29);
        }
        else {
          discfill("Kstarp", xE, 1./beta);
        }
      }
    }

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


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

    ///@}


    /// @name Histograms
    ///@{
    BinnedHistoPtr<int> _h_sig_rho,_h_sig_Kstar0;
    BinnedHistoPtr<int> _h_mult_rho,_h_mult_Kstar0;
    map<string,BinnedHistoPtr<string>> _h;
    map<string,YODA::Axis<double>> _axes;
    map<string,vector<string>> _edges;
    ///@}


  };


  RIVET_DECLARE_PLUGIN(HRS_1989_I276948);

}