Rivet analyses

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

namespace Rivet {


  /// @brief Analysis of kinematic distributions of $f_1$ meson decays
  class MC_DECAY_F1 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(MC_DECAY_F1);


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

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

      // Initialise and register projections
      declare(UnstableParticles(), "UFS");

      // eta pi0 pi0 mode
      book(_h_eta0_etapi0, "eta0_etapi0", 70, 0.66, 1.36);
      book(_h_eta0_pi0pi0, "eta0_pi0pi0", 80, 0.2, 1.0);
      book(_h_eta0_etapi0pi0, "eta0_etapi0pi0", 70, 1.0, 1.7);
      // eta pi+pi- mode
      book(_h_eta1_etapip, "eta1_etapip", 70, 0.66, 1.36);
      book(_h_eta1_etapim, "eta1_etapim", 70, 0.66, 1.36);
      book(_h_eta1_pippim, "eta1_pippim", 80, 0.2, 1.0);
      book(_h_eta1_etapippim, "eta1_etapippim", 70, 1.0, 1.7);
      // pi+pi-2pi0
      book(_h_4pi0_pi0pi0, "4pi0_pi0pi0", 80, 0.2, 1.0);
      book(_h_4pi0_pippi0, "4pi0_pippi0", 80, 0.2, 1.0);
      book(_h_4pi0_pimpi0, "4pi0_pimpi0", 80, 0.2, 1.0);
      book(_h_4pi0_pippim, "4pi0_pippim", 80, 0.2, 1.0);
      book(_h_4pi0_pippimpi0, "4pi0_pippimpi0", 100, 0.4, 1.4);
      book(_h_4pi0_pippi0pi0, "4pi0_pippi0pi0", 100, 0.4, 1.4);
      book(_h_4pi0_pimpi0pi0, "4pi0_pimpi0pi0", 100, 0.4, 1.4);
      book(_h_4pi0_4pi, "4pi0_4pi", 70, 1.0, 1.7);
      // 2pi+ 2pi- mode
      book(_h_4pi1_pippip, "4pi1_pippip", 80, 0.2, 1.0);
      book(_h_4pi1_pimpim, "4pi1_pimpim", 80, 0.2, 1.0);
      book(_h_4pi1_pippim, "4pi1_pippim", 80, 0.2, 1.0);
      book(_h_4pi1_pimpimpip, "4pi1_pimpimpip", 100, 0.4, 1.4);
      book(_h_4pi1_pippippim, "4pi1_pippippim", 100, 0.4, 1.4);
      book(_h_4pi1_4pi, "4pi1_4pi", 70, 1.0, 1.7);
    }

    void findDecayProducts(const Particle& mother,
                           unsigned int& nstable,
                           Particles& pip,
                           Particles& pim,
                           Particles& pi0,
                           Particles& eta) {
      for (const Particle& p : mother.children()) {
        int id = p.pid();
        if (id == PID::ETA) {
          eta.push_back(p);
          ++nstable;
        }
        else if (id == PID::PIPLUS) {
          pip.push_back(p);
          ++nstable;
        }
        else if (id == PID::PIMINUS) {
          pim.push_back(p);
          ++nstable;
        }
        else if (id == PID::PI0) {
          pi0.push_back(p);
          ++nstable;
        }
        else if (id == PID::K0S || id == PID::K0L || id == PID::KPLUS || id == PID::KMINUS)
          ++nstable;
        else if (!p.children().empty()) {
          findDecayProducts(p, nstable, pip, pim, pi0, eta);
        }
        else
          ++nstable;
      }
    }

    /// Perform the per-event analysis
    void analyze(const Event& event) {
      // Loop over f_1 mesons
      for (const Particle& f1 : apply<UnstableParticles>(event, "UFS").particles(Cuts::abspid == 20223)) {
        unsigned int nstable(0);
        Particles pip, pim, pi0, eta;
        findDecayProducts(f1, nstable, pip, pim, pi0, eta);
        // pi+ pi- pi0 pi0
        if (nstable == 4 && pip.size() == 1 && pim.size() == 1 && pi0.size() == 2) {
          _h_4pi0_pi0pi0->fill((pi0[0].momentum() + pi0[1].momentum()).mass(), 1.);
          _h_4pi0_pippi0->fill((pip[0].momentum() + pi0[0].momentum()).mass(), 1.);
          _h_4pi0_pippi0->fill((pip[0].momentum() + pi0[1].momentum()).mass(), 1.);
          _h_4pi0_pimpi0->fill((pim[0].momentum() + pi0[0].momentum()).mass(), 1.);
          _h_4pi0_pimpi0->fill((pim[0].momentum() + pi0[1].momentum()).mass(), 1.);
          _h_4pi0_pippim->fill((pip[0].momentum() + pim[0].momentum()).mass(), 1.);
          _h_4pi0_pippimpi0->fill((pip[0].momentum() + pim[0].momentum() + pi0[0].momentum()).mass(), 1.);
          _h_4pi0_pippimpi0->fill((pip[0].momentum() + pim[0].momentum() + pi0[1].momentum()).mass(), 1.);
          _h_4pi0_pippi0pi0->fill((pi0[0].momentum() + pi0[1].momentum() + pip[0].momentum()).mass(), 1.);
          _h_4pi0_pimpi0pi0->fill((pi0[0].momentum() + pi0[1].momentum() + pim[0].momentum()).mass(), 1.);
          _h_4pi0_4pi->fill(
              (pi0[0].momentum() + pi0[1].momentum() + pim[0].momentum() + pip[0].momentum()).mass(), 1.);
        }
        else if (nstable == 4 && pip.size() == 2 && pim.size() == 2) {
          _h_4pi1_pippip->fill((pip[0].momentum() + pip[1].momentum()).mass(), 1.);
          _h_4pi1_pimpim->fill((pim[0].momentum() + pim[1].momentum()).mass(), 1.);
          _h_4pi1_pippim->fill((pip[0].momentum() + pim[0].momentum()).mass(), 1.);
          _h_4pi1_pippim->fill((pip[0].momentum() + pim[1].momentum()).mass(), 1.);
          _h_4pi1_pippim->fill((pip[1].momentum() + pim[0].momentum()).mass(), 1.);
          _h_4pi1_pippim->fill((pip[1].momentum() + pim[1].momentum()).mass(), 1.);
          _h_4pi1_pimpimpip->fill((pim[0].momentum() + pim[1].momentum() + pip[0].momentum()).mass(), 1.);
          _h_4pi1_pimpimpip->fill((pim[0].momentum() + pim[1].momentum() + pip[1].momentum()).mass(), 1.);
          _h_4pi1_pippippim->fill((pip[0].momentum() + pip[1].momentum() + pim[0].momentum()).mass(), 1.);
          _h_4pi1_pippippim->fill((pip[0].momentum() + pip[1].momentum() + pim[1].momentum()).mass(), 1.);
          _h_4pi1_4pi->fill(
              (pip[0].momentum() + pip[1].momentum() + pim[0].momentum() + pim[1].momentum()).mass(), 1.);
        }
        else if (nstable == 3 && eta.size() == 1 && pip.size() == 1 && pim.size() == 1) {
          _h_eta1_etapip->fill((eta[0].momentum() + pip[0].momentum()).mass(), 1.);
          _h_eta1_etapim->fill((eta[0].momentum() + pim[0].momentum()).mass(), 1.);
          _h_eta1_pippim->fill((pim[0].momentum() + pip[0].momentum()).mass(), 1.);
          _h_eta1_etapippim->fill((eta[0].momentum() + pim[0].momentum() + pip[0].momentum()).mass(), 1.);
        }
        else if (nstable == 3 && eta.size() == 1 && pi0.size() == 2) {
          _h_eta0_etapi0->fill((eta[0].momentum() + pi0[0].momentum()).mass(), 1.);
          _h_eta0_etapi0->fill((eta[0].momentum() + pi0[1].momentum()).mass(), 1.);
          _h_eta0_pi0pi0->fill((pi0[0].momentum() + pi0[1].momentum()).mass(), 1.);
          _h_eta0_etapi0pi0->fill((eta[0].momentum() + pi0[0].momentum() + pi0[1].momentum()).mass(), 1.);
        }
      }
    }


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

      normalize(_h_eta0_etapi0);
      normalize(_h_eta0_pi0pi0);
      normalize(_h_eta0_etapi0pi0);

      normalize(_h_eta1_etapip);
      normalize(_h_eta1_etapim);
      normalize(_h_eta1_pippim);
      normalize(_h_eta1_etapippim);

      normalize(_h_4pi0_pi0pi0);
      normalize(_h_4pi0_pippi0);
      normalize(_h_4pi0_pimpi0);
      normalize(_h_4pi0_pippim);
      normalize(_h_4pi0_pippimpi0);
      normalize(_h_4pi0_pippi0pi0);
      normalize(_h_4pi0_pimpi0pi0);
      normalize(_h_4pi0_4pi);

      normalize(_h_4pi1_pippip);
      normalize(_h_4pi1_pimpim);
      normalize(_h_4pi1_pippim);
      normalize(_h_4pi1_pimpimpip);
      normalize(_h_4pi1_pippippim);
      normalize(_h_4pi1_4pi);
    }

    /// @}


    // @name Histograms
    /// @{
    Histo1DPtr _h_eta0_etapi0, _h_eta0_pi0pi0, _h_eta0_etapi0pi0;
    Histo1DPtr _h_eta1_etapip, _h_eta1_etapim, _h_eta1_pippim, _h_eta1_etapippim;
    Histo1DPtr _h_4pi0_pi0pi0, _h_4pi0_pippi0, _h_4pi0_pimpi0, _h_4pi0_pippim, _h_4pi0_pippimpi0,
        _h_4pi0_pippi0pi0, _h_4pi0_pimpi0pi0, _h_4pi0_4pi;
    Histo1DPtr _h_4pi1_pippip, _h_4pi1_pimpim, _h_4pi1_pippim, _h_4pi1_pimpimpip, _h_4pi1_pippippim,
        _h_4pi1_4pi;
    /// @}
  };


  RIVET_DECLARE_PLUGIN(MC_DECAY_F1);

}