Rivet analyses

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

namespace Rivet {


  /// @brief Add a short analysis description here
  class TOPAZ_1995_I381900 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(TOPAZ_1995_I381900);


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

    /// Book histograms and initialise projections before the run
    void init() {

      // Initialise and register projections
      declare(Beam(), "Beams");
      declare(UnstableParticles(), "UFS");
      const ChargedFinalState cfs;
      declare(cfs, "CFS");

      // Book histograms
      book(_h["charged"], 1, 1, 1);
      book(_h["pi"]     , 2, 1, 1);
      book(_h["Kp"]     , 2, 1, 2);
      book(_h["proton"] , 2, 1, 3);
      book(_h["K0"]     , 3, 1, 1);
      book(_wSum,"TMP/wSum");

      _axes["charged"] = YODA::Axis<double>(22, 0.6, 5.0);
      _axes["pi"] = YODA::Axis<double>{1.17, 1.47, 1.77, 2.07, 2.37, 2.545, 2.695, 2.85,3.1,3.315,
                                       3.46, 3.56, 3.69, 3.865, 4.05, 4.28, 4.565, 4.77, 4.89};
      _axes["K0"] = YODA::Axis<double>{1.5, 2.1, 2.4, 2.7, 3.0, 3.3, 3.6, 4.5};
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      if (_edges.empty()) {
        for (const auto& item : _h) {
          _edges[item.first] = item.second->xEdges();
        }
      }
      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
      int nch = cfs.particles().size();
      if (nch<5)  vetoEvent;
      // Get beams and average beam momentum
      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
      const double meanBeamMom = ( beams.first.p3().mod() + beams.second.p3().mod() ) / 2.0;
      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
      _wSum->fill();

      // neutral kaons
      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
      for (const Particle & p : ufs.particles(Cuts::pid==130 || Cuts::pid==310)) {
        const double xi = -log(p.p3().mod()/meanBeamMom);
        fillhist("K0", "K0", xi);
      }
      // charged particles
      for (const Particle& p : cfs.particles()) {
        double xi = -log(p.p3().mod()/meanBeamMom);
        fillhist("charged", "charged", xi);
        int id = abs(p.pid());
        if (id==211) {
          fillhist("pi", "pi", xi);
        }
        else if (id==321) {
          fillhist("Kp", "pi", xi);
        }
        else if (id==2212) {
          fillhist("proton", "pi", xi);
        }
      }
    }

    void fillhist(const string& label_hist, const string& label_axis, const double value) {
      string edge = "OTHER";
      size_t idx = _axes[label_axis].index(value);
      // binning for pi, Kp and proton same but different masked bins (pion only 1 masked bin 9)
      if ( label_axis=="pi" ) {
        if (idx==9) idx=0;
        else if (idx>9) --idx;
      }
      if (idx && idx <= _edges[label_axis].size())  edge = _edges[label_axis][idx-1];
      _h[label_hist]->fill(edge);
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      scale(_h, 1./ *_wSum);
      for( auto & hist : _h) {
        for(auto & b: hist.second->bins()) {
          size_t idx = b.index();
          if (hist.first=="pi" || _axes.find(hist.first)==_axes.end()) {
            auto it = std::find(_edges["pi"].begin(), _edges["pi"].end(), _edges[hist.first][idx-1]);
            idx = std::distance(_edges["pi"].begin(), it)+1;
            if(idx>=9) ++idx;
            b.scaleW(1./_axes["pi"].width(idx));
          }
          else
            b.scaleW(1./_axes[hist.first].width(idx));
        }
      }
    }

    /// @}


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


  };


  RIVET_DECLARE_PLUGIN(TOPAZ_1995_I381900);

}