Rivet analyses

// -*- 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_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);

}