1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Add a short analysis description here
  9  class BELLE_2013_I1238273 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2013_I1238273);
 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_q2_B0bar_pi     ,1, 1, 1);
 27      book(_h_q2_B0bar_rho    ,3, 1, 1);
 28      book(_h_q2_Bminus_pi    ,2, 1, 1);
 29      book(_h_q2_Bminus_rho   ,4, 1, 1);
 30      book(_h_q2_Bminus_omega ,5, 1, 1);
 31
 32    }
 33
 34    // Calculate the Q2 using mother and daugher meson
 35    double q2(const Particle& B, int mesonID) {
 36      FourMomentum q = B.mom() - filter_select(B.children(), Cuts::pid==mesonID)[0];
 37      return q*q;
 38    }
 39
 40    // Check for explicit decay into pdgids
 41    bool isSemileptonicDecay(const Particle& mother, vector<int> ids) {
 42      // Trivial check to ignore any other decays but the one in question modulo photons
 43      const Particles children = mother.children(Cuts::pid!=PID::PHOTON);
 44      if (children.size()!=ids.size()) return false;
 45      // Check for the explicit decay
 46      return all(ids, [&](int i){return count(children, hasPID(i))==1;});
 47    }
 48
 49    /// Perform the per-event analysis
 50    void analyze(const Event& event) {
 51      // Loop over B0bar Mesons
 52      for(const Particle& p : apply<UnstableParticles>(event, "UFS").particles(Cuts::pid==PID::B0BAR)) {
 53        if (isSemileptonicDecay(p, {PID::PIPLUS, PID::ELECTRON, PID::NU_EBAR}) ||
 54            isSemileptonicDecay(p, {PID::PIPLUS, PID::MUON,     PID::NU_MUBAR})) {
 55            _h_q2_B0bar_pi->fill(q2(p, PID::PIPLUS));
 56        }
 57        if (isSemileptonicDecay(p, {PID::RHOPLUS, PID::ELECTRON, PID::NU_EBAR}) ||
 58            isSemileptonicDecay(p, {PID::RHOPLUS, PID::MUON,     PID::NU_MUBAR})) {
 59            _h_q2_B0bar_rho->fill(q2(p, PID::RHOPLUS));
 60        }
 61      }
 62      // Loop over B- Mesons
 63      for(const Particle& p : apply<UnstableParticles>(event, "UFS").particles(Cuts::pid==PID::BMINUS)) {
 64        if (isSemileptonicDecay(p, {PID::PI0, PID::ELECTRON, PID::NU_EBAR}) ||
 65            isSemileptonicDecay(p, {PID::PI0, PID::MUON,     PID::NU_MUBAR})) {
 66            _h_q2_Bminus_pi->fill(q2(p, PID::PI0));
 67        }
 68        if (isSemileptonicDecay(p, {PID::RHO0, PID::ELECTRON, PID::NU_EBAR}) ||
 69            isSemileptonicDecay(p, {PID::RHO0, PID::MUON,    PID::NU_MUBAR})) {
 70            _h_q2_Bminus_rho->fill(q2(p,PID::RHO0));
 71        }
 72        if (isSemileptonicDecay(p, {PID::OMEGA, PID::ELECTRON, PID::NU_EBAR}) ||
 73            isSemileptonicDecay(p, {PID::OMEGA, PID::MUON,     PID::NU_MUBAR})) {
 74            _h_q2_Bminus_omega->fill(q2(p, PID::OMEGA));
 75        }
 76      }
 77    }
 78
 79
 80    /// Normalise histograms etc., after the run
 81    void finalize() {
 82
 83      normalize(_h_q2_B0bar_pi    , 298.8);  // normalize to BF*dQ2
 84      normalize(_h_q2_B0bar_rho   , 1304.8); // normalize to BF*dQ2
 85      normalize(_h_q2_Bminus_pi   , 324.8);  // normalize to BF*dQ2
 86      normalize(_h_q2_Bminus_rho  , 367.0);  // normalize to BF*dQ2
 87      normalize(_h_q2_Bminus_omega, 793.1);  // normalize to BF*dQ2
 88
 89    }
 90
 91    //@}
 92
 93
 94  private:
 95
 96
 97    /// @name Histograms
 98    //@{
 99    Histo1DPtr _h_q2_B0bar_pi    ;
100    Histo1DPtr _h_q2_B0bar_rho   ;
101    Histo1DPtr _h_q2_Bminus_pi   ;
102    Histo1DPtr _h_q2_Bminus_rho  ;
103    Histo1DPtr _h_q2_Bminus_omega;
104    //@}
105
106
107  };
108
109
110  // The hook for the plugin system
111  RIVET_DECLARE_PLUGIN(BELLE_2013_I1238273);
112
113
114}