Rivet analyses
Identified Particle Spectra and rates in Υ(1S) decays and continuum at 10.49 GeV
Experiment: CLEO (CESR)
Inspire ID: 205668
Status: VALIDATED
Authors: - Peter Richardson
References: - Phys.Rev. D31 (1985) 2161
Beams: e+ e-
Beam energies: (4.7, 4.7); (5.2, 5.2)GeV
Run details: - e+e- > hadrons at Upslion1s and 10.49 GeV
Spectra and rates for π±, K±, π0, K0, Λ, Ξ−, ρ0, K*±, K*0 and ϕ production in Υ(1S) decays and continuum at 10.49 GeV.
Source
code:CLEO_1985_I205668.cc
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/UnstableParticles.hh"
namespace Rivet {
/// @brief Spectra in Upsilon(1S) decay and nearby continuum
class CLEO_1985_I205668 : public Analysis {
public:
/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(CLEO_1985_I205668);
/// @name Analysis methods
/// @{
/// Book histograms and initialise projections before the run
void init() {
// projections
declare(FinalState(), "FS");
declare(UnstableParticles(), "UFS");
// histos
book(_weightSum_cont, "TMP/weightSumcont");
book(_weightSum_Ups1, "TMP/weightSumUps1");
// multiplcities
for (size_t ix = 0; ix < 2; ++ix) {
for (size_t iy = 0; iy < 12; ++iy) {
book(_mult[ix][iy], "/TMP/MULT_" + toString(ix) + "_" + toString(iy));
}
}
// cont spectra
book(_cont["pip"], 1, 1, 1);
book(_cont["Kp"], 2, 1, 1);
book(_cont["p"], 3, 1, 1);
book(_cont["pi0"], 4, 1, 1);
book(_cont["K0"], 5, 1, 1);
book(_cont["lam"], 6, 1, 1);
book(_cont["xi"], 7, 1, 1);
book(_cont["rho"], 8, 1, 1);
book(_cont["Kstarp"], 9, 1, 1);
book(_cont["Kstar0"], 10, 1, 1);
book(_cont["phi"], 11, 1, 1);
// ups spectra
book(_ups1["pip"], 1, 1, 2);
book(_ups1["Kp"], 2, 1, 2);
book(_ups1["p"], 3, 1, 2);
book(_ups1["pi0"], 4, 1, 2);
book(_ups1["K0"], 5, 1, 2);
book(_ups1["lam"], 6, 1, 2);
book(_ups1["xi"], 7, 1, 2);
book(_ups1["rho"], 8, 1, 2);
book(_ups1["Kstarp"], 9, 1, 2);
book(_ups1["Kstar0"], 10, 1, 2);
book(_ups1["phi"], 11, 1, 2);
_axes[0]["pip"] = YODA::Axis<double>(
{0.05, 0.07, 0.09, 0.11, 0.13, 0.15, 0.17, 0.19, 0.48, 0.58, 0.68, 0.78, 0.98});
_axes[0]["Kp"] = YODA::Axis<double>({0.03, 0.09, 0.11, 0.13, 0.15, 0.17, 0.19});
_axes[0]["p"] = YODA::Axis<double>({0.06, 0.14, 0.155, 0.185, 0.215, 0.245, 0.275});
_axes[0]["pi0"] = YODA::Axis<double>({0.1, 0.2, 0.3, 0.4, 0.5});
_axes[0]["K0"] = YODA::Axis<double>(
{0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.9});
_axes[0]["lam"] = YODA::Axis<double>(
{0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.5, 0.65, 0.8, 0.95});
_axes[0]["xi"] = YODA::Axis<double>({0.2, 0.3, 0.4, 0.5, 0.6, 0.7});
_axes[0]["rho"] = YODA::Axis<double>({0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0});
_axes[0]["Kstarp"] = YODA::Axis<double>({0.06, 0.12, 0.24, 0.36, 0.48, 0.6});
_axes[0]["Kstar0"] = YODA::Axis<double>({0.0, 0.06, 0.12, 0.24, 0.36, 0.48});
_axes[0]["phi"] = YODA::Axis<double>({0.195, 0.385, 0.575, 0.945});
_axes[1]["pip"] = YODA::Axis<double>(
{0.05, 0.07, 0.09, 0.11, 0.13, 0.15, 0.17, 0.19, 0.48, 0.58, 0.68, 0.88});
_axes[1]["Kp"] = YODA::Axis<double>({0.02, 0.1, 0.11, 0.13, 0.15, 0.17, 0.19});
_axes[1]["p"] = _axes[0]["p"];
_axes[1]["pi0"] = _axes[0]["pi0"];
_axes[1]["K0"] = YODA::Axis<double>({0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6,
0.65, 0.7, 0.75, 0.8, 0.85, 0.9});
_axes[1]["lam"] = _axes[0]["lam"];
_axes[1]["xi"] = _axes[0]["xi"];
_axes[1]["rho"] = YODA::Axis<double>({0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7});
_axes[1]["Kstarp"] = _axes[0]["Kstarp"];
_axes[1]["Kstar0"] = YODA::Axis<double>({0.06, 0.12, 0.24, 0.36, 0.48});
_axes[1]["phi"] = YODA::Axis<double>({0.28, 0.36, 0.7, 1.0});
}
/// Recursively walk the decay tree to find decay products of @a p
void findDecayProducts(Particle mother, Particles& unstable) {
for (const Particle& p : mother.children()) {
const int id = p.abspid();
if (id == PID::PIPLUS || id == PID::KPLUS || id == PID::PROTON || id == PID::PI0 || id == PID::K0S
|| id == PID::K0L || id == PID::LAMBDA || id == PID::XIMINUS || id == PID::RHO0 || id == 323
|| id == 313 || id == 225 || id == PID::PHI) {
unstable.push_back(p);
}
if (!p.children().empty()) findDecayProducts(p, unstable);
}
}
/// Perform the per-event analysis
void analyze(const Event& event) {
if (_edges[0].empty()) {
for (const auto& item : _cont) {
_edges[0][item.first] = item.second->xEdges();
_edges[1][item.first] = _ups1[item.first]->xEdges();
}
}
// Find the upsilons
// First in unstable final state
const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
Particles upsilons = ufs.particles(Cuts::pid == 553);
// continuum
if (upsilons.empty()) {
_weightSum_cont->fill();
const FinalState& fs = apply<FinalState>(event, "FS");
// FS particles
for (const Particle& p : fs.particles()) {
int id = p.abspid();
double xp = 2. * p.p3().mod() / sqrtS();
if (id == PID::PIPLUS) {
discfill("pip", xp, 0);
_mult[1][0]->fill();
}
else if (id == PID::KPLUS) {
discfill("Kp", xp, 0);
_mult[1][1]->fill();
}
else if (id == PID::PROTON) {
discfill("p", xp, 0);
_mult[1][2]->fill();
}
}
// Unstable particles
for (const Particle& p : ufs.particles()) {
int id = p.abspid();
double xp = 2. * p.p3().mod() / sqrtS();
if (id == PID::PI0) {
discfill("pi0", xp, 0);
_mult[1][3]->fill();
}
else if (id == PID::K0S || id == PID::K0L) {
discfill("K0", xp, 0);
_mult[1][4]->fill();
}
else if (id == PID::LAMBDA) {
discfill("lam", xp, 0);
_mult[1][5]->fill();
}
else if (id == PID::XIMINUS) {
discfill("xi", xp, 0);
_mult[1][6]->fill();
}
else if (id == PID::RHO0) {
discfill("rho", xp, 0);
_mult[1][7]->fill();
}
else if (id == 323) {
discfill("Kstarp", xp, 0);
_mult[1][8]->fill();
}
else if (id == 313) {
discfill("Kstar0", xp, 0);
_mult[1][9]->fill();
}
else if (id == PID::PHI) {
discfill("phi", xp, 0);
_mult[1][10]->fill();
}
else if (id == 225) {
_mult[1][11]->fill();
}
}
}
else {
for (const Particle& ups : upsilons) {
_weightSum_Ups1->fill();
Particles unstable;
LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(ups.momentum().betaVec());
// Find the decay products we want
findDecayProducts(ups, unstable);
for (const Particle& p : unstable) {
int id = p.abspid();
double xp = 2. * boost.transform(p.momentum()).p3().mod() / ups.mass();
if (id == PID::PIPLUS) {
discfill("pip", xp, 1);
_mult[0][0]->fill();
}
else if (id == PID::KPLUS) {
discfill("Kp", xp, 1);
_mult[0][1]->fill();
}
else if (id == PID::PROTON) {
discfill("p", xp, 1);
_mult[0][2]->fill();
}
else if (id == PID::PI0) {
discfill("pi0", xp, 1);
_mult[0][3]->fill();
}
else if (id == PID::K0S || id == PID::K0L) {
discfill("K0", xp, 1);
_mult[0][4]->fill();
}
else if (id == PID::LAMBDA) {
discfill("lam", xp, 1);
_mult[0][5]->fill();
}
else if (id == PID::XIMINUS) {
discfill("xi", xp, 1);
_mult[0][6]->fill();
}
else if (id == PID::RHO0) {
discfill("rho", xp, 1);
_mult[0][7]->fill();
}
else if (id == 323) {
discfill("Kstarp", xp, 1);
_mult[0][8]->fill();
}
else if (id == 313) {
discfill("Kstar0", xp, 1);
_mult[0][9]->fill();
}
else if (id == PID::PHI) {
discfill("phi", xp, 1);
_mult[0][10]->fill();
}
else if (id == 225) {
_mult[0][11]->fill();
}
}
}
}
}
void discfill(const string& name, const double value, const size_t k) {
string edge = "OTHER";
size_t idx = _axes[k][name].index(value);
if (name == "pip") {
if (idx == 8)
idx = 0;
else if (idx > 8)
idx -= 1;
}
if (idx && idx <= _edges[k][name].size()) edge = _edges[k][name][idx - 1];
(k ? _ups1 : _cont)[name]->fill(edge);
}
/// Normalise histograms etc., after the run
void finalize() {
// multiplicities
const vector<CounterPtr> scales = {_weightSum_Ups1, _weightSum_cont};
for (size_t ix = 0; ix < 12; ++ix) {
BinnedEstimatePtr<string> est;
book(est, ix + 12, 1, 1);
for (size_t iy = 0; iy < 2; ++iy) {
if (scales[iy]->val() > 0.) {
unsigned int iz = iy == 0 ? 2 : 1;
scale(_mult[iy][ix], 1. / *scales[iy]);
est->bin(iz).set(_mult[iy][ix]->val(), _mult[iy][ix]->err());
}
}
}
// spectra
if (_weightSum_cont->val() > 0.) {
scale(_cont, 1. / *_weightSum_cont);
for (auto& hist : _cont) {
for (auto& b : hist.second->bins()) {
size_t idx = b.index();
if (hist.first == "pip" && idx >= 8) idx += 1;
b.scaleW(1. / _axes[0][hist.first].width(idx));
}
}
}
if (_weightSum_Ups1->val() > 0.) {
scale(_ups1, 1. / *_weightSum_Ups1);
for (auto& hist : _ups1) {
for (auto& b : hist.second->bins()) {
size_t idx = b.index();
if (hist.first == "pip" && idx >= 8) idx += 1;
b.scaleW(1. / _axes[1][hist.first].width(idx));
}
}
}
}
/// @}
/// @name Histograms
/// @{
map<string, BinnedHistoPtr<string>> _cont, _ups1;
map<string, YODA::Axis<double>> _axes[2];
map<string, vector<string>> _edges[2];
CounterPtr _weightSum_cont, _weightSum_Ups1;
CounterPtr _mult[2][12];
/// @}
};
RIVET_DECLARE_PLUGIN(CLEO_1985_I205668);
}