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 PLUTO_1984_I204487 : public Analysis {
10  public:
11
12    /// Constructor
13    RIVET_DEFAULT_ANALYSIS_CTOR(PLUTO_1984_I204487);
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      if (sqrtS()<0.36*GeV || sqrtS()>1.72*GeV)
25        throw Error("Invalid CMS energy for PLUTO_1984_I1260740");
26      book(_cPi, "/TMP/nPi");
27    }
28
29
30    /// Perform the per-event analysis
31    void analyze(const Event& event) {
32      Particles part = apply<FinalState>(event,"FS").particles();
33      if (part.size()!=2) vetoEvent;
34      double cTheta(0.);
35      bool foundP(false),foundM(false);
36      for (const Particle& p : part) {
37        if (p.pid()==PID::PIPLUS) {
38          foundP=true;
39          cTheta = abs(p.momentum().z()/p.momentum().p3().mod());
40        }
41        else if (p.pid()==PID::PIMINUS)
42          foundM=true;
43      }
44      if (!foundP || !foundM) vetoEvent;
45      if (cTheta<=0.2)    _cPi->fill();
46    }
47
48
49    /// Normalise histograms etc., after the run
50    void finalize() {
51      double fact = crossSection()/nanobarn/sumOfWeights();
52      double sigma = _cPi->val()*fact/0.2;
53      double error = _cPi->err()*fact/0.2;
54      Estimate1DPtr mult;
55      book(mult, 1, 1, 1);
56      for (auto& b : mult->bins()) {
57        if (inRange(sqrtS(), b.xMin(), b.xMax())) {
58          b.set(sigma, error);
59        }
60      }
61    }
62
63    ///@}
64
65
66    /// @name Histograms
67    ///@{
68    CounterPtr _cPi;
69    ///@}
70
71
72  };
73
74
75  RIVET_DECLARE_PLUGIN(PLUTO_1984_I204487);
76
77}