1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/Thrust.hh"
  5#include "Rivet/Projections/Beam.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief pi, kaon and proton spectra at 10.58 GeV
 11  class BELLE_2019_I1718551 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2019_I1718551);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23      // projections
 24      const FinalState fs;
 25      declare(fs, "FS");
 26      declare(Thrust(fs),"Thrust");
 27      declare(Beam(), "Beams");
 28
 29      for (size_t ix=0; ix<6; ++ix) {
 30
 31        if (ix == 1 || ix == 2) {
 32          book(_proton[ix], {0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45,
 33                             0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85});
 34        }
 35        else if (ix) {
 36          book(_proton[ix], {0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45,
 37                             0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.90});
 38        }
 39        else {
 40          book(_proton[ix], {0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.70});
 41        }
 42
 43        if (ix) {
 44          book(_pion[ix], {0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55,
 45                           0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0});
 46          book(_kaon[ix], {0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55,
 47                           0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0});
 48        }
 49        else {
 50          book(_pion[ix], {0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45,
 51                           0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85,0.90});
 52          book(_kaon[ix], {0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45,
 53                           0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85});
 54        }
 55        for (size_t iy=0; iy < _pion[ix]->numBins(); ++iy) {
 56          book(_pion[ix]->bin(iy+1), 1, ix+1, iy+1);
 57          book(_kaon[ix]->bin(iy+1), 2, ix+1, iy+1);
 58        }
 59        for (size_t iy=0; iy < _proton[ix]->numBins(); ++iy) {
 60          book(_proton[ix]->bin(iy+1), 3, ix+1, iy+1);
 61        }
 62      }
 63    }
 64
 65
 66    /// Perform the per-event analysis
 67    void analyze(const Event& event) {
 68      // Get beams and average beam momentum
 69      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 70      const double meanBeamMom = ( beams.first.p3().mod() + beams.second.p3().mod() ) / 2.0;
 71      // get the thrust
 72      const vector<double> tbins{0.7,0.8,0.85,0.9,0.95,1.0};
 73      const Thrust& thrust = apply<Thrust>(event, "Thrust");
 74      // find the thrust bin
 75      unsigned int ithrust=0;
 76      for (; ithrust<tbins.size(); ++ithrust) {
 77        if (thrust.thrust() <= tbins[ithrust]) break;
 78      }
 79      Vector3 a1 = thrust.thrustMajorAxis();
 80      Vector3 a2 = thrust.thrustMinorAxis();
 81      // loop over the charged hadrons
 82      const FinalState& fs = apply<FinalState>(event, "FS");
 83      for (const Particle& charged : fs.particles(Cuts::abspid==211 or Cuts::abspid==321 or Cuts::abspid==2212)) {
 84        double xE = charged.momentum().t()/meanBeamMom;
 85        double pT = sqrt(sqr(a1.dot(charged.momentum().p3()))+sqr(a2.dot(charged.momentum().p3())));
 86        if (charged.abspid()==211) {
 87          _pion[ithrust]->fill(xE, pT);
 88        }
 89        else if (charged.abspid()==321) {
 90          _kaon[ithrust]->fill(xE, pT);
 91        }
 92        else if (charged.abspid()==2212) {
 93          _proton[ithrust]->fill(xE, pT);
 94        }
 95      }
 96    }
 97
 98
 99    /// Normalise histograms etc., after the run
100    void finalize() {
101      // scale factor including bin width in x_E
102      const double fact = crossSection()/femtobarn/sumOfWeights()/0.05;
103      scale(_pion, fact);
104      scale(_kaon, fact);
105      scale(_proton, fact);
106    }
107
108    /// @}
109
110
111    /// @name Histograms
112    /// @{
113    Histo1DGroupPtr _pion[6], _kaon[6], _proton[6];
114    /// @}
115
116
117  };
118
119
120  RIVET_DECLARE_PLUGIN(BELLE_2019_I1718551);
121
122}