1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "fastjet/JadePlugin.hh"
  6
  7namespace fastjet {
  8
  9class P_scheme : public JetDefinition::Recombiner {
 10 public:
 11  std::string description() const {return "";}
 12  void recombine(const PseudoJet & pa, const PseudoJet & pb,
 13		 PseudoJet & pab) const {
 14    PseudoJet tmp = pa + pb;
 15    double E = sqrt(tmp.px()*tmp.px() + tmp.py()*tmp.py() + tmp.pz()*tmp.pz());
 16    pab.reset_momentum(tmp.px(), tmp.py(), tmp.pz(), E);
 17  }
 18  void preprocess(PseudoJet & p) const {
 19    double E = sqrt(p.px()*p.px() + p.py()*p.py() + p.pz()*p.pz());
 20    p.reset_momentum(p.px(), p.py(), p.pz(), E);
 21  }
 22  ~P_scheme() { }
 23};
 24
 25class E0_scheme : public JetDefinition::Recombiner {
 26 public:
 27  std::string description() const {return "";}
 28  void recombine(const PseudoJet & pa, const PseudoJet & pb,
 29		 PseudoJet & pab) const {
 30    PseudoJet tmp = pa + pb;
 31    double fact = tmp.E()/sqrt(tmp.px()*tmp.px()+tmp.py()*tmp.py()+tmp.pz()*tmp.pz());
 32    pab.reset_momentum(fact*tmp.px(), fact*tmp.py(), fact*tmp.pz(), tmp.E());
 33  }
 34  void preprocess(PseudoJet & p) const {
 35    double fact = p.E()/sqrt(p.px()*p.px()+p.py()*p.py()+p.pz()*p.pz());
 36
 37    p.reset_momentum(fact*p.px(), fact*p.py(), fact*p.pz(), p.E());
 38  }
 39  ~E0_scheme() { }
 40};
 41
 42}
 43
 44
 45namespace Rivet {
 46
 47
 48  /// @brief AMY jets at
 49  class AMY_1995_I406129 : public Analysis {
 50  public:
 51
 52    /// Constructor
 53    RIVET_DEFAULT_ANALYSIS_CTOR(AMY_1995_I406129);
 54
 55
 56    /// @name Analysis methods
 57    /// @{
 58
 59    /// Book histograms and initialise projections before the run
 60    void init() {
 61      // Initialise and register projections
 62      FinalState fs;
 63      declare(fs, "FS");
 64      // Book histograms
 65      book(_h["jade_P"],  2, 1, 1);
 66      book(_h["jade_E"],  3, 1, 1);
 67      book(_h["jade_E0"], 6, 1, 1);
 68      book(_h["durham"],  4, 1, 1);
 69    }
 70
 71
 72    /// Perform the per-event analysis
 73    void analyze(const Event& event) {
 74      Particles particles =  apply<FinalState>(event, "FS").particles();
 75      MSG_DEBUG("Num particles = " << particles.size());
 76      PseudoJets pjs;
 77      double mpi=.13957;
 78      for (const Particle & p : particles) {
 79        Vector3 mom = p.p3();
 80        double energy = p.E();
 81        if(PID::isCharged(p.pid())) {
 82          energy = sqrt(mom.mod2()+sqr(mpi));
 83        }
 84        else {
 85          double fact = energy/mom.mod();
 86          mom *=fact;
 87        }
 88        pjs.push_back(fastjet::PseudoJet(mom.x(),mom.y(),mom.z(),energy));
 89      }
 90      // durham
 91      fastjet::JetDefinition durDef(fastjet::ee_kt_algorithm, fastjet::E_scheme);
 92      fastjet::ClusterSequence durham(pjs,durDef);
 93      double y_23 = durham.exclusive_ymerge_max(2);
 94      string label = _h["durham"]->xEdges()[durhamAxis.index(y_23)-1];
 95      _h["durham"]->fill(label);
 96      // jade e-scheme
 97      fastjet::JetDefinition::Plugin *plugin = new fastjet::JadePlugin();
 98      fastjet::JetDefinition jadeEDef(plugin);
 99      jadeEDef.set_recombination_scheme(fastjet::E_scheme);
100      fastjet::ClusterSequence jadeE(pjs,jadeEDef);
101      y_23 = jadeE.exclusive_ymerge_max(2);
102      label = _h["jade_E"]->xEdges()[jadeAxis.index(y_23)-1];
103      _h["jade_E"]->fill(label);
104      // jade p-scheme
105      fastjet::P_scheme p_scheme;
106      fastjet::JetDefinition jadePDef(plugin);
107      jadePDef.set_recombiner(&p_scheme);
108      fastjet::ClusterSequence jadeP(pjs,jadePDef);
109      y_23 = jadeP.exclusive_ymerge_max(2);
110      label = _h["jade_P"]->xEdges()[jadeAxis.index(y_23)-1];
111      _h["jade_P"]->fill(label);
112      // jade E0-scheme
113      fastjet::E0_scheme e0_scheme;
114      fastjet::JetDefinition jadeE0Def(plugin);
115      jadeE0Def.set_recombiner(&e0_scheme);
116      fastjet::ClusterSequence jadeE0(pjs,jadeE0Def);
117      y_23 = jadeE0.exclusive_ymerge_max(2);
118      label = _h["jade_E0"]->xEdges()[jadeAxis.index(y_23)-1];
119      _h["jade_E0"]->fill(label);
120    }
121
122
123    /// Normalise histograms etc., after the run
124    void finalize() {
125      scale(_h, 1./sumOfWeights());
126    }
127
128    /// @}
129
130
131    /// @name Histograms
132    /// @{
133    map<string,BinnedHistoPtr<string>> _h;
134    YODA::Axis<double> jadeAxis{0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.0825, 0.1,
135                                0.12, 0.14, 0.16, 0.1825, 0.21, 0.24, 0.27, 0.3};
136    YODA::Axis<double> durhamAxis{0.0, 0.002, 0.004, 0.0065, 0.0115, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07,
137                                  0.0825, 0.1, 0.12, 0.14, 0.16, 0.1825, 0.21, 0.24, 0.27, 0.3};
138    /// @}
139
140
141  };
142
143
144  RIVET_DECLARE_PLUGIN(AMY_1995_I406129);
145
146
147}