Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/DecayedParticles.hh"
#include "Rivet/Projections/UnstableParticles.hh"

namespace Rivet {


  /// @brief B -> D* omega pi
  class BELLE_2015_I1369998 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2015_I1369998);


    /// @name Analysis methods
    /// @{

    /// Book histograms and initialise projections before the run
    void init() {
      UnstableParticles ufs = UnstableParticles(Cuts::abspid == 511);
      declare(ufs, "UFS");
      DecayedParticles B0(ufs);
      B0.addStable(413);
      B0.addStable(-413);
      B0.addStable(223);
      declare(B0, "B0");
      for (unsigned int ix = 0; ix < 4; ++ix)
        for (unsigned int iy = 0; iy < 6; ++iy) book(_h[ix][iy], 1 + ix, 1, 1 + iy);
    }

    void findChildren(const Particle& p,
                      Particles& pim,
                      Particles& pip,
                      Particles& pi0,
                      unsigned int& ncount) {
      for (const Particle& child : p.children()) {
        if (child.pid() == PID::PIPLUS) {
          pip.push_back(child);
          ncount += 1;
        }
        else if (child.pid() == PID::PIMINUS) {
          pim.push_back(child);
          ncount += 1;
        }
        else if (child.pid() == PID::PI0) {
          pi0.push_back(child);
          ncount += 1;
        }
        else if (child.children().empty()) {
          ncount += 1;
        }
        else
          findChildren(child, pim, pip, pi0, ncount);
      }
    }

    /// Perform the per-event analysis
    void analyze(const Event& event) {
      static const map<PdgId, unsigned int>& mode = {{413, 1}, {223, 1}, {-211, 1}};
      static const map<PdgId, unsigned int>& modeCC = {{-413, 1}, {223, 1}, {211, 1}};
      DecayedParticles B0 = apply<DecayedParticles>(event, "B0");
      // loop over particles
      for (unsigned int ix = 0; ix < B0.decaying().size(); ++ix) {
        int sign = 1;
        if (B0.decaying()[ix].pid() < 0 && B0.modeMatches(ix, 3, mode))
          sign = 1;
        else if (B0.decaying()[ix].pid() > 0 && B0.modeMatches(ix, 3, modeCC))
          sign = -1;
        else
          continue;
        const Particle& Dstar = B0.decayProducts()[ix].at(sign * 413)[0];
        const Particle& omega = B0.decayProducts()[ix].at(223)[0];
        const Particle& pim1 = B0.decayProducts()[ix].at(-sign * 211)[0];
        // mass hists, no cuts
        double mOmegaPi2 = (omega.momentum() + pim1.momentum()).mass2();
        _h[0][0]->fill(mOmegaPi2);
        double mDstarpi2 = (Dstar.momentum() + pim1.momentum()).mass2();
        _h[1][0]->fill(mDstarpi2);
        // check the no of decay products
        if (Dstar.children().size() != 2 || omega.children().size() != 3) continue;
        // find the children of the D* meson
        Particle D0, pip1;
        if (Dstar.children()[0].pid() == sign * 211 && Dstar.children()[1].pid() == sign * 421) {
          pip1 = Dstar.children()[0];
          D0 = Dstar.children()[1];
        }
        else if (Dstar.children()[1].pid() == sign * 211 && Dstar.children()[0].pid() == sign * 421) {
          pip1 = Dstar.children()[1];
          D0 = Dstar.children()[0];
        }
        else
          continue;
        // children of the omega
        unsigned int ncount = 0;
        Particles pip, pim, pi0;
        findChildren(omega, pim, pip, pi0, ncount);
        if (ncount != 3 || !(pim.size() == 1 && pip.size() == 1 && pi0.size() == 1)) continue;
        // first bottom to the B frame
        LorentzTransform boostB = LorentzTransform::mkFrameTransformFromBeta(
            B0.decaying()[ix].momentum().betaVec());
        FourMomentum pOmega = boostB.transform(omega.momentum());
        FourMomentum pDstar = boostB.transform(Dstar.momentum());
        FourMomentum pD = boostB.transform(D0.momentum());
        FourMomentum ppim1 = boostB.transform(pim1.momentum());
        FourMomentum ppim2 = boostB.transform(pim[0].momentum());
        FourMomentum ppip1 = boostB.transform(pip1.momentum());
        FourMomentum ppip2 = boostB.transform(pip[0].momentum());
        // ---------------------- First set of angles --------------------------------------
        // first the angles for D* (pi omega)
        LorentzTransform boostD = LorentzTransform::mkFrameTransformFromBeta(pDstar.betaVec());
        Vector3 axisD = boostD.transform(pD).p3().unit();
        Vector3 axispip1 = boostD.transform(ppip1).p3().unit();
        Vector3 axisDstar = (pOmega + ppim1).p3().unit();
        double cBeta1 = axisDstar.dot(axisD);
        _h[0][3]->fill(cBeta1);
        LorentzTransform boostWpi = LorentzTransform::mkFrameTransformFromBeta((pOmega + ppim1).betaVec());
        FourMomentum pOmega2 = boostWpi.transform(pOmega);
        Vector3 axisW = pOmega2.p3().unit();
        Vector3 axisWpi = (pOmega + ppim1).p3().unit();
        double cXi1 = axisWpi.dot(axisW);
        _h[0][1]->fill(cXi1);
        // now angle between the two planes
        Vector3 transW = axisW - cXi1 * axisWpi;
        Vector3 transD = axisD - cBeta1 * axisDstar;
        double psi1 = atan2(transW.cross(transD).dot(axisDstar), transW.dot(transD));
        _h[0][5]->fill(psi1);
        // normal to omega decay plane
        LorentzTransform boostW = LorentzTransform::mkFrameTransformFromBeta(pOmega2.betaVec());
        FourMomentum ppim3 = boostW.transform(boostWpi.transform(ppim2));
        FourMomentum ppip3 = boostW.transform(boostWpi.transform(ppip2));
        Vector3 nW = ppim3.p3().cross(ppip3.p3()).unit();
        // boost B decay products to omega rest frame
        FourMomentum pOmegaPi = boostW.transform(boostWpi.transform(pOmega + ppim1));
        FourMomentum pDstar2 = boostW.transform(boostWpi.transform(pDstar));
        Vector3 axisWpi2 = pOmegaPi.p3().unit();
        double cTheta1 = axisWpi2.dot(nW);
        transW = nW - cTheta1 * axisWpi2;
        transD = pDstar2.p3().unit() - pDstar2.p3().unit().dot(axisWpi2) * axisWpi2;
        double phi1 = atan2(transW.cross(transD).dot(axisWpi2), transW.dot(transD));
        _h[0][2]->fill(cTheta1);
        _h[0][4]->fill(phi1);
        // ---------------------- Second set of angles --------------------------------------
        // boost to D* pi frame
        LorentzTransform boostDpi = LorentzTransform::mkFrameTransformFromBeta((pDstar + ppim1).betaVec());
        pDstar2 = boostDpi.transform(pDstar);
        pOmega2 = boostDpi.transform(pOmega);
        axisW = pOmega2.p3().unit();
        axisDstar = pDstar2.p3().unit();
        double cXi2 = axisW.dot(axisDstar);
        _h[1][1]->fill(cXi2);
        // boost to D* rest frame
        LorentzTransform boostDstar = LorentzTransform::mkFrameTransformFromBeta(pDstar2.betaVec());
        axisW = boostDstar.transform(pOmega2).p3().unit();
        Vector3 axisDSpi = boostDstar.transform(boostDpi.transform(pDstar + ppim1)).p3().unit();
        axisD = boostDstar.transform(boostDpi.transform(pD)).p3().unit();
        double cBeta2 = axisD.dot(axisDSpi);
        _h[1][3]->fill(cBeta2);
        transW = axisW - axisW.dot(axisDSpi) * axisDSpi;
        transD = axisD - cBeta2 * axisDSpi;
        double psi2 = atan2(transW.cross(transD).dot(axisDSpi), transW.dot(transD));
        _h[1][5]->fill(psi2);
        // boost to omega frame
        boostW = LorentzTransform::mkFrameTransformFromBeta(pOmega.betaVec());
        ppim3 = boostW.transform(ppim2);
        ppip3 = boostW.transform(ppip2);
        nW = ppim3.p3().cross(ppip3.p3()).unit();
        axisDSpi = boostW.transform(pDstar + ppim1).p3().unit();
        axisDstar = boostW.transform(pDstar).p3().unit();
        double cTheta2 = axisDSpi.dot(nW);
        _h[1][2]->fill(cTheta2);
        transW = nW - cTheta2 * axisDSpi;
        transD = axisDstar.unit() - axisDstar.dot(axisDSpi) * axisDSpi;
        double phi2 = atan2(transW.cross(transD).dot(axisDSpi), transW.dot(transD));
        _h[1][4]->fill(psi2);
        // restricted plots
        if (abs(cTheta1) > .5) {
          _h[2][0]->fill(mOmegaPi2);
        }
        else {
          _h[2][1]->fill(mOmegaPi2);
          _h[2][3]->fill(cBeta1);
          _h[2][5]->fill(psi1);
          _h[3][1]->fill(mDstarpi2);
          _h[3][3]->fill(cTheta2);
          _h[3][5]->fill(phi2);
        }
        if (cXi2 > -.4) {
          _h[2][2]->fill(cBeta1);
          _h[2][4]->fill(psi1);
          _h[3][0]->fill(mDstarpi2);
          _h[3][2]->fill(cTheta2);
          _h[3][4]->fill(phi2);
        }
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      for (unsigned int ix = 0; ix < 4; ++ix)
        for (unsigned int iy = 0; iy < 6; ++iy) normalize(_h[ix][iy], 1., false);
    }

    /// @}


    /// @name Histograms
    /// @{
    Histo1DPtr _h[4][6];
    /// @}
  };


  RIVET_DECLARE_PLUGIN(BELLE_2015_I1369998);

}