Rivet analyses
title: PLUTO_1981_I165122
Ratio of the cross section for the production of $K^0_S$ to that for $\mu^+\mu^-$ between $3.6$ and $31.6$ GeV
Experiment: PLUTO (DORIS)
Inspire ID: 165122
Status: VALIDATED
Authors: - Peter Richardson
References: - Phys.Lett. B104 (1981) 79-83, 1981
Beams: e- e+
Beam energies: ANY
Run details: - e+ e- to hadrons and e+ e- to mu+ mu- (for normalization)
Ratio of the cross section for the production of $K^0_S$ to that for $\mu^+\mu^-$ between $3.6$ and $31.6$ GeV. The average number of $K^0$ per hadronic event is also provided for some energies, together with the kaon spectrum at 9.4 and 31.6 GeV. N.B. The point at $9.45\to9.456$ GeV is from the $\Upsilon(1S)$ resonance.
Source code:PLUTO_1981_I165122.cc
```c++ // -- C++ --
include "Rivet/Analysis.hh"
include "Rivet/Projections/Beam.hh"
include "Rivet/Projections/FinalState.hh"
include "Rivet/Projections/UnstableParticles.hh"
namespace Rivet {
/// @brief kaon production at low energies class PLUTO_1981_I165122 : public Analysis { public:
/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(PLUTO_1981_I165122);
/// @name Analysis methods
/// @{
/// Book histograms and initialise projections before the run
void init() {
// Initialise and register projections
declare(Beam(), "Beams");
declare(FinalState(), "FS");
declare(UnstableParticles(), "UFS");
// // Book histograms
book(_c_hadrons, "/TMP/sigma_hadrons", refData<YODA::BinnedEstimate<string>>(1, 1, 1));
for (size_t ix = 0; ix < 2; ++ix) {
book(_c_muons[ix], "/TMP/sigma_muons_" + toString(ix),
refData<YODA::BinnedEstimate<string>>(1 + 2 * ix, 1, 1));
book(_c_kaons[ix], "/TMP/sigma_kaons_" + toString(ix),
refData<YODA::BinnedEstimate<string>>(1 + 2 * ix, 1, 1));
for (const string& en : _c_muons[ix].binning().edges<0>()) {
const size_t idx = en.find("-");
if (idx != string::npos) {
const double emin = stod(en.substr(0, idx));
const double emax = stod(en.substr(idx + 1, string::npos));
if (inRange(sqrtS() / GeV, emin, emax)) {
_sqs[ix] = en;
break;
}
}
else {
const double eval = stod(en) * GeV;
if (isCompatibleWithSqrtS(eval)) {
_sqs[ix] = en;
break;
}
}
}
}
raiseBeamErrorIf(_sqs[0].empty() && _sqs[1].empty());
size_t ih = 3;
for (const double eval : {30.0, 9.4, 9.456}) {
const string label = toString(round(eval / MeV));
if (isCompatibleWithSqrtS(eval)) _sqs[2] = label;
book(_h[label], ++ih, 1, 1);
}
book(_c_hadronsY, "TMP/nUps");
}
/// Recursively walk the decay tree to find decay products of @a p
void findDecayProducts(const Particle& mother, Particles& kaons, Particles& stable) const {
for (const Particle& p : mother.children()) {
const int id = p.pid();
if (id == 130 || id == 310) {
kaons += p;
}
if (id == 111 or p.children().empty())
stable += p;
else {
findDecayProducts(p, kaons, stable);
}
}
}
/// Perform the per-event analysis
void analyze(const Event& event) {
// Get beams and average beam momentum
const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
const double meanBeamMom = 0.5 * (beams.first.p3().mod() + beams.second.p3().mod());
MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
// Find the Upsilons among the unstables
const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
Particles upsilons = ufs.particles(Cuts::pid == 553);
// Continuum
if (upsilons.empty()) {
MSG_DEBUG("No Upsilons found => continuum event");
// final state particles
const FinalState& fs = apply<FinalState>(event, "FS");
map<long, int> nCount;
int ntotal(0);
for (const Particle& p : fs.particles()) {
nCount[p.pid()] += 1;
++ntotal;
}
if (nCount[-13] == 1 and nCount[13] == 1 && ntotal == 2 + nCount[22]) {
// mu+mu- + photons
for (size_t ix = 0; ix < 2; ++ix) {
if (!_sqs[ix].empty()) _c_muons[ix]->fill(_sqs[ix]);
}
}
else { // everything else
if (!_sqs[0].empty()) _c_hadrons->fill(_sqs[0]);
}
// unstable particles
for (const Particle& p : ufs.particles(Cuts::pid == 130 or Cuts::pid == 310)) {
if (!_sqs[2].empty()) {
const double xp = p.p3().mod() / meanBeamMom;
_h[_sqs[2]]->fill(xp);
}
for (size_t ix = 0; ix < 2; ++ix) {
if (!_sqs[ix].empty()) _c_kaons[ix]->fill(_sqs[ix]);
}
}
}
else {
MSG_DEBUG("Upsilons found => resonance event");
for (const Particle& ups : upsilons) {
_c_hadronsY->fill();
Particles kaons, stable;
// Find the decay products we want
findDecayProducts(ups, kaons, stable);
// boost to rest frame (if required)
LorentzTransform cms_boost;
if (ups.p3().mod() > 1 * MeV) {
cms_boost = LorentzTransform::mkFrameTransformFromBeta(ups.mom().betaVec());
}
const double mass = ups.mass();
map<long, int> nCount;
int ntotal(0);
for (const Particle& p : stable) {
nCount[p.pid()] += 1;
++ntotal;
}
for (const Particle& kaon : kaons) {
const FourMomentum p2 = cms_boost.transform(kaon.mom());
const double xp = 2. * p2.p3().mod() / mass;
_h["9456"s]->fill(xp);
_c_kaons[0]->fill("9.45 - 9.466"s);
_c_kaons[1]->fill("9.45 - 9.456"s);
}
// mu+mu- + photons
if (nCount[-13] == 1 and nCount[13] == 1 && ntotal == 2 + nCount[22]) {
for (size_t ix = 0; ix < 2; ++ix) {
_c_muons[ix]->fill("9.45 - 9.466"s);
}
}
// everything else
else {
_c_hadrons->fill("9.45 - 9.466"s);
}
}
}
}
/// Normalise histograms etc., after the run
void finalize() {
BinnedEstimatePtr<string> ratio;
book(ratio, 1, 1, 1);
divide(_c_kaons[0], _c_muons[0], ratio);
book(ratio, 2, 1, 1);
divide(_c_kaons[0], _c_hadrons, ratio);
book(ratio, 3, 1, 1);
divide(_c_kaons[1], _c_muons[1], ratio);
// normalize the spectra if required
for (auto& item : _h) {
if (item.first == "9456"s) {
if (_c_hadronsY->val()) scale(item.second, 1.0 / *_c_hadronsY);
}
else {
// beam energy not well defined during merging!
scale(item.second, sqr(sqrtS()) * crossSection() / microbarn / sumOfWeights());
}
}
}
/// @}
/// @name Histograms
/// @{
BinnedHistoPtr<string> _c_hadrons, _c_muons[2], _c_kaons[2];
map<string, Histo1DPtr> _h;
CounterPtr _c_hadronsY;
string _sqs[3];
/// @}
};
RIVET_DECLARE_PLUGIN(PLUTO_1981_I165122);
} ```