rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

OPAL_1993_S2692198

Measurement of photon production at LEP 1
Experiment: OPAL (LEP Run 1)
Inspire ID: 343181
Status: UNVALIDATED
Authors:
  • Peter Richardson
References: 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_S2692198.cc
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/Beam.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
#include "Rivet/Projections/FastJets.hh"
#include "fastjet/JadePlugin.hh"

namespace Rivet {


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

    /// Constructor
    OPAL_1993_S2692198()
      : Analysis("OPAL_1993_S2692198")
    {    }


    /// @name Analysis methods
    //@{

    void analyze(const Event& e) {
      // Get event weight for histo filling
      const double weight = e.weight();

      // Extract the photons
      Particles photons;
      Particles nonPhotons;
      FourMomentum ptotal;
      const FinalState& fs = apply<FinalState>(e, "FS");
      foreach (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
      foreach (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 = 0; j < _nPhotonDurham->numBins(); ++j) {
          bool accept(true);
          double ycut = _nPhotonDurham->bin(j).xMid(); ///< @todo Should this be xMin?
          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]));
            double cost = pjet.p3().unit().dot(pgamma.p3().unit());
            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(ycut, weight*_nPhotonDurham->bin(j).xWidth());
          size_t njet = min(size_t(4), exclusive_jets.size()) - 1;
          if (j < _nPhotonJetDurham[njet]->numBins()) {
            _nPhotonJetDurham[njet]->fillBin(j, weight*_nPhotonJetDurham[njet]->bin(j).xWidth());
          }
        }
        // Run the jet clustering JADE
        fastjet::ClusterSequence clust_seq2(input_particles, jade_def);
        for (size_t j = 0; j < _nPhotonJade->numBins(); ++j) {
          bool accept(true);
          double ycut = _nPhotonJade->bin(j).xMid(); ///< @todo Should this be xMin?
          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(ycut, weight*_nPhotonJade->bin(j).xWidth());
          size_t njet = min(size_t(4), exclusive_jets.size()) - 1;
          if (j < _nPhotonJetJade[njet]->numBins()) {
            _nPhotonJetJade[njet]->fillBin(j, weight*_nPhotonJetJade[njet]->bin(j).xWidth());
          }
        }
        // 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
      _nPhotonJade   = bookHisto1D(1, 1, 1);
      _nPhotonDurham = bookHisto1D(2, 1, 1);
      for (size_t ix = 0; ix < 4; ++ix) {
        _nPhotonJetJade  [ix] = bookHisto1D(3 , 1, 1+ix);
        _nPhotonJetDurham[ix] = bookHisto1D(4 , 1, 1+ix);
      }
    }


    /// Finalize
    void finalize() {
      const double fact = 1000/sumOfWeights();
      scale(_nPhotonJade, fact);
      scale(_nPhotonDurham, fact);
      for (size_t ix = 0; ix < 4; ++ix) {
        scale(_nPhotonJetJade  [ix],fact);
        scale(_nPhotonJetDurham[ix],fact);
      }
    }

    //@}

  private:

    Histo1DPtr _nPhotonJade;
    Histo1DPtr _nPhotonDurham;
    Histo1DPtr _nPhotonJetJade  [4];
    Histo1DPtr _nPhotonJetDurham[4];

  };



  // The hook for the plugin system
  DECLARE_RIVET_PLUGIN(OPAL_1993_S2692198);

}