Rivet analyses

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

namespace Rivet {


  /// @brief J/psi from chi_c decays
  class CDF_1997_I440446 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(CDF_1997_I440446);


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

    /// Book histograms and initialise projections before the run
    void init() {
      declare(UnstableParticles(), "UFS");
      const YODA::BinnedEstimate<int>& ref = refData<YODA::BinnedEstimate<int>>(1,1,1);
      book(_h_total[0], "TMP/h_chi", ref);
      book(_h_total[1], "TMP/h_psi", ref);
      for (unsigned int ix=0; ix<3; ++ix) {
        book(_h_psi[ix], 2+ix, 1, 1);
      }
      _axis = YODA::Axis<double>({5.,5.5,6.,6.5,7.,8.,9.,10.,12.,14.,17.,20.});
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      if (_edges.empty())  _edges = _h_psi[0]->xEdges();

      // Final state of unstable particles to get particle spectra
      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
      // first J/psi for denominator
      for (const Particle& p : ufs.particles(Cuts::pid==443)) {
        if (p.fromBottom()) continue;
        const double abseta = p.abseta();
        const double xp = p.perp();
        if (xp<4. || abseta>0.6) continue;
        _h_total[1]->fill(1800);
        // from those for higher charmonium
        Particle parent = p.parents()[0];
        if ( parent.pid()==100443 || parent.pid()==20443 || parent.pid()==445) continue;
        _h_psi[0]->fill(disc(xp));
      }
      // chi_1 and chi_2 for numerator
      for (const Particle& p : ufs.particles(Cuts::pid==20443 || Cuts::pid==100443 || Cuts::pid==445)) {
        if(p.fromBottom()) continue;
        Particle Jpsi;
        bool found(false);
        for (const Particle & child : p.children()) {
          if (child.pid()==443) {
            found = true;
            Jpsi=child;
          }
        }
        if(!found) continue;
        double abseta=Jpsi.abseta();
        double xp = Jpsi.perp();
        if (xp<4. || abseta>0.6) continue;
        if (p.pid()==100443) {
          _h_psi[2]->fill(disc(xp));
        }
        else {
          _h_psi[1]->fill(disc(xp));
          _h_total[0]->fill(1800);
        }
      }
    }

    string disc(const double value) {
      size_t idx = _axis.index(value);
      if (0 < idx && idx <= _axis.numBins())  return _edges[idx-1];
      return "OTHER"s;
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      // normalisation factor
      const double br = 0.05961;
      const double factor = br*crossSection()/nanobarn/sumOfWeights();
      scale(_h_psi,factor);
      for(unsigned int ix=0;ix<3;++ix) {
        for(auto & b : _h_psi[ix]->bins())
          b.scaleW(1./_axis.width(b.index()));
      }
      BinnedEstimatePtr<int> tmp;
      book(tmp, 1, 1, 1);
      efficiency(_h_total[0], _h_total[1], tmp);
    }

    /// @}


    /// @name Histograms
    /// @{
    BinnedHistoPtr<int> _h_total[2];
    BinnedHistoPtr<string> _h_psi[3];
    YODA::Axis<double> _axis;
    vector<string> _edges;
    /// @}


  };


  RIVET_DECLARE_PLUGIN(CDF_1997_I440446);

}