Rivet analyses

// -*- 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);

}