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 BABAR_2008_I769107 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2008_I769107);
 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(), "UFS");
 23
 24      // Book histograms
 25      book(_h_spectrum, 1, 1, 1);
 26      book(_nBottom, "TMP/BottomCounter");
 27    }
 28
 29    void findDecayProducts(const Particle& mother,
 30                           unsigned int& nK0, unsigned int& nKp, unsigned int& nKm) {
 31      for (const Particle & p : mother.children()) {
 32        int id = p.pid();
 33        if ( id == PID::KPLUS ) {
 34          ++nKp;
 35        }
 36        else if (id == PID::KMINUS ) {
 37          ++nKm;
 38        }
 39        else if (id == PID::K0S) {
 40          ++nK0;
 41        }
 42        else if (id == PID::PI0 || id == PID::PIPLUS || id == PID::PIMINUS) {
 43	  continue;
 44        }
 45        else if ( !p.children().empty() ) {
 46          findDecayProducts(p, nK0, nKp, nKm);
 47        }
 48        else
 49	  continue;
 50      }
 51    }
 52
 53    /// Perform the per-event analysis
 54    void analyze(const Event& event) {
 55
 56      // Loop over bottoms
 57      for (const Particle& bottom : apply<UnstableParticles>(event, "UFS").particles(Cuts::abspid==521 ||
 58										     Cuts::abspid==511)) {
 59	// remove mixing entries etc
 60	for (const Particle & child : bottom.children())
 61          if (child.abspid() == 511 || child.pid()==bottom.pid() ) continue;
 62        _nBottom->fill();
 63	FourMomentum pgamma(0.,0.,0.,0.);
 64	unsigned int ngamma = 0;
 65        for (const Particle & child : bottom.children()) {
 66	  if (child.pid() == PID::PHOTON) {
 67            ngamma += 1;
 68            pgamma += child.momentum();
 69          }
 70	}
 71	if (ngamma != 1) continue;
 72
 73	unsigned int nK0(0),nKp(0),nKm(0);
 74	findDecayProducts(bottom, nK0, nKp, nKm);
 75	unsigned int nk = nKp-nKm+nK0;
 76	if (nk % 2 == 1) {
 77	  const LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(bottom.momentum().betaVec());
 78	  double eGamma = boost.transform(pgamma).E();
 79	  _h_spectrum->fill(eGamma);
 80	}
 81      }
 82    }
 83
 84
 85    /// Normalise histograms etc., after the run
 86    void finalize() {
 87      scale(_h_spectrum, 1e4/_nBottom->sumW());
 88    }
 89
 90    /// @}
 91
 92
 93    /// @name Histograms
 94    /// @{
 95    Histo1DPtr _h_spectrum;
 96    CounterPtr _nBottom;
 97    /// @}
 98
 99
100  };
101
102
103  RIVET_DECLARE_PLUGIN(BABAR_2008_I769107);
104
105}