Rivet analyses
Measurements of event shapes by DELPHI, above and below $m_Z$
Experiment: DELPHI (LEP)
Inspire ID: 620250
Status: VALIDATED
Authors: - Peter Richardson
References: - Eur.Phys.J. C29 (2003) 285-312
Beams: e+ e-
Beam energies: (22.5, 22.5); (33.0, 33.0); (38.0, 38.0); (91.5, 91.5); (94.5, 94.5); (96.0, 96.0); (98.0, 98.0); (100.0, 100.0); (101.0, 101.0); (102.5, 102.5); (103.5, 103.5)GeV
Run details: - Hadronic Z decay events generated below the Z pole.
Measurement of a wide range of event shapes by DELPHI at energies below the Z pole using radiative events and above $m_Z$ from LEP2. This analyses allows the energy dependence of simulations to be studied. Only the distributions and not the means are implemented.
Source code:DELPHI_2003_I620250.cc
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/Beam.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/Hemispheres.hh"
#include "Rivet/Projections/ParisiTensor.hh"
#include "Rivet/Projections/Sphericity.hh"
#include "Rivet/Projections/Thrust.hh"
namespace Rivet {
/// @brief DELPHI event shapes below the Z pole
class DELPHI_2003_I620250 : public Analysis {
public:
/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(DELPHI_2003_I620250);
/// @name Analysis methods
/// @{
/// Book histograms and initialise projections before the run
void init() {
// Initialise and register projections.
declare(Beam(), "Beams");
const FinalState fs;
declare(fs, "FS");
const Thrust thrust(fs);
declare(thrust, "Thrust");
declare(Sphericity(fs), "Sphericity");
declare(ParisiTensor(fs), "Parisi");
declare(Hemispheres(thrust), "Hemispheres");
// Histogram booking offset numbers.
size_t ih = 0, iy = 1;
for (double eVal : allowedEnergies()) {
const string en = toString(round(eVal));
if (isCompatibleWithSqrtS(eVal)) {
_sqs = en;
skipBin = eVal < 200.;
isDisc = eVal >= 100.;
}
book(_c[en], "_sumW_" + en);
if (eVal < 100.) {
book(_h[en + "thrust"], 1, 1, iy);
book(_h[en + "major"], 2, 1, iy);
book(_h[en + "minor"], 3, 1, iy);
book(_h[en + "sphericity"], 4, 1, iy);
book(_h[en + "planarity"], 5, 1, iy);
book(_h[en + "oblateness"], 6, 1, iy);
book(_h[en + "heavy_jet_mass"], 7, 1, iy);
book(_h[en + "light_jet_mass"], 9, 1, iy);
book(_h[en + "diff_jet_mass"], 10, 1, iy);
book(_h[en + "total_jet_mass"], 11, 1, iy);
book(_h[en + "heavy_jet_mass_E"], 8, 1, iy);
book(_h[en + "total_jet_mass_E"], 12, 1, iy);
book(_h[en + "wide_broading"], 13, 1, iy);
book(_h[en + "narrow_broading"], 14, 1, iy);
book(_h[en + "total_broading"], 15, 1, iy);
book(_h[en + "diff_broading"], 16, 1, iy);
book(_h[en + "CParam"], 17, 1, iy);
if (iy == 3) iy = 0;
}
else {
book(_d[en]["rap"], 30 + ih, 1, iy);
book(_d[en]["xi"], 32 + ih, 1, iy);
book(_d[en]["pTIn"], 34 + ih, 1, iy);
book(_d[en]["pTOut"], 36 + ih, 1, iy);
book(_d[en]["thrust"], 38 + ih, 1, iy);
book(_d[en]["major"], 40 + ih, 1, iy);
book(_d[en]["minor"], 42 + ih, 1, iy);
book(_d[en]["oblateness"], 44 + ih, 1, iy);
book(_d[en]["wide_broading"], 46 + ih, 1, iy);
book(_d[en]["total_broading"], 48 + ih, 1, iy);
book(_d[en]["diff_broading"], 50 + ih, 1, iy);
book(_d[en]["CParam"], 52 + ih, 1, iy);
book(_d[en]["DParam"], 54 + ih, 1, iy);
book(_d[en]["heavy_jet_mass"], 56 + ih, 1, iy);
book(_d[en]["heavy_jet_mass_P"], 58 + ih, 1, iy);
book(_d[en]["heavy_jet_mass_E"], 60 + ih, 1, iy);
book(_d[en]["light_jet_mass"], 62 + ih, 1, iy);
book(_d[en]["diff_jet_mass"], 64 + ih, 1, iy);
book(_d[en]["sphericity"], 66 + ih, 1, iy);
book(_d[en]["planarity"], 68 + ih, 1, iy);
book(_d[en]["aplanarity"], 70 + ih, 1, iy);
}
if (iy == 4) {
++ih;
iy = 0;
}
++iy;
}
raiseBeamErrorIf(_sqs.empty());
_axis["rap"] = YODA::Axis<double>({0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.5});
_axis["xi"] = YODA::Axis<double>(
{0.0, 0.4, 0.8, 1.2, 1.6, 2.0, 2.4, 2.8, 3.2, 3.6, 4.0, 4.4, 4.8, 5.2, 5.6, 6.0, 6.4});
_axis["pTIn"] = YODA::Axis<double>({0.0, 0.1, 0.4, 0.65, 0.9, 1.1, 1.4, 2.0, 3.0, 4.0, 6.0, 8.0, 12.0});
_axis["pTOut"] = YODA::Axis<double>({0.0, 0.2, 0.4, 0.6, 0.85, 1.2, 1.6, 2.0, 3.0});
_axis["thrust"] = YODA::Axis<double>({0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09,
0.10, 0.12, 0.14, 0.16, 0.18, 0.20, 0.24, 0.28, 0.32, 0.36});
_axis["major"] = YODA::Axis<double>({0.0, 0.02, 0.04, 0.05, 0.06, 0.07, 0.08, 0.1, 0.12, 0.14, 0.16,
0.2, 0.24, 0.28, 0.32, 0.36, 0.4, 0.44, 0.48, 0.52, 0.56, 0.6});
_axis["minor"] = YODA::Axis<double>(
{0.0, 0.02, 0.04, 0.05, 0.06, 0.07, 0.08, 0.1, 0.12, 0.14, 0.16, 0.2, 0.24, 0.28, 0.32});
_axis["oblateness"] = YODA::Axis<double>({0.0, 0.02, 0.04, 0.06, 0.08, 0.1, 0.12, 0.14, 0.16, 0.18,
0.20, 0.24, 0.28, 0.32, 0.36, 0.4, 0.44});
_axis["wide_broading"] = YODA::Axis<double>(
{0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.1, 0.12, 0.14, 0.17, 0.20, 0.24, 0.28});
_axis["total_broading"] = YODA::Axis<double>({0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07,
0.08, 0.09, 0.1, 0.11, 0.13, 0.15, 0.17, 0.19,
0.21, 0.24, 0.27, 0.3, 0.33, 0.36});
_axis["diff_broading"] = YODA::Axis<double>({0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09,
0.1, 0.12, 0.14, 0.16, 0.18, 0.2, 0.24, 0.28});
_axis["CParam"] = YODA::Axis<double>({0.0, 0.04, 0.08, 0.12, 0.16, 0.2, 0.24, 0.28,
0.32, 0.36, 0.4, 0.44, 0.48, 0.52, 0.56, 0.6,
0.64, 0.68, 0.72, 0.76, 0.8, 0.84, 0.88});
_axis["DParam"] = YODA::Axis<double>({0.00, 0.02, 0.04, 0.06, 0.08, 0.10, 0.12, 0.14, 0.16, 0.20, 0.24,
0.28, 0.32, 0.36, 0.40, 0.44, 0.48, 0.54});
_axis["heavy_jet_mass"] = YODA::Axis<double>(
{0.00, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08, 0.10, 0.12, 0.14, 0.16, 0.20, 0.24, 0.28, 0.32});
_axis["heavy_jet_mass_P"] = YODA::Axis<double>(
{0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08, 0.1, 0.12, 0.14, 0.16, 0.2, 0.24, 0.28, 0.32});
_axis["heavy_jet_mass_E"] = YODA::Axis<double>(
{0.00, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08, 0.10, 0.12, 0.14, 0.16, 0.20, 0.24, 0.28, 0.32});
_axis["light_jet_mass"] = YODA::Axis<double>({0.0, 0.01, 0.02, 0.03, 0.04, 0.05});
_axis["diff_jet_mass"] = YODA::Axis<double>(
{0.0, 0.01, 0.02, 0.03, 0.04, 0.06, 0.08, 0.12, 0.16, 0.2, 0.25, 0.3});
_axis["sphericity"] = YODA::Axis<double>({0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08, 0.1, 0.12,
0.16, 0.2, 0.25, 0.3, 0.35, 0.4, 0.5, 0.6});
_axis["planarity"] = YODA::Axis<double>(
{0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08, 0.1, 0.12, 0.16, 0.2, 0.25, 0.3, 0.35, 0.4});
_axis["aplanarity"] = YODA::Axis<double>(
{0.0, 0.004, 0.01, 0.016, 0.02, 0.03, 0.04, 0.06, 0.08, 0.1, 0.12, 0.16});
}
/// Perform the per-event analysis
void analyze(const Event& event) {
if (isDisc && _edges.empty()) {
for (const auto& item : _axis) {
_edges[item.first] = _d[_sqs][item.first]->xEdges();
}
}
_c[_sqs]->fill();
// Get beams and average beam momentum
const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
const double meanBeamMom = (beams.first.p3().mod() + beams.second.p3().mod()) / 2.0;
MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
const Thrust& thrust = apply<Thrust>(event, "Thrust");
// sphericity related
const Sphericity& sphericity = apply<Sphericity>(event, "Sphericity");
// hemisphere related
const Hemispheres& hemi = apply<Hemispheres>(event, "Hemispheres");
smartfill("thrust", 1. - thrust.thrust());
smartfill("major", thrust.thrustMajor());
smartfill("minor", thrust.thrustMinor());
smartfill("oblateness", thrust.oblateness());
smartfill("sphericity", sphericity.sphericity());
smartfill("planarity", sphericity.planarity());
if (isDisc) smartfill("aplanarity", sphericity.aplanarity());
smartfill("heavy_jet_mass", hemi.scaledM2high());
smartfill("light_jet_mass", hemi.scaledM2low());
smartfill("diff_jet_mass", hemi.scaledM2diff());
smartfill("wide_broading", hemi.Bmax());
if (!isDisc) smartfill("narrow_broading", hemi.Bmin());
smartfill("total_broading", hemi.Bsum());
smartfill("diff_broading", hemi.Bdiff());
if (!isDisc) smartfill("total_jet_mass", hemi.scaledM2low() + hemi.scaledM2high());
// E and p scheme jet masses
Vector3 axis = thrust.thrustAxis();
FourMomentum p4WithE, p4AgainstE;
FourMomentum p4WithP, p4AgainstP;
double Evis(0);
for (const Particle& p : apply<FinalState>(event, "FS").particles()) {
Vector3 p3 = p.momentum().vector3().unitVec();
const double E = p.momentum().E();
Evis += E;
p3 = E * p3;
const double p3Para = dot(p3, axis);
FourMomentum p4E(E, p3.x(), p3.y(), p3.z());
FourMomentum p4P(p.p3().mod(), p.p3().x(), p.p3().y(), p.p3().z());
if (p3Para > 0) {
p4WithE += p4E;
p4WithP += p4P;
}
else if (p3Para < 0) {
p4AgainstE += p4E;
p4AgainstP += p4P;
}
else {
MSG_WARNING("Particle split between hemispheres");
p4WithE += 0.5 * p4E;
p4AgainstE += 0.5 * p4E;
p4WithP += 0.5 * p4P;
p4AgainstP += 0.5 * p4P;
}
}
// E scheme
const double mass2With_E = p4WithE.mass2() / sqr(Evis);
const double mass2Against_E = p4AgainstE.mass2() / sqr(Evis);
// fill the histograms
smartfill("heavy_jet_mass_E", max(mass2With_E, mass2Against_E));
if (!isDisc) smartfill("total_jet_mass_E", mass2With_E + mass2Against_E);
// pscheme
const double mass2With_P = p4WithP.mass2() / sqr(Evis);
const double mass2Against_P = p4AgainstP.mass2() / sqr(Evis);
// fill the histograms
if (isDisc) smartfill("heavy_jet_mass_P", max(mass2With_P, mass2Against_P));
MSG_DEBUG("Calculating Parisi params");
const ParisiTensor& parisi = apply<ParisiTensor>(event, "Parisi");
smartfill("CParam", parisi.C());
if (isDisc) smartfill("DParam", parisi.D());
// single particle distributions
const FinalState& fs = apply<FinalState>(event, "FS");
if (isDisc) {
for (const Particle& p : fs.particles()) {
if (!PID::isCharged(p.pid())) continue;
// Get momentum and energy of each particle.
const Vector3 mom3 = p.p3();
const double energy = p.E();
// Scaled momenta.
const double mom = mom3.mod();
const double scaledMom = mom / meanBeamMom;
const double logInvScaledMom = -std::log(scaledMom);
smartfill("xi", logInvScaledMom);
// Get momenta components w.r.t. thrust and sphericity.
const double momT = dot(thrust.thrustAxis(), mom3);
const double pTinT = dot(mom3, thrust.thrustMajorAxis());
const double pToutT = dot(mom3, thrust.thrustMinorAxis());
smartfill("pTIn", fabs(pTinT / GeV));
smartfill("pTOut", fabs(pToutT / GeV));
// Calculate rapidities w.r.t. thrust and sphericity.
const double rapidityT = 0.5 * std::log((energy + momT) / (energy - momT));
smartfill("rap", fabs(rapidityT));
MSG_TRACE(fabs(rapidityT) << " " << scaledMom / GeV);
}
}
}
void smartfill(const string& tag, const double value) {
if (isDisc) {
size_t idx = _axis[tag].index(value);
// skip masked bin in wide broadening
if (tag == "wide_broading" && skipBin) {
if (idx == 8)
idx = 0;
else if (idx > 8)
--idx;
}
if (idx && idx <= _edges[tag].size()) {
_d[_sqs][tag]->fill(_edges[tag][idx - 1]);
}
else {
_d[_sqs][tag]->fill(string("OTHER"));
}
}
else {
_h[_sqs + tag]->fill(value);
}
}
/// Normalise histograms etc., after the run
void finalize() {
normalize(_h);
scale(_c, crossSectionPerEvent());
scale(_d, crossSectionPerEvent());
for (double eVal : allowedEnergies()) {
const string en = toString(round(eVal));
for (auto& item : _d[en]) {
if (item.first == "rap" || item.first == "xi" || item.first == "pTIn" || item.first == "pTOut") {
if (!isZero(_c[en]->sumW())) scale(item.second, 1.0 / _c[en]->sumW());
}
else
normalize(item.second);
for (auto& b : item.second->bins()) {
size_t idx = b.index();
// skip masked bin in wide broadening
if (skipBin && item.first == "wide_broading") {
if (idx >= 8) ++idx;
}
b.scaleW(1. / _axis[item.first].width(idx));
}
}
}
}
/// @}
/// @name Histograms
/// @{
map<string, map<string, BinnedHistoPtr<string>>> _d;
map<string, YODA::Axis<double>> _axis;
map<string, vector<string>> _edges;
map<string, Histo1DPtr> _h;
map<string, CounterPtr> _c;
bool isDisc, skipBin;
string _sqs = "";
/// @}
};
RIVET_DECLARE_PLUGIN(DELPHI_2003_I620250);
}