Rivet analyses

title: MC_DECAY_ONIUM_PIPI

Analysis of kinematics in $\mathcal{O}\to\mathcal{O}^\prime\pi\pi$ decays

Experiment: ()

Status: VALIDATED

Authors: - Peter Richardson

References: none listed

Beams: * *

Beam energies: ANY

Run details: - Any type of process producing bottom or charmonium mesons

Analysis of the kinematics in the decay of bottom and charmonium resonances to lighter resonances and $\pi\pi$

Source code:MC_DECAY_ONIUM_PIPI.cc

```c++ // -- C++ --

include "Rivet/Analysis.hh"

include "Rivet/Projections/UnstableParticles.hh"

namespace Rivet {

class MC_DECAY_ONIUM_PIPI : public Analysis { public:

/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(MC_DECAY_ONIUM_PIPI);


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

/// Book histograms and initialise projections before the run
void init() {

  // Initialise and register projections
  declare(UnstableParticles(), "UFS");
  // psi 2S
  bookHistos(100443, 443, 0.6);
  // psi(3770)
  bookHistos(30443, 443, 0.7);
  // Upsilon (4S)
  bookHistos(300553, 553, 1.2);
  bookHistos(300553, 100553, 0.6);
  // Upsilon (3S)
  bookHistos(200553, 553, 0.9);
  bookHistos(200553, 100553, 0.4);
  // Upsilon (2S)
  bookHistos(100553, 553, 0.6);
}

void bookHistos(int id1, int id2, double deltaM) {
  double twompi = 0.378;
  _incoming.push_back(id1);
  _outgoing.push_back(id2);
  std::ostringstream title;
  title << "h_" << id1 << "_" << id2 << "_";
  _mpipi.push_back(make_pair(Histo1DPtr(), Histo1DPtr()));
  book(_mpipi.back().first, title.str() + "mpippim", 100, twompi / GeV, deltaM / GeV);
  book(_mpipi.back().second, title.str() + "mpi0pi0", 100, twompi / GeV, deltaM / GeV);
  _hel.push_back(make_pair(Histo1DPtr(), Histo1DPtr()));
  book(_hel.back().first, title.str() + "hpippim", 100, -1., 1.);
  book(_hel.back().second, title.str() + "hpi0pi0", 100, 0., 1.);
}

void findDecayProducts(const Particle& mother,
                       unsigned int& nstable,
                       Particles& pip,
                       Particles& pim,
                       Particles& pi0,
                       Particles& onium) {
  for (const Particle& p : mother.children()) {
    int id = p.pid();
    if (id == PID::PIMINUS) {
      pim.push_back(p);
      ++nstable;
    }
    else if (id == PID::PIPLUS) {
      pip.push_back(p);
      ++nstable;
    }
    else if (id == PID::PI0) {
      pi0.push_back(p);
      ++nstable;
    }
    else if (abs(id) % 1000 == 443 || abs(id) % 1000 == 553) {
      onium.push_back(p);
      ++nstable;
    }
    else if (!p.children().empty()) {
      findDecayProducts(p, nstable, pip, pim, pi0, onium);
    }
    else
      ++nstable;
  }
}

/// Perform the per-event analysis
void analyze(const Event& event) {
  // loop over unstable particles
  for (const Particle& vMeson : apply<UnstableParticles>(event, "UFS").particles()) {
    int id = vMeson.pid();
    if (id % 1000 != 443 && id % 1000 != 553) continue;
    unsigned int nstable(0);
    Particles pip, pim, pi0, onium;
    findDecayProducts(vMeson, nstable, pip, pim, pi0, onium);
    // check for onium
    if (onium.size() != 1 || nstable != 3) continue;
    // check for pipi
    if (!((pip.size() == 1 && pim.size() == 1) || pi0.size() == 2)) continue;
    // check if histos already made
    unsigned int iloc = 0;
    bool found(false);
    while (!found && iloc < _incoming.size()) {
      if (_incoming[iloc] == vMeson.pid() && _outgoing[iloc] == onium[0].pid())
        found = true;
      else
        ++iloc;
    }
    // if histos not made, make them
    if (!found) {
      MSG_WARNING("MC_DECAY_ONIUM_PIPI" << vMeson.pid() << " " << onium[0].pid() << " "
                                        << vMeson.mass() - onium[0].mass() << "\n");
      continue;
    }
    // boost to rest frame of the pion pair
    FourMomentum q = vMeson.momentum() - onium[0].momentum();
    LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(q.betaVec());
    FourMomentum qp = onium[0].momentum();
    FourMomentum ppi = pip.size() == 1 ? pip[0].momentum() : pi0[0].momentum();
    qp = boost.transform(qp);
    ppi = boost.transform(ppi);
    double cX = -ppi.p3().unit().dot(qp.p3().unit());
    if (pi0.size() == 2) {
      _mpipi[iloc].second->fill(q.mass());
      _hel[iloc].second->fill(abs(cX));
    }
    else {
      _mpipi[iloc].first->fill(q.mass());
      _hel[iloc].first->fill(cX);
    }
  }
}


/// Normalise histograms etc., after the run
void finalize() {

  // normalize to unity
  for (unsigned int ix = 0; ix < _mpipi.size(); ++ix) {
    normalize(_mpipi[ix].first);
    normalize(_mpipi[ix].second);
    normalize(_hel[ix].first);
    normalize(_hel[ix].second);
  }
}

/// @}

/**
 *  Incoming onium states
 */
vector<long> _incoming;

/**
 *  Outgoing onium states
 */
vector<long> _outgoing;

/**
 *  Histograms for the \f$\pi^+\pi^-\f$ masses
 */
vector<pair<Histo1DPtr, Histo1DPtr>> _mpipi;

/**
 *  Histmgrams for the helicity angles
 */
vector<pair<Histo1DPtr, Histo1DPtr>> _hel;

};

RIVET_DECLARE_PLUGIN(MC_DECAY_ONIUM_PIPI);

} ```