Rivet analyses referenceMC_HFJETSMonte Carlo validation analysis to study heavy flavour productionExperiment: () Status: VALIDATED Authors:
Beams: * * Beam energies: ANY Run details:
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 // set clustering radius from input option
23 const double R = getOption<double>("R", 0.6);
24
25 // set clustering algorithm from input option
26 JetAlg clusterAlgo;
27 const string algoopt = getOption("ALGO", "ANTIKT");
28 if ( algoopt == "KT" ) {
29 clusterAlgo = JetAlg::KT;
30 } else if ( algoopt == "CA" ) {
31 clusterAlgo = JetAlg::CA;
32 } else if ( algoopt == "ANTIKT" ) {
33 clusterAlgo = JetAlg::ANTIKT;
34 } else {
35 MSG_WARNING("Unknown jet clustering algorithm option " + algoopt + ". Defaulting to anti-kT");
36 clusterAlgo = JetAlg::ANTIKT;
37 }
38
39 FastJets fj(FinalState(Cuts::abseta < 5), clusterAlgo, R);
40 fj.useInvisibles();
41 declare(fj, "Jets");
42 declare(HeavyHadrons(Cuts::abseta < 5 && Cuts::pT > 500*MeV), "BCHadrons");
43
44 book(_h_ptCJetLead ,"ptCJetLead", linspace(5, 0, 20, false) + logspace(25, 20, 200));
45 book(_h_ptCHadrLead ,"ptCHadrLead", linspace(5, 0, 10, false) + logspace(25, 10, 200));
46 book(_h_ptFracC ,"ptfracC", 50, 0, 1.5);
47 book(_h_eFracC ,"efracC", 50, 0, 1.5);
48
49 book(_h_ptBJetLead ,"ptBJetLead", linspace(5, 0, 20, false) + logspace(25, 20, 200));
50 book(_h_ptBHadrLead ,"ptBHadrLead", linspace(5, 0, 10, false) + logspace(25, 10, 200));
51 book(_h_ptFracB ,"ptfracB", 50, 0, 1.5);
52 book(_h_eFracB ,"efracB", 50, 0, 1.5);
53 }
54
55
56 /// Perform the per-event analysis
57 void analyze(const Event& event) {
58
59 // Get jets and heavy hadrons
60 const Jets& jets = apply<JetFinder>(event, "Jets").jetsByPt();
61 const Particles bhadrons = sortByPt(apply<HeavyHadrons>(event, "BCHadrons").bHadrons());
62 const Particles chadrons = sortByPt(apply<HeavyHadrons>(event, "BCHadrons").cHadrons());
63 MSG_DEBUG("# b hadrons = " << bhadrons.size() << ", # c hadrons = " << chadrons.size());
64
65 // Loop over jets and use ghost-tag info
66 for (const Jet& j : jets) {
67 bool gotLeadingB = false, gotLeadingC = false;
68 // b-tag testing
69 if (!gotLeadingB && j.bTagged(Cuts::pT > 500*MeV)) {
70 gotLeadingB = true;
71 Particle bhad = sortByPt(j.bTags(Cuts::pT > 500*MeV))[0];
72 _h_ptBJetLead->fill(j.pT()/GeV);
73 _h_ptBHadrLead->fill(bhad.pT()/GeV);
74 _h_ptFracB->fill(bhad.pT() / j.pT());
75 _h_eFracB->fill(bhad.E() / j.E());
76 continue;
77 }
78 // c-tag testing
79 if (!gotLeadingC && j.cTagged(Cuts::pT > 500*MeV) && !j.bTagged(Cuts::pT > 500*MeV)) {
80 gotLeadingC = true;
81 Particle chad = sortByPt(j.cTags(Cuts::pT > 500*MeV))[0];
82 _h_ptCJetLead->fill(j.pT()/GeV);
83 _h_ptCHadrLead->fill(chad.pT()/GeV);
84 _h_ptFracC->fill(chad.pT() / j.pT());
85 _h_eFracC->fill(chad.E() / j.E());
86 }
87 // Escape early if we've found both the leading b and c jets
88 if (gotLeadingB && gotLeadingC) break;
89 }
90
91 // // Tag the leading b and c jets with a deltaR < 0.3 match
92 // // b-tagged jet are excluded from also being considered as c-tagged
93 // MSG_DEBUG("Getting b/c-tags");
94 // const double MAX_DR = 0.3;
95 // bool gotLeadingB = false, gotLeadingC = false;
96 // for (const Jet& j : jets) {
97 // if (!gotLeadingB) {
98 // FourMomentum leadBJet, leadBHadr;
99 // double dRmin = MAX_DR;
100 // for (const Particle& b : bhadrons) {
101 // const double dRcand = min(dRmin, deltaR(j, b));
102 // if (dRcand < dRmin) {
103 // dRmin = dRcand;
104 // leadBJet = j.momentum();
105 // leadBHadr = b.momentum();
106 // MSG_DEBUG("New closest b-hadron jet tag candidate: dR = " << dRmin
107 // << " for jet pT = " << j.pT()/GeV << " GeV, "
108 // << " b hadron pT = " << b.pT()/GeV << " GeV, PID = " << b.pid());
109 // }
110 // }
111 // if (dRmin < MAX_DR) {
112 // // A jet has been tagged, so fill the histos and break the loop
113 // _h_ptBJetLead->fill(leadBJet.pT()/GeV);
114 // _h_ptBHadrLead->fill(leadBHadr.pT()/GeV);
115 // _h_ptFracB->fill(leadBHadr.pT() / leadBJet.pT());
116 // _h_eFracB->fill(leadBHadr.E() / leadBJet.E());
117 // gotLeadingB = true;
118 // continue; // escape this loop iteration so the same jet isn't c-tagged
119 // }
120 // }
121 // if (!gotLeadingC) {
122 // FourMomentum leadCJet, leadCHadr;
123 // double dRmin = MAX_DR;
124 // for (const Particle& c : chadrons) {
125 // const double dRcand = min(dRmin, deltaR(j, c));
126 // if (dRcand < dRmin) {
127 // dRmin = dRcand;
128 // leadCJet = j.momentum();
129 // leadCHadr = c.momentum();
130 // MSG_DEBUG("New closest c-hadron jet tag candidate: dR = " << dRmin
131 // << " for jet pT = " << j.pT()/GeV << " GeV, "
132 // << " c hadron pT = " << c.pT()/GeV << " GeV, PID = " << c.pid());
133 // }
134 // }
135 // if (dRmin < MAX_DR) {
136 // // A jet has been tagged, so fill the histos and break the loop
137 // _h_ptCJetLead->fill(leadCJet.pT()/GeV);
138 // _h_ptCHadrLead->fill(leadCHadr.pT()/GeV);
139 // _h_ptFracC->fill(leadCHadr.pT() / leadCJet.pT());
140 // _h_eFracC->fill(leadCHadr.E() / leadCJet.E());
141 // gotLeadingB = true;
142 // }
143 // }
144 // // If we've found both a leading b and a leading c jet, break the loop over jets
145 // if (gotLeadingB && gotLeadingC) break;
146 // }
147
148 }
149
150
151 /// Normalise histograms etc., after the run
152 void finalize() {
153 normalize({_h_ptCJetLead, _h_ptCHadrLead, _h_ptBJetLead, _h_ptBHadrLead,
154 _h_ptFracC, _h_eFracC, _h_ptFracB, _h_eFracB});
155 }
156
157
158 /// @name Histograms
159 /// @{
160 Histo1DPtr _h_ptCJetLead, _h_ptCHadrLead, _h_ptFracC, _h_eFracC;
161 Histo1DPtr _h_ptBJetLead, _h_ptBHadrLead, _h_ptFracB, _h_eFracB;
162 /// @}
163
164
165 };
166
167
168
169 RIVET_DECLARE_PLUGIN(MC_HFJETS);
170
171}
|