Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/DileptonFinder.hh"
#include "Rivet/Projections/DirectFinalState.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/MissingMomentum.hh"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/JetDefinition.hh"
#include "fastjet/Selector.hh" // For creating Filter for jet Projector

namespace Rivet {

  /// @brief Measurement of charged hadron production in Z-tagged jets at 8 TeV
  class LHCB_2019_I1730448 : public Analysis {

    /////////////////////////////////////////////////////////////////////////////////////////////////
    // Analysis parameters

    // Particle reconstruction
    const double ETA_MIN_PARTICLES = 2.;  // Min pseudorapidity of final-state particles
    const double ETA_MAX_PARTICLES = 4.5; // Max "

    // Jet reconstruction
    const double JET_R = 0.5;        // Jet resolution parameter for anti-kT
    const double ETA_MIN_JETS = 2.5; // Min rapidity of constructed Z-jets
    const double ETA_MAX_JETS = 4.;  // Max "
    const double PT_MIN_JETS = 20.;  // Min jet transverse momentum
    const double PT_MAX_JETS = 100.; // Max "

    // Z Boson
    const double MASS_MIN_ZBOSONS = 60.;  // Min invariant mass of Z bosons
    const double MASS_MAX_ZBOSONS = 120.; // Max "
    const double PT_MIN_ZBOSONS = 0.;     // Min Z boson transverse momentum
    const double PT_MAX_ZBOSONS = 100.;   // Max "

    // Charged Hadron
    const double P_MIN_HADRONS = 4.;    // Min momentum of charged hadrons
    const double P_MAX_HADRONS = 1000.; // Max "
    const double PT_MIN_HADRONS = 0.25; // Min transverse momentum of charged hadrons


  public:

    // Attaches constituent PID to a PseudoJet
    class JetInfo : public fastjet::PseudoJet::UserInfoBase {
    public:

      JetInfo(const int& _pid)
          : pid(_pid) { };
      int get_pid() const {
        return this->pid;
      }

    private:

      int pid = -999;
    };

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(LHCB_2019_I1730448);


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

    /////////////////////////////////////////////////////////////////////////////////////////////////
    /// Book histograms and initialise projections before the run
    void init() {
      // Initialise and register projections

      // The basic final-state projection:
      // all final-state particles within the given eta acceptance
      Cut particle_selector = Cuts::etaIn(ETA_MIN_PARTICLES, ETA_MAX_PARTICLES);
      const FinalState fs(particle_selector);
      declare(fs, "fs_particles");

      Cut Z_boson_selector = Cuts::ptIn(PT_MIN_ZBOSONS, PT_MAX_ZBOSONS)
          & Cuts::massIn(MASS_MIN_ZBOSONS, MASS_MAX_ZBOSONS) & particle_selector;

      // Find *decayed* Z bosons from dimuon children
      // pT > 20 GeV according to JHEP 08 (2015) 039;  DOI:10.1007/JHEP08%282015%29039 referenced indirectly in paper
      DileptonFinder zfinder(fs, 91.2 * GeV, -1., (Cuts::abspid == PID::MUON) && Cuts::pT > 20. * GeV,
                             Z_boson_selector);
      declare(zfinder, "ZFinder");

      // Jet pT bin edges
      const vector<double> PT_BIN_EDGES = {20., 30., 50., 100.};

      // Book histograms
      book(_h_ufld_z, PT_BIN_EDGES);
      book(_h_ufld_jt, PT_BIN_EDGES);
      book(_h_ufld_r, PT_BIN_EDGES);
      for (size_t i = 0; i < _h_ufld_z->numBins(); i++) {
        book(_h_ufld_z->bin(i + 1), 1, 1, i + 1);  // HEPData Table 1
        book(_h_ufld_jt->bin(i + 1), 2, 1, i + 1); // HEPData Table 2
        book(_h_ufld_r->bin(i + 1), 3, 1, i + 1);  // HEPData Table 3
      };
    }


