1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// Four-jet angles using Durham algorithm
 10  class OPAL_2001_I552446 : public Analysis {
 11  public:
 12
 13    RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_2001_I552446);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21
 22      // Initialise and register projections
 23      const FinalState fs;
 24      declare(fs, "FS");
 25      declare(FastJets(fs, JetAlg::DURHAM, 0.7), "Jets");
 26
 27      // Book histograms here
 28      book(_h_BZ      ,3, 1, 1);
 29      book(_h_KSW     ,4, 1, 1);
 30      book(_h_NR      ,5, 1, 1);
 31      book(_h_ALPHA34 ,6, 1, 1);
 32    }
 33
 34
 35    /// Perform the per-event analysis
 36    void analyze(const Event& event) {
 37
 38      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 39      if (apply<FinalState>(event, "FS").particles().size() < 2) {
 40        vetoEvent;
 41      }
 42
 43      const FastJets& fastjets = apply<FastJets>(event, "Jets");
 44      if (fastjets.clusterSeq()) {
 45        vector<fastjet::PseudoJet> jets;
 46        for (const fastjet::PseudoJet& jet :
 47                 fastjet::sorted_by_E(fastjets.clusterSeq()->exclusive_jets_ycut(0.008))) {
 48          if (jet.E()>3.0*GeV) jets.push_back(jet);
 49        }
 50        if (jets.size() == 4) {
 51          // Prevent nan-fill due to division by zero in calc_BZ
 52          double bz = fabs(calc_BZ(jets));
 53          if (!std::isnan(bz)) _h_BZ->fill(bz);
 54          _h_KSW->fill(calc_KSW(jets));
 55          _h_NR->fill(fabs(calc_NR(jets)));
 56          _h_ALPHA34->fill(calc_ALPHA34(jets));
 57        }
 58      }
 59
 60    }
 61
 62
 63    /// Normalise histograms etc., after the run
 64    void finalize() {
 65      normalize(_h_BZ);
 66      normalize(_h_KSW);
 67      normalize(_h_NR);
 68      normalize(_h_ALPHA34);
 69    }
 70
 71    /// @}
 72
 73
 74  private:
 75
 76    /// @name Jet angle calculator functions
 77    /// @{
 78
 79    /// @todo Use Jet or FourMomentum interface rather than PseudoJet
 80    /// @todo Move to utils?
 81    double calc_BZ(const vector<fastjet::PseudoJet>& jets) {
 82      assert(jets.size() == 4);
 83      Vector3 p12 = cross( momentum3(jets[0]), momentum3(jets[1]));
 84      Vector3 p34 = cross( momentum3(jets[2]), momentum3(jets[3]));
 85      return dot(p12,p34) / (p12.mod()*p34.mod());
 86    }
 87
 88
 89    /// @todo Use Jet or FourMomentum interface rather than PseudoJet
 90    /// @todo Move to utils?
 91    double calc_KSW(const vector<fastjet::PseudoJet>& jets) {
 92      assert(jets.size() == 4);
 93      Vector3 p13 = cross( momentum3(jets[0]), momentum3(jets[2]));
 94      Vector3 p24 = cross( momentum3(jets[1]), momentum3(jets[3]));
 95      Vector3 p14 = cross( momentum3(jets[0]), momentum3(jets[3]));
 96      Vector3 p23 = cross( momentum3(jets[1]), momentum3(jets[2]));
 97      return cos (0.5*( acos (dot(p14,p23) / (p14.mod()*p23.mod())) +
 98                        acos (dot(p13,p24) / (p13.mod()*p24.mod())) ));
 99    }
100
101
102    /// @todo Use Jet or FourMomentum interface rather than PseudoJet
103    /// @todo Move to utils?
104    double calc_NR(const vector<fastjet::PseudoJet>& jets) {
105      assert(jets.size() == 4);
106      Vector3 p12 = momentum3(jets[0]) - momentum3(jets[1]);
107      Vector3 p34 = momentum3(jets[2]) - momentum3(jets[3]);
108      return dot(p12,p34) / (p12.mod()*p34.mod());
109    }
110
111    /// @todo Use Jet or FourMomentum interface rather than PseudoJet
112    /// @todo Move to utils?
113    double calc_ALPHA34(const vector<fastjet::PseudoJet>& jets) {
114      assert(jets.size() == 4);
115      Vector3 p3 = momentum3(jets[2]);
116      Vector3 p4 = momentum3(jets[3]);
117      return dot(p3,p4) / (p3.mod()*p4.mod());
118    }
119
120    /// @}
121
122
123    /// @name Histograms
124    /// @{
125    Histo1DPtr _h_BZ;
126    Histo1DPtr _h_KSW;
127    Histo1DPtr _h_NR;
128    Histo1DPtr _h_ALPHA34;
129    /// @}
130
131  };
132
133
134
135  RIVET_DECLARE_ALIASED_PLUGIN(OPAL_2001_I552446, OPAL_2001_S4553896);
136
137}