1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief B -> Xs gamma
  9  class BELLE_2009_I825222 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2009_I825222);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      // Initialise and register projections
 22      declare(UnstableParticles(Cuts::abspid==521 || Cuts::abspid==511), "UFS");
 23      // Book histograms
 24      book(_h_br, 1, 1, 1);
 25      book(_p_E,  1, 1, 2);
 26      book(_p_E2,"TMP/E2",refData(1,1,3));
 27      book(_nBottom, "TMP/BottomCounter");
 28    }
 29
 30    void findDecayProducts(const Particle& mother,
 31                           unsigned int& nK0, unsigned int& nKp,
 32			   unsigned int& nKm) {
 33      for (const Particle & p : mother.children()) {
 34        int id = p.pid();
 35        if ( id == PID::KPLUS )      ++nKp;
 36        else if (id == PID::KMINUS ) ++nKm;
 37        else if (id == PID::K0S)     ++nK0;
 38        else if (id == PID::PI0 || id == PID::PIPLUS || id == PID::PIMINUS) {
 39	  continue;
 40        }
 41        else if ( !p.children().empty() ) {
 42          findDecayProducts(p, nK0, nKp, nKm);
 43        }
 44      }
 45    }
 46
 47    /// Perform the per-event analysis
 48    void analyze(const Event& event) {
 49      // Loop over bottoms
 50      for (const Particle& bottom : apply<UnstableParticles>(event, "UFS").particles()) {
 51	// remove mixing entries etc
 52	if(bottom.children()[0].abspid()==bottom.abspid()) continue;
 53        _nBottom->fill();
 54	FourMomentum pgamma(0.,0.,0.,0.);
 55	unsigned int ngamma = 0;
 56        for (const Particle & child : bottom.children()) {
 57	  if (child.pid() == PID::PHOTON) {
 58            ngamma += 1;
 59            pgamma += child.momentum();
 60          }
 61	}
 62	if (ngamma != 1) continue;
 63        unsigned int nK0(0),nKp(0),nKm(0);
 64        FourMomentum p_tot(0,0,0,0);
 65        findDecayProducts(bottom, nK0, nKp, nKm);
 66        unsigned int nk = nKp-nKm+nK0;
 67        if (nk % 2 == 1) {
 68          const LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(bottom.momentum().betaVec());
 69          double eGamma = boost.transform(pgamma).E();
 70          for (const auto& bin : _h_br->bins()) {
 71            if(eGamma>bin.xMin()) {
 72              _h_br->fill(bin.xMid());
 73              _p_E ->fill(bin.xMid(),eGamma);
 74              _p_E2->fill(bin.xMid(),sqr(eGamma));
 75            }
 76          }
 77        }
 78      }
 79    }
 80
 81
 82    /// Normalise histograms etc., after the run
 83    void finalize() {
 84      // 1e4 for br ormalization and 0.1 for bin width
 85      scale(_h_br, 1e3/_nBottom->sumW());
 86      // dispersion
 87      Estimate1DPtr dispersion;
 88      book(dispersion,1,1,3);
 89      for (auto& b : dispersion->bins()) {
 90        const auto& bE  = _p_E->bin(b.index());
 91        const auto& bE2 = _p_E2->bin(b.index());
 92        const double val = bE2.xMean()-sqr(bE.xMean());
 93        const double err = val*sqrt(sqr(bE2.xStdErr()/bE2.xMean())+4.*sqr(bE.xStdErr()/bE.xMean()));
 94        b.set(val,err);
 95      }
 96    }
 97
 98    /// @}
 99
100
101    /// @name Histograms
102    /// @{
103    Histo1DPtr _h_br;
104    Profile1DPtr _p_E,_p_E2;
105    CounterPtr _nBottom;
106    /// @}
107
108
109  };
110
111
112  RIVET_DECLARE_PLUGIN(BELLE_2009_I825222);
113
114}