1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/FinalState.hh"
 4
 5namespace Rivet {
 6
 7
 8  /// @brief gamma gamma -> p pbar
 9  class TPC_1987_I246557 : public Analysis {
10  public:
11
12    /// Constructor
13    RIVET_DEFAULT_ANALYSIS_CTOR(TPC_1987_I246557);
14
15
16    /// @name Analysis methods
17    /// @{
18
19    /// Book histograms and initialise projections before the run
20    void init() {
21      // Initialise and register projections
22      declare(FinalState(), "FS");
23      // book histos
24      if (inRange(sqrtS()/GeV,2.,2.8)) {
25        if (inRange(sqrtS()/GeV,2.1,2.4)) {
26          book(_h_cTheta,3,1,1);
27        }
28        else if(sqrtS()<=2.8) {
29          book(_h_cTheta,3,1,2);
30        }
31        book(_cProton, 1, 1, 1);
32      }
33      else {
34        throw Error("Invalid CMS energy for TPC_1987_I246557");
35      }
36
37      // Determine which sqrtS bin needs to be filled
38      _edge = "OTHER"s;
39      for (const string& edge : _cProton.binning().edges<0>()) {
40        if (isCompatibleWithSqrtS(std::stod(edge)*GeV)) {
41          _edge = edge;
42          break;
43        }
44      }
45    }
46
47
48    /// Perform the per-event analysis
49    void analyze(const Event& event) {
50      Particles part = apply<FinalState>(event,"FS").particles();
51      if (part.size()!=2) vetoEvent;
52      double cTheta(0.);
53      bool foundP(false), foundM(false);
54      for (const Particle& p : part) {
55      if (p.pid()==PID::PROTON) {
56        foundP=true;
57        cTheta = abs(p.momentum().z()/p.momentum().p3().mod());
58      }
59      else if (p.pid()==PID::ANTIPROTON)
60        foundM=true;
61      }
62      if (!foundP || !foundM) vetoEvent;
63      if (cTheta<=0.6)    _cProton->fill(_edge);
64      if (_h_cTheta) _h_cTheta->fill(cTheta);
65    }
66
67
68    /// Normalise histograms etc., after the run
69    void finalize() {
70      const double fact = crossSection()/nanobarn/sumOfWeights();
71      if (_h_cTheta) scale(_h_cTheta,fact);
72      scale(_cProton, fact);
73    }
74
75    /// @}
76
77
78    /// @name Histograms
79    /// @{
80    Histo1DPtr _h_cTheta;
81    BinnedHistoPtr<string> _cProton;
82    string _edge;
83    /// @}
84
85
86  };
87
88
89  RIVET_DECLARE_PLUGIN(TPC_1987_I246557);
90
91}