Rivet analyses

title: OPAL_1993_I343181

Measurement of photon production at LEP 1

Experiment: OPAL (LEP Run 1)

Inspire ID: 343181

Status: UNVALIDATED

Authors: - Peter Richardson

References: - Z.Phys.C58:405-418,1993 - DOI: 10.1007/BF01557697

Beams: e+ e-

Beam energies: (45.6, 45.6)GeV

Run details: - $e^+ e^- \to$ jet jet (+ photons)

Measurement of the production of photons in $e^+ e^-\to q \bar q$ events at LEP 1.

Source code:OPAL_1993_I343181.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 "fastjet/JadePlugin.hh"

namespace Rivet {

/// @brief OPAL photon production /// /// @author Peter Richardson class OPAL_1993_I343181 : public Analysis { public:

RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_1993_I343181);


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

void analyze(const Event& e) {
  // Extract the photons
  Particles photons;
  Particles nonPhotons;
  FourMomentum ptotal;
  const FinalState& fs = apply<FinalState>(e, "FS");
  for (const Particle& p : fs.particles()) {
    ptotal += p.momentum();
    if (p.pid() == PID::PHOTON) {
      photons.push_back(p);
    }
    else {
      nonPhotons.push_back(p);
    }
  }
  // No photon return but still count for normalisation
  if (photons.empty()) return;
  // Definition of the Durham algorithm
  fastjet::JetDefinition durham_def(fastjet::ee_kt_algorithm, fastjet::E_scheme, fastjet::Best);
  // Definition of the JADE algorithm
  fastjet::JadePlugin jade;
  fastjet::JetDefinition jade_def = fastjet::JetDefinition(&jade);
  // Now for the weird jet algorithm
  double evis = ptotal.mass();
  vector<fastjet::PseudoJet> input_particles;
  // Pseudo-jets from the non photons
  for (const Particle& p : nonPhotons) {
    const FourMomentum p4 = p.momentum();
    input_particles.push_back(fastjet::PseudoJet(p4.px(), p4.py(), p4.pz(), p4.E()));
  }
  // Pseudo-jets from all bar the first photon
  for (size_t ix = 1; ix < photons.size(); ++ix) {
    const FourMomentum p4 = photons[ix].momentum();
    input_particles.push_back(fastjet::PseudoJet(p4.px(), p4.py(), p4.pz(), p4.E()));
  }
  // Now loop over the photons
  for (size_t ix = 0; ix < photons.size(); ++ix) {
    FourMomentum pgamma = photons[ix].momentum();
    // Run the jet clustering DURHAM
    fastjet::ClusterSequence clust_seq(input_particles, durham_def);
    // Cluster the jets
    for (size_t j = 1; j < _nPhotonDurham->numBinsX() + 1; ++j) {
      bool accept(true);
      const string edge = _nPhotonDurham->bin(j).xEdge();
      const double ycut = std::stod(edge);
      const double dcut = sqr(evis) * ycut;
      vector<fastjet::PseudoJet> exclusive_jets = sorted_by_E(clust_seq.exclusive_jets(dcut));
      for (size_t iy = 0; iy < exclusive_jets.size(); ++iy) {
        FourMomentum pjet(momentum(exclusive_jets[iy]));
        const double cost = pjet.p3().unit().dot(pgamma.p3().unit());
        const double ygamma = 2 * min(sqr(pjet.E() / evis), sqr(pgamma.E() / evis)) * (1 - cost);
        if (ygamma < ycut) {
          accept = false;
          break;
        }
      }
      if (!accept) continue;
      _nPhotonDurham->fill(edge);
      size_t njet = min(size_t(4), exclusive_jets.size()) - 1;
      if (j < _nPhotonJetDurham[njet]->numBins() + 1) {
        const auto& b = _nPhotonJetDurham[njet]->bin(j);
        _nPhotonJetDurham[njet]->fill(b.xEdge());
      }
    }
    // Run the jet clustering JADE
    fastjet::ClusterSequence clust_seq2(input_particles, jade_def);
    for (size_t j = 1; j < _nPhotonJade->numBinsX() + 1; ++j) {
      bool accept(true);
      const string edge = _nPhotonJade->bin(j).xEdge();
      const double ycut = std::stod(edge);
      const double dcut = sqr(evis) * ycut;
      vector<fastjet::PseudoJet> exclusive_jets = sorted_by_E(clust_seq2.exclusive_jets(dcut));
      for (size_t iy = 0; iy < exclusive_jets.size(); ++iy) {
        FourMomentum pjet(momentum(exclusive_jets[iy]));
        double cost = pjet.p3().unit().dot(pgamma.p3().unit());
        double ygamma = 2. * pjet.E() * pgamma.E() / sqr(evis) * (1. - cost);
        if (ygamma < ycut) {
          accept = false;
          break;
        }
      }
      if (!accept) continue;
      /// @todo Really want to use a "bar graph" here (i.e. ignore bin width)
      _nPhotonJade->fill(edge);
      size_t njet = min(size_t(4), exclusive_jets.size()) - 1;
      if (j < _nPhotonJetJade[njet]->numBins() + 1) {
        const auto& b = _nPhotonJetJade[njet]->bin(j);
        _nPhotonJetJade[njet]->fill(b.xEdge());
      }
    }
    // Add this photon back in for the next iteration of the loop
    if (ix + 1 != photons.size()) {
      input_particles[nonPhotons.size() + ix] = fastjet::PseudoJet(pgamma.px(), pgamma.py(), pgamma.pz(),
                                                                   pgamma.E());
    }
  }
}


void init() {
  // Projections
  declare(FinalState(), "FS");

  // Book datasets
  book(_nPhotonJade, 1, 1, 1);
  book(_nPhotonDurham, 2, 1, 1);
  for (size_t ix = 0; ix < 4; ++ix) {
    book(_nPhotonJetJade[ix], 3, 1, 1 + ix);
    book(_nPhotonJetDurham[ix], 4, 1, 1 + ix);
  }
}


/// Finalize
void finalize() {
  const double fact = 1000 / sumOfWeights();
  scale(_nPhotonJade, fact);
  scale(_nPhotonDurham, fact);
  scale(_nPhotonJetJade, fact);
  scale(_nPhotonJetDurham, fact);
}

/// @}

private:

BinnedHistoPtr<string> _nPhotonJade;
BinnedHistoPtr<string> _nPhotonDurham;
BinnedHistoPtr<string> _nPhotonJetJade[4];
BinnedHistoPtr<string> _nPhotonJetDurham[4];

};

RIVET_DECLARE_ALIASED_PLUGIN(OPAL_1993_I343181, OPAL_1993_S2692198);

} ```

Aliases: - OPAL_1993_S2692198