1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/Thrust.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6#include "Rivet/Projections/Hemispheres.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief Thrust, heavy jet mass, and y3 at 58 GeV
 12  class TOPAZ_1993_I361661 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(TOPAZ_1993_I361661);
 17
 18
 19    /// @name Analysis methods
 20    /// @{
 21
 22    /// Book histograms and initialise projections before the run
 23    void init() {
 24
 25      const FinalState fs;
 26      declare(fs, "FS");
 27      declare(FastJets(fs, JetAlg::DURHAM, 0.7), "DurhamJets");
 28      const Thrust thrust(fs);
 29      declare(thrust, "Thrust");
 30      declare(Hemispheres(thrust), "Hemispheres");
 31
 32      // Book histograms
 33      book(_h["thrust"], 1, 1, 1);
 34      book(_h["rho"]   , 2, 1, 1);
 35      book(_h["y23"]   , 3, 1, 1);
 36
 37      _axes["thrust"] = YODA::Axis<double>(22, 0.8, 5.2);
 38      _axes["rho"] = YODA::Axis<double>(21, 1.0, 5.2);
 39      _axes["y23"] = YODA::Axis<double>(18, 1.2, 8.4);
 40
 41    }
 42
 43
 44    /// Perform the per-event analysis
 45    void analyze(const Event& event) {
 46      if (_edges.empty()) {
 47        for (const auto& item : _h) {
 48          _edges[item.first] = item.second->xEdges();
 49        }
 50      }
 51      // First, veto on leptonic events by requiring at least 4 charged FS particles
 52      const FinalState& fs = apply<FinalState>(event, "FS");
 53      const size_t numParticles = fs.particles().size();
 54      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 55      if (numParticles < 2) {
 56        MSG_DEBUG("Failed leptonic event cut");
 57        vetoEvent;
 58      }
 59      // thrust
 60      const Thrust& thrust = apply<Thrust>(event, "Thrust");
 61      fillhist("thrust", -log(1.-thrust.thrust()));
 62      // jet mass
 63      const Hemispheres& hemi = apply<Hemispheres>(event, "Hemispheres");
 64      fillhist("rho", -log(hemi.scaledM2high()));
 65      // Jets
 66      const FastJets& durjet = apply<FastJets>(event, "DurhamJets");
 67      if (numParticles>=3) {
 68        if (durjet.clusterSeq()) fillhist("y23", -log(durjet.clusterSeq()->exclusive_ymerge_max(2)));
 69      }
 70    }
 71
 72    void fillhist(const string& label, const double value) {
 73      string edge = "OTHER";
 74      const size_t idx = _axes[label].index(value);
 75      if (idx && idx <= _edges[label].size()) {
 76        edge = _edges[label][idx-1];
 77      }
 78      _h[label]->fill(edge);
 79    }
 80
 81
 82    /// Normalise histograms etc., after the run
 83    void finalize() {
 84      normalize(_h);
 85      for( auto & hist : _h) {
 86        for(auto & b: hist.second->bins()) {
 87          const size_t idx = b.index();
 88          b.scaleW(1./_axes[hist.first].width(idx));
 89        }
 90      }
 91    }
 92
 93    /// @}
 94
 95
 96    /// @name Histograms
 97    /// @{
 98    map<string, BinnedHistoPtr<string>> _h;
 99    map<string, YODA::Axis<double>> _axes;
100    map<string, vector<string>> _edges;
101    /// @}
102
103  };
104
105
106  RIVET_DECLARE_PLUGIN(TOPAZ_1993_I361661);
107
108
109}