Rivet analyses

title: BELLE_2009_I825222

$B\to X_s\gamma$ with different photon energy cuts

Experiment: BELLE (KEKB)

Inspire ID: 825222

Status: VALIDATED NOHEPDATA SINGLEWEIGHT

Authors: - Peter Richardson

References: - Phys.Rev.Lett. 103 (2009) 241801

Beams: * *

Beam energies: ANY

Run details: - Any process producing B+ and B0 mesons, originally Upsilon(4S) decay

Measurement of the branching ratios, average photon energies and photon energy dispersion for $B\to X_s\gamma$ with different photon energy cuts.

Source code:BELLE_2009_I825222.cc

```c++ // -- C++ --

include "Rivet/Analysis.hh"

include "Rivet/Projections/UnstableParticles.hh"

namespace Rivet {

/// @brief B -> Xs gamma class BELLE_2009_I825222 : public Analysis { public:

/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2009_I825222);


/// @name Analysis methods
/// @{

/// Book histograms and initialise projections before the run
void init() {
  // Initialise and register projections
  declare(UnstableParticles(Cuts::abspid == 521 || Cuts::abspid == 511), "UFS");
  // Book histograms
  book(_h_br, 1, 1, 1);
  book(_p_E, 1, 1, 2);
  book(_p_E2, "TMP/E2", refData<YODA::BinnedEstimate<string>>(1, 1, 3));
  book(_nBottom, "TMP/BottomCounter");
}

void findDecayProducts(const Particle& mother, unsigned int& nK0, unsigned int& nKp, unsigned int& nKm) {
  for (const Particle& p : mother.children()) {
    int id = p.pid();
    if (id == PID::KPLUS)
      ++nKp;
    else if (id == PID::KMINUS)
      ++nKm;
    else if (id == PID::K0S)
      ++nK0;
    else if (id == PID::PI0 || id == PID::PIPLUS || id == PID::PIMINUS) {
      continue;
    }
    else if (!p.children().empty()) {
      findDecayProducts(p, nK0, nKp, nKm);
    }
  }
}

/// Perform the per-event analysis
void analyze(const Event& event) {
  if (_edges.empty()) {
    _edges = _h_br->xEdges();
    for (const string& en : _edges) _eCut.push_back(std::stod(en) * GeV);
  }
  // Loop over bottoms
  for (const Particle& bottom : apply<UnstableParticles>(event, "UFS").particles()) {
    // remove mixing entries etc
    if (bottom.children()[0].abspid() == bottom.abspid()) continue;
    _nBottom->fill();
    FourMomentum pgamma(0., 0., 0., 0.);
    unsigned int ngamma = 0;
    for (const Particle& child : bottom.children()) {
      if (child.pid() == PID::PHOTON) {
        ngamma += 1;
        pgamma += child.momentum();
      }
    }
    if (ngamma != 1) continue;
    unsigned int nK0(0), nKp(0), nKm(0);
    FourMomentum p_tot(0, 0, 0, 0);
    findDecayProducts(bottom, nK0, nKp, nKm);
    unsigned int nk = nKp - nKm + nK0;
    if (nk % 2 == 1) {
      const LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(
          bottom.momentum().betaVec());
      double eGamma = boost.transform(pgamma).E();
      for (unsigned int ix = 0; ix < _eCut.size(); ++ix) {
        if (eGamma > _eCut[ix]) {
          _h_br->fill(_edges[ix]);
          _p_E->fill(_edges[ix], eGamma);
          _p_E2->fill(_edges[ix], sqr(eGamma));
        }
      }
    }
  }
}


/// Normalise histograms etc., after the run
void finalize() {
  // 1e4 for br ormalization
  scale(_h_br, 1e4 / _nBottom->sumW());
  // dispersion
  BinnedEstimatePtr<string> dispersion;
  book(dispersion, 1, 1, 3);
  for (auto& b : dispersion->bins()) {
    const auto& bE = _p_E->bin(b.index());
    const auto& bE2 = _p_E2->bin(b.index());
    const double val = bE2.mean(2) - sqr(bE.mean(2));
    const double err = val * sqrt(sqr(bE2.stdErr(2) / bE2.mean(2)) + 4. * sqr(bE.stdErr(2) / bE.mean(2)));
    b.set(val, err);
  }
}

/// @}


/// @name Histograms
/// @{
BinnedHistoPtr<string> _h_br;
BinnedProfilePtr<string> _p_E, _p_E2;
CounterPtr _nBottom;
vector<string> _edges;
vector<double> _eCut;
/// @}

};

RIVET_DECLARE_PLUGIN(BELLE_2009_I825222);

} ```