Rivet analyses

Measurement of b-jet mass with and without grooming

Experiment: LHCB (LHC)

Inspire ID: 2922449

Status: VALIDATED

Authors: - Ezra D. Lesser

References: - Phys. Lett. B 869 (2025) 139854 - DOI:10.1016/J.PhysLetB.2025.139854 - arXiv: 2505.11955 - Expt page: LHCb-PAPER-2025-009 - CERN-EP-2025-097

Beams: p+ p+

Beam energies: (6500.0, 6500.0)GeV

Run details: - Hard QCD events from proton-proton interactions at 13 TeV centre-of-mass energy in which B± meson decay is suppressed since they are fully reconstructed from the J/ψK± decay channel before the b-jet mass measurement is performed.

The LHCb Collaboration presents measurements of the b-jet mass both with and without grooming using proton-proton data collected at a center-of-mass energy of $\sqrt{s} = 13$ TeV. B± mesons are fully reconstructed from the J/ψK± decay channel. Jets with transverse momentum pT, jet > 10 GeV are reconstructed within the rapidity range 2.5 < yjet < 4.0 using the anti-kT algorithm with resolution paparmeter R = 0.5. Groomed jets are produced using Soft Drop with zcut = 0.1 and β = 0. Both groomed and ungroomed jets are further evaluated using the Winner-Take-All flavor algorithm, and figures are produced both with and without WTA flavor tagging.

Source code:LHCB_2025_I2922449.cc

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Analyses/LHCbCommon.hh"
#include "Rivet/Projections/FastJets.hh"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/JetDefinition.hh"
#include "fastjet/contrib/LundGenerator.hh" // For Soft Drop grooming

namespace Rivet {


  /// @brief Measurements of the $b$-jet mass both with and without grooming in proton-proton collision events at 13 TeV
  class LHCB_2025_I2922449 : public Analysis {
    /////////////////////////////////////////////////////////////////////////////////////////////////
    // Analysis parameters

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

    // Jet reconstruction
    const double JET_R = 0.5;             // Jet resolution parameter for anti-kT
    const double RAPIDITY_MIN_JETS = 2.5; // Min rapidity of constructed B-jets
    const double RAPIDITY_MAX_JETS = 4.0; // Max "
    const double PT_MIN_JETS = 10.;       // Min jet transverse momentum
    const double PT_MAX_JETS = 100.;      // Max "

    // HF reconstruction
    const vector<PdgId> CH_B_MESON_PID = {PID::BPLUS, PID::BMINUS};

    // WTA reclustering
    const double SMALL_NUMBER = 1e-4;
    const double& MAX_WTA_DISTANCE = SMALL_NUMBER;
    const fastjet::JetDefinition WTA_JET_DEF = fastjet::JetDefinition(fastjet::cambridge_algorithm,
                                                                      fastjet::JetDefinition::max_allowable_R,
                                                                      fastjet::WTA_pt_scheme);

    // Soft Drop grooming
    const double SD_ZCUT = 0.1;
    const double SD_BETA = 0.;
    const fastjet::contrib::LundGenerator B_JET_LG = fastjet::contrib::LundGenerator(
        fastjet::cambridge_algorithm);

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

  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(LHCB_2025_I2922449);


    /// @name Analysis methods
    /// @{
    /////////////////////////////////////////////////////////////////////////////////////////////////
    /// Book histograms and initialise projections before the run
    void init() {

      // Initialise and register projections
      // HVRecoFinalState projection selects also targeted charged B mesons
      Cut particle_selector = Cuts::etaIn(ETA_MIN_PARTICLES, ETA_MAX_PARTICLES);
      const LHCb::HVRecoFinalState hvfs(particle_selector, CH_B_MESON_PID, particle_selector);
      declare(hvfs, "hvfs");

      // The final-state particles declared above are clustered using FastJet with
      // the anti-kT algorithm and a jet-radius parameter 0.5
      // muons are included and neutrinos are excluded from the clustering
      FastJets jetfs(hvfs, JetAlg::ANTIKT, JET_R, JetMuons::ALL, JetInvisibles::NONE);
      declare(jetfs, "jets");

      // Book histograms

      book(_hh, PT_BIN_EDGES);        // HEPData Table 12-17
      book(_hh_gr, PT_BIN_EDGES);     // HEPData Table 0-5
      book(_hh_WTA, PT_BIN_EDGES);    // HEPData Table 36-41
      book(_hh_gr_WTA, PT_BIN_EDGES); // HEPData Table 18-23
      for (size_t i = 0; i < _hh->numBins(); ++i) {
        book(_hh->bin(i + 1), 13 + i, 1, 1);
        book(_hh_gr->bin(i + 1), 1 + i, 1, 1);
        book(_hh_WTA->bin(i + 1), 37 + i, 1, 1);
        book(_hh_gr_WTA->bin(i + 1), 19 + i, 1, 1);
      }

      return;
    }


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

      // retrieve reconstructed charged B mesons from the HVRecoFinalState projection
      Particles B_mesons = apply<LHCb::HVRecoFinalState>(event, "hvfs").hfHadrons();

