1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/UnstableParticles.hh"
 4#include "Rivet/Projections/ChargedFinalState.hh"
 5
 6namespace Rivet {
 7
 8
 9  /// @brief Xi- and Sigma*+/- spectra
10  class HRS_1987_I246162 : public Analysis {
11  public:
12
13    /// Constructor
14    RIVET_DEFAULT_ANALYSIS_CTOR(HRS_1987_I246162);
15
16
17    /// @name Analysis methods
18    /// @{
19
20    /// Book histograms and initialise projections before the run
21    void init() {
22      // projects
23      declare(UnstableParticles(Cuts::abspid==3114 || Cuts::abspid==3224 ||
24                                Cuts::abspid==3312 || Cuts::abspid==3122), "UFS");
25      const ChargedFinalState cfs;
26      declare(cfs, "CFS");
27      // histos
28      for (unsigned int ix=0; ix<2; ++ix) {
29        for(unsigned int iy=0;iy<2;++iy)
30          book(_h_total[ix][iy],1+ix,1,1+iy);
31        book(_h_x[ix],4,1,1+ix);
32        book(_h_ratio[ix],3,1,1+ix);
33      }
34      book(_c,"TMP/nLam");
35    }
36
37
38    /// Perform the per-event analysis
39    void analyze(const Event& event) {
40      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
41      const size_t numParticles = cfs.particles().size();
42
43      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
44      if (numParticles < 2) {
45        MSG_DEBUG("Failed leptonic event cut");
46        vetoEvent;
47      }
48      MSG_DEBUG("Passed leptonic event cut");
49      for(const Particle & p : apply<UnstableParticles>(event,"UFS").particles()) {
50        if(p.abspid()==3122) {
51          _c->fill();
52          continue;
53        }
54      	double xE = 2.*p.E()/sqrtS();
55      	Vector3 mom3 = p.p3();
56        const double energy = p.E();
57      	double modp = mom3.mod();
58      	double beta = modp/energy;
59        unsigned int iloc = p.abspid()==3312 ? 1 : 0;
60        _h_ratio[iloc]->fill(29);
61        if(iloc==0&&xE>0.1&&xE<0.8) 
62          _h_total[iloc][0]->fill(29);
63        else if(iloc==1&&xE>0.05&&xE<0.75)
64          _h_total[iloc][0]->fill(29);
65        _h_total[iloc][1]->fill(29);
66        _h_x[iloc]->fill(xE,1./beta);
67      }
68    }
69
70
71    /// Normalise histograms etc., after the run
72    void finalize() {
73      for (unsigned int ix=0; ix<2; ++ix) {
74        scale(_h_total[ix],1./sumOfWeights());
75        scale(_h_x[ix], sqr(sqrtS())*crossSection()/nanobarn/sumOfWeights());
76        scale(_h_ratio[ix], 1./ *_c);
77      }
78    }
79
80    /// @}
81
82
83    /// @name Histograms
84    /// @{
85    BinnedHistoPtr<int> _h_total[2][2],_h_ratio[2];
86    Histo1DPtr _h_x[2];
87    CounterPtr _c;
88    /// @}
89
90
91  };
92
93
94  RIVET_DECLARE_PLUGIN(HRS_1987_I246162);
95
96}