Rivet analyses
Jet rates and event shapes at LEP I+II
Experiment: L3 (LEP I+II)
Inspire ID: 652683
Status: VALIDATED
Authors: - Adil Jueid
References: - Phys.Rept. 399 (2004) 71-174 - 10.1016/j.physrep.2004.07.002 - arXiv: hep-ex/0406049
Beams: e+ e-
Beam energies: (20.7, 20.7); (27.6, 27.6); (32.7, 32.7); (37.9, 37.9); (41.1, 41.1); (42.5, 42.5); (45.6, 45.6); (65.0, 65.0); (68.0, 68.0); (80.7, 80.7); (86.2, 86.2); (91.4, 91.4); (94.3, 94.3); (97.2, 97.2); (100.1, 100.1); (103.1, 103.1)GeV
Run details: - e+e− → qq̄ at several energies.
Event shapes, charged particle multiplicity and scaled momentum distributions, flavour separated at the Z-boson peak. The analysis also includes non-flavour separated event shapes, charged particle multiplicity and scaled momentum distributions at events both above and below the Z pole. There are also jet distributions in the paper which are not currently implemented.
Source
code:L3_2004_I652683.cc
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/Beam.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
#include "Rivet/Projections/Thrust.hh"
#include "Rivet/Projections/ParisiTensor.hh"
#include "Rivet/Projections/Hemispheres.hh"
#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT
#include "Rivet/Projections/InitialQuarks.hh"
#include "fastjet/EECambridgePlugin.hh"
namespace Rivet {
/// Jet rates and event shapes at LEP I+II
class L3_2004_I652683 : public Analysis {
public:
/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(L3_2004_I652683);
/// Book histograms and initialise projections before the run
void init() {
// Projections to use
const FinalState FS;
declare(FS, "FS");
declare(Beam(), "beams");
const ChargedFinalState CFS;
declare(CFS, "CFS");
const Thrust thrust(FS);
declare(thrust, "thrust");
declare(ParisiTensor(FS), "Parisi");
declare(Hemispheres(thrust), "Hemispheres");
declare(InitialQuarks(), "initialquarks");
FastJets jadeJets = FastJets(FS, JetAlg::JADE, -1, JetMuons::ALL, JetInvisibles::DECAY);
FastJets durhamJets = FastJets(FS, JetAlg::DURHAM, 0.7, JetMuons::ALL, JetInvisibles::DECAY);
FastJets cambridgeJets = FastJets(FS, JetAlg::CAM, 0.7, JetMuons::ALL, JetInvisibles::DECAY);
declare(jadeJets, "JadeJets");
declare(durhamJets, "DurhamJets");
// Book the histograms
size_t ih = 0, iy = 1, jh = 0, jy = 1;
for (double eVal : allowedEnergies()) {
const string en = toString(round(eVal/MeV));
if (isCompatibleWithSqrtS(eVal)) _sqs = en;
book(_w[en], "_sumW_"+en);
if (eVal < 90.) {
book(_n[en+"thrust"], 21+ih, 1, iy);
book(_n[en+"rho"], 26+ih, 1, iy);
book(_n[en+"B_T"], 31+ih, 1, iy);
book(_n[en+"B_W"], 36+ih, 1, iy);
}
else if (eVal > 120.) {
book(_n[en+"thrust"], 23+ih, 1, iy);
book(_n[en+"rho"], 28+ih, 1, iy);
book(_n[en+"B_T"], 33+ih, 1, iy);
book(_n[en+"B_W"], 38+ih, 1, iy);
book(_n[en+"C"], 41+ih, 1, iy);
book(_n[en+"D"], 44+ih, 1, iy);
book(_n[en+"xi"], 66+ih, 1, iy);
// and the jets
size_t j = 3*ih+iy;
book(_s[en+"y_2_JADE"], j, 1, 1);
book(_s[en+"y_3_JADE"], j, 1, 2);
book(_s[en+"y_4_JADE"], j, 1, 3);
book(_s[en+"y_5_JADE"], j, 1, 4);
book(_s[en+"y_2_Durham"], 9+j, 1, 1);
book(_s[en+"y_3_Durham"], 9+j, 1, 2);
book(_s[en+"y_4_Durham"], 9+j, 1, 3);
book(_s[en+"y_5_Durham"], 9+j, 1, 4);
if (j==8 || j==9) {
book(_s[en+"y_2_Cambridge"], 11+j, 1, 1);
book(_s[en+"y_3_Cambridge"], 11+j, 1, 2);
book(_s[en+"y_4_Cambridge"], 11+j, 1, 3);
book(_s[en+"y_5_Cambridge"], 11+j, 1, 4);
}
book(_d[en+"mult"], 60+jh, 1, jy);
if (jy==3) {
++jh; jy = 0;
}
++jy;
}
else {
book(_c["l"], "_sumW_udsc");
book(_c["b"], "_sumW_b");
book(_h["l"]["Thrust"], 47, 1, 1);
book(_h["b"]["Thrust"], 47, 1, 2);
book(_h["l"]["heavyJetmass"], 48, 1, 1);
book(_h["b"]["heavyJetmass"], 48, 1, 2);
book(_h["l"]["totalJetbroad"], 49, 1, 1);
book(_h["b"]["totalJetbroad"], 49, 1, 2);
book(_h["l"]["wideJetbroad"], 50, 1, 1);
book(_h["b"]["wideJetbroad"], 50, 1, 2);
book(_h["l"]["Cparameter"], 51, 1, 1);
book(_h["b"]["Cparameter"], 51, 1, 2);
book(_h["l"]["Dparameter"], 52, 1, 1);
book(_h["b"]["Dparameter"], 52, 1, 2);
book(_h["N"]["scaledMomentum"], 65, 1, 1);
book(_h["Nl"]["scaledMomentum"], 65, 1, 2);
book(_h["Nb"]["scaledMomentum"], 65, 1, 3);
book(_d["N"], 59, 1, 1);
book(_d["Nl"], 59, 1, 2);
book(_d["Nb"], 59, 1, 3);
ih = iy = 0;
}
if (iy==3) {
++ih; iy = 0;
}
++iy;
}
raiseBeamErrorIf(_sqs.empty());
}
/// Perform the per-event analysis
void analyze(const Event& event) {
_w[_sqs]->fill();
// Get beam average momentum
const ParticlePair& beams = apply<Beam>(event, "beams").beams();
const double beamMomentum = 0.5*(beams.first.p3().mod() + beams.second.p3().mod());
// InitialQuarks projection to have udsc events separated from b events
/// @todo Yuck!!! Eliminate when possible...
int iflav = _sqs.size() == 6? 6 : 0;
// only need the flavour at Z pole
if (_sqs == "91200"s) {
int flavour = 0;
const InitialQuarks& iqf = apply<InitialQuarks>(event, "initialquarks");
Particles quarks;
if (iqf.particles().size() == 2) {
flavour = iqf.particles().front().abspid();
quarks = iqf.particles();
}
else {
map<int, Particle> quarkmap;
for (const Particle& p : iqf.particles()) {
if (quarkmap.find(p.pid()) == quarkmap.end()) quarkmap[p.pid()] = p;
else if (quarkmap[p.pid()].E() < p.E()) quarkmap[p.pid()] = p;
}
double max_energy = 0.;
for (int i = 1; i <= 5; ++i) {
double energy = 0.;
if (quarkmap.find(i) != quarkmap.end()) {
energy += quarkmap[ i].E();
}
if (quarkmap.find(-i) != quarkmap.end()) {
energy += quarkmap[-i].E();
}
if (energy > max_energy) {
flavour = i;
}
}
if (quarkmap.find(flavour) != quarkmap.end()) {
quarks.push_back(quarkmap[flavour]);
}
if (quarkmap.find(-flavour) != quarkmap.end()) {
quarks.push_back(quarkmap[-flavour]);
}
}
// Flavour label
/// @todo Change to a bool?
if (flavour == PID::DQUARK || flavour == PID::UQUARK ||
flavour == PID::SQUARK || flavour == PID::CQUARK) {
iflav = 1;
_c["l"]->fill();
}
else if (flavour == PID::BQUARK) {
iflav = 5;
_c["b"]->fill();
}
}
// Charged multiplicity
const FinalState& cfs = apply<FinalState>(event, "CFS");
if (_sqs == "91200"s) {
_d["N"]->fill(cfs.size());
if (iflav == 1) _d["Nl"]->fill(cfs.size());
else if (iflav == 5) _d["Nb"]->fill(cfs.size());
}
else if (iflav) _d[_sqs+"mult"]->fill(cfs.size());
// Scaled momentum
const Particles& chparticles = cfs.particlesByPt();
for (const Particle& p : chparticles) {
const Vector3 momentum3 = p.p3();
const double mom = momentum3.mod();
const double scaledMom = mom/beamMomentum;
const double logScaledMom = std::log(scaledMom);
if (_sqs == "91200"s) {
_h["N"]["scaledMomentum"]->fill(-logScaledMom);
if (iflav == 1) _h["Nl"]["scaledMomentum"]->fill(-logScaledMom);
else if (iflav == 5) _h["Nb"]["scaledMomentum"]->fill(-logScaledMom);
}
else if (iflav) _n[_sqs+"xi"]->fill(-logScaledMom);
}
// Thrust
const Thrust& thrust = apply<Thrust>(event, "thrust");
if (_sqs == "91200"s) {
if (iflav == 1) _h["l"]["Thrust"]->fill(thrust.thrust());
else if (iflav == 5) _h["b"]["Thrust"]->fill(thrust.thrust());
}
else {
_n[_sqs+"thrust"]->fill(1.-thrust.thrust());
}
// C and D Parisi parameters
const ParisiTensor& parisi = apply<ParisiTensor>(event, "Parisi");
if (_sqs == "91200"s) {
if (iflav == 1) {
_h["l"]["Cparameter"]->fill(parisi.C());
_h["l"]["Dparameter"]->fill(parisi.D());
} else if (iflav == 5) {
_h["b"]["Cparameter"]->fill(parisi.C());
_h["b"]["Dparameter"]->fill(parisi.D());
}
}
else if (iflav) {
_n[_sqs+"C"]->fill(parisi.C());
_n[_sqs+"D"]->fill(parisi.D());
}
// The hemisphere variables
const Hemispheres& hemisphere = apply<Hemispheres>(event, "Hemispheres");
if (_sqs == "91200"s) {
if (iflav == 1) {
_h["l"]["heavyJetmass"]->fill(hemisphere.scaledM2high());
_h["l"]["totalJetbroad"]->fill(hemisphere.Bsum());
_h["l"]["wideJetbroad"]->fill(hemisphere.Bmax());
} else if (iflav == 5) {
_h["b"]["heavyJetmass"]->fill(hemisphere.scaledM2high());
_h["b"]["totalJetbroad"]->fill(hemisphere.Bsum());
_h["b"]["wideJetbroad"]->fill(hemisphere.Bmax());
}
}
else {
_n[_sqs+"rho"]->fill(hemisphere.scaledM2high());
_n[_sqs+"B_T"]->fill(hemisphere.Bsum());
_n[_sqs+"B_W"]->fill(hemisphere.Bmax());
}
// jade jet rates
if (_s.count(_sqs+"y_2_JADE")) {
const FastJets& jadejet = apply<FastJets>(event, "JadeJets");
if (jadejet.clusterSeq()) {
const double y_23 = jadejet.clusterSeq()->exclusive_ymerge_max(2);
const double y_34 = jadejet.clusterSeq()->exclusive_ymerge_max(3);
const double y_45 = jadejet.clusterSeq()->exclusive_ymerge_max(4);
const double y_56 = jadejet.clusterSeq()->exclusive_ymerge_max(5);
for (auto& b : _s[_sqs+"y_2_JADE"]->bins()) {
const double ycut = stod(b.xEdge());
if (y_23 < ycut) _s[_sqs+"y_2_JADE"]->fill(b.xEdge());
}
for (auto& b : _s[_sqs+"y_3_JADE"]->bins()) {
const double ycut = stod(b.xEdge());
if (y_34 < ycut && y_23 > ycut) {
_s[_sqs+"y_3_JADE"]->fill(b.xEdge());
}
}
for (auto& b : _s[_sqs+"y_4_JADE"]->bins()) {
const double ycut = stod(b.xEdge());
if (y_45 < ycut && y_34 > ycut) {
_s[_sqs+"y_4_JADE"]->fill(b.xEdge());
}
}
for (auto& b : _s[_sqs+"y_5_JADE"]->bins()) {
const double ycut = stod(b.xEdge());
if (y_56 < ycut && y_45 > ycut) {
_s[_sqs+"y_5_JADE"]->fill(b.xEdge());
}
}
}
}
// Durham jet rates
if (_s.count(_sqs+"y_2_Durham")) {
const FastJets& durhamjet = apply<FastJets>(event, "DurhamJets");
if (durhamjet.clusterSeq()) {
const double y_23 = durhamjet.clusterSeq()->exclusive_ymerge_max(2);
const double y_34 = durhamjet.clusterSeq()->exclusive_ymerge_max(3);
const double y_45 = durhamjet.clusterSeq()->exclusive_ymerge_max(4);
const double y_56 = durhamjet.clusterSeq()->exclusive_ymerge_max(5);
for (auto& b : _s[_sqs+"y_2_Durham"]->bins()) {
const double ycut = stod(b.xEdge());
if (y_23 < ycut) _s[_sqs+"y_2_Durham"]->fill(b.xEdge());
}
for (auto& b : _s[_sqs+"y_3_Durham"]->bins()) {
const double ycut = stod(b.xEdge());
if (y_34 < ycut && y_23 > ycut) {
_s[_sqs+"y_3_Durham"]->fill(b.xEdge());
}
}
for (auto& b : _s[_sqs+"y_4_Durham"]->bins()) {
const double ycut = stod(b.xEdge());
if (y_45 < ycut && y_34 > ycut) {
_s[_sqs+"y_4_Durham"]->fill(b.xEdge());
}
}
for (auto& b : _s[_sqs+"y_5_Durham"]->bins()) {
const double ycut = stod(b.xEdge());
if (y_56 < ycut && y_45 > ycut) {
_s[_sqs+"y_5_Durham"]->fill(b.xEdge());
}
}
}
}
// Cambridge
if (_s.count(_sqs+"y_2_Cambridge")) {
PseudoJets pjs;
const FinalState& fs = apply<FinalState>(event, "FS");
for (size_t i = 0; i < fs.particles().size(); ++i) {
fastjet::PseudoJet pj = fs.particles()[i];
pjs.push_back(pj);
}
for (size_t i = 1; i < _s[_sqs+"y_2_Cambridge"]->numBins()+1; ++i) {
const auto& b = _s[_sqs+"y_2_Cambridge"]->bin(i);
const double ycut = stod(b.xEdge());
fastjet::EECambridgePlugin plugin(ycut);
fastjet::JetDefinition jdef(&plugin);
fastjet::ClusterSequence cseq(pjs, jdef);
size_t njet = cseq.inclusive_jets().size();
if (njet==2) {
_s[_sqs+"y_2_Cambridge"]->fill(b.xEdge());
}
else if (njet==3) {
if (i < _s[_sqs+"y_3_Cambridge"]->numBins()+1) {
_s[_sqs+"y_3_Cambridge"]->fill(b.xEdge());
}
}
else if(njet==4) {
if (i < _s[_sqs+"y_4_Cambridge"]->numBins()+1) {
_s[_sqs+"y_4_Cambridge"]->fill(b.xEdge());
}
}
else if(njet==5) {
if (i < _s[_sqs+"y_5_Cambridge"]->numBins()+1) {
_s[_sqs+"y_5_Cambridge"]->fill(b.xEdge());
}
}
}
}
}
/// Normalise histograms etc., after the run
void finalize() {
// Z pole plots
scale(_c, crossSectionPerEvent());
scale(_w, crossSectionPerEvent());
scale(_h, crossSectionPerEvent());
scale(_d, crossSectionPerEvent());
scale(_n, crossSectionPerEvent());
scale(_s, crossSectionPerEvent());
for (const auto& item : _c) {
const double c = item.second->sumW();
if (c) scale(_h[item.first], 1.0/c);
}
for (const auto& item : _d) {
const double c = item.second->sumW();
if (c) scale(_h[item.first]["scaledMomentum"], 1.0/c);
}
normalize(_d, 1.0, true);
for (auto& item : _n) {
size_t pos = item.first.find("xi");
if (pos != string::npos) {
const double w = _w[item.first.substr(0,pos)]->sumW();
if (!isZero(w)) scale(item.second, 1.0/w);
}
else normalize(item.second);
}
// the jets
for (auto& item : _s) {
const double w = _w[item.first.substr(0,6)]->sumW();
if (!isZero(w)) scale(item.second, 1.0/w);
}
}
/// @name Histograms
/// @{
map<string, CounterPtr> _c, _w;
map<string, Histo1DPtr> _n;
map<string, map<string, Histo1DPtr>> _h;
map<string, BinnedHistoPtr<string>> _s;
map<string, BinnedHistoPtr<int>> _d;
string _sqs = "";
/// @}
};
RIVET_DECLARE_PLUGIN(L3_2004_I652683);
}