Rivet analyses

title: LHCB_2019_I1730448

Measurement of charged hadron production in Z-tagged jets at 8 TeV

Experiment: LHCB (LHC)

Inspire ID: 1730448

Status: VALIDATED

Authors: - Alice Xu

References: - Phys. Rev. Lett. 123 (2019) 232001 - DOI:10.1103/PhysRevLett.123.232001 - arXiv: 1904.08878 - Expt page: LHCb-PAPER-2019-012

Beams: p+ p+

Beam energies: (4000.0, 4000.0)GeV

Run details: - Proton-proton interactions at 8 TeV centre-of-mass energy.

The LHCb Collaboration presents measurements of the longitudinal, transverse, and radial distributions of charged hadrons in $Z$-jets using proton-proton data collected at a center-of-mass energy of $\sqrt{s} = 8$TeV. $Z$ bosons are fully reconstructed from the $\mu^- \mu^+$ decay channel. Jets with transverse momentum $p_\mathrm{T,jet} > 20$ GeV are reconstructed within the rapidity range $2.5 < y_\mathrm{jet} < 4.0$ using the anti-$k_\mathrm{T}$ algorithm with resolution paparmeter $R = 0.5$.

Source code:LHCB_2019_I1730448.cc

```c++ // -- 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);

} ```