1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief B -> X_s gamma spectrum
  9  class BELLE_2015_I1330289 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2015_I1330289);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21
 22      // Initialise and register projections
 23      declare(UnstableParticles(), "UFS");
 24
 25      // Book histograms
 26      book(_h_spectrum, 1, 1, 2);
 27      book(_nBottom, "TMP/BottomCounter");
 28    }
 29
 30
 31    void findDecayProducts(const Particle& mother,
 32                           unsigned int& nK0, unsigned int& nKp, unsigned int& nKm,
 33                           FourMomentum& ptot) {
 34      for (const Particle & p : mother.children()) {
 35        int id = p.pid();
 36        if ( id == PID::KPLUS ) {
 37          ++nKp;
 38          ptot += p.momentum();
 39        }
 40        else if (id == PID::KMINUS ) {
 41          ++nKm;
 42          ptot += p.momentum();
 43        }
 44        else if (id == PID::K0S) {
 45          ++nK0;
 46          ptot += p.momentum();
 47        }
 48        else if (id == PID::PI0 || id == PID::PIPLUS || id == PID::PIMINUS) {
 49          ptot += p.momentum();
 50        }
 51        else if ( !p.children().empty() ) {
 52          findDecayProducts(p, nK0, nKp, nKm, ptot);
 53        }
 54        else
 55          ptot += p.momentum();
 56      }
 57    }
 58
 59
 60    /// Perform the per-event analysis
 61    void analyze(const Event& event) {
 62
 63      // Loop over bottoms
 64      for (const Particle& bottom : apply<UnstableParticles>(event, "UFS").particles()) {
 65        if (bottom.pid() != -521 && bottom.pid() != -511) continue;
 66        FourMomentum pgamma(0.,0.,0.,0.);
 67        unsigned int ngamma = 0;
 68        bool fs = true;
 69        for (const Particle & child : bottom.children()) {
 70          if (child.pid() == bottom.pid()) {
 71            fs = false;
 72            break;
 73          }
 74          else if (child.pid() == PID::PHOTON) {
 75            ngamma += 1;
 76            pgamma += child.momentum();
 77          }
 78        }
 79        if (!fs) continue;
 80        _nBottom->fill();
 81        if (ngamma != 1) continue;
 82        unsigned int nK0(0),nKp(0),nKm(0);
 83        FourMomentum p_tot(0,0,0,0);
 84        findDecayProducts(bottom, nK0, nKp, nKm, p_tot);
 85        unsigned int nk = nKp-nKm+nK0;
 86        if (nk % 2 == 1) {
 87          p_tot -= pgamma;
 88          _h_spectrum->fill(p_tot.mass()/GeV);
 89        }
 90      }
 91
 92    }
 93
 94
 95    /// Normalise histograms etc., after the run
 96    void finalize() {
 97      scale(_h_spectrum, 1e6/_nBottom->sumW());
 98      // multiply by the bin width
 99      for (size_t ix = 0; ix < _h_spectrum->numBins(); ++ix) {
100        _h_spectrum->bin(ix+1).scaleW(_h_spectrum->bin(ix+1).xWidth());
101      }
102    }
103
104    /// @}
105
106
107    /// @name Histograms
108    /// @{
109    Histo1DPtr _h_spectrum;
110    CounterPtr _nBottom;
111    /// @}
112
113  };
114
115
116
117  RIVET_DECLARE_PLUGIN(BELLE_2015_I1330289);
118
119}