Rivet analyses

Event shape distributions and moments in e+e -> hadrons at 91–209 GeV

Experiment: OPAL (LEP 1 & 2)

Inspire ID: 669402

Status: VALIDATED

Authors: - Andy Buckley

References: - Eur.Phys.J.C40:287-316,2005 - arXiv: hep-ex/0503051

Beams: e+ e-

Beam energies: (45.6, 45.6); (66.5, 66.5); (88.5, 88.5); (98.5, 98.5)GeV

Run details: - Hadronic e+e events at 4 representative energies (91, 133, 177, 197). Runs need to have ISR suppressed, since the analysis was done using a cut of $\sqrt{s} - \sqrt{s_\text{reco}} < 1\,\text{GeV}$. Particles with a lifetime  > 3 ⋅ 10−10 s are considered to be stable.

Measurement of e+e event shape variable distributions and their 1st to 5th moments in LEP running from the Z pole to the highest LEP 2 energy of 209 GeV.

Source code:OPAL_2004_I669402.cc

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/Beam.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
#include "Rivet/Projections/Sphericity.hh"
#include "Rivet/Projections/Thrust.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/ParisiTensor.hh"
#include "Rivet/Projections/Hemispheres.hh"
#include <cmath>

namespace Rivet {


  /// @brief OPAL event shapes and moments at 91, 133, 177, and 197 GeV
  ///
  /// @author Andy Buckley
  class OPAL_2004_I669402 : public Analysis {
  public:

    RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_2004_I669402);


    /// @name Analysis methods
    /// @{

    void init() {
      // Projections
      declare(Beam(), "Beams");
      const FinalState fs;
      declare(fs, "FS");
      const ChargedFinalState cfs;
      declare(cfs, "CFS");
      declare(FastJets(fs, JetAlg::DURHAM, 0.7), "DurhamJets");
      declare(Sphericity(fs), "Sphericity");
      declare(ParisiTensor(fs), "Parisi");
      const Thrust thrust(fs);
      declare(thrust, "Thrust");
      declare(Hemispheres(thrust), "Hemispheres");

      size_t ih = 0;
      for (double eVal : allowedEnergies()) {

        const string en = toString(round(eVal));
        if (isCompatibleWithSqrtS(eVal))  _sqs = en;

        book(_h[en]["1MinusT"],     1, 1, ih+1);
        book(_h[en]["HemiMassH"],   2, 1, ih+1);
        book(_h[en]["CParam"],      3, 1, ih+1);
        book(_h[en]["HemiBroadT"],  4, 1, ih+1);
        book(_h[en]["HemiBroadW"],  5, 1, ih+1);
        book(_h[en]["TMajor"],      7, 1, ih+1);
        book(_h[en]["TMinor"],      8, 1, ih+1);
        book(_h[en]["Aplanarity"],  9, 1, ih+1);
        book(_h[en]["Sphericity"], 10, 1, ih+1);
        book(_h[en]["Oblateness"], 11, 1, ih+1);
        book(_h[en]["HemiMassL"],  12, 1, ih+1);
        book(_h[en]["HemiBroadN"], 13, 1, ih+1);
        book(_h[en]["DParam"],     14, 1, ih+1);
        book(_h["jet"+en]["Y23Durham"], 6, 1, ih+1);

        book(_d[en]["1MinusTMom"],    15, 1, ih+1);
        book(_d[en]["HemiMassHMom"],  16, 1, ih+1);
        book(_d[en]["CParamMom"],     17, 1, ih+1);
        book(_d[en]["HemiBroadTMom"], 18, 1, ih+1);
        book(_d[en]["HemiBroadWMom"], 19, 1, ih+1);
        book(_d[en]["TMajorMom"]    , 21, 1, ih+1);
        book(_d[en]["TMinorMom"]    , 22, 1, ih+1);
        book(_d[en]["SphericityMom"], 23, 1, ih+1);
        book(_d[en]["OblatenessMom"], 24, 1, ih+1);
        book(_d[en]["HemiMassLMom"] , 25, 1, ih+1);
        book(_d[en]["HemiBroadNMom"], 26, 1, ih+1);
        book(_d["jet"+en]["Y23DurhamMom"], 20, 1, ih+1);

        book(_c[en], "_sumWTrack2_"+en);
        book(_c["jet"+en], "_sumWJet3_"+en);
        ++ih;
      }
      raiseBeamErrorIf(_sqs.empty());
    }


