1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief B -> K* psi(2S)
  9  class CLEO_2001_I535113 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(CLEO_2001_I535113);
 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      UnstableParticles ufs = UnstableParticles(Cuts::abspid==511 || Cuts::abspid==521);
 23      declare(ufs, "UFS");
 24      // histos
 25      for (unsigned int ix=0; ix<3; ++ix) {
 26        book(_h[ix],1,1,1+ix);
 27      }
 28    }
 29
 30    bool isKstar(int pid) const {
 31      return pid==313 || pid==323;
 32    }
 33
 34    bool isK(int pid) const {
 35      return pid==130 || pid==310 || pid==311 || pid==321;
 36    }
 37
 38    bool isPi(int pid) const {
 39      return pid==211 || pid==111;
 40    }
 41
 42    /// Perform the per-event analysis
 43    void analyze(const Event& event) {
 44      UnstableParticles ufs = apply<UnstableParticles>(event, "UFS");
 45      for (const Particle& B : ufs.particles()) {
 46        if (B.children().size()!=2) continue;
 47        Particle onium,Kstar;
 48        if (B.children()[0].pid()==100443 &&
 49           isKstar(B.children()[1].abspid())) {
 50          onium = B.children()[0];
 51          Kstar = B.children()[1];
 52        }
 53        else if (B.children()[1].pid()==100443 &&
 54          isKstar(B.children()[0].abspid())) {
 55          onium = B.children()[1];
 56          Kstar = B.children()[0];
 57        }
 58        else {
 59          continue;
 60        }
 61        // find Kstar decay products
 62        Particle K;
 63        if (isK (Kstar.children()[0].abspid()) &&
 64            isPi(Kstar.children()[1].abspid())) {
 65          K = Kstar.children()[0];
 66        }
 67        else if(isK (Kstar.children()[1].abspid()) &&
 68          isPi(Kstar.children()[0].abspid())) {
 69          K = Kstar.children()[1];
 70        }
 71        else {
 72          continue;
 73        }
 74        unsigned int imode=0;
 75        if (B.abspid()==511 && K.abspid()==321) imode=2;
 76        else if(B.abspid()==521) {
 77          if(K.abspid()==321) imode=1;
 78          else                imode=0;
 79        }
 80        else {
 81          continue;
 82        }
 83        // boost to B rest frame
 84        LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(B.mom().betaVec());
 85        FourMomentum pOnium = boost.transform(onium.mom());
 86        FourMomentum pKstar = boost.transform(Kstar.mom());
 87        FourMomentum pK     = boost.transform(K    .mom());
 88        // axes
 89        Vector3 axisX = pOnium.p3().unit();
 90        // kaon helicity angle
 91        LorentzTransform boostK = LorentzTransform::mkFrameTransformFromBeta(pKstar.betaVec());
 92        double cosK = -axisX.dot(boostK.transform(pK).p3().unit());
 93        _h[imode]->fill(cosK);
 94      }
 95    }
 96
 97
 98    /// Normalise histograms etc., after the run
 99    void finalize() {
100      normalize(_h, 1.0, false);
101    }
102
103    /// @}
104
105
106    /// @name Histograms
107    /// @{
108    Histo1DPtr _h[3];
109    /// @}
110
111
112  };
113
114
115  RIVET_DECLARE_PLUGIN(CLEO_2001_I535113);
116
117}