1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/UnstableParticles.hh"
 4#include "Rivet/Projections/FinalState.hh"
 5#include "Rivet/Projections/ChargedFinalState.hh"
 6#include "Rivet/Projections/Thrust.hh"
 7
 8namespace Rivet {
 9
10
11  /// @brief Lambda0 production at 29 GeV
12  class MARKII_1985_I209198 : public Analysis {
13  public:
14
15    /// Constructor
16    RIVET_DEFAULT_ANALYSIS_CTOR(MARKII_1985_I209198);
17
18
19    /// @name Analysis methods
20    ///@{
21
22    /// Book histograms and initialise projections before the run
23    void init() {
24      // Initialise and register projections
25      declare(UnstableParticles(), "UFS");
26      const ChargedFinalState cfs;
27      declare(cfs, "CFS");
28      const FinalState fs;
29      declare(cfs, "FS");
30      declare(Thrust(fs), "Thrust");
31      //Histograms
32      book(_h_spect, 2, 1, 1);
33      book(_h_pT   , 3, 1, 1);
34    }
35
36
37    /// Perform the per-event analysis
38    void analyze(const Event& event) {
39      if (_edges.empty())  _edges = _h_spect->xEdges();
40      // 5 charged particles
41      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
42      if(cfs.particles().size()<5) vetoEvent;
43      // thrust
44      const Thrust& thrust = apply<Thrust>(event, "Thrust");
45      // lambdas
46      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
47      for (const Particle& p : ufs.particles(Cuts::abspid==3122)) {
48        const double xp = 2.*p.E()/sqrtS();
49      	Vector3 mom3 = p.p3();
50        const double beta = mom3.mod() / p.E();
51        const double pTin  = dot(mom3, thrust.thrustMajorAxis());
52        const double pTout = dot(mom3, thrust.thrustMinorAxis());
53        const double pT2 = sqr(pTin)+sqr(pTout);
54        const size_t idx = _axis.index(xp);
55        string edge = "OTHER";
56        if (idx && idx <= _axis.numBins())  edge = _edges[idx-1];
57        _h_spect->fill(edge, 1./beta);
58        _h_pT   ->fill(pT2);
59      }
60    }
61
62
63    /// Normalise histograms etc., after the run
64    void finalize() {
65      scale(_h_spect, sqr(sqrtS())*crossSection()/nanobarn/sumOfWeights());
66      normalize(_h_pT);
67      for(auto & b: _h_spect->bins()) {
68        const size_t idx = b.index();
69        b.scaleW(1./_axis.width(idx));
70      }
71    }
72
73    ///@}
74
75
76    /// @name Histograms
77    ///@{
78    BinnedHistoPtr<string> _h_spect;
79    Histo1DPtr _h_pT;
80    vector<string> _edges;
81    YODA::Axis<double> _axis{0.084,0.094,0.104,0.124,0.157,0.190,0.223,0.256,
82                             0.289,0.322,0.355,0.423,0.491,0.559,0.627,0.695};
83    ///@}
84
85
86  };
87
88
89  RIVET_DECLARE_PLUGIN(MARKII_1985_I209198);
90
91}