Rivet analyses

title: MC_DECAY_OMEGAPHIA1_3PION

MC analysis of $\omega$, $\phi$ and $a_1\to3\pi$ decays

Experiment: ()

Status: VALIDATED

Authors: - Peter Richardson

References: none listed

Beams: * *

Beam energies: ANY

Run details: - Any process producing a1, omega or phi mesons

Analysis of the mass distributions and Dalitz plots in $a_1\to3\pi$ and $\omega,\phi\to\pi^+\pi^-\pi^0$ decays. Based on old Herwig++ internal analysis.

Source code:MC_DECAY_OMEGAPHIA1_3PION.cc

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

include "Rivet/Analysis.hh"

include "Rivet/Projections/UnstableParticles.hh"

namespace Rivet {

class MC_DECAY_OMEGAPHIA1_3PION : public Analysis { public:

/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(MC_DECAY_OMEGAPHIA1_3PION);


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

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

  // Initialise and register projections
  declare(UnstableParticles(), "UFS");
  // book histograms a_1
  // Histograms for a_10 -> pi0pi0pi0
  book(_hist0, "hist0", 200, 0.2, 1.5);
  // // dalitz plot
  book(_dalitz0, "dalitz0", 50, 0.2, 1.5, 50, 0.2, 1.5);
  // Histograms for a_1+ -> pi0pi0pi+
  // Mass of the pi0pi0 pair
  book(_hist1A, "hist1A", 200, 0.2, 1.5);
  // Mass of the pi0pi+ pair
  book(_hist1B, "hist1B", 200, 0.2, 1.5);
  // dalitz plot
  book(_dalitz1, "dalitz1", 50, 0.2, 1.5, 50, 0.2, 1.5);
  // Histograms for a_10 -> pi+pi-pi0
  // Mass of the pi+pi- pair
  book(_hist2A, "hist2A", 200, 0.2, 1.5);
  // Mass of the pi+pi0 pair
  book(_hist2B, "hist2B", 200, 0.2, 1.5);
  // Mass of the pi-pi0 pair
  book(_hist2C, "hist2C", 200, 0.2, 1.5);
  // dalitz plot
  book(_dalitz2, "dalitz2", 50, 0.2, 1.5, 50, 0.2, 1.5);
  //  Histograms for a_1+ -> pi+pi+pi-
  // Mass of the pi+pi+ pair
  book(_hist3A, "hist3A", 200, 0.2, 1.5);
  // Mass of the pi+pi- pair
  book(_hist3B, "hist3B", 200, 0.2, 1.5);
  // dalitz plot
  book(_dalitz3, "dalitz3", 50, 0.2, 1.5, 50, 0.2, 1.5);

