Rivet analyses
MC analysis of ω, ϕ and a1 → 3π 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 a1 → 3π and ω, ϕ → π+π−π0 decays. Based on old Herwig++ internal analysis.
Source
code:MC_DECAY_OMEGAPHIA1_3PION.cc
// -*- 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);
}