Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/UnstableParticles.hh"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/JetDefinition.hh"


namespace Rivet {


  /// @brief J/psi production in jets at 13 TeV
  class LHCB_2017_I1509507 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(LHCB_2017_I1509507);


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

    /// Book histograms and initialise projections before the run
    void init() {
      declare(UnstableParticles(), "UFS");
      declare(FinalState(), "FS");
      for (unsigned int ix = 0; ix < 2; ++ix) book(_h_frag[ix], 1, 1, ix + 1);
    }

    /// Perform the per-event analysis
    void analyze(const Event& event) {
      // first see if we have any J/psi in the region
      Particles Jpsi;
      for (const Particle& p :
           apply<UnstableParticles>(event, "UFS").particles(Cuts::pid == 443 && Cuts::etaIn(2., 4.5))) {
        if (p.children().size() != 2) continue;
        bool found = true;
        for (const Particle& child : p.children()) {
          double eta = child.eta();
          if (child.abspid() != 13 || eta < 2. || eta > 4.5 || child.perp() < 0.5 || child.p3().mod() < 5.) {
            found = false;
            break;
          }
        }
        if (found) Jpsi.push_back(p);
      }
      // no jpsi veto
      if (Jpsi.empty()) vetoEvent;
      // now get the final-state particles for the jets
      const Particles& fsParticles = apply<FinalState>(event, "FS").particles();
      vector<PseudoJet> particles;
      particles.reserve(fsParticles.size());
      // fs for fastjet omitting any J/psi decay products
      for (const Particle& p : fsParticles) {
        // skip muons and neutrinos (seems standard for LHCb)
        if (p.abspid() == 13 || p.abspid() == 12 || p.abspid() == 14 or p.abspid() == 16) continue;
        // skip anything coming from the decay of one of the jpsis
        Particle parent = p;
        while (!parent.parents().empty()) {
          if (parent.pid() == 443) break;
          parent = parent.parents()[0];
        }
        bool match = parent.pid() == 443;
        if (match) {
          match = false;
          for (const Particle& psi : Jpsi) {
            match = fuzzyEquals(parent.momentum(), psi.momentum());
            if (match) break;
          }
        }
        if (!match) {
          PseudoJet j = p.pseudojet();
          j.set_user_index(0);
          particles.push_back(j);
        }
      }
      // add the jpsis to the particles for fastjet
      for (const Particle& p : Jpsi) {
        PseudoJet j = p.pseudojet();
        j.set_user_index(p.fromBottom() + 1);
        particles.push_back(j);
      }
      JetDefinition jet_def(fastjet::antikt_algorithm, 0.5);
      fastjet::ClusterSequence clu = ClusterSequence(particles, jet_def);
      vector<PseudoJet> jets = clu.inclusive_jets();
      for (const PseudoJet& jet : jets) {
        // pt and eta cut
        if (jet.perp() < 20. || jet.eta() < 2. || jet.eta() > 4.) continue;
        // loop over constituents and find jpsi
        for (const PseudoJet& sub : jet.constituents()) {
          if (sub.user_index() == 0) continue;
          double z = sub.perp() / jet.perp();
          _h_frag[sub.user_index() - 1]->fill(z);
        }
      }
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      for (unsigned int ix = 0; ix < 2; ++ix) normalize(_h_frag[ix], 1., false);
    }

    /// @}


    /// @name Histograms
    /// @{
    Histo1DPtr _h_frag[2];
    /// @}
  };


  RIVET_DECLARE_PLUGIN(LHCB_2017_I1509507);

}