Rivet analyses
title: DELPHI_2001_I526164
Charged and Identified particle spectra in $q\bar{q}$ and $WW$ events
Experiment: DELPHI (LEP)
Inspire ID: 526164
Status: VALIDATED
Authors: - Peter Richardson
References: - Eur.Phys.J.C 18 (2000) 203-228, 2000
Beams: e+ e-
Beam energies: (66.5, 66.5); (80.5, 80.5); (86.0, 86.0); (91.5, 91.5); (94.5, 94.5); (100.0, 100.0)GeV
Run details: - e+ e- -> q qbar amd W+W- events
Measurement of charged and identified particle production in both $q\bar{q}$ and $W^+W^-$ events at LEP. Need $q\bar{q}$ at all the energies and $W^+W^-$ at 183 and 189 GeV.
Source code:DELPHI_2001_I526164.cc
```c++ // -- C++ --
include "Rivet/Analysis.hh"
include "Rivet/Projections/FinalState.hh"
include "Rivet/Projections/Thrust.hh"
include "Rivet/Projections/UnstableParticles.hh"
namespace Rivet {
/// @brief DELPHI W decay analysis class DELPHI_2001_I526164 : public Analysis { public:
/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(DELPHI_2001_I526164);
/// @name Analysis methods
///@{
/// Book histograms and initialise projections before the run
void init() {
// projections
declare(UnstableParticles(), "UFS");
declare(FinalState(), "FS");
// Histogram booking
for (size_t ih = 1; ih < 7; ++ih) {
book(_e[ih], (ih == 1) ? 1 : 2, 1, (ih == 1) ? 1 : ih - 1);
if (ih < 3) book(_e[99 + ih], 3, 1, 2 * ih - 1);
}
for (double eVal : allowedEnergies()) {
const string en = toString(round(eVal));
if (isCompatibleWithSqrtS(eVal)) _sqs = en;
for (size_t ih = 0; ih < 7; ++ih) {
book(_c[en + toString(ih)], "_n_qq_"s + toString(ih) + "_"s + en);
}
if (en == "183"s) {
book(_q[en]["K0"], 16, 1, 1);
book(_q[en]["Lam"], 16, 1, 3);
book(_d[en + "p"], 8, 1, 1);
book(_d[en + "xi"], 10, 1, 3);
book(_d[en + "pT"], 12, 1, 3);
for (size_t ih = 0; ih < 2; ++ih) {
const string iw = toString(ih);
book(_h[en + iw]["p_charged"], 7, 1, 2 - ih);
book(_h[en + iw]["p_chargedB"], "/TMP/h_7_1_" + toString(2 - ih),
_d[en + "p"].binning().edges<0>());
book(_h[en + iw]["xi_charged"], 10, 1, 2 - ih);
book(_h[en + iw]["xi_chargedB"], "/TMP/h_10_1_" + toString(2 - ih),
_d[en + "xi"].binning().edges<0>());
book(_h[en + iw]["pT_charged"], 12, 1, 2 - ih);
book(_h[en + iw]["pT_chargedB"], "/TMP/h_12_1_" + toString(2 - ih),
_d[en + "pT"].binning().edges<0>());
}
}
else if (en == "189"s) {
book(_q[en]["K0"], 16, 1, 2);
book(_q[en]["Lam"], 16, 1, 4);
book(_d[en + "p"], 6, 1, 1);
book(_d[en + "xi"], 9, 1, 3);
book(_d[en + "pT"], 11, 1, 3);
for (size_t ih = 0; ih < 2; ++ih) {
const string iw = toString(ih);
book(_h[en + iw]["p_charged"], 5, 1, 2 - ih);
book(_h[en + iw]["p_chargedB"], "/TMP/h_5_1_" + toString(2 - ih),
_d[en + "p"].binning().edges<0>());
book(_h[en + iw]["xi_charged"], 9, 1, 2 - ih);
book(_h[en + iw]["xi_chargedB"], "/TMP/h_9_1_" + toString(2 - ih),
_d[en + "xi"].binning().edges<0>());
book(_h[en + iw]["pT_charged"], 11, 1, 2 - ih);
book(_h[en + iw]["pT_chargedB"], "/TMP/h_11_1_" + toString(2 - ih),
_d[en + "pT"].binning().edges<0>());
for (size_t iy = 0; iy < 4; ++iy) {
book(_h[en + iw]["xi_ident" + toString(iy)], 13 + ih, 1, 1 + iy);
if (!ih) book(_h[en + iw]["xi_ident2" + toString(iy)], 15, 1, 1 + iy);
if (iy) book(_e[10 + 2 * iy - ih], 4, 1, 2 * iy - ih);
}
}
}
for (size_t ix = 0; ix < 2; ++ix) {
for (size_t iy = 0; iy < 5; ++iy) {
const string suff = toString(ix) + toString(iy);
book(_c[en + suff], "_n_WW_"s + suff + "_"s + en);
}
}
}
raiseBeamErrorIf(_sqs.empty());
isWW = _sqs == "183"s || _sqs == "189"s;
}
/// Perform the per-event analysis
void analyze(const Event& event) {
const double mW = 80.379 * GeV;
multimap<Particle, Particles, ParticleOrdering> Wbosons;
// loop over final state particles to find W's
const Particles& finalState = apply<FinalState>(event, "FS").particles();
// loop over FS particles
for (const Particle& p : finalState) {
Particle parent = p;
while (!parent.parents().empty()) {
parent = parent.parents()[0];
if (parent.abspid() == 24 && parent.mass() > 20.) break;
}
if (parent.abspid() != 24) continue;
// find those which came from W
bool found = false;
for (auto& W : Wbosons) {
// W already in list add particle to its decay products
if (fuzzyEquals(W.first.mom(), parent.mom())) {
W.second.push_back(p);
found = true;
break;
}
}
if (!found) {
// check W not child
bool Wchild = false;
for (const Particle& child : parent.children()) {
if (child.abspid() == 24) {
Wchild = true;
break;
}
}
// add to list
if (!Wchild) {
Particles temp = {p};
Wbosons.insert(make_pair(parent, temp));
}
}
}
// no W's => q qbar event
if (Wbosons.empty()) {
_c[_sqs + "0"s]->fill();
UnstableParticles ufs = apply<UnstableParticles>(event, "UFS");
for (const Particle& p : ufs.particles()) {
if (p.abspid() == PID::PIPLUS) {
_c[_sqs + "1"s]->fill();
_c[_sqs + "2"s]->fill();
}
else if (p.abspid() == PID::KPLUS) {
_c[_sqs + "1"s]->fill();
_c[_sqs + "3"s]->fill();
}
else if (p.abspid() == PID::PROTON) {
_c[_sqs + "1"s]->fill();
_c[_sqs + "5"s]->fill();
}
else if (p.abspid() == 130 || p.abspid() == 310) {
_c[_sqs + "4"s]->fill();
if (isWW) _q[_sqs]["K0"]->fill(-log(2. * p.p3().mod() / sqrtS()));
}
else if (p.abspid() == PID::LAMBDA) {
_c[_sqs + "6"s]->fill();
if (isWW) _q[_sqs]["Lam"]->fill(-log(2. * p.p3().mod() / sqrtS()));
}
}
}
else if (Wbosons.size() == 2) {
bool leptonic[2] = {false, false};
unsigned int iboson = 0;
for (auto& W : Wbosons) {
for (const Particle& child : W.first.children()) {
// find lepton bosons
if (child.abspid() == 11 || child.abspid() == 13) {
leptonic[iboson] = true;
break;
}
// veto events with W-> tau decays
else if (child.abspid() == 15)
vetoEvent;
}
++iboson;
}
// only fully hadronic or semi-leptonic events
if (leptonic[0] && leptonic[1]) vetoEvent;
// 0 SL 1 hadronic
const string iw = (leptonic[0] || leptonic[1]) ? "0"s : "1"s;
_c[_sqs + iw + "0"s]->fill();
Particles particles;
iboson = 0;
for (auto& W : Wbosons) {
if (!leptonic[iboson]) {
particles.insert(particles.end(), W.second.begin(), W.second.end());
}
++iboson;
}
// calculate thrust
Thrust thrust;
thrust.calc(particles);
// type of event
// loop over particles and fill histos
for (const Particle& p : particles) {
if (!PID::isCharged(p.pid())) continue;
// Get momentum of each particle.
const Vector3 mom3 = p.p3();
// Scaled momenta.
const double mom = mom3.mod();
const double xiW = -log(2. * mom / mW);
const double xiL = -log(2. * mom / sqrtS());
// Get momenta components w.r.t. thrust
const double pTinT = dot(mom3, thrust.thrustMajorAxis());
const double pToutT = dot(mom3, thrust.thrustMinorAxis());
double pT = sqrt(sqr(pTinT) + sqr(pToutT));
// fill charged particle hists
if (isWW) {
_h[_sqs + iw]["p_charged"s]->fill(mom);
_h[_sqs + iw]["xi_charged"s]->fill(xiL);
_h[_sqs + iw]["pT_charged"s]->fill(pT / GeV);
_h[_sqs + iw]["p_chargedB"s]->fill(mom);
_h[_sqs + iw]["xi_chargedB"s]->fill(xiL);
_h[_sqs + iw]["pT_chargedB"s]->fill(pT / GeV);
}
// and identified particles
_c[_sqs + iw + "1"s]->fill();
if (p.abspid() == PID::PIPLUS)
_c[_sqs + iw + "2"s]->fill();
else if (p.abspid() == PID::KPLUS)
_c[_sqs + iw + "3"s]->fill();
else if (p.abspid() == PID::PROTON)
_c[_sqs + iw + "4"s]->fill();
if (_sqs == "189"s) {
_h[_sqs + iw]["xi_ident0"]->fill(xiW);
if (p.abspid() == PID::PIPLUS)
_h[_sqs + iw]["xi_ident1"]->fill(xiW);
else if (p.abspid() == PID::KPLUS)
_h[_sqs + iw]["xi_ident2"]->fill(xiW);
else if (p.abspid() == PID::PROTON)
_h[_sqs + iw]["xi_ident3"]->fill(xiW);
}
}
// boosted specta in W rest frame
if (iw == "0"s && _sqs == "189"s) {
iboson = 0;
for (auto& W : Wbosons) {
++iboson;
if (leptonic[iboson - 1]) continue;
// boost to rest frame
LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(W.first.momentum().betaVec());
FourMomentum psum;
for (const Particle& p : W.second) psum += p.mom();
// spectra
for (const Particle& p : W.second) {
if (!PID::isCharged(p.pid())) continue;
// Get momentum of each particle.
FourMomentum phad = boost.transform(p.mom());
const Vector3 mom3 = phad.p3();
// Scaled momenta.
const double mom = mom3.mod();
const double scaledMom = 2. * mom / mW;
const double xi = -log(scaledMom);
// /identified particle spectra
_h[_sqs + iw]["xi_ident20"]->fill(xi);
if (p.abspid() == PID::PIPLUS)
_h[_sqs + iw]["xi_ident21"]->fill(xi);
else if (p.abspid() == PID::KPLUS)
_h[_sqs + iw]["xi_ident22"]->fill(xi);
else if (p.abspid() == PID::PROTON)
_h[_sqs + iw]["xi_ident23"]->fill(xi);
}
}
}
}
}
/// Normalise histograms etc., after the run
void finalize() {
// qq
scale(_h, crossSectionPerEvent());
scale(_q, crossSectionPerEvent());
scale(_c, crossSectionPerEvent());
for (double eVal : allowedEnergies()) {
const int edge = round(eVal);
const string en = toString(edge);
if (_c[en + "0"s]->effNumEntries()) {
scale(_q[en], 1.0 / *_c[en + "0"s]);
for (size_t ih = 1; ih < 7; ++ih) {
if (!_c[en + toString(ih)]->effNumEntries()) continue;
// hist for axis
for (auto& b : _e[ih]->bins()) {
if (edge == b.xEdge()) {
b = *_c[en + toString(ih)] / *_c[en + "0"s];
}
}
}
}
if (edge != 183 && edge != 189) continue;
for (size_t ih = 0; ih < 2; ++ih) {
const string iw = toString(ih);
if (!_c[en + iw + "0"s]->effNumEntries()) continue;
const double sf = 1.0 / _c[en + iw + "0"s]->val();
scale(_h[en + iw], sf);
// charged multis
Estimate0D ratio = *_c[en + toString(iw) + "1"s] / *_c[en + toString(iw) + "0"s];
_e[100 + ih]->binAt(edge).set(ratio.val(), ratio.errPos());
if (en != "189"s) continue;
for (size_t iy = 1; iy < 4; ++iy) {
Estimate0D R = *_c[en + toString(iw) + toString(iy + 1)] / *_c[en + toString(iw) + "0"s];
for (auto& b : _e[10 + 2 * iy - ih]->bins()) {
b.set(R.val(), R.errPos());
}
}
}
// difference histos
for (const string& obs : {"p"s, "xi"s, "pT"s}) {
YODA::Estimate1D htemp = _h[en + "0"][obs + "_chargedB"]->mkEstimate();
htemp.scale(2.0);
const string hpath = _d[en + obs]->path();
*_d[en + obs] = *_h[en + "1"][obs + "_chargedB"] - htemp;
_d[en + obs]->setPath(hpath);
}
}
}
struct ParticleOrdering {
bool operator()(const Particle& p1, const Particle& p2) const {
return p1.pid() > p2.pid();
}
};
///@}
/// @name Histograms
///@{
// qqbar
map<string, map<string, Histo1DPtr>> _h, _q;
map<size_t, BinnedEstimatePtr<int>> _e;
map<string, Estimate1DPtr> _d;
map<string, CounterPtr> _c;
string _sqs = "";
bool isWW;
///@}
};
RIVET_DECLARE_PLUGIN(DELPHI_2001_I526164);
} ```