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 BELLE_2005_I677625 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2005_I677625);
 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.025,4.)) {
 25        book(_sigmaProton,"TMP/nProton",refData(1,1,1));
 26        if (inRange(sqrtS()/GeV, 2.075, 2.9) || inRange(sqrtS()/GeV, 3.1, 4.)) {
 27          double sMin=2.075, sMax=2.1, step=0.1;
 28          unsigned int ihist=3,iy=1;
 29          while (sMin<4.) {
 30            if (inRange(sqrtS()/GeV, sMin, sMax)) break;
 31            sMin=sMax;
 32            sMax+=step;
 33            iy+=1;
 34            if (iy==4) {
 35              ihist+=1;
 36              iy=1;
 37            }
 38            if (fuzzyEquals(2.9, sMin)) {
 39              sMin=3.1;
 40              sMax=3.5;
 41              step=0.5;
 42            }
 43          }
 44          book(_h_cTheta[0],ihist,1,iy);
 45        }
 46        if (inRange(sqrtS()/GeV, 2.075, 2.5)) {
 47          book(_h_cTheta[1],2,1,1);
 48        }
 49        else if (inRange(sqrtS()/GeV, 2.5, 3.)) {
 50          book(_h_cTheta[1],2,1,2);
 51        }
 52        else if (inRange(sqrtS()/GeV, 3., 4.)) {
 53          book(_h_cTheta[1],2,1,3);
 54        }
 55      }
 56      else {
 57        throw Error("Invalid CMS energy for BELLE_2005_I677625");
 58      }
 59    }
 60
 61
 62    /// Perform the per-event analysis
 63    void analyze(const Event& event) {
 64      Particles part = apply<FinalState>(event,"FS").particles();
 65      if (part.size()!=2) vetoEvent;
 66      double cTheta(0.);
 67      bool foundP(false),foundM(false);
 68      for (const Particle & p : part) {
 69        if (p.pid()==PID::PROTON) {
 70          foundP=true;
 71          cTheta = abs(p.mom().z()/p.mom().p3().mod());
 72        }
 73        else if (p.pid()==PID::ANTIPROTON) {
 74          foundM=true;
 75        }
 76      }
 77      if (!foundP || !foundM) vetoEvent;
 78      if (cTheta<=0.6)    _sigmaProton->fill(sqrtS());
 79      if (_h_cTheta[0] ) _h_cTheta[0]->fill(cTheta);
 80      if (_h_cTheta[1] ) _h_cTheta[1]->fill(cTheta);
 81    }
 82
 83
 84    /// Normalise histograms etc., after the run
 85    void finalize() {
 86      const double fact = crossSection()/nanobarn/sumOfWeights();
 87      for (unsigned int ix=0; ix<2; ++ix) {
 88        if (_h_cTheta[ix] ) scale(_h_cTheta[ix], fact);
 89      }
 90      scale(_sigmaProton,fact);
 91      Estimate1DPtr tmp;
 92      book(tmp,1,1,1);
 93      barchart(_sigmaProton,tmp);
 94    }
 95
 96    /// @}
 97
 98
 99    /// @name Histograms
100    /// @{
101    Histo1DPtr _h_cTheta[2];
102    Histo1DPtr _sigmaProton;
103    /// @}
104
105
106  };
107
108
109  RIVET_DECLARE_PLUGIN(BELLE_2005_I677625);
110
111}