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, "/TMP/nPi");
 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 )  scale(_h_cTheta, fact);
 72      if (_h_cTheta2)  scale(_h_cTheta2, fact);
 73      scale(_cPi, fact);
 74    }
 75
 76    /// @}
 77
 78    string map2string(const double theta, const size_t alt) const {
 79      const size_t idx = thetaAxis.index(theta) - 1;
 80      if (idx < thetaEdges[alt].size())  return thetaEdges[alt][idx];
 81      return "OTHER";
 82    }
 83
 84    /// @name Histograms
 85    /// @{
 86    BinnedHistoPtr<string> _cPi, _h_cTheta, _h_cTheta2;
 87    vector<string> thetaEdges[2];
 88    YODA::Axis<double> thetaAxis;
 89    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,
 90                               1.025, 1.05, 1.075, 1.1, 1.125, 1.15, 1.175, 1.2, 1.225, 1.25, 1.275, 1.3,
 91                               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};
 92    string ecmEdge = "";
 93    /// @}
 94
 95  };
 96
 97
 98  RIVET_DECLARE_PLUGIN(CELLO_1992_I345437);
 99
100}