Rivet analyses

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

namespace Rivet {


  /// @brief MC Semi-leptonic decay analysis
  class MC_DECAY_SEMILEPTONIC : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(MC_DECAY_SEMILEPTONIC);


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

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

      // Initialise and register projections
      declare(UnstableParticles(), "UFS");
      // B decays
      bookHistos(511, -413, -11, 5.3);
      bookHistos(-511, 413, 11, 5.3);
      bookHistos(511, -413, -13, 5.3);
      bookHistos(-511, 413, 13, 5.3);
      bookHistos(521, -423, -11, 5.3);
      bookHistos(-521, 423, 11, 5.3);
      bookHistos(521, -423, -13, 5.3);
      bookHistos(-521, 423, 13, 5.3);
      bookHistos(511, -411, -11, 5.3);
      bookHistos(-511, 411, 11, 5.3);
      bookHistos(511, -411, -13, 5.3);
      bookHistos(-511, 411, 13, 5.3);
      bookHistos(521, -421, -11, 5.3);
      bookHistos(-521, 421, 11, 5.3);
      bookHistos(521, -421, -13, 5.3);
      bookHistos(-521, 421, 13, 5.3);
      bookHistos(511, -415, -11, 5.3);
      bookHistos(-511, 415, 11, 5.3);
      bookHistos(511, -415, -13, 5.3);
      bookHistos(-511, 415, 13, 5.3);
      bookHistos(521, -425, -11, 5.3);
      bookHistos(-521, 425, 11, 5.3);
      bookHistos(521, -425, -13, 5.3);
      bookHistos(-521, 425, 13, 5.3);
      bookHistos(511, -10411, -11, 5.3);
      bookHistos(-511, 10411, 11, 5.3);
      bookHistos(511, -10411, -13, 5.3);
      bookHistos(-511, 10411, 13, 5.3);
      bookHistos(521, -10421, -11, 5.3);
      bookHistos(-521, 10421, 11, 5.3);
      bookHistos(521, -10421, -13, 5.3);
      bookHistos(-521, 10421, 13, 5.3);
      bookHistos(511, -10413, -11, 5.3);
      bookHistos(-511, 10413, 11, 5.3);
      bookHistos(511, -10413, -13, 5.3);
      bookHistos(-511, 10413, 13, 5.3);
      bookHistos(521, -10423, -11, 5.3);
      bookHistos(-521, 10423, 11, 5.3);
      bookHistos(521, -10423, -13, 5.3);
      bookHistos(-521, 10423, 13, 5.3);
      bookHistos(511, -20413, -11, 5.3);
      bookHistos(-511, 20413, 11, 5.3);
      bookHistos(511, -20413, -13, 5.3);
      bookHistos(-511, 20413, 13, 5.3);
      bookHistos(521, -20423, -11, 5.3);
      bookHistos(-521, 20423, 11, 5.3);
      bookHistos(521, -20423, -13, 5.3);
      bookHistos(-521, 20423, 13, 5.3);

      bookHistos(511, -213, -11, 5.3);
      bookHistos(-511, 213, 11, 5.3);
      bookHistos(511, -213, -13, 5.3);
      bookHistos(-511, 213, 13, 5.3);
      bookHistos(521, 113, -11, 5.3);
      bookHistos(-521, 113, 11, 5.3);
      bookHistos(521, 113, -13, 5.3);
      bookHistos(-521, 113, 13, 5.3);
      bookHistos(521, 223, -11, 5.3);
      bookHistos(-521, 223, 11, 5.3);
      bookHistos(521, 223, -13, 5.3);
      bookHistos(-521, 223, 13, 5.3);
      bookHistos(511, -211, -11, 5.3);
      bookHistos(-511, 211, 11, 5.3);
      bookHistos(511, -211, -13, 5.3);
      bookHistos(-511, 211, 13, 5.3);
      bookHistos(521, 111, -11, 5.3);
      bookHistos(-521, 111, 11, 5.3);
      bookHistos(521, 111, -13, 5.3);
      bookHistos(-521, 111, 13, 5.3);
      bookHistos(521, 221, -11, 5.3);
      bookHistos(-521, 221, 11, 5.3);
      bookHistos(521, 221, -13, 5.3);
      bookHistos(-521, 221, 13, 5.3);
      bookHistos(521, 331, -11, 5.3);
      bookHistos(-521, 331, 11, 5.3);
      bookHistos(521, 331, -13, 5.3);
      bookHistos(-521, 331, 13, 5.3);

