1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/Beam.hh"
 4#include "Rivet/Projections/Thrust.hh"
 5#include "Rivet/Projections/UnstableParticles.hh"
 6#include "Rivet/Projections/ChargedFinalState.hh"
 7
 8namespace Rivet {
 9
10
11  /// @brief Lambda production at 29 GeV
12  class HRS_1986_I17781 : public Analysis {
13  public:
14
15    /// Constructor
16    RIVET_DEFAULT_ANALYSIS_CTOR(HRS_1986_I17781);
17
18
19    /// @name Analysis methods
20    ///@{
21
22    /// Book histograms and initialise projections before the run
23    void init() {
24      declare(UnstableParticles(), "UFS");
25      const ChargedFinalState cfs;
26      declare(cfs, "CFS");
27      declare(Thrust(cfs), "Thrust");
28      book(_h_spect,1,1,1);
29      book(_h_rap  ,2,1,1);
30      book(_h_mult ,3,1,1);
31    }
32
33
34    /// Perform the per-event analysis
35    void analyze(const Event& event) {
36      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
37      const size_t numParticles = cfs.particles().size();
38
39      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
40      if (numParticles < 2) {
41        MSG_DEBUG("Failed leptonic event cut");
42        vetoEvent;
43      }
44      MSG_DEBUG("Passed leptonic event cut");
45      // Thrust
46      const Thrust& thrust = apply<Thrust>(event, "Thrust");
47      UnstableParticles ufs = apply<UnstableParticles>(event,"UFS");
48      for (const Particle & p : ufs.particles(Cuts::abspid==3122)) {
49      	double xE = 2.*p.E()/sqrtS();
50      	Vector3 mom3 = p.p3();
51        const double energy = p.E();
52      	double modp = mom3.mod();
53      	double beta = modp/energy;
54        const double momT = dot(thrust.thrustAxis(), mom3);
55        const double rapidityT = 0.5 * std::log((energy + momT) / (energy - momT));
56        _h_spect->fill(xE,1./beta);
57        _h_rap->fill(abs(rapidityT));
58        _h_mult->fill(29);
59      }
60    }
61
62
63    /// Normalise histograms etc., after the run
64    void finalize() {
65      scale( _h_spect, sqr(sqrtS())*crossSection()/nanobarn/sumOfWeights());
66      scale( _h_rap  , 1./sumOfWeights());
67      scale( _h_mult , 1./sumOfWeights());
68    }
69
70    ///@}
71
72
73    /// @name Histograms
74    ///@{
75    Histo1DPtr _h_spect,_h_rap;
76    BinnedHistoPtr<int> _h_mult;
77    ///@}
78
79
80  };
81
82
83  RIVET_DECLARE_PLUGIN(HRS_1986_I17781);
84
85}