1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief J/psi from chi_c decays
  9  class CDF_1997_I440446 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(CDF_1997_I440446);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      declare(UnstableParticles(), "UFS");
 22      const YODA::BinnedEstimate<int>& ref = refData<YODA::BinnedEstimate<int>>(1,1,1);
 23      book(_h_total[0], "TMP/h_chi", ref);
 24      book(_h_total[1], "TMP/h_psi", ref);
 25      for (unsigned int ix=0; ix<3; ++ix) {
 26        book(_h_psi[ix], 2+ix, 1, 1);
 27      }
 28      _axis = YODA::Axis<double>({5.,5.5,6.,6.5,7.,8.,9.,10.,12.,14.,17.,20.});
 29    }
 30
 31
 32    /// Perform the per-event analysis
 33    void analyze(const Event& event) {
 34      if (_edges.empty())  _edges = _h_psi[0]->xEdges();
 35
 36      // Final state of unstable particles to get particle spectra
 37      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 38      // first J/psi for denominator
 39      for (const Particle& p : ufs.particles(Cuts::pid==443)) {
 40        if (p.fromBottom()) continue;
 41        const double abseta = p.abseta();
 42        const double xp = p.perp();
 43        if (xp<4. || abseta>0.6) continue;
 44        _h_total[1]->fill(1800);
 45        // from those for higher charmonium
 46        Particle parent = p.parents()[0];
 47        if ( parent.pid()==100443 || parent.pid()==20443 || parent.pid()==445) continue;
 48        _h_psi[0]->fill(disc(xp));
 49      }
 50      // chi_1 and chi_2 for numerator
 51      for (const Particle& p : ufs.particles(Cuts::pid==20443 || Cuts::pid==100443 || Cuts::pid==445)) {
 52        if(p.fromBottom()) continue;
 53        Particle Jpsi;
 54        bool found(false);
 55        for (const Particle & child : p.children()) {
 56          if (child.pid()==443) {
 57            found = true;
 58            Jpsi=child;
 59          }
 60        }
 61        if(!found) continue;
 62        double abseta=Jpsi.abseta();
 63        double xp = Jpsi.perp();
 64        if (xp<4. || abseta>0.6) continue;
 65        if (p.pid()==100443) {
 66          _h_psi[2]->fill(disc(xp));
 67        }
 68        else {
 69          _h_psi[1]->fill(disc(xp));
 70          _h_total[0]->fill(1800);
 71        }
 72      }
 73    }
 74
 75    string disc(const double value) {
 76      size_t idx = _axis.index(value);
 77      if (0 < idx && idx <= _axis.numBins())  return _edges[idx-1];
 78      return "OTHER"s;
 79    }
 80
 81
 82    /// Normalise histograms etc., after the run
 83    void finalize() {
 84      // normalisation factor
 85      const double br = 0.05961;
 86      const double factor = br*crossSection()/nanobarn/sumOfWeights();
 87      scale(_h_psi,factor);
 88      for(unsigned int ix=0;ix<3;++ix) {
 89        for(auto & b : _h_psi[ix]->bins())
 90          b.scaleW(1./_axis.width(b.index()));
 91      }
 92      BinnedEstimatePtr<int> tmp;
 93      book(tmp, 1, 1, 1);
 94      efficiency(_h_total[0], _h_total[1], tmp);
 95    }
 96
 97    /// @}
 98
 99
100    /// @name Histograms
101    /// @{
102    BinnedHistoPtr<int> _h_total[2];
103    BinnedHistoPtr<string> _h_psi[3];
104    YODA::Axis<double> _axis;
105    vector<string> _edges;
106    /// @}
107
108
109  };
110
111
112  RIVET_DECLARE_PLUGIN(CDF_1997_I440446);
113
114}