1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Thrust.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/Hemispheres.hh"
  6#include "Rivet/Projections/ChargedFinalState.hh"
  7
  8namespace Rivet {
  9
 10
 11  class JADE_1998_I447560 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(JADE_1998_I447560);
 16
 17
 18    /// Book histograms and initialise projections before the run
 19    void init() {
 20      const ChargedFinalState cfs(Cuts::pT > 0.1*GeV);
 21      declare(cfs, "CFS");
 22      declare(FastJets(cfs, JetAlg::DURHAM, 0.7), "DurhamJets");
 23
 24      // Thrust
 25      const Thrust thrust(cfs);
 26      declare(thrust, "Thrust");
 27      declare(Hemispheres(thrust), "Hemispheres");
 28
 29      // Histos
 30      int offset = 0;
 31      switch ( int(sqrtS()/GeV) ) {
 32
 33        case 44:
 34          offset = 0;
 35          book(_h["thrust"], 2+offset, 1, 1);
 36          book(_h["MH"],     3+offset, 1, 1);
 37          book(_h["BT"],     4+offset, 1, 1);
 38          book(_h["BW"],     5+offset, 1, 1);
 39          book(_h["y23"], 10, 1, 1);
 40          break;
 41        case 35:
 42          offset = 4;
 43          book(_h["thrust"], 2+offset, 1, 1);
 44          book(_h["MH"],     3+offset, 1, 1);
 45          book(_h["BT"],     4+offset, 1, 1);
 46          book(_h["BW"],     5+offset, 1, 1);
 47          book(_h["y23"], 11, 1, 1);
 48          break;
 49        case 22:
 50          book(_h["y23"], 12, 1, 1);
 51          break;
 52      }
 53
 54      _axes["thrust"] = YODA::Axis<double>{0.0, 0.02, 0.04, 0.06, 0.08, 0.1, 0.12, 0.14, 0.16,
 55                                           0.18, 0.2, 0.23, 0.26, 0.3, 0.35, 0.4, 0.5};
 56      _axes["MH"] = YODA::Axis<double>{0.0, 0.06, 0.1, 0.14, 0.18, 0.22, 0.26, 0.3,
 57                                       0.34, 0.38, 0.42, 0.46, 0.5, 0.55, 0.6};
 58      _axes["BT"] = YODA::Axis<double>{0.0, 0.03, 0.06, 0.08, 0.1, 0.12, 0.14, 0.16,
 59                                       0.18, 0.2, 0.22, 0.24, 0.27, 0.3, 0.35, 0.4};
 60      _axes["BW"] = YODA::Axis<double>{0.0, 0.02, 0.04, 0.06, 0.08, 0.1, 0.12,
 61                                       0.14, 0.16, 0.18, 0.2, 0.23, 0.26, 0.3};
 62      _axes["y23"] = YODA::Axis<double>{0.0, 0.001, 0.002, 0.005, 0.01, 0.02, 0.03, 0.04,
 63                                        0.05, 0.06, 0.08, 0.1, 0.13, 0.16, 0.2, 0.25, 0.4};
 64    }
 65
 66
 67    /// Perform the per-event analysis
 68    void analyze(const Event& event) {
 69
 70      if (_edges.empty()) {
 71        for (const auto& item : _h) {
 72          _edges[item.first] = item.second->xEdges();
 73        }
 74      }
 75
 76      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
 77
 78      // JADE hadronic event selection
 79      if (cfs.particles().size() < 3 ) vetoEvent;
 80
 81      const Thrust& thrust = apply<Thrust>(event, "Thrust");
 82      const Vector3 & thrustAxis = thrust.thrustAxis ();
 83      double theta = thrustAxis.theta();
 84      if ( fabs(cos(theta)) >= 0.8 ) {
 85        MSG_DEBUG("Failed thrust angle cut: " << fabs(cos(theta)));
 86        vetoEvent;
 87      }
 88      /// @todo Evis, pmiss, pbal
 89
 90      const Hemispheres& hemi = apply<Hemispheres>(event, "Hemispheres");
 91      const FastJets& durjet = apply<FastJets>(event, "DurhamJets");
 92
 93      const double y23 = durjet.clusterSeq()->exclusive_ymerge_max(2);
 94
 95      // Make sure we don't run into a segfault by trying to fill non-existing histos
 96      const int s = int(sqrtS()/GeV);
 97      if (s == 44 || s == 35) {
 98        discfill("thrust", 1. - thrust.thrust());
 99        discfill("MH", sqrt(hemi.scaledM2high()));
100        discfill("BT", hemi.Bsum());
101        discfill("BW", hemi.Bmax());
102      }
103      discfill("y23", y23);
104    }
105
106    /// Normalise histograms etc., after the run
107    void finalize() {
108      normalize(_h);
109      for (auto& item : _h) {
110        for (auto& b : item.second->bins()) {
111          b.scaleW(_axes[item.first].width(b.index()));
112        }
113      }
114    }
115
116    void discfill(const string& name, const double value) {
117      string edge = "OTHER";
118      const size_t idx = _axes[name].index(value);
119      if (idx && idx <= _axes[name].numBins())  edge = _edges[name][idx-1];
120      _h[name]->fill(edge);
121    }
122
123
124  private:
125
126    map<string, BinnedHistoPtr<string>> _h;
127    map<string, vector<string>> _edges;
128    map<string, YODA::Axis<double>> _axes;
129
130  };
131
132
133  RIVET_DECLARE_ALIASED_PLUGIN(JADE_1998_I447560, JADE_1998_S3612880);
134
135}