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