Rivet analyses

title: L3_2004_I652683

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^- \to q \bar{q}$ 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++ // -- C++ --

include "Rivet/Analysis.hh"

include "Rivet/Projections/Beam.hh"

include "Rivet/Projections/ChargedFinalState.hh"

include "Rivet/Projections/FastJets.hh"

include "Rivet/Projections/FinalState.hh"

include "Rivet/Projections/Hemispheres.hh"

include "Rivet/Projections/ParisiTensor.hh"

include "Rivet/Projections/Thrust.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);

} ```