  // Book histograms omega/phi
  for (unsigned int ix = 0; ix < 2; ++ix) {
    double mmax = ix == 0 ? 0.8 : 1.0;
    std::ostringstream title1;
    title1 << "xhist_" << ix + 1;
    _h_xhist.push_back(Histo1DPtr());
    book(_h_xhist.back(), title1.str(), 200, -300., 300.);
    std::ostringstream title2;
    title2 << "yhist_" << ix + 1;
    _h_yhist.push_back(Histo1DPtr());
    book(_h_yhist.back(), title2.str(), 200, 0., 400.);
    std::ostringstream title3;
    title3 << "mplus_" << ix + 1;
    _h_mplus.push_back(Histo1DPtr());
    book(_h_mplus.back(), title3.str(), 200, 200., mmax * 1000.);
    std::ostringstream title4;
    title4 << "mminus_" << ix + 1;
    _h_mminus.push_back(Histo1DPtr());
    book(_h_mminus.back(), title4.str(), 200, 200., mmax * 1000.);
    std::ostringstream title5;
    title5 << "m0_" << ix + 1;
    _h_m0.push_back(Histo1DPtr());
    book(_h_m0.back(), title5.str(), 200, 200., mmax * 1000.);
    std::ostringstream title6;
    title6 << "dalitz_" << ix + 1;
    _h_dalitz.push_back(Histo2DPtr());
    book(_h_dalitz.back(), title6.str(), 50, 0.2, mmax, 50, 0.2, mmax);
  }
}


void findDecayProducts(const Particle& mother,
                       unsigned int& nstable,
                       Particles& pip,
                       Particles& pim,
                       Particles& pi0) {
  for (const Particle& p : mother.children()) {
    int id = p.pid();
    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 (!p.children().empty()) {
      findDecayProducts(p, nstable, pip, pim, pi0);
    }
    else
      ++nstable;
  }
}

/// Perform the per-event analysis
void analyze(const Event& event) {
  for (const Particle& meson : apply<UnstableParticles>(event, "UFS")
                                   .particles(Cuts::pid == PID::PHI || Cuts::pid == PID::OMEGA
                                              || Cuts::abspid == 20213 || Cuts::pid == 20113)) {
    unsigned int nstable(0);
    Particles pip, pim, pi0;
    findDecayProducts(meson, nstable, pip, pim, pi0);
    if (nstable != 3) continue;
    if (meson.pid() < 0) {
      swap(pim, pip);
    }
    if (meson.pid() == PID::PHI || meson.pid() == PID::OMEGA) {
      if (pip.size() != 1 || pim.size() != 1 || pi0.size() != 1) continue;
      unsigned int iloc = meson.pid() == PID::OMEGA ? 0 : 1;
      LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(meson.momentum().betaVec());
      FourMomentum pp = boost.transform(pip[0].momentum());
      FourMomentum pm = boost.transform(pim[0].momentum());
      FourMomentum p0 = boost.transform(pi0[0].momentum());
      double mp = (pp + p0).mass(), mm = (pm + pp).mass();
      _h_mplus[iloc]->fill(mp / MeV);
      _h_mminus[iloc]->fill((pm + p0).mass() / MeV);
      _h_m0[iloc]->fill(mm / MeV);
      double x = pp.t() - pm.t();
      double y = p0.t() - p0.mass();
      _h_xhist[iloc]->fill(x / MeV);
      _h_yhist[iloc]->fill(y / MeV);
      _h_dalitz[iloc]->fill(mp, mm);
    }
    else {
      // a_1+ -> pi+pi+pi-
      if (pip.size() == 2 && pim.size() == 1) {
        _hist3A->fill((pip[0].momentum() + pip[1].momentum()).mass());
        _hist3B->fill((pip[0].momentum() + pim[0].momentum()).mass());
        _hist3B->fill((pip[1].momentum() + pim[0].momentum()).mass());
        _dalitz3->fill((pip[0].momentum() + pim[0].momentum()).mass(),
                       (pip[1].momentum() + pim[0].momentum()).mass());
        _dalitz3->fill((pip[1].momentum() + pim[0].momentum()).mass(),
                       (pip[0].momentum() + pim[0].momentum()).mass());
      }
      // a_1+ -> pi0pi0pi+
      else if (pip.size() == 1 && pi0.size() == 2) {
        _hist1A->fill((pi0[0].momentum() + pi0[1].momentum()).mass());
        _hist1B->fill((pip[0].momentum() + pi0[0].momentum()).mass());
        _hist1B->fill((pip[0].momentum() + pi0[1].momentum()).mass());
        _dalitz1->fill((pip[0].momentum() + pi0[0].momentum()).mass(),
                       (pip[0].momentum() + pi0[1].momentum()).mass());
        _dalitz1->fill((pip[0].momentum() + pi0[1].momentum()).mass(),
                       (pip[0].momentum() + pi0[0].momentum()).mass());
      }
      // a_10 -> pi0pi0pi0
      else if (pi0.size() == 3) {
        _hist0->fill((pi0[0].momentum() + pi0[1].momentum()).mass());
        _hist0->fill((pi0[0].momentum() + pi0[2].momentum()).mass());
        _hist0->fill((pi0[1].momentum() + pi0[2].momentum()).mass());
        _dalitz0->fill((pi0[0].momentum() + pi0[1].momentum()).mass(),
                       (pi0[0].momentum() + pi0[2].momentum()).mass());
        _dalitz0->fill((pi0[0].momentum() + pi0[1].momentum()).mass(),
                       (pi0[1].momentum() + pi0[2].momentum()).mass());
        _dalitz0->fill((pi0[0].momentum() + pi0[2].momentum()).mass(),
                       (pi0[1].momentum() + pi0[2].momentum()).mass());
        _dalitz0->fill((pi0[0].momentum() + pi0[2].momentum()).mass(),
                       (pi0[0].momentum() + pi0[1].momentum()).mass());
        _dalitz0->fill((pi0[1].momentum() + pi0[2].momentum()).mass(),
                       (pi0[0].momentum() + pi0[1].momentum()).mass());
        _dalitz0->fill((pi0[1].momentum() + pi0[2].momentum()).mass(),
                       (pi0[0].momentum() + pi0[2].momentum()).mass());
      }
      // a_10 -> pi+pi-pi0
      else if (pi0.size() == 1 && pip.size() == 1 && pim.size() == 1) {
        _hist2A->fill((pim[0].momentum() + pip[0].momentum()).mass());
        _hist2B->fill((pip[0].momentum() + pi0[0].momentum()).mass());
        _hist2C->fill((pim[0].momentum() + pi0[0].momentum()).mass());
        _dalitz2->fill((pim[0].momentum() + pi0[0].momentum()).mass(),
                       (pip[0].momentum() + pi0[0].momentum()).mass());
      }
    }
  }
}


/// Normalise histograms etc., after the run
void finalize() {
  // a_1
  normalize(_hist0);
  normalize(_dalitz0);
  normalize(_hist1A);
  normalize(_hist1B);
  normalize(_dalitz1);
  normalize(_hist2A);
  normalize(_hist2B);
  normalize(_hist2C);
  normalize(_dalitz2);
  normalize(_hist3A);
  normalize(_hist3B);
  normalize(_dalitz3);
  // omega/phi
  for (unsigned int ix = 0; ix < 2; ++ix) {
    normalize(_h_xhist[ix]);
    normalize(_h_yhist[ix]);
    normalize(_h_mplus[ix]);
    normalize(_h_mminus[ix]);
    normalize(_h_m0[ix]);
    normalize(_h_dalitz[ix]);
  }
}

/// @}

/// @name Histograms a_1
/// @{
// Histograms for a_10 -> pi0pi0pi0
Histo1DPtr _hist0;
// dalitz plot
Histo2DPtr _dalitz0;
// Histograms for a_1+ -> pi0pi0pi+
// Mass of the pi0pi0 pair
Histo1DPtr _hist1A;
// Mass of the pi0pi+ pair
Histo1DPtr _hist1B;
// dalitz plot
Histo2DPtr _dalitz1;
// Histograms for a_10 -> pi+pi-pi0
// Mass of the pi+pi- pair
Histo1DPtr _hist2A;
// Mass of the pi+pi0 pair
Histo1DPtr _hist2B;
// Mass of the pi-pi0 pair
Histo1DPtr _hist2C;
// dalitz plot
Histo2DPtr _dalitz2;
//  Histograms for a_1+ -> pi+pi+pi-
// Mass of the pi+pi+ pair
Histo1DPtr _hist3A;
// Mass of the pi+pi- pair
Histo1DPtr _hist3B;
// dalitz plot
Histo2DPtr _dalitz3;
/// @}

/// @name Histograms omega/phi
/// @{
// Histogram for the x-values
vector<Histo1DPtr> _h_xhist;
// Histogram for the y-values
vector<Histo1DPtr> _h_yhist;
//  The mass of the \rho^+
vector<Histo1DPtr> _h_mplus;
//  The mass of the \rho^-
vector<Histo1DPtr> _h_mminus;
// The mass of the \rho^0
vector<Histo1DPtr> _h_m0;
// Dalitz plot
vector<Histo2DPtr> _h_dalitz;
/// @}

};

RIVET_DECLARE_PLUGIN(MC_DECAY_OMEGAPHIA1_3PION);

} ```