1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4#include "Rivet/Projections/DecayedParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief B+ -> p pbar K+, pi+
 10  class BELLE_2008_I754259 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2008_I754259);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // projections
 23      UnstableParticles ufs = UnstableParticles(Cuts::abspid==521);
 24      declare(ufs, "UFS");
 25      DecayedParticles BP(ufs);
 26      declare(BP, "BP");
 27      // histograms
 28      for (unsigned int ix=0; ix<2; ++ix) {
 29        book(_h_mass [ix],1+ix,1,1);
 30        book(_h_angle[ix],3,1,1+ix);
 31      }
 32      book(_h_angle2, {2.*.9382722, 2.0, 2.2, 2.4, 2.6, 2.85});
 33      for (unsigned int ix=0; ix<5; ++ix) {
 34        book(_h_angle2->bin(ix+1), 4, 1, 1+ix);
 35      }
 36      book(_c,"TMP/nB");
 37    }
 38
 39
 40    /// Perform the per-event analysis
 41    void analyze(const Event& event) {
 42      DecayedParticles BP = apply<DecayedParticles>(event, "BP");
 43      // loop over particles
 44      for (unsigned int ix=0; ix<BP.decaying().size(); ++ix) {
 45        int imode=0,sign=1;
 46        _c->fill();
 47        if (BP.modeMatches(ix,3,mode1)) {
 48          imode= 0;
 49          sign = 1;
 50        }
 51        else if (BP.modeMatches(ix,3,mode1CC)) {
 52          imode= 0;
 53          sign =-1;
 54        }
 55        else if (BP.modeMatches(ix,3,mode2)) {
 56          imode= 1;
 57          sign = 1;
 58        }
 59        else if (BP.modeMatches(ix,3,mode2CC)) {
 60          imode= 1;
 61          sign =-1;
 62        }
 63        else {
 64          continue;
 65        }
 66        const Particle& pp   = BP.decayProducts()[ix].at( sign*2212)[0];
 67       	const Particle& pbar = BP.decayProducts()[ix].at(-sign*2212)[0];
 68        const Particle & meson = BP.decayProducts()[ix].at( sign*(imode==0 ? 321 : 211))[0];
 69        FourMomentum pbaryon = pp.mom()+pbar.mom();
 70        double mass = pbaryon.mass();
 71        _h_mass[imode]->fill(mass);
 72        if (mass>2.85) continue;
 73        // boost to B rest frame
 74        LorentzTransform boost =
 75          LorentzTransform::mkFrameTransformFromBeta(BP.decaying()[ix]. mom().betaVec());
 76        FourMomentum pp1 = boost.transform(pbar.mom());
 77        pbaryon = boost.transform(pbaryon);
 78        Vector3 axis = boost.transform(meson.mom()).p3().unit();
 79        LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(pbaryon.betaVec());
 80        double cosp = axis.dot(boost1.transform(pp1).p3().unit());
 81        _h_angle[imode]->fill(cosp);
 82        if (imode==0) _h_angle2->fill(mass,cosp);
 83      }
 84    }
 85
 86
 87    /// Normalise histograms etc., after the run
 88    void finalize() {
 89      scale(_h_mass,  1e6/ *_c);
 90      scale(_h_angle, 1e6/ *_c);
 91      normalize(_h_angle2);
 92    }
 93
 94    /// @}
 95
 96
 97    /// @name Histograms
 98    /// @{
 99    Histo1DPtr _h_mass[2],_h_angle[2];
100    Histo1DGroupPtr _h_angle2;
101    CounterPtr _c;
102    const map<PdgId,unsigned int> mode1   = { { 2212,1}, {-2212,1}, { 321,1} };
103    const map<PdgId,unsigned int> mode1CC = { { 2212,1}, {-2212,1}, {-321,1} };
104    const map<PdgId,unsigned int> mode2   = { { 2212,1}, {-2212,1}, { 211,1} };
105    const map<PdgId,unsigned int> mode2CC = { { 2212,1}, {-2212,1}, {-211,1} };
106    /// @}
107
108
109  };
110
111
112  RIVET_DECLARE_PLUGIN(BELLE_2008_I754259);
113
114}