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_S3612880 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(JADE_1998_S3612880);
 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, FastJets::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
 55
 56    /// Perform the per-event analysis
 57    void analyze(const Event& event) {
 58      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
 59
 60      // JADE hadronic event selection
 61      if (cfs.particles().size() < 3 ) vetoEvent;
 62
 63      const Thrust& thrust = apply<Thrust>(event, "Thrust");
 64      const Vector3 & thrustAxis = thrust.thrustAxis ();
 65      double theta = thrustAxis.theta();
 66      if ( fabs(cos(theta)) >= 0.8 ) {
 67        MSG_DEBUG("Failed thrust angle cut: " << fabs(cos(theta)));
 68        vetoEvent;
 69      }
 70      /// @todo Evis, pmiss, pbal
 71
 72      const Hemispheres& hemi = apply<Hemispheres>(event, "Hemispheres");
 73      const FastJets& durjet = apply<FastJets>(event, "DurhamJets");
 74
 75      const double y23 = durjet.clusterSeq()->exclusive_ymerge_max(2);
 76
 77      // Make sure we don't run into a segfault by trying to fill non-existing histos
 78      int s = int(sqrtS()/GeV);
 79      if (s == 44 || s == 35) {
 80        _h_thrust->fill(1. - thrust.thrust());
 81        _h_MH->fill(sqrt(hemi.scaledM2high()));
 82        _h_BT->fill(hemi.Bsum());
 83        _h_BW->fill(hemi.Bmax());
 84      }
 85      _h_y23->fill(y23);
 86    }
 87
 88    /// Normalise histograms etc., after the run
 89    void finalize() {
 90      // Make sure we don't try to normalise non-existing histos
 91      const int s = int(sqrtS()/GeV);
 92      if (s == 44 || s == 35) {
 93        normalize(_h_thrust);
 94        normalize(_h_MH);
 95        normalize(_h_BT);
 96        normalize(_h_BW);
 97      }
 98      normalize(_h_y23);
 99    }
100
101    //@}
102
103
104  private:
105
106    Histo1DPtr _h_thrust, _h_MH, _h_BT, _h_BW, _h_y23;
107
108  };
109
110
111
112  RIVET_DECLARE_ALIASED_PLUGIN(JADE_1998_S3612880, JADE_1998_I447560);
113
114}