rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

MC_HFJETS

Monte Carlo validation analysis to study heavy flavour production
Experiment: ()
Status: VALIDATED
Authors:
  • Andy Buckley
No references listed
Beams: * *
Beam energies: ANY
Run details:
  • Run any events which will produce jets above 20 GeV. Of most interest for processes where c and b hadrons can be produced (either hard or soft) of course!

Plots to study fragmentation of heavy flavour hadrons in jets.

Source code: MC_HFJETS.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/UnstableParticles.hh"
  6#include "Rivet/Projections/PrimaryHadrons.hh"
  7#include "Rivet/Projections/HeavyHadrons.hh"
  8
  9namespace Rivet {
 10
 11
 12  class MC_HFJETS : public Analysis {
 13  public:
 14
 15    // Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(MC_HFJETS);
 17
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21
 22      FastJets fj(FinalState(Cuts::abseta < 5), FastJets::ANTIKT, 0.6);
 23      fj.useInvisibles();
 24      declare(fj, "Jets");
 25      declare(HeavyHadrons(Cuts::abseta < 5 && Cuts::pT > 500*MeV), "BCHadrons");
 26
 27      book(_h_ptCJetLead ,"ptCJetLead", linspace(5, 0, 20, false) + logspace(25, 20, 200));
 28      book(_h_ptCHadrLead ,"ptCHadrLead", linspace(5, 0, 10, false) + logspace(25, 10, 200));
 29      book(_h_ptFracC ,"ptfracC", 50, 0, 1.5);
 30      book(_h_eFracC ,"efracC", 50, 0, 1.5);
 31
 32      book(_h_ptBJetLead ,"ptBJetLead", linspace(5, 0, 20, false) + logspace(25, 20, 200));
 33      book(_h_ptBHadrLead ,"ptBHadrLead", linspace(5, 0, 10, false) + logspace(25, 10, 200));
 34      book(_h_ptFracB ,"ptfracB", 50, 0, 1.5);
 35      book(_h_eFracB ,"efracB", 50, 0, 1.5);
 36    }
 37
 38
 39    /// Perform the per-event analysis
 40    void analyze(const Event& event) {
 41
 42      // Get jets and heavy hadrons
 43      const Jets& jets = apply<JetAlg>(event, "Jets").jetsByPt();
 44      const Particles bhadrons = sortByPt(apply<HeavyHadrons>(event, "BCHadrons").bHadrons());
 45      const Particles chadrons = sortByPt(apply<HeavyHadrons>(event, "BCHadrons").cHadrons());
 46      MSG_DEBUG("# b hadrons = " << bhadrons.size() << ", # c hadrons = " << chadrons.size());
 47
 48      // Loop over jets and use ghost-tag info
 49      for (const Jet& j : jets) {
 50        bool gotLeadingB = false, gotLeadingC = false;
 51        // b-tag testing
 52        if (!gotLeadingB && j.bTagged(Cuts::pT > 500*MeV)) {
 53          gotLeadingB = true;
 54          Particle bhad = sortByPt(j.bTags(Cuts::pT > 500*MeV))[0];
 55          _h_ptBJetLead->fill(j.pT()/GeV);
 56          _h_ptBHadrLead->fill(bhad.pT()/GeV);
 57          _h_ptFracB->fill(bhad.pT() / j.pT());
 58          _h_eFracB->fill(bhad.E() / j.E());
 59          continue;
 60        }
 61        // c-tag testing
 62        if (!gotLeadingC && j.cTagged(Cuts::pT > 500*MeV) && !j.bTagged(Cuts::pT > 500*MeV)) {
 63          gotLeadingC = true;
 64          Particle chad = sortByPt(j.cTags(Cuts::pT > 500*MeV))[0];
 65          _h_ptCJetLead->fill(j.pT()/GeV);
 66          _h_ptCHadrLead->fill(chad.pT()/GeV);
 67          _h_ptFracC->fill(chad.pT() / j.pT());
 68          _h_eFracC->fill(chad.E() / j.E());
 69        }
 70        // Escape early if we've found both the leading b and c jets
 71        if (gotLeadingB && gotLeadingC) break;
 72      }
 73
 74      // // Tag the leading b and c jets with a deltaR < 0.3 match
 75      // // b-tagged jet are excluded from also being considered as c-tagged
 76      // MSG_DEBUG("Getting b/c-tags");
 77      // const double MAX_DR = 0.3;
 78      // bool gotLeadingB = false, gotLeadingC = false;
 79      // for (const Jet& j : jets) {
 80      //   if (!gotLeadingB) {
 81      //     FourMomentum leadBJet, leadBHadr;
 82      //     double dRmin = MAX_DR;
 83      //     for (const Particle& b : bhadrons) {
 84      //       const double dRcand = min(dRmin, deltaR(j, b));
 85      //       if (dRcand < dRmin) {
 86      //         dRmin = dRcand;
 87      //         leadBJet = j.momentum();
 88      //         leadBHadr = b.momentum();
 89      //         MSG_DEBUG("New closest b-hadron jet tag candidate: dR = " << dRmin
 90      //                   << " for jet pT = " << j.pT()/GeV << " GeV, "
 91      //                   << " b hadron pT = " << b.pT()/GeV << " GeV, PID = " << b.pid());
 92      //       }
 93      //     }
 94      //     if (dRmin < MAX_DR) {
 95      //       // A jet has been tagged, so fill the histos and break the loop
 96      //       _h_ptBJetLead->fill(leadBJet.pT()/GeV, weight);
 97      //       _h_ptBHadrLead->fill(leadBHadr.pT()/GeV, weight);
 98      //       _h_ptFracB->fill(leadBHadr.pT() / leadBJet.pT(), weight);
 99      //       _h_eFracB->fill(leadBHadr.E() / leadBJet.E(), weight);
100      //       gotLeadingB = true;
101      //       continue; // escape this loop iteration so the same jet isn't c-tagged
102      //     }
103      //   }
104      //   if (!gotLeadingC) {
105      //     FourMomentum leadCJet, leadCHadr;
106      //     double dRmin = MAX_DR;
107      //     for (const Particle& c : chadrons) {
108      //       const double dRcand = min(dRmin, deltaR(j, c));
109      //       if (dRcand < dRmin) {
110      //         dRmin = dRcand;
111      //         leadCJet = j.momentum();
112      //         leadCHadr = c.momentum();
113      //         MSG_DEBUG("New closest c-hadron jet tag candidate: dR = " << dRmin
114      //                   << " for jet pT = " << j.pT()/GeV << " GeV, "
115      //                   << " c hadron pT = " << c.pT()/GeV << " GeV, PID = " << c.pid());
116      //       }
117      //     }
118      //     if (dRmin < MAX_DR) {
119      //       // A jet has been tagged, so fill the histos and break the loop
120      //       _h_ptCJetLead->fill(leadCJet.pT()/GeV, weight);
121      //       _h_ptCHadrLead->fill(leadCHadr.pT()/GeV, weight);
122      //       _h_ptFracC->fill(leadCHadr.pT() / leadCJet.pT(), weight);
123      //       _h_eFracC->fill(leadCHadr.E() / leadCJet.E(), weight);
124      //       gotLeadingB = true;
125      //     }
126      //   }
127      //   // If we've found both a leading b and a leading c jet, break the loop over jets
128      //   if (gotLeadingB && gotLeadingC) break;
129      // }
130
131    }
132
133
134    /// Normalise histograms etc., after the run
135    void finalize() {
136      normalize({_h_ptCJetLead, _h_ptCHadrLead, _h_ptBJetLead, _h_ptBHadrLead,
137            _h_ptFracC, _h_eFracC, _h_ptFracB, _h_eFracB});
138    }
139
140
141    /// @name Histograms
142    ///@{
143    Histo1DPtr _h_ptCJetLead, _h_ptCHadrLead, _h_ptFracC, _h_eFracC;
144    Histo1DPtr _h_ptBJetLead, _h_ptBHadrLead, _h_ptFracB, _h_eFracB;
145    ///@}
146
147
148  };
149
150
151
152  // The hook for the plugin system
153  RIVET_DECLARE_PLUGIN(MC_HFJETS);
154
155}