1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// gamma gamma -> pi+pi-/K+ K- for energies between 2.4 and 4.1 GeV
  9  class BELLE_2005_I667712 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2005_I667712);
 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()<2.4*GeV || sqrtS()>4.1*GeV)
 25        throw Error("Invalid CMS energy for BELLE_2005_I667712");
 26      for (unsigned int ix=0;ix<7;++ix) {
 27        std::ostringstream title;
 28        title << "/TMP/nPi_" << ix;
 29        book(_cPi[ix], title.str());
 30      }
 31      for (unsigned int ix=0;ix<7;++ix) {
 32        std::ostringstream title;
 33        title << "/TMP/nK_" << ix;
 34        book(_cK[ix], title.str());
 35      }
 36    }
 37
 38
 39    /// Perform the per-event analysis
 40    void analyze(const Event& event) {
 41      Particles part = apply<FinalState>(event,"FS").particles();
 42      if (part.size() != 2) vetoEvent;
 43      if (part[0].pid() != -part[1].pid()) vetoEvent;
 44      double cTheta(0.);
 45      bool foundPi(false),foundK(false);
 46      for (const Particle& p : part) {
 47        if (p.pid() == PID::PIPLUS) {
 48          foundPi=true;
 49          cTheta = abs(p.momentum().z()/p.momentum().p3().mod());
 50        }
 51        else if (p.pid() == PID::KPLUS) {
 52          foundK=true;
 53          cTheta = abs(p.momentum().z()/p.momentum().p3().mod());
 54        }
 55      }
 56      if (!foundPi && !foundK) vetoEvent;
 57      int ibin = cTheta/0.1;
 58      if (ibin>5) vetoEvent;
 59      if (foundPi) {
 60        _cPi[   0  ]->fill();
 61        _cPi[ibin+1]->fill();
 62      }
 63      else if (foundK) {
 64        _cK[   0  ]->fill();
 65        _cK[ibin+1]->fill();
 66      }
 67    }
 68
 69
 70    /// Normalise histograms etc., after the run
 71    void finalize() {
 72      double fact = crossSection()/nanobarn/sumOfWeights();
 73      for (unsigned int ip=0;ip<2;++ip) {
 74        CounterPtr denom = ip==0 ? _cPi[0] : _cK[0];
 75        if (denom->numEntries()==0) continue;
 76        for (unsigned int ih=0;ih<7;++ih) {
 77          CounterPtr numer = ip==0 ? _cPi[ih] : _cK[ih];
 78          double sigma = numer->val()*fact;
 79          double error = numer->err()*fact;
 80          // bin width for 2d dist
 81          if (ih!=0) {
 82            sigma /=0.1;
 83            error /=0.1;
 84          }
 85          unsigned int ix=5+ip, iy=ih;
 86          if (ih==0) {
 87            ix=1;
 88            iy = ip==0 ? 2 : 1;
 89          }
 90          // ratio
 91          Estimate0D temp;
 92          Estimate1DPtr cross, ratio;
 93          book(cross, ix, 1, iy);
 94          if (ih!=0) {
 95            book(ratio,ix-2,1,iy);
 96            temp = *numer / *denom;
 97          }
 98          for (auto& b : cross->bins()) {
 99            if (inRange(sqrtS(), b.xMin(), b.xMax())) {
100              b.set(sigma, error);
101              if (ih!=0) {
102                ratio->bin(b.index()).set(temp.val()/0.1, temp.errPos()/0.1);
103              }
104            }
105          }
106        }
107      }
108
109    }
110
111    ///@}
112
113
114    /// @name Histograms
115    ///@{
116    CounterPtr _cPi[7],_cK[7];
117    ///@}
118
119  };
120
121
122
123  RIVET_DECLARE_PLUGIN(BELLE_2005_I667712);
124
125}