1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief B -> K* J/psi
  9  class BELLE_2005_I680703 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2005_I680703);
 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      // histograms
 25      for (unsigned int iz=0; iz<3; ++iz) {
 26        book(_h[iz],1,1,1+iz);
 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 jPsi,Kstar;
 48        if (B.children()[0].pid()==443 && isKstar(B.children()[1].abspid())) {
 49          jPsi = B.children()[0];
 50          Kstar = B.children()[1];
 51        }
 52        else if(B.children()[1].pid()==443 && isKstar(B.children()[0].abspid())) {
 53          jPsi = B.children()[1];
 54          Kstar = B.children()[0];
 55        }
 56        else {
 57          continue;
 58        }
 59        if (jPsi.children().size()!=2) continue;
 60        // find Kstar decay products
 61        Particle K;
 62        if (isK (Kstar.children()[0].abspid()) &&
 63            isPi(Kstar.children()[1].abspid())) {
 64          K = Kstar.children()[0];
 65        }
 66        else if (isK(Kstar.children()[1].abspid()) &&
 67                 isPi(Kstar.children()[0].abspid())) {
 68          K = Kstar.children()[1];
 69        }
 70        else {
 71          continue;
 72        }
 73        // find jPsi decay products
 74        Particle oDec;
 75        if (jPsi.children()[0].pid()==-13 && jPsi.children()[1].pid()== 13) {
 76          oDec = jPsi.children()[1];
 77        }
 78        else if (jPsi.children()[1].pid()==-13 && jPsi.children()[0].pid()== 13) {
 79          oDec = jPsi.children()[0];
 80        }
 81        else if (jPsi.children()[0].pid()==-11 && jPsi.children()[1].pid()== 11) {
 82          oDec = jPsi.children()[1];
 83        }
 84        else if (jPsi.children()[1].pid()==-11 && jPsi.children()[0].pid()== 11) {
 85          oDec = jPsi.children()[0];
 86        }
 87        else {
 88          continue;
 89        }
 90        // boost to B rest frame
 91        LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(B.mom().betaVec());
 92        FourMomentum pJPsi = boost.transform(jPsi.mom());
 93        FourMomentum pKstar = boost.transform(Kstar.mom());
 94        FourMomentum pK     = boost.transform(K    .mom());
 95        FourMomentum pOdec  = boost.transform(oDec .mom());
 96        // axes
 97        Vector3 axisX = pJPsi.p3().unit();
 98        Vector3 axisY = (pK.p3()-axisX.dot(pK.p3())*axisX).unit();
 99        Vector3 axisZ = axisX.cross(axisY).unit();
100        // kaon helicity angle
101        LorentzTransform boostK = LorentzTransform::mkFrameTransformFromBeta(pKstar.betaVec());
102        double cosK = -axisX.dot(boostK.transform(pK).p3().unit());
103        // transversality angles
104        LorentzTransform boostL = LorentzTransform::mkFrameTransformFromBeta(pJPsi.betaVec());
105        Vector3 axisL = boostL.transform(pOdec).p3().unit();
106        double cosL = axisL.dot(axisZ);
107        double phiL = atan2(axisL.dot(axisY),axisL.dot(axisX));
108        // fill hists
109        _h[0]->fill(cosL);
110        _h[1]->fill(phiL);
111        _h[2]->fill(cosK);
112      }
113    }
114
115
116    /// Normalise histograms etc., after the run
117    void finalize() {
118      normalize(_h, 1.0, false);
119    }
120
121    /// @}
122
123
124    /// @name Histograms
125    /// @{
126    Histo1DPtr _h[3];
127    /// @}
128
129
130  };
131
132
133  RIVET_DECLARE_PLUGIN(BELLE_2005_I680703);
134
135}