      // D decays
      bookHistos(411, -311, -11, 1.9);
      bookHistos(-411, 311, 11, 1.9);
      bookHistos(411, -311, -13, 1.9);
      bookHistos(-411, 311, 13, 1.9);
      bookHistos(421, -321, -11, 1.9);
      bookHistos(-421, 321, 11, 1.9);
      bookHistos(421, -321, -13, 1.9);
      bookHistos(-421, 321, 13, 1.9);
      bookHistos(411, -313, -11, 1.9);
      bookHistos(-411, 313, 11, 1.9);
      bookHistos(411, -313, -13, 1.9);
      bookHistos(-411, 313, 13, 1.9);
      bookHistos(421, -323, -11, 1.9);
      bookHistos(-421, 323, 11, 1.9);
      bookHistos(421, -323, -13, 1.9);
      bookHistos(-421, 323, 13, 1.9);

      bookHistos(411, -10313, -11, 1.9);
      bookHistos(-411, 10313, 11, 1.9);
      bookHistos(411, -10313, -13, 1.9);
      bookHistos(-411, 10313, 13, 1.9);
      bookHistos(421, -10323, -11, 1.9);
      bookHistos(-421, 10323, 11, 1.9);
      bookHistos(421, -10323, -13, 1.9);
      bookHistos(-421, 10323, 13, 1.9);

      bookHistos(411, -315, -11, 1.9);
      bookHistos(-411, 315, 11, 1.9);
      bookHistos(411, -315, -13, 1.9);
      bookHistos(-411, 315, 13, 1.9);
      bookHistos(421, -325, -11, 1.9);
      bookHistos(-421, 325, 11, 1.9);
      bookHistos(421, -325, -13, 1.9);
      bookHistos(-421, 325, 13, 1.9);

      bookHistos(411, 111, -11, 1.9);
      bookHistos(-411, 111, 11, 1.9);
      bookHistos(411, 111, -13, 1.9);
      bookHistos(-411, 111, 13, 1.9);
      bookHistos(411, 221, -11, 1.9);
      bookHistos(-411, 221, 11, 1.9);
      bookHistos(411, 221, -13, 1.9);
      bookHistos(-411, 221, 13, 1.9);
      bookHistos(411, 331, -11, 1.9);
      bookHistos(-411, 331, 11, 1.9);
      bookHistos(411, 331, -13, 1.9);
      bookHistos(-411, 331, 13, 1.9);
      bookHistos(421, -211, -11, 1.9);
      bookHistos(-421, 211, 11, 1.9);
      bookHistos(421, -211, -13, 1.9);
      bookHistos(-421, 211, 13, 1.9);

      bookHistos(411, 113, -11, 1.9);
      bookHistos(-411, 113, 11, 1.9);
      bookHistos(411, 113, -13, 1.9);
      bookHistos(-411, 113, 13, 1.9);
      bookHistos(411, 223, -11, 1.9);
      bookHistos(-411, 223, 11, 1.9);
      bookHistos(411, 223, -13, 1.9);
      bookHistos(-411, 223, 13, 1.9);
      bookHistos(421, -213, -11, 1.9);
      bookHistos(-421, 213, 11, 1.9);
      bookHistos(421, -213, -13, 1.9);
      bookHistos(-421, 213, 13, 1.9);

      // D_s decays
      bookHistos(431, 221, -11, 1.9);
      bookHistos(-431, 221, 11, 1.9);
      bookHistos(431, 221, -13, 1.9);
      bookHistos(-431, 221, 13, 1.9);
      bookHistos(431, 331, -11, 1.9);
      bookHistos(-431, 331, 11, 1.9);
      bookHistos(431, 331, -13, 1.9);
      bookHistos(-431, 331, 13, 1.9);
      bookHistos(431, 333, -11, 1.9);
      bookHistos(-431, 333, 11, 1.9);
      bookHistos(431, 333, -13, 1.9);
      bookHistos(-431, 333, 13, 1.9);

