rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2021_I1913061

b-fragmentation at 13 TeV
Experiment: ATLAS (LHC)
Inspire ID: 1913061
Status: VALIDATED
Authors:
  • Javier Llorente
References: Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • pp interactions at 13000 GeV. Particles with c*tau>10mm are stable.

The fragmentation properties of jets containing $b$-hadrons are studied using charged $B$-mesons in 139 fb$^{-1}$ of $pp$ collisions at $\sqrt{s} = 13$ TeV, recorded with the ATLAS detector at the LHC during the period from 2015 to 2018. The $B$ mesons are reconstructed using the decay of $B^{\pm}$ into $J/\psi K^{\pm}$, with the $J/\psi$ decaying into a pair of muons. Jets are reconstructed using the anti-$k_t$ algorithm with radius parameter $R=0.4$. The measurement determines the longitudinal and transverse momentum profiles of the reconstructed $B$-hadrons with respect to the axes of the jets to which they are geometrically associated. These distributions are measured in intervals of the jet transverse momentum, ranging from 50 GeV to above 100 GeV. The results are corrected for detector effects and compared with several Monte Carlo predictions using different parton shower and hadronisation models. The results for the longitudinal and transverse profiles provide useful inputs to improve the description of heavy-flavour fragmentation in jets.

Source code: ATLAS_2021_I1913061.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/PromptFinalState.hh"
  5#include "Rivet/Projections/LeptonFinder.hh"
  6#include "Rivet/Projections/VetoedFinalState.hh"
  7#include "Rivet/Projections/UnstableParticles.hh"
  8#include "Rivet/Projections/FastJets.hh"
  9
 10namespace Rivet {
 11
 12
 13  /// @brief b-fragmentation at 13 TeV
 14  class ATLAS_2021_I1913061 : public Analysis {
 15
 16  public:
 17
 18    /// Default constructor
 19    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2021_I1913061);
 20
 21    /// @name Analysis methods
 22    /// @{
 23
 24    void init() {
 25
 26      //Jet building: anti-kT R = 0.4, including muons and neutrinos.
 27      FinalState photons(Cuts::abspid == PID::PHOTON);
 28
 29      PromptFinalState bare_mu(Cuts::abspid == PID::MUON, TauDecaysAs::PROMPT);
 30      LeptonFinder all_dressed_mu(bare_mu, photons, 0.1, Cuts::abseta < 2.5);
 31
 32      PromptFinalState bare_el(Cuts::abspid == PID::ELECTRON, TauDecaysAs::PROMPT);
 33      LeptonFinder all_dressed_el(bare_el, photons, 0.1, Cuts::abseta < 2.5);
 34
 35      VetoedFinalState vfs(FinalState(Cuts::abseta < 4.5));
 36      vfs.addVetoOnThisFinalState(all_dressed_el);
 37      vfs.addVetoOnThisFinalState(all_dressed_mu);
 38
 39      FastJets jets(vfs, JetAlg::ANTIKT, 0.4, JetMuons::ALL, JetInvisibles::DECAY);
 40      declare(jets, "JETS");
 41
 42      //Charged B mesons
 43      UnstableParticles bpm_fs(Cuts::abspid == 521);
 44      declare(bpm_fs, "BPM_FS");
 45
 46      // Book histograms
 47      book(_h["zFrag_pt01"], 1,1,1);
 48      book(_h["ptRel_pt01"], 2,1,1);
 49      book(_h["zFrag_pt02"], 3,1,1);
 50      book(_h["ptRel_pt02"], 4,1,1);
 51      book(_h["zFrag_pt03"], 5,1,1);
 52      book(_h["ptRel_pt03"], 6,1,1);
 53      book(_p["zFrag"], 7,1,1);
 54      book(_p["ptRel"], 8,1,1);
 55    }
 56
 57
 58    /// Perform the per-event analysis
 59    void analyze(const Event& event) {
 60
 61      const Jets& jets = apply<FastJets>(event, "JETS").jetsByPt(Cuts::pT > 7*GeV);
 62      const Particles& bpmFS = apply<UnstableParticles>(event, "BPM_FS").particlesByPt();
 63
 64      //Preselect B mesons in J/psi K decay channel
 65      Particles goodHadrons;
 66
 67      for (const Particle& p : bpmFS) {
 68
 69        bool passKaon = 0;
 70        bool passMuon1 = 0;
 71        bool passMuon2 = 0;
 72
 73        const Particles& bDecays = p.children();
 74
 75        for (const Particle& s : bDecays){
 76
 77          if (s.abspid() == PID::KPLUS && s.pt() > 4*GeV) passKaon = 1;
 78          if (s.abspid() == PID::JPSI) {
 79
 80            const Particles& jpsiDecays = s.children();
 81
 82            for (const Particle& m : jpsiDecays){
 83
 84              if (m.pid() == PID::ANTIMUON && m.pt() > 6*GeV) passMuon1 = 1;
 85              if (m.pid() == PID::MUON && m.pt() > 6*GeV) passMuon2 = 1;
 86            }
 87          }
 88        }
 89
 90        if (passKaon && passMuon1 && passMuon2) goodHadrons += p;
 91      }
 92
 93      //Preselect jets passing kinematic cuts
 94      Jets goodJets;
 95      for (size_t i = 0; i < jets.size(); ++i) {
 96        if (jets[i].pT() <= 20*GeV) continue;
 97        if (jets[i].abseta() >= 2.1) continue;
 98        bool overlaps = false;
 99        for (size_t j = i + 1; j < jets.size(); ++j) {
100          if (jets[j].pT() > 20*GeV && deltaR(jets[i], jets[j]) < 0.8) {
101            overlaps = true; break;
102          }
103        }
104        if (!overlaps)  goodJets += jets[i];
105      }
106
107      //Associate the jets to the B hadrons.
108      for (const Jet& thisJet : goodJets) {
109
110        Particles jetHads;
111        for (const Particle& thisHad : goodHadrons) {
112          if (deltaR(thisJet, thisHad) < 0.4) jetHads += thisHad;
113        }
114
115        if (jetHads.size() == 1){
116          const Particle& jetHadron = jetHads[0];
117
118          //Out-of-cone correction for B-decay products
119          Vector3 jetVector = momentum3(thisJet); Vector3 jetVector0 = jetVector;
120          Vector3 hadronVector = momentum3(jetHadron);
121          Vector3 kaonVector; Vector3 muonVector1; Vector3 muonVector2;
122
123          const Particles& bDecays = jetHadron.children();
124          for (const Particle& s : bDecays){
125
126            if (s.abspid() == PID::KPLUS) kaonVector = momentum3(s);
127            if (s.abspid() == PID::JPSI) {
128
129              const Particles& jpsiDecays = s.children();
130
131              for (const Particle& m : jpsiDecays){
132
133                if (m.pid() == PID::ANTIMUON && m.pt() > 6*GeV) muonVector1 = momentum3(m);
134                if (m.pid() == PID::MUON && m.pt() > 6*GeV) muonVector2 = momentum3(m);
135              }
136            }
137          }
138
139          if (deltaR(muonVector1, jetVector0) > 0.4) jetVector += muonVector1;
140          if (deltaR(muonVector2, jetVector0) > 0.4) jetVector += muonVector2;
141          if (deltaR(kaonVector, jetVector0) > 0.4) jetVector += kaonVector;
142
143          if (jetVector.perp() <= 30.0*GeV) vetoEvent;
144          if (jetVector.abseta() >= 2.1)    vetoEvent;
145          if (deltaR(jetVector, hadronVector) >= 0.4)  vetoEvent;
146
147          //Longitudinal and transverse profiles
148          double zFrag = hadronVector.dot(jetVector)/jetVector.mod2();
149          double ptRel = (hadronVector.cross(jetVector)).mod()/jetVector.mod();
150
151          if (jetVector.perp() >= 50*GeV && jetVector.perp() < 70*GeV) {
152            if (zFrag <= 0.23) zFrag = 0.24;
153            if (zFrag >= 1.00) zFrag = 0.99;
154            if (ptRel >= 8.00) ptRel = 7.90;
155
156            _h["zFrag_pt01"]->fill(zFrag);
157            _h["ptRel_pt01"]->fill(ptRel);
158          }
159
160          if (jetVector.perp() >= 70*GeV && jetVector.perp() < 100*GeV) {
161            if (zFrag <= 0.23) zFrag = 0.24;
162            if (zFrag >= 1.00) zFrag = 0.99;
163            if (ptRel >= 10.0) ptRel = 9.90;
164
165            _h["zFrag_pt02"]->fill(zFrag);
166            _h["ptRel_pt02"]->fill(ptRel);
167          }
168
169          if (jetVector.perp() >= 100*GeV) {
170            if (zFrag <= 0.16) zFrag = 0.17;
171            if (zFrag >= 1.00) zFrag = 0.99;
172            if (ptRel >= 14.0) ptRel = 13.9;
173
174            _h["zFrag_pt03"]->fill(zFrag);
175            _h["ptRel_pt03"]->fill(ptRel);
176          }
177
178          double ptFill = jetVector.perp()/GeV;
179          if (ptFill >= 150.)  ptFill = 125.;
180
181          // For pTrel, we need to weight the mean by 1/binWidth
182          // (what is being plotted is the mean of the histogram, not the variable!)
183          // This doesn't apply to z since all bins have the same width.
184          double binWidth = 0.0;
185          double transBins[12] = { 0.0, 0.5, 1.0, 1.5, 2.2, 3.0, 4.0, 5.0, 6.5, 8.0, 10.0, 14.0 };
186
187          for (size_t iBin = 0; iBin < 11; ++iBin) {
188            if (ptRel >= transBins[iBin] && ptRel < transBins[iBin+1])  binWidth = transBins[iBin+1]-transBins[iBin];
189          }
190
191          _p["zFrag"]->fill(ptFill, zFrag);
192          _p["ptRel"]->fill(ptFill, ptRel, 1./binWidth);
193
194        }
195      }
196    }
197
198
199    void finalize() {
200      normalize(_h);
201    }
202
203
204  private:
205
206    //Histograms
207    map<string, Histo1DPtr> _h;
208    map<string, Profile1DPtr> _p;
209
210  };
211
212
213  RIVET_DECLARE_PLUGIN(ATLAS_2021_I1913061);
214
215}