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      _h_durham->fill(y_23);
 95      // jade e-scheme
 96      fastjet::JetDefinition::Plugin *plugin = new fastjet::JadePlugin();
 97      fastjet::JetDefinition jadeEDef(plugin);
 98      jadeEDef.set_recombination_scheme(fastjet::E_scheme);
 99      fastjet::ClusterSequence jadeE(pjs,jadeEDef);
100      y_23 = jadeE.exclusive_ymerge_max(2);
101      _h_jade_E->fill(y_23);
102      // jade p-scheme
103      fastjet::P_scheme p_scheme;
104      fastjet::JetDefinition jadePDef(plugin);
105      jadePDef.set_recombiner(&p_scheme);
106      fastjet::ClusterSequence jadeP(pjs,jadePDef);
107      y_23 = jadeP.exclusive_ymerge_max(2);
108      _h_jade_P->fill(y_23);
109      // jade E0-scheme
110      fastjet::E0_scheme e0_scheme;
111      fastjet::JetDefinition jadeE0Def(plugin);
112      jadeE0Def.set_recombiner(&e0_scheme);
113      fastjet::ClusterSequence jadeE0(pjs,jadeE0Def);
114      y_23 = jadeE0.exclusive_ymerge_max(2);
115      _h_jade_E0->fill(y_23);
116    }
117
118
119    /// Normalise histograms etc., after the run
120    void finalize() {
121      scale(_h_jade_E , 1./sumOfWeights());
122      scale(_h_jade_E0, 1./sumOfWeights());
123      scale(_h_jade_P , 1./sumOfWeights());
124      scale(_h_durham , 1./sumOfWeights());
125
126    }
127
128    //@}
129
130
131    /// @name Histograms
132    //@{
133    Histo1DPtr _h_jade_P, _h_jade_E, _h_jade_E0, _h_durham; 
134    //@}
135
136
137  };
138
139
140  // The hook for the plugin system
141  RIVET_DECLARE_PLUGIN(AMY_1995_I406129);
142
143
144}