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 CELLO_1992_I345437 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(CELLO_1992_I345437);
 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.75*GeV || sqrtS() > 2*GeV)
 25        throw Error("Invalid CMS energy for CELLO_1992_I345437");
 26      int ibin = (sqrtS()-0.70)/0.05;
 27      if (ibin > 0 && ibin < 19)
 28        book(_h_cTheta,2,1,ibin);
 29      if (inRange(sqrtS()/GeV,0.85,.95))
 30        book(_h_cTheta2,2,1,19);
 31      else if (inRange(sqrtS()/GeV,1.15,1.25))
 32        book(_h_cTheta2,2,1,20);
 33      else if (inRange(sqrtS()/GeV,1.25,1.35))
 34        book(_h_cTheta2,2,1,21);
 35
 36      book(_cPi, 1, 1, 1);
 37      _thetaAxis = YODA::Axis<double>(16, 0.0, 0.8);
 38    }
 39
 40
 41    /// Perform the per-event analysis
 42    void analyze(const Event& event) {
 43
 44      if (ecmEdge == "")  ecmEdge = _cPi->bin(ecmAxis.index(sqrtS()/GeV)).xEdge();
 45      if (_h_cTheta  && thetaEdges[0].empty())  thetaEdges[0] = _h_cTheta->xEdges();
 46      if (_h_cTheta2 && thetaEdges[1].empty())  thetaEdges[1] = _h_cTheta2->xEdges();
 47
 48      Particles part = apply<FinalState>(event,"FS").particles();
 49      if (part.size() != 2) vetoEvent;
 50      double cTheta(0.);
 51      bool foundP(false), foundM(false);
 52      for (const Particle& p : part) {
 53        if (p.pid()==PID::PIPLUS) {
 54          foundP=true;
 55          cTheta = abs(p.momentum().z()/p.momentum().p3().mod());
 56        }
 57        else if (p.pid()==PID::PIMINUS) {
 58          foundM=true;
 59        }
 60      }
 61      if (!foundP || !foundM)  vetoEvent;
 62      if (cTheta <= 0.6)  _cPi->fill(ecmEdge);
 63      if (_h_cTheta )  _h_cTheta ->fill(map2string(cTheta, 0));
 64      if (_h_cTheta2)  _h_cTheta2->fill(map2string(cTheta, 1));
 65    }
 66
 67
 68    /// Normalise histograms etc., after the run
 69    void finalize() {
 70      const double fact = crossSection()/nanobarn/sumOfWeights();
 71      if (_h_cTheta ) {
 72        scale(_h_cTheta, fact);
 73        for(auto & b : _h_cTheta->bins()) {
 74          const size_t idx = b.index();
 75          b.scaleW(1./_thetaAxis.width(idx));
 76        }
 77      }
 78      if (_h_cTheta2) {
 79        scale(_h_cTheta2, fact);
 80        for(auto & b : _h_cTheta2->bins()) {
 81          const size_t idx = b.index();
 82          b.scaleW(1./_thetaAxis.width(idx));
 83        }
 84      }
 85      scale(_cPi, fact);
 86    }
 87
 88    /// @}
 89
 90    string map2string(const double theta, const size_t alt) const {
 91      const size_t idx = _thetaAxis.index(theta) - 1;
 92      if (idx < thetaEdges[alt].size())  return thetaEdges[alt][idx];
 93      return "OTHER";
 94    }
 95
 96    /// @name Histograms
 97    /// @{
 98    BinnedHistoPtr<string> _cPi, _h_cTheta, _h_cTheta2;
 99    vector<string> thetaEdges[2];
100    YODA::Axis<double> _thetaAxis;
101    YODA::Axis<double> ecmAxis{0.75, 0.775, 0.8, 0.825, 0.85, 0.875, 0.9, 0.925, 0.95, 0.975, 1.0,
102                               1.025, 1.05, 1.075, 1.1, 1.125, 1.15, 1.175, 1.2, 1.225, 1.25, 1.275, 1.3,
103                               1.325, 1.35, 1.375, 1.4, 1.425, 1.45, 1.475, 1.51875, 1.6, 1.7, 1.8, 1.9, 2.0};
104    string ecmEdge = "";
105    /// @}
106
107  };
108
109
110  RIVET_DECLARE_PLUGIN(CELLO_1992_I345437);
111
112}