    /////////////////////////////////////////////////////////////////////////////////////////////////
    /// Perform the per-event analysis
    void analyze(const Event& event) {

      const FinalState& fs = apply<FinalState>(event, "fs_particles");

      if (getLog().isActive(Rivet::Log::DEBUG)) {
        int nmuons = 0;
        for (const Particle& p : fs.particles()) {
          if (abs(p.pid()) == 13) nmuons++;
        }
        MSG_DEBUG("fs size: " << fs.particles().size() << " | muons: " << nmuons << std::endl);
      }

      // The final-state particles declared above are clustered using FastJet with
      // the anti-kT algorithm and a jet-radius parameter 0.5
      Rivet::PseudoJets trimmed_particles;
      for (const Particle& fs_particle : fs.particles()) {
        fastjet::PseudoJet InfoJet = fs_particle.pseudojet();
        InfoJet.set_user_info_shared_ptr(
            fastjet::SharedPtr<fastjet::PseudoJet::UserInfoBase>(new JetInfo(fs_particle.pid())));
        trimmed_particles.push_back(InfoJet);
      }
      fastjet::JetDefinition jet_def(fastjet::antikt_algorithm, JET_R);
      fastjet::ClusterSequence cs(trimmed_particles, jet_def);
      Rivet::PseudoJets jets = sorted_by_pt(cs.inclusive_jets());

      // Retrieve *decayed* Z bosons and muon children
      const DileptonFinder& zfinder = apply<DileptonFinder>(event, "ZFinder");
      if (zfinder.bosons().empty()) vetoEvent;

      const Particles& Z_bosons = zfinder.bosons();
      const Particles& muons = zfinder.constituents(); // mu+ and mu-

      // Retrieve clustered jets, sorted by pT, with applied rapidity and pT cuts
      fastjet::Selector jet_selector = fastjet::SelectorPtRange(PT_MIN_JETS, PT_MAX_JETS)
          && fastjet::SelectorEtaRange(ETA_MIN_JETS, ETA_MAX_JETS);
      jets = jet_selector(jets);
      if (jets.empty()) {
        vetoEvent;
      }

      // Create a vector of jets that are back to back with a Z boson, with no muons in jet cone
      Jets Z_jets;
      for (const Rivet::Jet jet : jets) {
        for (const Particle& Z_boson : Z_bosons) {
          if ((Z_boson.pT() < PT_MIN_ZBOSONS) || (Z_boson.pT() > PT_MAX_ZBOSONS))
            continue; // Check Z pT - already in cuts!!!
          bool muon_in_jet = false;
          for (const Particle& muon : muons) {
            if (deltaR(muon, jet) < JET_R) {
              muon_in_jet = true;
              break;
            }
          }
          if (muon_in_jet) continue;
          if (deltaPhi(Z_boson, jet) > (7 * M_PI) / 8) { // Azimuthal cut
            Z_jets.push_back(jet);
            // Recording the number of jets in each bin, for scaling purposes
            if ((jet.pT() >= 20.) && (jet.pT() < 30.)) {
              num_jets_20_30 += 1.;
            }
            else if ((jet.pT() >= 30.) && (jet.pT() < 50.)) {
              num_jets_30_50 += 1.;
            }
            else if ((jet.pT() >= 50.) && (jet.pT() < 100.)) {
              num_jets_50_100 += 1.;
            }
            MSG_DEBUG("deltaphi: " << deltaPhi(Z_boson, jet) << std::endl
                                   << "jet eta: " << jet.eta() << std::endl
                                   << "jet pt: " << jet.pT() << std::endl
                                   << "N_jet <20-30>: " << num_jets_20_30 << std::endl
                                   << "N_jet <30-50>: " << num_jets_30_50 << std::endl
                                   << "N_jet <50-100>: " << num_jets_50_100 << std::endl);
          }
        }
      }
      if (Z_jets.empty()) {
        vetoEvent;
      }

      MSG_DEBUG(Z_jets.size() << " Z jets found." << std::endl);

      ////////////////////////////////////////////////////////////////////////////////////////

      // Loop over jets and apply operations & fill histograms
      for (const Rivet::Jet& Z_jet : Z_jets) {
        for (const fastjet::PseudoJet& constituent : Z_jet.pseudojet().constituents()) {
          const JetInfo& myinfo = constituent.user_info<JetInfo>();
          Rivet::FourMomentum cmom(constituent.E(), constituent.px(), constituent.py(), constituent.pz());
          // Fill histogram
          if (((abs(myinfo.get_pid()) == PID::PIPLUS) || (abs(myinfo.get_pid()) == PID::KPLUS)
               || (abs(myinfo.get_pid()) == PID::PROTON))
              && (constituent.modp() > P_MIN_HADRONS) && (constituent.modp() < P_MAX_HADRONS)
              && (constituent.pt() > PT_MIN_HADRONS)
              && (deltaR(cmom, Z_jet, Rivet::RapScheme::YRAP) < JET_R)) {
            MSG_DEBUG("Filling histograms for " << myinfo.get_pid() << " ..." << std::endl);
            double num_z = Z_jet.p3().dot(cmom.p3());
            double den_z = Z_jet.p2();
            double num_jt = (Z_jet.p3().cross(cmom.p3())).mod();
            double den_jt = Z_jet.p();
            double r = deltaR(cmom, Z_jet, Rivet::RapScheme::YRAP);
            double zj_pt = Z_jet.pT();
            _h_ufld_z->fill(zj_pt, num_z / den_z);    // Fill histogram with longitudinal momentum
            _h_ufld_jt->fill(zj_pt, num_jt / den_jt); // Fill histogram with transverse momentum
            _h_ufld_r->fill(zj_pt, r);                // Fill histogram with radial profile distribution
          };
        }
      }
    }


    /////////////////////////////////////////////////////////////////////////////////////////////////
    /// Scale histograms etc., after the run
    void finalize() {
      const vector<double> scaleFactors = {1., (fuzzyEquals(0., num_jets_20_30) ? 1. : 1. / num_jets_20_30),
                                           (fuzzyEquals(0., num_jets_30_50) ? 1. : 1. / num_jets_30_50),
                                           (fuzzyEquals(0., num_jets_50_100) ? 1. : 1. / num_jets_50_100),
                                           1.};
      scale(_h_ufld_z, scaleFactors);
      scale(_h_ufld_jt, scaleFactors);
      scale(_h_ufld_r, scaleFactors);
    }

    /// @}

  private:

    /// @name Histograms
    /// @{
    Histo1DGroupPtr _h_ufld_z;
    Histo1DGroupPtr _h_ufld_jt;
    Histo1DGroupPtr _h_ufld_r;
    /// @}

    // Jet reconstruction counters
    double num_jets_20_30 = 0.;  // Total number of Z-tagged jets from the 20<pT<30 GeV bin
    double num_jets_30_50 = 0.;  // " 30<pT<50 GeV bin
    double num_jets_50_100 = 0.; // " 50<pT<100 GeV bin
  };


  RIVET_DECLARE_PLUGIN(LHCB_2019_I1730448);

}