      if (B_mesons.empty()) {
        MSG_DEBUG("No charged B mesons found. Veto event...");
        vetoEvent;
      }

      // listing B mesons in DEBUG mode
      if (getLog().isActive(Log::Level::DEBUG)) {
        MSG_DEBUG("B meson(s) found with ");
        bool add_comma = false;
        for (const Particle& B_meson : B_mesons) {
          if (add_comma) MSG_DEBUG(", ");
          MSG_DEBUG("rapidity=" << B_meson.rap() << " and pT=" << B_meson.pt());
          add_comma = true;
        }
        MSG_DEBUG(std::endl);
      }

      // Retrieve clustered jets, sorted by pT, with applied rapidity and pT cuts
      Cut jet_selector = Cuts::rapIn(RAPIDITY_MIN_JETS, RAPIDITY_MAX_JETS)
          & Cuts::ptIn(PT_MIN_JETS, PT_MAX_JETS);
      Jets jets = apply<FastJets>(event, "jets").jetsByPt(jet_selector);
      if (jets.empty()) {
        MSG_DEBUG("Could not reconstruct jets in LHCb detector fiducial phase space" << std::endl);
        vetoEvent;
      }

      Jets B_jets;
      for (const Jet& jet : jets) {
        // Jet contains undecayed B+ or B-
        if (jet.containsParticleId(CH_B_MESON_PID)) B_jets.push_back(jet);
      }

      MSG_DEBUG("#" << B_jets.size() << " B jets found" << endl);
      if (B_jets.empty()) {
        vetoEvent;
      }
      ////////////////////////////////////////////////////////////////////////////////////////

      // Loop over jets and apply operations & fill histograms
      Jets B_jets_WTA;
      for (const Jet& B_jet : B_jets) {
        // Find the WTA axis of the B jet using C/A algorithm
        fastjet::ClusterSequence B_jet_WTA_cs(B_jet.pseudojet().constituents(), WTA_JET_DEF);
        Jet B_jet_WTA = fastjet::sorted_by_pt(B_jet_WTA_cs.inclusive_jets())[0];
        assert(B_jet.pseudojet().constituents().size()
               == B_jet_WTA.pseudojet()
                      .constituents()
                      .size()); // Make sure that reclustering preserves all particles

        // Check if charged B meson is on WTA axis (arXiv:2205.01117)
        bool WTA_tagged = false;
        for (const Particle& B_meson : B_mesons) {
          double WTA_distance = B_jet_WTA.pseudojet().delta_R(B_meson);
          if (WTA_distance < MAX_WTA_DISTANCE) {
            WTA_tagged = true;
            break;
          }
        }

        // Groom the jets with Soft Drop
        Jet B_jet_gr;
        for (fastjet::contrib::LundDeclustering& ld : B_JET_LG.result(B_jet)) {
          if (ld.z() > SD_ZCUT * std::pow(ld.Delta() / JET_R, SD_BETA)) { // SD condition
            B_jet_gr = ld.pair();
            break;
          }
        }
        // Check that charged B meson survives
        bool gr_tagged = false;
        for (const Particle& B_meson : B_mesons) {
          if (gr_tagged || !B_jet_gr.pseudojet().has_constituents()) break;
          for (const fastjet::PseudoJet& particle : B_jet_gr.pseudojet().constituents()) {
            if ((particle.delta_R(B_meson) < SMALL_NUMBER)
                && fuzzyEquals(particle.pt(), B_meson.pt(), SMALL_NUMBER)) {
              // Should be the same particle (could always compare barcodes?!)
              gr_tagged = true;
              break;
            }
          }
        }

        if (gr_tagged) {
          MSG_DEBUG("** B jet survived grooming" << endl);
        }
        else {
          MSG_DEBUG("** B jet did NOT survive grooming" << endl);
        }

        // Fill histograms
        MSG_DEBUG("Filling histograms..." << endl);
        double jpt = B_jet.pT() / GeV;
        double mpt = B_jet.pseudojet().m() / jpt;
        _hh->fill(jpt, mpt);
        if (WTA_tagged) _hh_WTA->fill(jpt, mpt);
        if (gr_tagged) {
          // Note: Groomed observable uses ungroomed jet pT
          double grmpt = B_jet_gr.pseudojet().m() / jpt;
          _hh_gr->fill(jpt, grmpt);
          if (WTA_tagged) _hh_gr_WTA->fill(jpt, grmpt);
        }
      }

      return;
    }


    /////////////////////////////////////////////////////////////////////////////////////////////////
    /// Normalise histograms etc., after the run
    void finalize() {

      normalize(_hh_WTA);
      normalize(_hh_gr);
      normalize(_hh);
      normalize(_hh_gr_WTA);

      return;
    }

    /// @}

    /////////////////////////////////////////////////////////////////////////////////////////////////
    /// @name Histograms
    /// @{
    Histo1DGroupPtr _hh, _hh_gr, _hh_WTA, _hh_gr_WTA;
    /// @}
  };


  RIVET_DECLARE_PLUGIN(LHCB_2025_I2922449);

}