1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/FinalState.hh"
 4
 5
 6namespace Rivet {
 7
 8
 9  /// @brief e+ e- > K+K-, KS0 KL0 and pi+pi-pi0
10  class SND_2001_I533574 : public Analysis {
11  public:
12
13    /// Constructor
14    RIVET_DEFAULT_ANALYSIS_CTOR(SND_2001_I533574);
15
16
17    /// @name Analysis methods
18    /// @{
19
20    /// Book histograms and initialise projections before the run
21    void init() {
22
23      // Initialise and register projections
24      declare(FinalState(), "FS");
25      for (unsigned int ix=0;ix<2;++ix) {
26        for (unsigned int iy=0;iy<4;++iy)
27          book(_sigma[ix][iy], 1+ix, 1, 1+iy);
28        
29        for (const string& en : _sigma[ix][3].binning().edges<0>()) {
30          double end = std::stod(en)*MeV;
31          if(isCompatibleWithSqrtS(end)) {
32            _ecms[ix] = en;
33            break;
34          }
35        }
36      }
37      if(_ecms[0].empty() && _ecms[1].empty())
38        MSG_ERROR("Beam energy incompatible with analysis.");
39    }
40    
41    
42    /// Perform the per-event analysis
43    void analyze(const Event& event) {
44      const FinalState& fs = apply<FinalState>(event, "FS");
45
46      map<long,int> nCount;
47      int ntotal(0);
48      for (const Particle& p : fs.particles()) {
49	nCount[p.pid()] += 1;
50	++ntotal;
51      }
52      if(ntotal==2) {
53	if(nCount[321]==1 && nCount[-321]==1) {
54          for(unsigned int ix=0;ix<2;++ix)
55            _sigma[ix][0]->fill(_ecms[ix]);
56        }
57        else if(nCount[130]==1 && nCount[310]==1) {
58          for(unsigned int ix=0;ix<2;++ix) {
59            _sigma[ix][1]->fill(_ecms[ix]);
60            _sigma[ix][2]->fill(_ecms[ix]);
61          }
62        }
63      }
64      else if(ntotal==3 && nCount[211] == 1 && nCount[-211] == 1 && nCount[111] == 1) {
65        for(unsigned int ix=0;ix<2;++ix)
66          _sigma[ix][3]->fill(_ecms[ix]);
67      }
68
69    }
70
71
72    /// Normalise histograms etc., after the run
73    void finalize() {
74      double fact = crossSection()/ sumOfWeights() /nanobarn;
75      for (unsigned int ix=0;ix<2;++ix) {
76        for (unsigned int iy=0;iy<4;++iy) {
77          scale(_sigma[ix][iy],fact);
78        }
79      }
80    }
81
82    /// @}
83
84
85    /// @name Histograms
86    /// @{
87    BinnedHistoPtr<string> _sigma[2][4];
88    string _ecms[2];
89    /// @}
90
91
92  };
93
94
95  RIVET_DECLARE_PLUGIN(SND_2001_I533574);
96
97
98}