Rivet analyses

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

namespace Rivet {


  /// @brief Upsilon(1S) -> gamma pi+pi- K+K-
  class BABAR_2018_I1667191 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2018_I1667191);


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

    /// Book histograms and initialise projections before the run
    void init() {
      // projections
      UnstableParticles ufs = UnstableParticles(Cuts::abspid == 200553 or Cuts::abspid == 100553);
      declare(ufs, "UFS");
      DecayedParticles UPS(ufs);
      UPS.addStable(553);
      declare(UPS, "UPS");
      // histograms
      for (unsigned int ix = 0; ix < 2; ++ix) book(_h_mass[ix], 1 + ix, 1, 1);
      book(_h_pi, 3, 1, 1);
      for (unsigned int ix = 0; ix < 3; ++ix)
        for (unsigned int iy = 0; iy < 2; ++iy) book(_h_angle[ix][iy], 4 + ix, 1, 1 + iy);
    }


    /// Recursively walk the decay tree to find decay products of @a p
    void findDecayProducts(Particle mother,
                           Particles& gamma,
                           Particles& pip,
                           Particles& pim,
                           Particles& Kp,
                           Particles& Km,
                           unsigned int& nstable) {
      for (const Particle& p : mother.children()) {
        if (p.pid() == 211)
          pip.push_back(p);
        else if (p.pid() == -211)
          pim.push_back(p);
        else if (p.pid() == 321)
          Kp.push_back(p);
        else if (p.pid() == -321)
          Km.push_back(p);
        else if (p.pid() == 22)
          gamma.push_back(p);
        else if (p.children().empty())
          nstable += 1;
        else {
          findDecayProducts(p, gamma, pip, pim, Kp, Km, nstable);
        }
      }
    }

    /// Perform the per-event analysis
    void analyze(const Event& event) {
      static const map<PdgId, unsigned int>& mode = {{553, 1}, {211, 1}, {-211, 1}};
      DecayedParticles UPS = apply<DecayedParticles>(event, "UPS");
      // loop over particles
      for (unsigned int ix = 0; ix < UPS.decaying().size(); ++ix) {
        // check pi+pi- upslion(1S) decay mode
        if (!UPS.modeMatches(ix, 3, mode)) continue;
        const Particle& ups1 = UPS.decayProducts()[ix].at(553)[0];
        const Particle& pips = UPS.decayProducts()[ix].at(211)[0];
        const Particle& pims = UPS.decayProducts()[ix].at(-211)[0];
        // boost to rest frame
        LorentzTransform boost;
        if (UPS.decaying()[ix].p3().mod() > 1 * MeV)
          boost = LorentzTransform::mkFrameTransformFromBeta(UPS.decaying()[ix].momentum().betaVec());
        FourMomentum ppipi = boost.transform(pips.momentum() + pims.momentum());
        Vector3 axis1 = ppipi.p3().unit();
        LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(ppipi);
        FourMomentum ppi = boost1.transform(boost.transform(pips.momentum()));
        _h_pi->fill(abs(ppi.p3().unit().dot(axis1)));
        unsigned int nstable = 0;
        Particles gamma, pip, pim, Kp, Km;
        findDecayProducts(ups1, gamma, pip, pim, Kp, Km, nstable);
        if (gamma.size() != 1 || nstable != 0) continue;
        // gamma pi+pi-
        if (Kp.empty() && Km.empty() && pip.size() == 1 && pim.size() == 1) {
          ppipi = pip[0].momentum() + pim[0].momentum();
          double mpipi = ppipi.mass();
          _h_mass[0]->fill(mpipi);
          axis1 *= -1;
          LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(
              boost.transform(ups1.momentum()).betaVec());
          FourMomentum pgamma = boost2.transform(boost.transform(gamma[0].momentum()));
          Vector3 axis2 = pgamma.p3().unit();
          double cGamma = axis2.dot(axis1);
          ppipi = boost2.transform(boost1.transform(ppipi));
          LorentzTransform boost3 = LorentzTransform::mkFrameTransformFromBeta(ppipi.betaVec());
          Vector3 axis3 = boost3.transform(boost2.transform(boost1.transform(pip[0].momentum()))).p3().unit();
          double cH = axis3.dot(axis2);
          int iloc = -1;
          if (mpipi > 0.6 && mpipi < 1.)
            iloc = 0;
          else if (mpipi > 1.092 && mpipi < 1.46)
            iloc = 1;
          if (iloc >= 0) {
            _h_angle[iloc][0]->fill(cGamma);
            _h_angle[iloc][1]->fill(cH);
          }
        }
        // gamma K+K-
        else if (pip.empty() && pim.empty() && Kp.size() == 1 && Km.size() == 1) {
          FourMomentum pKK = Kp[0].momentum() + Km[0].momentum();
          double mKK = pKK.mass();
          _h_mass[1]->fill(mKK);
          axis1 *= -1;
          LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(
              boost.transform(ups1.momentum()).betaVec());
          FourMomentum pgamma = boost2.transform(boost.transform(gamma[0].momentum()));
          Vector3 axis2 = pgamma.p3().unit();
          double cGamma = axis2.dot(axis1);
          pKK = boost2.transform(boost1.transform(pKK));
          LorentzTransform boost3 = LorentzTransform::mkFrameTransformFromBeta(pKK.betaVec());
          Vector3 axis3 = boost3.transform(boost2.transform(boost1.transform(Kp[0].momentum()))).p3().unit();
          double cH = axis3.dot(axis2);
          if (mKK > 1.424 && mKK < 1.62) {
            _h_angle[2][0]->fill(cGamma);
            _h_angle[2][1]->fill(cH);
          }
        }
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      for (unsigned int ix = 0; ix < 2; ++ix) normalize(_h_mass[ix], 1., false);
      normalize(_h_pi, 1., false);
      for (unsigned int ix = 0; ix < 3; ++ix)
        for (unsigned int iy = 0; iy < 2; ++iy) normalize(_h_angle[ix][iy], 1., false);
    }

    /// @}


    /// @name Histograms
    /// @{
    Histo1DPtr _h_mass[2];
    Histo1DPtr _h_pi;
    Histo1DPtr _h_angle[3][2];
    /// @}
  };


  RIVET_DECLARE_PLUGIN(BABAR_2018_I1667191);

}