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      size_t idx = durhamAxis.index(y_23);
 95      string label = idx>0 && idx <= _h["durham"]->xEdges().size() ?  _h["durham"]->xEdges()[idx-1] : "OTHER";
 96      _h["durham"]->fill(label);
 97      // jade e-scheme
 98      fastjet::JetDefinition::Plugin *plugin = new fastjet::JadePlugin();
 99      fastjet::JetDefinition jadeEDef(plugin);
100      jadeEDef.set_recombination_scheme(fastjet::E_scheme);
101      fastjet::ClusterSequence jadeE(pjs,jadeEDef);
102      y_23 = jadeE.exclusive_ymerge_max(2);
103      idx = jadeAxis.index(y_23);
104      label = idx>0 && idx <= _h["jade_E"]->xEdges().size() ?  _h["jade_E"]->xEdges()[idx-1] : "OTHER";
105      _h["jade_E"]->fill(label);
106      // jade p-scheme
107      fastjet::P_scheme p_scheme;
108      fastjet::JetDefinition jadePDef(plugin);
109      jadePDef.set_recombiner(&p_scheme);
110      fastjet::ClusterSequence jadeP(pjs,jadePDef);
111      y_23 = jadeP.exclusive_ymerge_max(2);
112      idx = jadeAxis.index(y_23);
113      label = idx>0 && idx <= _h["jade_P"]->xEdges().size() ?  _h["jade_P"]->xEdges()[idx-1] : "OTHER";
114      _h["jade_P"]->fill(label);
115      // jade E0-scheme
116      fastjet::E0_scheme e0_scheme;
117      fastjet::JetDefinition jadeE0Def(plugin);
118      jadeE0Def.set_recombiner(&e0_scheme);
119      fastjet::ClusterSequence jadeE0(pjs,jadeE0Def);
120      y_23 = jadeE0.exclusive_ymerge_max(2);
121      idx = jadeAxis.index(y_23);
122      label = idx>0 && idx <= _h["jade_E0"]->xEdges().size() ?  _h["jade_E0"]->xEdges()[idx-1] : "OTHER";
123      _h["jade_E0"]->fill(label);
124    }
125
126
127    /// Normalise histograms etc., after the run
128    void finalize() {
129      scale(_h, 1./sumOfWeights());
130      for( auto & hist : _h) {
131        for(auto & b: hist.second->bins()) {
132          const size_t idx = b.index();
133          b.scaleW(hist.first=="durham" ? 1./durhamAxis.width(idx) : 1./jadeAxis.width(idx));
134        }
135      }
136    }
137
138    /// @}
139
140
141    /// @name Histograms
142    /// @{
143    map<string,BinnedHistoPtr<string>> _h;
144    YODA::Axis<double> jadeAxis{0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.0825, 0.1,
145                                0.12, 0.14, 0.16, 0.1825, 0.21, 0.24, 0.27, 0.3};
146    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,
147                                  0.0825, 0.1, 0.12, 0.14, 0.16, 0.1825, 0.21, 0.24, 0.27, 0.3};
148    /// @}
149
150
151  };
152
153
154  RIVET_DECLARE_PLUGIN(AMY_1995_I406129);
155
156
157}