1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/FinalState.hh"
 4
 5namespace Rivet {
 6
 7
 8  /// @brief gamma gamma -> pi+pi-
 9  class BELLE_2007_I749358 : public Analysis {
10  public:
11
12    /// Constructor
13    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2007_I749358);
14
15
16    /// @name Analysis methods
17    ///@{
18
19    /// Book histograms and initialise projections before the run
20    void init() {
21      // Final state
22      declare(FinalState(),"FS");
23      // check CMS energy in range
24      _energyAxis = YODA::Axis<double>(140, 0.8, 1.5);
25      _thetaAxis = YODA::Axis<double>(12, 0.0, 0.6);
26      if (sqrtS()<0.8*GeV || sqrtS()>1.5*GeV)
27        throw Error("Invalid CMS energy for BELLE_2007_I749358");
28      eIndex = _energyAxis.index(sqrtS()/GeV) - 1;
29
30      // bin for the angle plots
31      int ibin = (sqrtS()-0.8)/0.005 + 2;
32      book(_h_cTheta,ibin,1,1);
33      book(_cPi, 1, 1, 1);
34    }
35
36
37    /// Perform the per-event analysis
38    void analyze(const Event& event) {
39      if (_EcmEdges.empty())  _EcmEdges = _cPi->xEdges();
40      if (_thetaEdges.empty())  _thetaEdges = _h_cTheta->xEdges();
41
42      Particles part = apply<FinalState>(event,"FS").particles();
43      if (part.size()!=2) vetoEvent;
44      double cTheta(0.);
45      bool foundP(false),foundM(false);
46      for (const Particle& p : part) {
47        if (p.pid()==PID::PIPLUS) {
48          foundP=true;
49          cTheta = abs(p.momentum().z()/p.momentum().p3().mod());
50        }
51        else if (p.pid()==PID::PIMINUS)
52          foundM=true;
53      }
54      if (!foundP || !foundM)  vetoEvent;
55      if (cTheta<=0.6)  _cPi->fill(_EcmEdges[eIndex]);
56      if (_h_cTheta ) {
57        size_t idx = _thetaAxis.index(cTheta) - 1;
58        _h_cTheta->fill(idx < _thetaEdges.size() ? _thetaEdges[idx] : "OTHER"s);
59      }
60    }
61
62
63    /// Normalise histograms etc., after the run
64    void finalize() {
65      const double fact = crossSection()/nanobarn/sumOfWeights();
66      if (_h_cTheta ) {
67        scale(_h_cTheta ,fact);
68        for(auto & b : _h_cTheta->bins()) {
69          const size_t idx = b.index();
70          b.scaleW(1./_thetaAxis.width(idx));
71        }
72      }
73      scale(_cPi, fact);
74    }
75
76    ///@}
77
78
79    /// @name Histograms
80    ///@{
81    BinnedHistoPtr<string> _h_cTheta, _cPi;
82    YODA::Axis<double> _energyAxis;
83    YODA::Axis<double> _thetaAxis;
84    vector<string> _EcmEdges, _thetaEdges;
85    size_t eIndex;
86    ///@}
87
88
89  };
90
91
92  RIVET_DECLARE_PLUGIN(BELLE_2007_I749358);
93
94}