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