    void analyze(const Event& event) {
      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
      const FinalState& cfs = apply<FinalState>(event, "CFS");
      if (cfs.size() < 2) vetoEvent;

      _c[_sqs]->fill();

      // Thrusts
      const Thrust& thrust = apply<Thrust>(event, "Thrust");
      _h[_sqs]["1MinusT"]->fill(1-thrust.thrust());
      _h[_sqs]["TMajor"]->fill(thrust.thrustMajor());
      _h[_sqs]["TMinor"]->fill(thrust.thrustMinor());
      _h[_sqs]["Oblateness"]->fill(thrust.oblateness());
      for (int n = 1; n <= 5; ++n) {
        _d[_sqs]["1MinusTMom"]->fill(n, intpow(1-thrust.thrust(), n));
        _d[_sqs]["TMajorMom"]->fill(n, intpow(thrust.thrustMajor(), n));
        _d[_sqs]["TMinorMom"]->fill(n, intpow(thrust.thrustMinor(), n));
        _d[_sqs]["OblatenessMom"]->fill(n, intpow(thrust.oblateness(), n));
      }

      // Jets
      const FastJets& durjet = apply<FastJets>(event, "DurhamJets");
      if (durjet.clusterSeq()) {
        _c["jet"+_sqs]->fill();
        const double y23 = durjet.clusterSeq()->exclusive_ymerge_max(2);
        if (y23>0.0) {
          _h["jet"+_sqs]["Y23Durham"]->fill(y23);
          for (int n = 1; n <= 5; ++n) {
            _d["jet"+_sqs]["Y23DurhamMom"]->fill(n, intpow(y23, n));
          }
        }
      }

      // Sphericities
      const Sphericity& sphericity = apply<Sphericity>(event, "Sphericity");
      const double sph = sphericity.sphericity();
      const double apl = sphericity.aplanarity();
      _h[_sqs]["Sphericity"]->fill(sph);
      _h[_sqs]["Aplanarity"]->fill(apl);
      for (int n = 1; n <= 5; ++n) {
        _d[_sqs]["SphericityMom"]->fill(n, intpow(sph, n));
      }

      // C & D params
      const ParisiTensor& parisi = apply<ParisiTensor>(event, "Parisi");
      const double cparam = parisi.C();
      const double dparam = parisi.D();
      _h[_sqs]["CParam"]->fill(cparam);
      _h[_sqs]["DParam"]->fill(dparam);
      for (int n = 1; n <= 5; ++n) {
        _d[_sqs]["CParamMom"]->fill(n, intpow(cparam, n));
      }

      // Hemispheres
      const Hemispheres& hemi = apply<Hemispheres>(event, "Hemispheres");
      // The paper says that M_H/L are scaled by sqrt(s), but scaling by E_vis is the way that fits the data...
      const double hemi_mh = hemi.scaledMhigh();
      const double hemi_ml = hemi.scaledMlow();
      /// @todo This shouldn't be necessary... what's going on? Memory corruption suspected :(
      // if (std::isnan(hemi_ml)) {
      //   MSG_ERROR("NaN in HemiL! Event = " << numEvents());
      //   MSG_ERROR(hemi.M2low() << ", " << hemi.E2vis());
      // }
      if (!std::isnan(hemi_mh) && !std::isnan(hemi_ml)) {
        const double hemi_bmax = hemi.Bmax();
        const double hemi_bmin = hemi.Bmin();
        const double hemi_bsum = hemi.Bsum();
        _h[_sqs]["HemiMassH"]->fill(hemi_mh);
        _h[_sqs]["HemiMassL"]->fill(hemi_ml);
        _h[_sqs]["HemiBroadW"]->fill(hemi_bmax);
        _h[_sqs]["HemiBroadN"]->fill(hemi_bmin);
        _h[_sqs]["HemiBroadT"]->fill(hemi_bsum);
        for (int n = 1; n <= 5; ++n) {
          // if (std::isnan(pow(hemi_ml, n))) MSG_ERROR("NaN in HemiL moment! Event = " << numEvents());
          _d[_sqs]["HemiMassHMom"]->fill(n, intpow(hemi_mh, n));
          _d[_sqs]["HemiMassLMom"]->fill(n, intpow(hemi_ml, n));
          _d[_sqs]["HemiBroadWMom"]->fill(n, intpow(hemi_bmax, n));
          _d[_sqs]["HemiBroadNMom"]->fill(n, intpow(hemi_bmin, n));
          _d[_sqs]["HemiBroadTMom"]->fill(n, intpow(hemi_bsum, n));
        }
      }
    }


    void finalize() {
      scale(_h, crossSectionPerEvent());
      scale(_d, crossSectionPerEvent());
      scale(_c, crossSectionPerEvent());
      for (const auto& item : _c) {
        if (!item.second->val())  continue;
        scale(_h[item.first], 1.0 / *item.second);
        scale(_d[item.first], 1.0 / *item.second);
      }
    }

    /// @}


  private:

    /// Counters of event weights passing the cuts
    /// @{
    map<string,CounterPtr> _c;
    map<string,map<string,Histo1DPtr>> _h;
    map<string,map<string,BinnedHistoPtr<int>>> _d;
    string _sqs = "";
    /// @}

  };



  RIVET_DECLARE_ALIASED_PLUGIN(OPAL_2004_I669402, OPAL_2004_S6132243);

}

Aliases: - OPAL_2004_S6132243