rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

MC_HFDECAYS

Monte Carlo validation observables for heavy-flavour decays
Experiment: ()
Status: VALIDATED
Authors:
  • Francesco Curcio
  • Christian Gutschow
No references listed
Beams: * *
Beam energies: ANY
Run details:
  • Anything with heavy-flavour jets

Jets with $p_\perp>25$ GeV are constructed with an anti-$k_\perp$ jet finder with $R=0.4$ and projected onto many different observables sensitive to the decays of $b$- and $c$-hadrons in the heavy-flavour-tagged jet.

Source code: MC_HFDECAYS.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/FastJets.hh"
  5#include "Rivet/Projections/HeavyHadrons.hh"
  6#include "Rivet/Tools/BinnedHistogram.hh"
  7
  8namespace Rivet {
  9
 10
 11    class MC_HFDECAYS : public Analysis {
 12    public:
 13
 14      /// Constructor
 15      RIVET_DEFAULT_ANALYSIS_CTOR(MC_HFDECAYS);
 16
 17      const string whoDis(const int pid) const {
 18        switch (pid) {
 19          case PID::B0:           return "B0";
 20          case PID::BPLUS:        return "BPLUS";
 21          case PID::B0S:          return "B0S";
 22          case PID::LAMBDAB:      return "LAMBDAB";
 23          case PID::D0:           return "D0";
 24          case PID::DPLUS:        return "DPLUS";
 25          case PID::DSPLUS:       return "DSPLUS";
 26          case PID::LAMBDACPLUS:  return "LAMBDACPLUS";
 27          default:                return "";
 28        }
 29      }
 30
 31      void addBinned(const string &name, const YODA::HistoBin1D &ptbin, const int nbins, const double lo, const double hi) {
 32        string suff = "_" + to_string(int(ptbin.xMin())) + "_" + to_string(int(ptbin.xMax()));
 33        { Histo1DPtr tmp; _b[name].add(ptbin.xMin(), ptbin.xMax(), book(tmp, name+suff, nbins, lo, hi)); }
 34      }
 35
 36      /// Book histograms and initialise projections before the run
 37      void init() {
 38
 39        declare(HeavyHadrons(),"HA");
 40
 41        FastJets jetpro(FinalState(), FastJets::ANTIKT, 0.4, JetAlg::Muons::DECAY, JetAlg::Invisibles::DECAY);
 42        declare(jetpro, "Jets");
 43
 44        // bar charts
 45        book(_h["bar_b_jet_width"], "width_B_jet", 7, 0., 0.3);
 46        book(_h["bar_c_jet_width"], "width_C_jet", 7, 0., 0.3);
 47        // profiles
 48        book(_p["b_jet_rho"], "avg_rho_B_jet", 10, 0., 0.4);
 49        book(_p["c_jet_rho"], "avg_rho_C_jet", 10, 0., 0.4);
 50        book(_p["b_jet_psi"], "avg_psi_B_jet", 10, 0., 0.4);
 51        book(_p["c_jet_psi"], "avg_psi_C_jet", 10, 0., 0.4);
 52        // histograms
 53        book(_h["b_frac"], "prod_frac_B", 11, 0.5,11.5);
 54        book(_h["c_frac"], "prod_frac_C",  7, 0.5, 7.5);
 55        book(_h["b_jet_frag"], "frag_B_jet", 50, 0., 1.1);
 56        book(_h["c_jet_frag"], "frag_C_jet", 50, 0., 1.1);
 57        book(_h["b_jet_pT"], "pT_B_jet", 25, 25., 1025.);
 58        book(_h["c_jet_pT"], "pT_C_jet", 25, 25., 1025.);
 59        book(_h["b_jet_pThad"], "pThad_B_jet", 10, 0., 200.);
 60        book(_h["c_jet_pThad"], "pThad_C_jet", 10, 0., 200.);
 61        book(_h["b_jet_high_pT"], "high_pT_frag_B_jet", 46, 0., 1.);
 62        book(_h["c_jet_high_pT"], "high_pT_frag_C_jet", 46, 0., 1.);
 63        book(_h["b_jet_high_pT_1H"], "high_pT_frag_B_jet_1H", 46, 0., 1.);
 64        book(_h["c_jet_high_pT_1H"], "high_pT_frag_C_jet_1H", 46, 0., 1.);
 65
 66
 67        book(_h["ch_B0"         ], "B0_charged_mult",          30, 0.5, 30.5);
 68        book(_h["ch_BPLUS"      ], "BPLUS_charged_mult",       30, 0.5, 30.5);
 69        book(_h["ch_B0S"        ], "B0S_charged_mult",         30, 0.5, 30.5);
 70        book(_h["ch_LAMBDAB"    ], "LAMBDAB_charged_mult",     30, 0.5, 30.5);
 71        book(_h["ch_D0"         ], "D0_charged_mult",          16, 0.5, 16.5);
 72        book(_h["ch_DPLUS"      ], "DPLUS_charged_multh",      16, 0.5, 16.5);
 73        book(_h["ch_DSPLUS"     ], "DSPLUS_charged_mult",      16, 0.5, 16.5);
 74        book(_h["ch_LAMBDACPLUS"], "LAMBDACPLUS_charged_mult", 16, 0.5, 16.5);
 75
 76        book(_h["st_B0"         ], "B0_stable_mult",          40, 0.5, 40.5);
 77        book(_h["st_BPLUS"      ], "BPLUS_stable_mult",       40, 0.5, 40.5);
 78        book(_h["st_B0S"        ], "B0S_stable_mult",         40, 0.5, 40.5);
 79        book(_h["st_LAMBDAB"    ], "LAMBDAB_stable_mult",     40, 0.5, 40.5);
 80        book(_h["st_D0"         ], "D0_stable_mult",          20, 0.5, 20.5);
 81        book(_h["st_DPLUS"      ], "DPLUS_stable_mult",       20, 0.5, 20.5);
 82        book(_h["st_DSPLUS"     ], "DSPLUS_stable_mult",      20, 0.5, 20.5);
 83        book(_h["st_LAMBDACPLUS"], "LAMBDACPLUS_stable_mult", 20, 0.5, 20.5);
 84
 85        book(_h["pt_B0"         ], "B0_pT",          10, 25., 425.);
 86        book(_h["pt_BPLUS"      ], "BPLUS_pT",       10, 25., 425.);
 87        book(_h["pt_B0S"        ], "B0S_pT",         10, 25., 425.);
 88        book(_h["pt_LAMBDAB"    ], "LAMBDAB_pT",     10, 25., 425.);
 89        book(_h["pt_D0"         ], "D0_pT",          10, 25., 425.);
 90        book(_h["pt_DPLUS"      ], "DPLUS_pT",       10, 25., 425.);
 91        book(_h["pt_DSPLUS"     ], "DSPLUS_pT",      10, 25., 425.);
 92        book(_h["pt_LAMBDACPLUS"], "LAMBDACPLUS_pT", 10, 25., 425.);
 93
 94        book(_h["b_jet_ch_mult"], "charged_mult_B_jets", 40, 0.5, 40.5);
 95        book(_h["c_jet_ch_mult"], "charged_mult_C_jets", 40, 0.5, 40.5);
 96
 97        book(_h["b_jet_l_pTrel"], "lepton_pTrel_B_jets", 8, 0., 15.);
 98        book(_h["c_jet_l_pTrel"], "lepton_pTrel_C_jets", 8, 0., 10.);
 99
100        book(_h["b_jet_l_pT"], "lepton_pT_B_jets", 10, 0., 100.);
101        book(_h["c_jet_l_pT"], "lepton_pT_C_jets", 10, 0., 100.);
102
103        // double-differentials
104        ptaxis = YODA::Histo1DAxis({25, 30, 50, 70, 100, 150, 300, 500, 1000});
105        for (size_t i = 0; i < ptaxis.numBins(); ++i) {
106          string suff = to_string(int(ptaxis.bin(i).xMin())) + "_" + to_string(int(ptaxis.bin(i).xMax()));
107          book(_p["avg_B_jet_rho_"+to_string(i)], "avg_B_jet_rho_"+suff, 10, 0., 0.4);
108          book(_p["avg_C_jet_rho_"+to_string(i)], "avg_C_jet_rho_"+suff, 10, 0., 0.4);
109
110          addBinned("avg_B_jet_ch_mult", ptaxis.bin(i), 40, 0.5, 40.5);
111          addBinned("avg_C_jet_ch_mult", ptaxis.bin(i), 40, 0.5, 40.5);
112          if (i == 0) {
113            addBinned("avg_B_jet_l_pTrel", ptaxis.bin(i), 2, 0., 4.);
114            addBinned("avg_C_jet_l_pTrel", ptaxis.bin(i), 2, 0., 3.);
115          }
116          else if (i <= 2) {
117            addBinned("avg_B_jet_l_pTrel", ptaxis.bin(i), 4, 0., 8.);
118            addBinned("avg_C_jet_l_pTrel", ptaxis.bin(i), 5, 0., 6.);
119          }
120          else {
121            addBinned("avg_B_jet_l_pTrel", ptaxis.bin(i), 8, 0., 15.);
122            addBinned("avg_C_jet_l_pTrel", ptaxis.bin(i), 8, 0., 10.);
123          }
124        }
125      }
126
127      double p_annulus(const Jet &jet, const double a, const double b) const {
128        // calculate the total momentum inside an annulus with a <= R < b
129        return sum(filter_select(jet.particles(), [&](const Particle &p) {
130          const double dr = deltaR(p, jet);
131          return (dr < b && dr >= a);
132        }), Kin::pT, 0.)/GeV;
133      }
134
135      void count_mult(const Particle& p) {
136        unsigned int nst = p.stableDescendants().size() + 0.5;
137        unsigned int nch = p.stableDescendants(Cuts::charge != 0).size() + 0.5;
138        _h["st_"+whoDis(p.pid())]->fill(nst);
139        _h["ch_"+whoDis(p.pid())]->fill(nch);
140      }
141
142      double pTrel (const Jet& jet, const Particle& p) const {
143        return (p.p3().cross(jet.p3())).mod()/(jet.p3().mod());
144      }
145
146
147      /// Perform the per-event analysis
148      void analyze(const Event& event) {
149
150        const HeavyHadrons &ha = apply<HeavyHadrons>(event,"HA");
151
152        if (ha.bHadrons().empty() && ha.cHadrons().empty()) vetoEvent;
153
154        for (const Particle &p : ha.bHadrons()) {
155          const string name = "pt_" + whoDis(p.pid());
156          if (p.pid() == PID::B0) {
157            count_mult(p);
158            _h["b_frac"]->fill(1);
159            _h[name]->fill(p.pT()/GeV);
160          }
161          else if (p.pid() == PID::BPLUS) {
162            count_mult(p);
163            _h["b_frac"]->fill(2);
164            _h[name]->fill(p.pT()/GeV);
165          }
166          else if (p.pid() == PID::B0S) {
167            count_mult(p);
168            _h["b_frac"]->fill(3);
169            _h[name]->fill(p.pT()/GeV);
170          }
171          else if (p.pid() == PID::BCPLUS)   _h["b_frac"]->fill(4);
172          else if (p.pid() == PID::LAMBDAB) {
173            count_mult(p);
174            _h["b_frac"]->fill(5);
175            _h[name]->fill(p.pT()/GeV);
176          }
177          else if (p.pid() == PID::XIBMINUS)     _h["b_frac"]->fill(6);
178          else if (p.pid() == PID::XI0B)         _h["b_frac"]->fill(7);
179          else if (p.pid() == PID::OMEGABMINUS)  _h["b_frac"]->fill(8);
180          else if (p.pid() == PID::SIGMABMINUS)  _h["b_frac"]->fill(9);
181          else if (p.pid() == PID::SIGMAB)       _h["b_frac"]->fill(10);
182          else if (p.pid() == PID::SIGMABPLUS)   _h["b_frac"]->fill(11);
183        }
184
185        for (const Particle &p : ha.cHadrons()) {
186          const string name = "pt_" + whoDis(p.pid());
187          if (p.pid() == PID::DPLUS) {
188            count_mult(p);
189            _h["c_frac"]->fill(1);
190            _h[name]->fill(p.pT()/GeV);
191          }
192          else if (p.pid() == PID::D0) {
193            count_mult(p);
194            _h["c_frac"]->fill(2);
195            _h[name]->fill(p.pT()/GeV);
196          }
197          else if (p.pid() == PID::DSPLUS) {
198            count_mult(p);
199            _h["c_frac"]->fill(3);
200            _h[name]->fill(p.pT()/GeV);
201          }
202          else if (p.pid() == PID::LAMBDACPLUS) {
203            count_mult(p);
204            _h["c_frac"]->fill(4);
205            _h[name]->fill(p.pT()/GeV);
206          }
207          else if (p.pid() == PID::XICPLUS)  _h["c_frac"]->fill(5);
208          else if (p.pid() == PID::XI0C)     _h["c_frac"]->fill(6);
209          else if (p.pid() == PID::OMEGA0C)  _h["c_frac"]->fill(7);
210        }
211
212        Jets jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 25.*GeV && Cuts::absrap < 2.5);
213        if (jets.empty()) vetoEvent;
214
215        auto r_bins = _p["b_jet_rho"]->xEdges();
216        const double dr = r_bins[1]-r_bins[0]; //dr is equal for all bins
217        for (const Jet& thisJet : jets) {
218          double p_0_R = p_annulus(thisJet, 0., 0.4);
219          if (fuzzyEquals(p_0_R, 0., 1e-5))  continue;
220          Particles bjets = thisJet.bTags(Cuts::pT > 5*GeV);
221          ifilter_select(bjets, deltaRLess(thisJet, 0.3));
222          for (const Particle &thisB : bjets) {
223            _h["b_jet_pThad"]->fill(thisB.pT()/GeV);
224            double z = thisJet.p3().dot(thisB.p3())/thisJet.p2();
225            _h["b_jet_frag"]->fill(z);
226            if (inRange(thisJet.pT(), 500*GeV, 1000*GeV)) {
227              _h["b_jet_high_pT"]->fill(z);
228              if (bjets.size() == 1)  _h["b_jet_high_pT_1H"]->fill(z);
229            }
230            for (size_t i = 0; i < r_bins.size()-1; ++i) {
231              double r = 0.5*(r_bins[i]+r_bins[i+1]);
232              double rho = p_annulus(thisJet, r_bins[i], r_bins[i+1])/p_0_R/dr;
233              double psi = p_annulus(thisJet, 0, r_bins[i+1])/p_0_R;
234              _p["b_jet_rho"]->fill(r, rho);
235              _p["b_jet_psi"]->fill(r, psi);
236              size_t ptb = ptaxis.binIndexAt(thisJet.pT()/GeV); // index in {0, nBins - 1}
237              if (ptb < ptaxis.numBins()) _p["avg_B_jet_rho_"+to_string(ptb)]->fill(r, rho);
238            }
239          }
240
241          Particles cjets = thisJet.cTags(Cuts::pT > 5*GeV);
242          ifilter_select(cjets, deltaRLess(thisJet, 0.3));
243          if (bjets.empty()) {
244            for (const Particle& thisC : cjets) {
245              _h["c_jet_pThad"]->fill(thisC.pT()/GeV);
246              double z = thisJet.p3().dot(thisC.p3())/thisJet.p2();
247              _h["c_jet_frag"]->fill(z);
248              if (inRange(thisJet.pT(), 500*GeV, 1000*GeV)) {
249                _h["c_jet_high_pT"]->fill(z);
250                if (cjets.size() == 1) {
251                  _h["c_jet_high_pT_1H"]->fill(z);
252                }
253              }
254              for (size_t i = 0; i < r_bins.size()-1; ++i) {
255                double r = 0.5*(r_bins[i]+r_bins[i+1]);
256                double rho = p_annulus(thisJet, r_bins[i], r_bins[i+1])/p_0_R/dr;
257                double psi = p_annulus(thisJet, 0, r_bins[i+1])/p_0_R;
258                _p["c_jet_rho"]->fill(r, rho);
259                _p["c_jet_psi"]->fill(r, psi);
260                size_t ptb = ptaxis.binIndexAt(thisJet.pT()/GeV); // index in {0, nBins - 1}
261                if (ptb < ptaxis.numBins()) _p["avg_C_jet_rho_"+to_string(ptb)]->fill(r, rho);
262              }
263            }
264          }
265
266          if (bjets.size()) {
267            double W_num = 0., W_den = 0.;
268            long N_charged = 0;
269            for (const Particle& pp : thisJet.particles()) {
270              if(pp.isStable()) {
271                W_num += pp.pT()*deltaR(thisJet,pp);
272                W_den += pp.pT();
273                if (pp.isCharged())  ++N_charged;
274              }
275              if(pp.isLepton()) {
276                _h["b_jet_l_pT"]->fill(pp.pT()/GeV);
277                _h["b_jet_l_pTrel"]->fill(pTrel(thisJet,pp)/GeV);
278                _b["avg_B_jet_l_pTrel"].fill(thisJet.pT()/GeV, pTrel(thisJet,pp)/GeV);
279              }
280            }
281            if (W_den)  _h["bar_b_jet_width"]->fill(W_num/W_den);
282            _h["b_jet_pT"]->fill(thisJet.pT()/GeV);
283            if (N_charged) {
284              _h["b_jet_ch_mult"]->fill(N_charged);
285              _b["avg_B_jet_ch_mult"].fill(thisJet.pT()/GeV, N_charged);
286            }
287          }
288          else if(cjets.size()) {
289            double W_num = 0., W_den = 0.;
290            long N_charged = 0;
291            for (const Particle& pp : thisJet.particles()) {
292              if(pp.isStable()) {
293                W_num += pp.pT()*deltaR(thisJet,pp);
294                W_den += pp.pT();
295                if(pp.isCharged()) ++N_charged;
296              }
297              if(pp.isLepton()) {
298                _h["c_jet_l_pT"]->fill(pp.pT()/GeV);
299                _h["c_jet_l_pTrel"]->fill(pTrel(thisJet,pp)/GeV);
300                _b["avg_C_jet_l_pTrel"].fill(thisJet.pT()/GeV, pTrel(thisJet,pp)/GeV);
301              }
302            }
303            if (W_den)  _h["bar_c_jet_width"]->fill(W_num/W_den);
304            _h["c_jet_pT"]->fill(thisJet.pT()/GeV);
305            if (N_charged) {
306              _h["c_jet_ch_mult"]->fill(N_charged);
307              _b["avg_C_jet_ch_mult"].fill(thisJet.pT()/GeV, N_charged);
308            }
309          }
310        }
311      }
312
313      /// Normalise histograms etc., after the run
314      void finalize() {
315        for (const auto &hit : _h) {
316          double sf = 1.0;
317          if (hit.first.find("bar_") != string::npos)  sf = (hit.second->xMax()-hit.second->xMin())/hit.second->numBins();
318          normalize(hit.second, sf);
319        }
320        for (const auto &bit : _b) {
321          for (const auto& hist : bit.second.histos()) { normalize(hist); }
322        }
323      }
324
325    private:
326      map<string, Histo1DPtr> _h;
327      map<string, Profile1DPtr> _p;
328      map<string, BinnedHistogram> _b;
329      YODA::Histo1DAxis ptaxis;
330
331  };
332
333  RIVET_DECLARE_PLUGIN(MC_HFDECAYS);
334}