1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/ChargedFinalState.hh"
  6#include "Rivet/Projections/FastJets.hh"
  7#include "fastjet/JadePlugin.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief OPAL photon production
 13  ///
 14  /// @author Peter Richardson
 15  class OPAL_1993_I343181 : public Analysis {
 16  public:
 17
 18    RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_1993_I343181);
 19
 20
 21    /// @name Analysis methods
 22    /// @{
 23
 24    void analyze(const Event& e) {
 25      // Extract the photons
 26      Particles photons;
 27      Particles nonPhotons;
 28      FourMomentum ptotal;
 29      const FinalState& fs = apply<FinalState>(e, "FS");
 30      for (const Particle& p : fs.particles()) {
 31        ptotal+= p.momentum();
 32        if (p.pid() == PID::PHOTON) {
 33          photons.push_back(p);
 34        } else {
 35          nonPhotons.push_back(p);
 36        }
 37      }
 38      // No photon return but still count for normalisation
 39      if (photons.empty()) return;
 40      // Definition of the Durham algorithm
 41      fastjet::JetDefinition durham_def(fastjet::ee_kt_algorithm, fastjet::E_scheme, fastjet::Best);
 42      // Definition of the JADE algorithm
 43      fastjet::JadePlugin jade;
 44      fastjet::JetDefinition jade_def = fastjet::JetDefinition(&jade);
 45      // Now for the weird jet algorithm
 46      double evis = ptotal.mass();
 47      vector<fastjet::PseudoJet> input_particles;
 48      // Pseudo-jets from the non photons
 49      for (const Particle& p :  nonPhotons) {
 50        const FourMomentum p4 = p.momentum();
 51        input_particles.push_back(fastjet::PseudoJet(p4.px(), p4.py(), p4.pz(), p4.E()));
 52      }
 53      // Pseudo-jets from all bar the first photon
 54      for (size_t ix = 1; ix < photons.size(); ++ix) {
 55        const FourMomentum p4 = photons[ix].momentum();
 56        input_particles.push_back(fastjet::PseudoJet(p4.px(), p4.py(), p4.pz(), p4.E()));
 57      }
 58      // Now loop over the photons
 59      for (size_t ix = 0; ix < photons.size(); ++ix) {
 60        FourMomentum pgamma = photons[ix].momentum();
 61        // Run the jet clustering DURHAM
 62        fastjet::ClusterSequence clust_seq(input_particles, durham_def);
 63        // Cluster the jets
 64        for (size_t j = 1; j < _nPhotonDurham->numBinsX()+1; ++j) {
 65          bool accept(true);
 66          const string edge = _nPhotonDurham->bin(j).xEdge();
 67          const double ycut = std::stod(edge);
 68          const double dcut = sqr(evis)*ycut;
 69          vector<fastjet::PseudoJet> exclusive_jets = sorted_by_E(clust_seq.exclusive_jets(dcut));
 70          for (size_t iy = 0; iy < exclusive_jets.size(); ++iy) {
 71            FourMomentum pjet(momentum(exclusive_jets[iy]));
 72            const double cost = pjet.p3().unit().dot(pgamma.p3().unit());
 73            const double ygamma = 2 * min(sqr(pjet.E()/evis), sqr(pgamma.E()/evis)) * (1 - cost);
 74            if (ygamma < ycut) {
 75              accept = false;
 76              break;
 77            }
 78          }
 79          if (!accept) continue;
 80          _nPhotonDurham->fill(edge);
 81          size_t njet = min(size_t(4), exclusive_jets.size()) - 1;
 82          if (j < _nPhotonJetDurham[njet]->numBins()+1) {
 83            const auto& b = _nPhotonJetDurham[njet]->bin(j);
 84            _nPhotonJetDurham[njet]->fill(b.xEdge());
 85          }
 86        }
 87        // Run the jet clustering JADE
 88        fastjet::ClusterSequence clust_seq2(input_particles, jade_def);
 89        for (size_t j = 1; j < _nPhotonJade->numBinsX()+1; ++j) {
 90          bool accept(true);
 91          const string edge = _nPhotonJade->bin(j).xEdge();
 92          const double ycut = std::stod(edge);
 93          const double dcut = sqr(evis)*ycut;
 94          vector<fastjet::PseudoJet> exclusive_jets = sorted_by_E(clust_seq2.exclusive_jets(dcut));
 95          for (size_t iy = 0; iy < exclusive_jets.size(); ++iy) {
 96            FourMomentum pjet(momentum(exclusive_jets[iy]));
 97            double cost = pjet.p3().unit().dot(pgamma.p3().unit());
 98            double ygamma = 2.*pjet.E()*pgamma.E()/sqr(evis)*(1.-cost);
 99            if (ygamma < ycut) {
100              accept = false;
101              break;
102            }
103          }
104          if (!accept) continue;
105          /// @todo Really want to use a "bar graph" here (i.e. ignore bin width)
106          _nPhotonJade->fill(edge);
107          size_t njet = min(size_t(4), exclusive_jets.size()) - 1;
108          if (j < _nPhotonJetJade[njet]->numBins()+1) {
109            const auto& b = _nPhotonJetJade[njet]->bin(j);
110            _nPhotonJetJade[njet]->fill(b.xEdge());
111          }
112        }
113        // Add this photon back in for the next iteration of the loop
114        if (ix+1 != photons.size()) {
115          input_particles[nonPhotons.size()+ix] = fastjet::PseudoJet(pgamma.px(), pgamma.py(), pgamma.pz(), pgamma.E());
116        }
117      }
118    }
119
120
121    void init() {
122      // Projections
123      declare(FinalState(), "FS");
124
125      // Book datasets
126      book(_nPhotonJade, 1, 1, 1);
127      book(_nPhotonDurham, 2, 1, 1);
128      for (size_t ix = 0; ix < 4; ++ix) {
129        book(_nPhotonJetJade[ix], 3, 1, 1+ix);
130        book(_nPhotonJetDurham[ix], 4, 1, 1+ix);
131      }
132    }
133
134
135    /// Finalize
136    void finalize() {
137      const double fact = 1000/sumOfWeights();
138      scale(_nPhotonJade, fact);
139      scale(_nPhotonDurham, fact);
140      scale(_nPhotonJetJade, fact);
141      scale(_nPhotonJetDurham, fact);
142    }
143
144    /// @}
145
146
147  private:
148
149    BinnedHistoPtr<string> _nPhotonJade;
150    BinnedHistoPtr<string> _nPhotonDurham;
151    BinnedHistoPtr<string> _nPhotonJetJade[4];
152    BinnedHistoPtr<string> _nPhotonJetDurham[4];
153
154  };
155
156
157
158  RIVET_DECLARE_ALIASED_PLUGIN(OPAL_1993_I343181, OPAL_1993_S2692198);
159
160}