      bookHistos(431, 311, -11, 1.9);
      bookHistos(-431, -311, 11, 1.9);
      bookHistos(431, 311, -13, 1.9);
      bookHistos(-431, -311, 13, 1.9);
      bookHistos(431, 313, -11, 1.9);
      bookHistos(-431, -313, 11, 1.9);
      bookHistos(431, 313, -13, 1.9);
      bookHistos(-431, -313, 13, 1.9);
    }

    void bookHistos(int id1, int id2, int ilep, double mass) {
      _incoming.push_back(id1);
      _outgoing.push_back(id2);
      _outgoingL.push_back(ilep);
      std::ostringstream title;
      title << "h_" << abs(id1);
      if (id1 > 0)
        title << "p";
      else
        title << "m";
      title << "_" << abs(id2);
      if (id2 > 0)
        title << "p";
      else
        title << "m";
      title << "_" << abs(ilep);
      if (ilep > 0)
        title << "p";
      else
        title << "m";
      title << "_";
      _energy.push_back(Histo1DPtr());
      book(_energy.back(), title.str() + "energy", 200, 0.0, 0.5 * mass / MeV);
      _scale.push_back(Histo1DPtr());
      book(_scale.back(), title.str() + "scale", 200, 0.0, mass / MeV);
    }

    void findDecayProducts(const Particle& mother,
                           unsigned int& nstable,
                           Particles& lp,
                           Particles& lm,
                           Particles& nu,
                           Particles& nub,
                           Particles& out) {
      for (const Particle& p : mother.children()) {
        int id = p.pid();
        if (id == PID::EMINUS || id == PID::MUON) {
          lm.push_back(p);
          ++nstable;
        }
        else if (id == PID::EPLUS || id == PID::ANTIMUON) {
          lp.push_back(p);
          ++nstable;
        }
        else if (id == PID::NU_E || id == PID::NU_EBAR) {
          nu.push_back(p);
          ++nstable;
        }
        else if (id == PID::NU_MU || id == PID::NU_MUBAR) {
          nub.push_back(p);
          ++nstable;
        }
        else if (PID::isMeson(id)) {
          out.push_back(p);
          ++nstable;
        }
        else if (!p.children().empty()) {
          findDecayProducts(p, nstable, lp, lm, nu, nub, out);
        }
        else
          ++nstable;
      }
    }

    /// Perform the per-event analysis
    void analyze(const Event& event) {
      // loop over unstable particles
      for (const Particle& meson : apply<UnstableParticles>(event, "UFS").particles()) {
        int id = meson.pid();
        // spin 0 mesons
        if (!PID::isMeson(id)) continue;
        if (abs(id) % 10 != 1) continue;
        unsigned int nstable(0);
        Particles lp, lm, nu, nub, out;
        findDecayProducts(meson, nstable, lp, lm, nu, nub, out);
        if (nstable != 3 || out.size() != 1) continue;
        int ilep = 0;
        FourMomentum plep, pmnu = out[0].momentum();
        double me2(0.);
        if (lp.size() == 1 && nu.size() == 1 && out.size() == 1) {
          if (nu[0].pid() != -lp[0].pid() + 1) continue;
          ilep = lp[0].pid();
          plep = nu[0].momentum() + lp[0].momentum();
          pmnu += nu[0].momentum();
          me2 = lp[0].mass2();
        }
        else if (lm.size() == 1 && nub.size() == 1 && out.size() == 1) {
          if (nub[0].pid() != -lm[0].pid() - 1) continue;
          ilep = lm[0].pid();
          plep = nub[0].momentum() + lm[0].momentum();
          pmnu += nub[0].momentum();
          me2 = lm[0].mass2();
        }
        else
          continue;
        // check if histos already exist
        unsigned int iloc = 0;
        bool found(false);
        while (!found && iloc < _incoming.size()) {
          if (_incoming[iloc] == id && _outgoing[iloc] == out[0].pid() && ilep == _outgoingL[iloc])
            found = true;
          else
            ++iloc;
        }
        if (!found) {
          MSG_WARNING("MC_DECAY_SEMILEPTONIC" << id << " " << out[0].pid() << " " << " " << ilep << " "
                                              << meson.mass() << "\n");
          continue;
        }
        // add the results to the histogram
        _scale[iloc]->fill(plep.mass() / MeV);
        double ee = 0.5 / meson.mass() * (meson.mass2() - pmnu.mass2() + me2);
        _energy[iloc]->fill(ee / MeV);
      }
    }

    /// Normalise histograms etc., after the run
    void finalize() {
      for (unsigned int ix = 0; ix < _energy.size(); ++ix) {
        normalize(_energy);
        normalize(_scale);
      }
    }

    /// @}


    /// @name Histograms
    /// @{
    /**
     *  PDG codes of the decaying mesons
     */
    vector<long> _incoming;

    /**
     *  PDG codes of the decay products
     */
    vector<long> _outgoing;

    /**
     *  Identidies of the leptons
     */
    vector<long> _outgoingL;

    /**
     *  Histograms
     */
    /// @{
    /**
     *  The lepton energy
     */
    vector<Histo1DPtr> _energy;

    /**
     *  The \f$q\f$ value
     */
    vector<Histo1DPtr> _scale;
    /// @}
  };


  RIVET_DECLARE_PLUGIN(MC_DECAY_SEMILEPTONIC);

}