Rivet analyses referenceMC_LEADJETUEUnderlying event in leading jet events, extended to LHCExperiment: () Status: VALIDATED Authors:
Beams: * * Beam energies: ANY Run details:
Rick Field's measurement of the underlying event in leading jet events, extended to the LHC. As usual, the leading jet of the defines an azimuthal toward/transverse/away decomposition, in this case the event is accepted within $|\eta| < 2$, as in the CDF 2008 version of the analysis. Since this isn't the Tevatron, I've chosen to use $k_\perp$ rather than midpoint jets. Source code: MC_LEADJETUE.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/ChargedFinalState.hh"
5#include "Rivet/Projections/FastJets.hh"
6
7namespace Rivet {
8
9
10 /// @brief MC validation analysis for underlying event in jet events
11 /// @author Andy Buckley
12 class MC_LEADJETUE : public Analysis {
13 public:
14
15 /// Constructor
16 MC_LEADJETUE()
17 : Analysis("MC_LEADJETUE")
18 { }
19
20
21 /// @name Analysis methods
22 /// @{
23
24 // Book histograms
25 void init() {
26 // Final state for the jet finding
27 const FinalState fsj((Cuts::etaIn(-4.0, 4.0)));
28 declare(fsj, "FSJ");
29 declare(FastJets(fsj, JetAlg::KT, 0.7), "Jets");
30
31 // Charged final state for the distributions
32 const ChargedFinalState cfs((Cuts::etaIn(-1.0, 1.0) && Cuts::pT >= 0.5*GeV));
33 declare(cfs, "CFS");
34
35 const double maxpt1 = 500.0;
36 book(_hist_pnchg ,"trans-nchg", 50, 0.0, maxpt1);
37 book(_hist_pmaxnchg ,"trans-maxnchg", 50, 0.0, maxpt1);
38 book(_hist_pminnchg ,"trans-minnchg", 50, 0.0, maxpt1);
39 book(_hist_pcptsum ,"trans-ptsum", 50, 0.0, maxpt1);
40 book(_hist_pmaxcptsum ,"trans-maxptsum", 50, 0.0, maxpt1);
41 book(_hist_pmincptsum ,"trans-minptsum", 50, 0.0, maxpt1);
42 book(_hist_pcptave ,"trans-ptavg", 50, 0.0, maxpt1);
43 }
44
45
46 // Do the analysis
47 void analyze(const Event& e) {
48
49 const FinalState& fsj = apply<FinalState>(e, "FSJ");
50 if (fsj.particles().empty()) {
51 MSG_DEBUG("Failed multiplicity cut");
52 vetoEvent;
53 }
54
55 const FastJets& jetpro = apply<FastJets>(e, "Jets");
56 const Jets jets = jetpro.jetsByPt();
57 MSG_DEBUG("Jet multiplicity = " << jets.size());
58
59 // Require the leading jet to be within |eta| < 2
60 if (jets.size() < 1 || fabs(jets[0].eta()) > 2) {
61 MSG_DEBUG("Failed jet cut");
62 vetoEvent;
63 }
64
65 const double jetphi = jets[0].phi();
66 const double jetpT = jets[0].pT();
67 MSG_DEBUG("Leading jet: pT = " << jetpT/GeV << " GeV"
68 << ", eta = " << jets[0].eta()
69 << ", phi = " << jetphi);
70
71 // Get the final states to work with for filling the distributions
72 const FinalState& cfs = apply<ChargedFinalState>(e, "CFS");
73
74 //size_t numOverall(0), numToward(0), numAway(0);
75 size_t numTrans1(0), numTrans2(0);
76 double ptSumTrans1(0.0), ptSumTrans2(0.0);
77 double ptMaxOverall(0.0), ptMaxToward(0.0), ptMaxTrans1(0.0), ptMaxTrans2(0.0), ptMaxAway(0.0);
78
79 // Calculate all the charged stuff
80 for (const Particle& p : cfs.particles()) {
81 const double dPhi = deltaPhi(p.phi(), jetphi);
82 const double pT = p.pT();
83 const double phi = p.phi();
84 const double rotatedphi = phi - jetphi;
85
86 //ptSumOverall += pT;
87 //++numOverall;
88 if (pT > ptMaxOverall) ptMaxOverall = pT;
89
90 if (dPhi < PI/3.0) {
91 //ptSumToward += pT;
92 //++numToward;
93 if (pT > ptMaxToward) ptMaxToward = pT;
94 }
95 else if (dPhi < 2*PI/3.0) {
96 if (rotatedphi <= PI) {
97 ptSumTrans1 += pT;
98 ++numTrans1;
99 if (pT > ptMaxTrans1) ptMaxTrans1 = pT;
100 } else {
101 ptSumTrans2 += pT;
102 ++numTrans2;
103 if (pT > ptMaxTrans2) ptMaxTrans2 = pT;
104 }
105 }
106 else {
107 //ptSumAway += pT;
108 //++numAway;
109 if (pT > ptMaxAway) ptMaxAway = pT;
110 }
111 }
112
113
114 // Fill the histograms
115 //_hist_tnchg->fill(jetpT/GeV, numToward/(4*PI/3));
116 _hist_pnchg->fill(jetpT/GeV, (numTrans1+numTrans2)/(4*PI/3));
117 _hist_pmaxnchg->fill(jetpT/GeV, (numTrans1>numTrans2 ? numTrans1 : numTrans2)/(2*PI/3));
118 _hist_pminnchg->fill(jetpT/GeV, (numTrans1<numTrans2 ? numTrans1 : numTrans2)/(2*PI/3));
119 //_hist_pdifnchg->fill(jetpT/GeV, abs(numTrans1-numTrans2)/(2*PI/3));
120 //_hist_anchg->fill(jetpT/GeV, numAway/(4*PI/3));
121
122 //_hist_tcptsum->fill(jetpT/GeV, ptSumToward/GeV/(4*PI/3));
123 _hist_pcptsum->fill(jetpT/GeV, (ptSumTrans1+ptSumTrans2)/GeV/(4*PI/3));
124 _hist_pmaxcptsum->fill(jetpT/GeV, (ptSumTrans1>ptSumTrans2 ? ptSumTrans1 : ptSumTrans2)/GeV/(2*PI/3));
125 _hist_pmincptsum->fill(jetpT/GeV, (ptSumTrans1<ptSumTrans2 ? ptSumTrans1 : ptSumTrans2)/GeV/(2*PI/3));
126 //_hist_pdifcptsum->fill(jetpT/GeV, fabs(ptSumTrans1-ptSumTrans2)/GeV/(2*PI/3));
127 //_hist_acptsum->fill(jetpT/GeV, ptSumAway/GeV/(4*PI/3));
128
129 //if (numToward > 0) {
130 // _hist_tcptave->fill(jetpT/GeV, ptSumToward/GeV/numToward);
131 // _hist_tcptmax->fill(jetpT/GeV, ptMaxToward/GeV);
132 //}
133 if ((numTrans1+numTrans2) > 0) {
134 _hist_pcptave->fill(jetpT/GeV, (ptSumTrans1+ptSumTrans2)/GeV/(numTrans1+numTrans2));
135 //_hist_pcptmax->fill(jetpT/GeV, (ptMaxTrans1 > ptMaxTrans2 ? ptMaxTrans1 : ptMaxTrans2)/GeV);
136 }
137 //if (numAway > 0) {
138 // _hist_acptave->fill(jetpT/GeV, ptSumAway/GeV/numAway);
139 // _hist_acptmax->fill(jetpT/GeV, ptMaxAway/GeV);
140 //}
141 }
142
143
144 void finalize() {
145 //
146 }
147
148
149 private:
150
151 Profile1DPtr _hist_pnchg;
152 Profile1DPtr _hist_pmaxnchg;
153 Profile1DPtr _hist_pminnchg;
154 Profile1DPtr _hist_pcptsum;
155 Profile1DPtr _hist_pmaxcptsum;
156 Profile1DPtr _hist_pmincptsum;
157 Profile1DPtr _hist_pcptave;
158
159 };
160
161
162
163 RIVET_DECLARE_PLUGIN(MC_LEADJETUE);
164
165}
|