Rivet analyses referenceCDF_2010_I849042CDF Run 2 underlying event in Drell-YanExperiment: CDF (Tevatron Run 2) Inspire ID: 849042 Status: VALIDATED Authors:
Beam energies: (980.0, 980.0) GeV Run details:
Deepak Kar and Rick Field's measurement of the underlying event in Drell-Yan events. $Z -> ee$ and $Z -> \mu\mu$ events are selected using a $Z$ mass window cut between 70 and 110 GeV. ``Toward'', ``away'' and ``transverse'' regions are defined in the same way as in the original (2001) CDF underlying event analysis. The reconstructed $Z$ defines the $\phi$ direction of the toward region. The leptons are ignored after the $Z$ has been reconstructed. Thus the region most sensitive to the underlying event is the toward region (the recoil jet is boosted into the away region). Source code: CDF_2010_I849042.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/ChargedFinalState.hh"
4#include "Rivet/Projections/ChargedLeptons.hh"
5#include "Rivet/Projections/FastJets.hh"
6
7namespace Rivet {
8
9
10 /// @brief CDF Run II underlying event in Drell-Yan
11 ///
12 /// @author Hendrik Hoeth
13 ///
14 /// Measurement of the underlying event in Drell-Yan
15 /// \f$ Z/\gamma^* \to e^+ e^- \f$ and
16 /// \f$ Z/\gamma^* \to \mu^+ \mu^- \f$ events. The reconstructed
17 /// Z defines the \f$ \phi \f$ orientation. A Z mass window cut is applied.
18 ///
19 /// @par Run conditions
20 ///
21 /// @arg \f$ \sqrt{s} = \f$ 1960 GeV
22 /// @arg produce Drell-Yan events
23 /// @arg Set particles with c*tau > 10 mm stable
24 /// @arg Z decay mode: Z -> e+e- and Z -> mu+mu-
25 /// @arg gamma decay mode: gamma -> e+e- and gamma -> mu+mu-
26 /// @arg minimum invariant mass of the fermion pair coming from the Z/gamma: 70 GeV
27 class CDF_2010_I849042 : public Analysis {
28 public:
29
30 RIVET_DEFAULT_ANALYSIS_CTOR(CDF_2010_I849042);
31
32
33 /// @name Analysis methods
34 /// @{
35
36 void init() {
37
38 _mode = 0;
39 if ( getOption("MODE") == "DY" ) _mode = 1;
40 else if ( getOption("MODE") == "QCD" ) _mode = 2;
41
42 // Set up projections
43 const ChargedFinalState cfs(Cuts::abseta < 1.0 && Cuts::pT >= 0.5*GeV);
44 const ChargedFinalState clfs(Cuts::abseta < 1.0 && Cuts::pT >= 20*GeV);
45 declare(cfs, "CFS");
46 declare(ChargedLeptons(clfs), "CL");
47
48 // Final state for the jet finding
49 const FinalState fsj(Cuts::abseta < 4.0);
50 declare(fsj, "FSJ");
51 declare(FastJets(fsj, JetAlg::CDFMIDPOINT, 0.7), "MidpointJets");
52
53 // Book histograms
54 if (_mode == 0 || _mode == 1) {
55 book(_p["tnchg"], 1, 1, 1);
56 book(_p["pnchg"], 1, 1, 2);
57 book(_p["anchg"], 1, 1, 3);
58 book(_p["pmaxnchg"], 2, 1, 1);
59 book(_p["pminnchg"], 2, 1, 2);
60 book(_p["pdifnchg"], 2, 1, 3);
61 book(_p["tcptsum"], 3, 1, 1);
62 book(_p["pcptsum"], 3, 1, 2);
63 book(_p["acptsum"], 3, 1, 3);
64 book(_p["pmaxcptsum"], 4, 1, 1);
65 book(_p["pmincptsum"], 4, 1, 2);
66 book(_p["pdifcptsum"], 4, 1, 3);
67 book(_p["tcptave"], 5, 1, 1);
68 book(_p["pcptave"], 5, 1, 2);
69 book(_p["tcptmax"], 6, 1, 1);
70 book(_p["pcptmax"], 6, 1, 2);
71 book(_p["zptvsnchg"], 7, 1, 1);
72 book(_p["cptavevsnchg"], 8, 1, 1);
73 book(_p["cptavevsnchgsmallzpt"], 9, 1, 1);
74 }
75
76 if (_mode == 0 || _mode == 2) {
77 book(_p["tnchg"], 10, 1, 1);
78 book(_p["pnchg"], 10, 1, 2);
79 book(_p["anchg"], 10, 1, 3);
80 book(_p["pmaxnchg"], 11, 1, 1);
81 book(_p["pminnchg"], 11, 1, 2);
82 book(_p["pdifnchg"], 11, 1, 3);
83 book(_p["tcptsum"], 12, 1, 1);
84 book(_p["pcptsum"], 12, 1, 2);
85 book(_p["acptsum"], 12, 1, 3);
86 book(_p["pmaxcptsum"], 13, 1, 1);
87 book(_p["pmincptsum"], 13, 1, 2);
88 book(_p["pdifcptsum"], 13, 1, 3);
89 book(_p["pcptave"], 14, 1, 1);
90 book(_p["pcptmax"], 15, 1, 1);
91 }
92 }
93
94
95 /// Do the analysis
96 void analyze(const Event& event) {
97
98 if (_mode == 0 || _mode == 1) doDYanalysis(event);
99 if (_mode == 0 || _mode == 2) doQCDanalysis(event);
100
101 }
102
103
104 void doDYanalysis(const Event& e) {
105
106 const FinalState& fs = apply<FinalState>(e, "CFS");
107 const size_t numParticles = fs.particles().size();
108
109 // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
110 if (numParticles < 1) {
111 MSG_DEBUG("Failed multiplicity cut");
112 vetoEvent;
113 }
114
115 // Get the leptons
116 const Particles& leptons = apply<ChargedLeptons>(e, "CL").chargedLeptons();
117
118 // We want exactly two leptons of the same flavour.
119 MSG_DEBUG("lepton multiplicity = " << leptons.size());
120 if (leptons.size() != 2 || leptons[0].pid() != -leptons[1].pid() ) vetoEvent;
121
122 // Lepton pT > 20 GeV
123 if (leptons[0].pT()/GeV <= 20 || leptons[1].pT()/GeV <= 20) vetoEvent;
124
125 // Lepton pair should have an invariant mass between 70 and 110 and |eta| < 6
126 const FourMomentum dilepton = leptons[0].momentum() + leptons[1].momentum();
127 if (!inRange(dilepton.mass()/GeV, 70., 110.) || fabs(dilepton.eta()) >= 6) vetoEvent;
128 MSG_DEBUG("Dilepton mass = " << dilepton.mass()/GeV << " GeV");
129 MSG_DEBUG("Dilepton pT = " << dilepton.pT()/GeV << " GeV");
130
131 // Calculate the observables
132 size_t numToward(0), numAway(0);
133 long int numTrans1(0), numTrans2(0);
134 double ptSumToward(0.0), ptSumTrans1(0.0), ptSumTrans2(0.0), ptSumAway(0.0);
135 double ptMaxToward(0.0), ptMaxTrans1(0.0), ptMaxTrans2(0.0), ptMaxAway(0.0);
136 const double phiZ = dilepton.azimuthalAngle();
137 const double pTZ = dilepton.pT();
138 /// @todo Replace with for
139 for (Particles::const_iterator p = fs.particles().begin(); p != fs.particles().end(); ++p) {
140 // Don't use the leptons
141 /// @todo Replace with PID::isLepton
142 if (abs(p->pid()) < 20) continue;
143
144 const double dPhi = deltaPhi(p->momentum().phi(), phiZ);
145 const double pT = p->pT();
146 double rotatedphi = p->momentum().phi() - phiZ;
147 while (rotatedphi < 0) rotatedphi += 2*PI;
148
149 if (dPhi < PI/3.0) {
150 ptSumToward += pT;
151 ++numToward;
152 if (pT > ptMaxToward)
153 ptMaxToward = pT;
154 } else if (dPhi < 2*PI/3.0) {
155 if (rotatedphi <= PI) {
156 ptSumTrans1 += pT;
157 ++numTrans1;
158 if (pT > ptMaxTrans1)
159 ptMaxTrans1 = pT;
160 }
161 else {
162 ptSumTrans2 += pT;
163 ++numTrans2;
164 if (pT > ptMaxTrans2)
165 ptMaxTrans2 = pT;
166 }
167 } else {
168 ptSumAway += pT;
169 ++numAway;
170 if (pT > ptMaxAway)
171 ptMaxAway = pT;
172 }
173 // We need to subtract the two leptons from the number of particles to get the correct multiplicity
174 _p["cptavevsnchg"]->fill(numParticles-2, pT);
175 if (pTZ < 10)
176 _p["cptavevsnchgsmallzpt"]->fill(numParticles-2, pT);
177 }
178
179 // Fill the histograms
180 _p["tnchg"]->fill(pTZ, numToward/(4*PI/3));
181 _p["pnchg"]->fill(pTZ, (numTrans1+numTrans2)/(4*PI/3));
182 _p["pmaxnchg"]->fill(pTZ, (numTrans1>numTrans2 ? numTrans1 : numTrans2)/(2*PI/3));
183 _p["pminnchg"]->fill(pTZ, (numTrans1<numTrans2 ? numTrans1 : numTrans2)/(2*PI/3));
184 _p["pdifnchg"]->fill(pTZ, abs(numTrans1-numTrans2)/(2*PI/3));
185 _p["anchg"]->fill(pTZ, numAway/(4*PI/3));
186
187 _p["tcptsum"]->fill(pTZ, ptSumToward/(4*PI/3));
188 _p["pcptsum"]->fill(pTZ, (ptSumTrans1+ptSumTrans2)/(4*PI/3));
189 _p["pmaxcptsum"]->fill(pTZ, (ptSumTrans1>ptSumTrans2 ? ptSumTrans1 : ptSumTrans2)/(2*PI/3));
190 _p["pmincptsum"]->fill(pTZ, (ptSumTrans1<ptSumTrans2 ? ptSumTrans1 : ptSumTrans2)/(2*PI/3));
191 _p["pdifcptsum"]->fill(pTZ, fabs(ptSumTrans1-ptSumTrans2)/(2*PI/3));
192 _p["acptsum"]->fill(pTZ, ptSumAway/(4*PI/3));
193
194 if (numToward > 0) {
195 _p["tcptave"]->fill(pTZ, ptSumToward/numToward);
196 _p["tcptmax"]->fill(pTZ, ptMaxToward);
197 }
198 if ((numTrans1+numTrans2) > 0) {
199 _p["pcptave"]->fill(pTZ, (ptSumTrans1+ptSumTrans2)/(numTrans1+numTrans2));
200 _p["pcptmax"]->fill(pTZ, (ptMaxTrans1 > ptMaxTrans2 ? ptMaxTrans1 : ptMaxTrans2));
201 }
202
203 // We need to subtract the two leptons from the number of particles to get the correct multiplicity
204 _p["zptvsnchg"]->fill(numParticles-2, pTZ);
205 }
206
207
208 void doQCDanalysis(const Event& e) {
209
210 const FinalState& fsj = apply<FinalState>(e, "FSJ");
211 if (fsj.particles().size() < 1) {
212 MSG_DEBUG("Failed multiplicity cut");
213 vetoEvent;
214 }
215
216 const Jets& jets = apply<FastJets>(e, "MidpointJets").jetsByPt();
217 MSG_DEBUG("Jet multiplicity = " << jets.size());
218
219 // We require the leading jet to be within |eta|<2
220 if (jets.size() < 1 || fabs(jets[0].eta()) >= 2) {
221 MSG_DEBUG("Failed leading jet cut");
222 vetoEvent;
223 }
224
225 const double jetphi = jets[0].phi();
226 const double jeteta = jets[0].eta();
227 const double jetpT = jets[0].pT();
228 MSG_DEBUG("Leading jet: pT = " << jetpT
229 << ", eta = " << jeteta << ", phi = " << jetphi);
230
231 // Get the final states to work with for filling the distributions
232 const FinalState& cfs = apply<ChargedFinalState>(e, "CFS");
233
234 size_t numToward(0), numAway(0);
235 long int numTrans1(0), numTrans2(0);
236 double ptSumToward(0.0), ptSumTrans1(0.0), ptSumTrans2(0.0), ptSumAway(0.0);
237 double ptMaxOverall(0.0), ptMaxToward(0.0), ptMaxTrans1(0.0), ptMaxTrans2(0.0), ptMaxAway(0.0);
238
239 // Calculate all the charged stuff
240 for (const Particle& p : cfs.particles()) {
241 const double dPhi = deltaPhi(p.phi(), jetphi);
242 const double pT = p.pT();
243 const double phi = p.phi();
244 double rotatedphi = phi - jetphi;
245 while (rotatedphi < 0) rotatedphi += 2*PI;
246
247 if (pT > ptMaxOverall) {
248 ptMaxOverall = pT;
249 }
250
251 if (dPhi < PI/3.0) {
252 ptSumToward += pT;
253 ++numToward;
254 if (pT > ptMaxToward) ptMaxToward = pT;
255 }
256 else if (dPhi < 2*PI/3.0) {
257 if (rotatedphi <= PI) {
258 ptSumTrans1 += pT;
259 ++numTrans1;
260 if (pT > ptMaxTrans1) ptMaxTrans1 = pT;
261 } else {
262 ptSumTrans2 += pT;
263 ++numTrans2;
264 if (pT > ptMaxTrans2) ptMaxTrans2 = pT;
265 }
266 }
267 else {
268 ptSumAway += pT;
269 ++numAway;
270 if (pT > ptMaxAway) ptMaxAway = pT;
271 }
272 } // end charged particle loop
273
274 // Fill the histograms
275 _p["tnchg"]->fill(jetpT/GeV, numToward/(4*PI/3));
276 _p["pnchg"]->fill(jetpT/GeV, (numTrans1+numTrans2)/(4*PI/3));
277 _p["pmaxnchg"]->fill(jetpT/GeV, (numTrans1>numTrans2 ? numTrans1 : numTrans2)/(2*PI/3));
278 _p["pminnchg"]->fill(jetpT/GeV, (numTrans1<numTrans2 ? numTrans1 : numTrans2)/(2*PI/3));
279 _p["pdifnchg"]->fill(jetpT/GeV, abs(numTrans1-numTrans2)/(2*PI/3));
280 _p["anchg"]->fill(jetpT/GeV, numAway/(4*PI/3));
281
282 _p["tcptsum"]->fill(jetpT/GeV, ptSumToward/GeV/(4*PI/3));
283 _p["pcptsum"]->fill(jetpT/GeV, (ptSumTrans1+ptSumTrans2)/GeV/(4*PI/3));
284 _p["pmaxcptsum"]->fill(jetpT/GeV, (ptSumTrans1>ptSumTrans2 ? ptSumTrans1 : ptSumTrans2)/GeV/(2*PI/3));
285 _p["pmincptsum"]->fill(jetpT/GeV, (ptSumTrans1<ptSumTrans2 ? ptSumTrans1 : ptSumTrans2)/GeV/(2*PI/3));
286 _p["pdifcptsum"]->fill(jetpT/GeV, fabs(ptSumTrans1-ptSumTrans2)/GeV/(2*PI/3));
287 _p["acptsum"]->fill(jetpT/GeV, ptSumAway/GeV/(4*PI/3));
288
289 if ((numTrans1+numTrans2) > 0) {
290 _p["pcptave"]->fill(jetpT/GeV, (ptSumTrans1+ptSumTrans2)/GeV/(numTrans1+numTrans2));
291 _p["pcptmax"]->fill(jetpT/GeV, (ptMaxTrans1 > ptMaxTrans2 ? ptMaxTrans1 : ptMaxTrans2)/GeV);
292 }
293 }
294
295
296 // void finalize() { }
297
298 /// @}
299
300
301 private:
302
303 size_t _mode;
304 map<string, Profile1DPtr> _p;
305
306 };
307
308
309
310 RIVET_DECLARE_ALIASED_PLUGIN(CDF_2010_I849042, CDF_2010_S8591881);
311
312}
|