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_S2692198 : public Analysis {
 16  public:
 17
 18    RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_1993_S2692198);
 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 = 0; j < _nPhotonDurham->numBins(); ++j) {
 65          bool accept(true);
 66          double ycut = _nPhotonDurham->bin(j).xMid(); ///< @todo Should this be xMin?
 67          double dcut = sqr(evis)*ycut;
 68          vector<fastjet::PseudoJet> exclusive_jets = sorted_by_E(clust_seq.exclusive_jets(dcut));
 69          for (size_t iy = 0; iy < exclusive_jets.size(); ++iy) {
 70            FourMomentum pjet(momentum(exclusive_jets[iy]));
 71            double cost = pjet.p3().unit().dot(pgamma.p3().unit());
 72            double ygamma = 2 * min(sqr(pjet.E()/evis), sqr(pgamma.E()/evis)) * (1 - cost);
 73            if (ygamma < ycut) {
 74              accept = false;
 75              break;
 76            }
 77          }
 78          if (!accept) continue;
 79          _nPhotonDurham->fill(ycut, _nPhotonDurham->bin(j).xWidth());
 80          size_t njet = min(size_t(4), exclusive_jets.size()) - 1;
 81          if (j < _nPhotonJetDurham[njet]->numBins()) {
 82            _nPhotonJetDurham[njet]->fillBin(j, _nPhotonJetDurham[njet]->bin(j).xWidth());
 83          }
 84        }
 85        // Run the jet clustering JADE
 86        fastjet::ClusterSequence clust_seq2(input_particles, jade_def);
 87        for (size_t j = 0; j < _nPhotonJade->numBins(); ++j) {
 88          bool accept(true);
 89          double ycut = _nPhotonJade->bin(j).xMid(); ///< @todo Should this be xMin?
 90          double dcut = sqr(evis)*ycut;
 91          vector<fastjet::PseudoJet> exclusive_jets = sorted_by_E(clust_seq2.exclusive_jets(dcut));
 92          for (size_t iy = 0; iy < exclusive_jets.size(); ++iy) {
 93            FourMomentum pjet(momentum(exclusive_jets[iy]));
 94            double cost = pjet.p3().unit().dot(pgamma.p3().unit());
 95            double ygamma = 2.*pjet.E()*pgamma.E()/sqr(evis)*(1.-cost);
 96            if (ygamma < ycut) {
 97              accept = false;
 98              break;
 99            }
100          }
101          if (!accept) continue;
102          /// @todo Really want to use a "bar graph" here (i.e. ignore bin width)
103          _nPhotonJade->fill(ycut, _nPhotonJade->bin(j).xWidth());
104          size_t njet = min(size_t(4), exclusive_jets.size()) - 1;
105          if (j < _nPhotonJetJade[njet]->numBins()) {
106            _nPhotonJetJade[njet]->fillBin(j, _nPhotonJetJade[njet]->bin(j).xWidth());
107          }
108        }
109        // Add this photon back in for the next iteration of the loop
110        if (ix+1 != photons.size()) {
111          input_particles[nonPhotons.size()+ix] = fastjet::PseudoJet(pgamma.px(), pgamma.py(), pgamma.pz(), pgamma.E());
112        }
113      }
114    }
115
116
117    void init() {
118      // Projections
119      declare(FinalState(), "FS");
120
121      // Book datasets
122      book(_nPhotonJade   ,1, 1, 1);
123      book(_nPhotonDurham ,2, 1, 1);
124      for (size_t ix = 0; ix < 4; ++ix) {
125        book(_nPhotonJetJade  [ix] ,3 , 1, 1+ix);
126        book(_nPhotonJetDurham[ix] ,4 , 1, 1+ix);
127      }
128    }
129
130
131    /// Finalize
132    void finalize() {
133      const double fact = 1000/sumOfWeights();
134      scale(_nPhotonJade, fact);
135      scale(_nPhotonDurham, fact);
136      for (size_t ix = 0; ix < 4; ++ix) {
137        scale(_nPhotonJetJade  [ix],fact);
138        scale(_nPhotonJetDurham[ix],fact);
139      }
140    }
141
142    /// @}
143
144
145  private:
146
147    Histo1DPtr _nPhotonJade;
148    Histo1DPtr _nPhotonDurham;
149    Histo1DPtr _nPhotonJetJade  [4];
150    Histo1DPtr _nPhotonJetDurham[4];
151
152  };
153
154
155
156  RIVET_DECLARE_ALIASED_PLUGIN(OPAL_1993_S2692198, OPAL_1993_I343181);
157
158}