1 | // Analysis.cc is a part of the PYTHIA event generator. |
---|
2 | // Copyright (C) 2012 Torbjorn Sjostrand. |
---|
3 | // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details. |
---|
4 | // Please respect the MCnet Guidelines, see GUIDELINES for details. |
---|
5 | |
---|
6 | // Function definitions (not found in the header) for the |
---|
7 | // Sphericity, Thrust, ClusJet, CellJet and SlowJet classes. |
---|
8 | |
---|
9 | #include "Analysis.h" |
---|
10 | |
---|
11 | namespace Pythia8 { |
---|
12 | |
---|
13 | //========================================================================== |
---|
14 | |
---|
15 | // Sphericity class. |
---|
16 | // This class finds sphericity-related properties of an event. |
---|
17 | |
---|
18 | //-------------------------------------------------------------------------- |
---|
19 | |
---|
20 | // Constants: could be changed here if desired, but normally should not. |
---|
21 | // These are of technical nature, as described for each. |
---|
22 | |
---|
23 | // Minimum number of particles to perform study. |
---|
24 | const int Sphericity::NSTUDYMIN = 2; |
---|
25 | |
---|
26 | // Maximum number of times that an error warning will be printed. |
---|
27 | const int Sphericity::TIMESTOPRINT = 1; |
---|
28 | |
---|
29 | // Assign mimimum squared momentum in weight to avoid division by zero. |
---|
30 | const double Sphericity::P2MIN = 1e-20; |
---|
31 | |
---|
32 | // Second eigenvalue not too low or not possible to find eigenvectors. |
---|
33 | const double Sphericity::EIGENVALUEMIN = 1e-10; |
---|
34 | |
---|
35 | //-------------------------------------------------------------------------- |
---|
36 | |
---|
37 | // Analyze event. |
---|
38 | |
---|
39 | bool Sphericity::analyze(const Event& event, ostream& os) { |
---|
40 | |
---|
41 | // Initial values, tensor and counters zero. |
---|
42 | eVal1 = eVal2 = eVal3 = 0.; |
---|
43 | eVec1 = eVec2 = eVec3 = 0.; |
---|
44 | double tt[4][4]; |
---|
45 | for (int j = 1; j < 4; ++j) |
---|
46 | for (int k = j; k < 4; ++k) tt[j][k] = 0.; |
---|
47 | int nStudy = 0; |
---|
48 | double denom = 0.; |
---|
49 | |
---|
50 | // Loop over desired particles in the event. |
---|
51 | for (int i = 0; i < event.size(); ++i) |
---|
52 | if (event[i].isFinal()) { |
---|
53 | if (select > 2 && event[i].isNeutral() ) continue; |
---|
54 | if (select == 2 && !event[i].isVisible() ) continue; |
---|
55 | ++nStudy; |
---|
56 | |
---|
57 | // Calculate matrix to be diagonalized. Special cases for speed. |
---|
58 | double pNow[4]; |
---|
59 | pNow[1] = event[i].px(); |
---|
60 | pNow[2] = event[i].py(); |
---|
61 | pNow[3] = event[i].pz(); |
---|
62 | double p2Now = pNow[1]*pNow[1] + pNow[2]*pNow[2] + pNow[3]*pNow[3]; |
---|
63 | double pWeight = 1.; |
---|
64 | if (powerInt == 1) pWeight = 1. / sqrt(max(P2MIN, p2Now)); |
---|
65 | else if (powerInt == 0) pWeight = pow( max(P2MIN, p2Now), powerMod); |
---|
66 | for (int j = 1; j < 4; ++j) |
---|
67 | for (int k = j; k < 4; ++k) tt[j][k] += pWeight * pNow[j] * pNow[k]; |
---|
68 | denom += pWeight * p2Now; |
---|
69 | } |
---|
70 | |
---|
71 | // Very low multiplicities (0 or 1) not considered. |
---|
72 | if (nStudy < NSTUDYMIN) { |
---|
73 | if (nFew < TIMESTOPRINT) os << " PYTHIA Error in " << |
---|
74 | "Sphericity::analyze: too few particles" << endl; |
---|
75 | ++nFew; |
---|
76 | return false; |
---|
77 | } |
---|
78 | |
---|
79 | // Normalize tensor to trace = 1. |
---|
80 | for (int j = 1; j < 4; ++j) |
---|
81 | for (int k = j; k < 4; ++k) tt[j][k] /= denom; |
---|
82 | |
---|
83 | // Find eigenvalues to matrix (third degree equation). |
---|
84 | double qCoef = ( tt[1][1] * tt[2][2] + tt[1][1] * tt[3][3] |
---|
85 | + tt[2][2] * tt[3][3] - pow2(tt[1][2]) - pow2(tt[1][3]) |
---|
86 | - pow2(tt[2][3]) ) / 3. - 1./9.; |
---|
87 | double qCoefRt = sqrt( -qCoef); |
---|
88 | double rCoef = -0.5 * ( qCoef + 1./9. + tt[1][1] * pow2(tt[2][3]) |
---|
89 | + tt[2][2] * pow2(tt[1][3]) + tt[3][3] * pow2(tt[1][2]) |
---|
90 | - tt[1][1] * tt[2][2] * tt[3][3] ) |
---|
91 | + tt[1][2] * tt[1][3] * tt[2][3] + 1./27.; |
---|
92 | double pTemp = max( min( rCoef / pow3(qCoefRt), 1.), -1.); |
---|
93 | double pCoef = cos( acos(pTemp) / 3.); |
---|
94 | double pCoefRt = sqrt( 3. * (1. - pow2(pCoef)) ); |
---|
95 | eVal1 = 1./3. + qCoefRt * max( 2. * pCoef, pCoefRt - pCoef); |
---|
96 | eVal3 = 1./3. + qCoefRt * min( 2. * pCoef, -pCoefRt - pCoef); |
---|
97 | eVal2 = 1. - eVal1 - eVal3; |
---|
98 | |
---|
99 | // Begin find first and last eigenvector. |
---|
100 | for (int iVal = 0; iVal < 2; ++iVal) { |
---|
101 | double eVal = (iVal == 0) ? eVal1 : eVal3; |
---|
102 | |
---|
103 | // If all particles are back-to-back then only first axis meaningful. |
---|
104 | if (iVal > 1 && eVal2 < EIGENVALUEMIN) { |
---|
105 | if (nBack < TIMESTOPRINT) os << " PYTHIA Error in " |
---|
106 | "Sphericity::analyze: particles too back-to-back" << endl; |
---|
107 | ++nBack; |
---|
108 | return false; |
---|
109 | } |
---|
110 | |
---|
111 | // Set up matrix to diagonalize. |
---|
112 | double dd[4][4]; |
---|
113 | for (int j = 1; j < 4; ++j) { |
---|
114 | dd[j][j] = tt[j][j] - eVal; |
---|
115 | for (int k = j + 1; k < 4; ++k) { |
---|
116 | dd[j][k] = tt[j][k]; |
---|
117 | dd[k][j] = tt[j][k]; |
---|
118 | } |
---|
119 | } |
---|
120 | |
---|
121 | // Find largest = pivotal element in matrix. |
---|
122 | int jMax = 0; |
---|
123 | int kMax = 0; |
---|
124 | double ddMax = 0.; |
---|
125 | for (int j = 1; j < 4; ++j) |
---|
126 | for (int k = 1; k < 4; ++k) |
---|
127 | if (abs(dd[j][k]) > ddMax) { |
---|
128 | jMax = j; |
---|
129 | kMax = k; |
---|
130 | ddMax = abs(dd[j][k]); |
---|
131 | } |
---|
132 | |
---|
133 | // Subtract one row from the other two; find new largest element. |
---|
134 | int jMax2 = 0; |
---|
135 | ddMax = 0.; |
---|
136 | for (int j = 1; j < 4; ++j) |
---|
137 | if ( j != jMax) { |
---|
138 | double pivot = dd[j][kMax] / dd[jMax][kMax]; |
---|
139 | for (int k = 1; k < 4; ++k) { |
---|
140 | dd[j][k] -= pivot * dd[jMax][k]; |
---|
141 | if (abs(dd[j][k]) > ddMax) { |
---|
142 | jMax2 = j; |
---|
143 | ddMax = abs(dd[j][k]); |
---|
144 | } |
---|
145 | } |
---|
146 | } |
---|
147 | |
---|
148 | // Construct eigenvector. Normalize to unit length; sign irrelevant. |
---|
149 | int k1 = kMax + 1; if (k1 > 3) k1 -= 3; |
---|
150 | int k2 = kMax + 2; if (k2 > 3) k2 -= 3; |
---|
151 | double eVec[4]; |
---|
152 | eVec[k1] = -dd[jMax2][k2]; |
---|
153 | eVec[k2] = dd[jMax2][k1]; |
---|
154 | eVec[kMax] = (dd[jMax][k1] * dd[jMax2][k2] |
---|
155 | - dd[jMax][k2] * dd[jMax2][k1]) / dd[jMax][kMax]; |
---|
156 | double length = sqrt( pow2(eVec[1]) + pow2(eVec[2]) |
---|
157 | + pow2(eVec[3]) ); |
---|
158 | |
---|
159 | // Store eigenvectors. |
---|
160 | if (iVal == 0) eVec1 = Vec4( eVec[1] / length, |
---|
161 | eVec[2] / length, eVec[3] / length, 0.); |
---|
162 | else eVec3 = Vec4( eVec[1] / length, |
---|
163 | eVec[2] / length, eVec[3] / length, 0.); |
---|
164 | } |
---|
165 | |
---|
166 | // Middle eigenvector is orthogonal to the other two; sign irrelevant. |
---|
167 | eVec2 = cross3( eVec1, eVec3); |
---|
168 | |
---|
169 | // Done. |
---|
170 | return true; |
---|
171 | |
---|
172 | } |
---|
173 | |
---|
174 | //-------------------------------------------------------------------------- |
---|
175 | |
---|
176 | // Provide a listing of the info. |
---|
177 | |
---|
178 | void Sphericity::list(ostream& os) const { |
---|
179 | |
---|
180 | // Header. |
---|
181 | os << "\n -------- PYTHIA Sphericity Listing -------- \n"; |
---|
182 | if (powerInt !=2) os << " Nonstandard momentum power = " |
---|
183 | << fixed << setprecision(3) << setw(6) << power << "\n"; |
---|
184 | os << "\n no lambda e_x e_y e_z \n"; |
---|
185 | |
---|
186 | // The three eigenvalues and eigenvectors. |
---|
187 | os << setprecision(5); |
---|
188 | os << " 1" << setw(11) << eVal1 << setw(11) << eVec1.px() |
---|
189 | << setw(10) << eVec1.py() << setw(10) << eVec1.pz() << "\n"; |
---|
190 | os << " 2" << setw(11) << eVal2 << setw(11) << eVec2.px() |
---|
191 | << setw(10) << eVec2.py() << setw(10) << eVec2.pz() << "\n"; |
---|
192 | os << " 3" << setw(11) << eVal3 << setw(11) << eVec3.px() |
---|
193 | << setw(10) << eVec3.py() << setw(10) << eVec3.pz() << "\n"; |
---|
194 | |
---|
195 | // Listing finished. |
---|
196 | os << "\n -------- End PYTHIA Sphericity Listing ----" << endl; |
---|
197 | |
---|
198 | } |
---|
199 | |
---|
200 | |
---|
201 | //========================================================================== |
---|
202 | |
---|
203 | // Thrust class. |
---|
204 | // This class finds thrust-related properties of an event. |
---|
205 | |
---|
206 | //-------------------------------------------------------------------------- |
---|
207 | |
---|
208 | // Constants: could be changed here if desired, but normally should not. |
---|
209 | // These are of technical nature, as described for each. |
---|
210 | |
---|
211 | // Minimum number of particles to perform study. |
---|
212 | const int Thrust::NSTUDYMIN = 2; |
---|
213 | |
---|
214 | // Maximum number of times that an error warning will be printed. |
---|
215 | const int Thrust::TIMESTOPRINT = 1; |
---|
216 | |
---|
217 | // Major not too low or not possible to find major axis. |
---|
218 | const double Thrust::MAJORMIN = 1e-10; |
---|
219 | |
---|
220 | //-------------------------------------------------------------------------- |
---|
221 | |
---|
222 | // Analyze event. |
---|
223 | |
---|
224 | bool Thrust::analyze(const Event& event, ostream& os) { |
---|
225 | |
---|
226 | // Initial values and counters zero. |
---|
227 | eVal1 = eVal2 = eVal3 = 0.; |
---|
228 | eVec1 = eVec2 = eVec3 = 0.; |
---|
229 | int nStudy = 0; |
---|
230 | vector<Vec4> pOrder; |
---|
231 | Vec4 pSum, nRef, pPart, pFull, pMax; |
---|
232 | |
---|
233 | // Loop over desired particles in the event. |
---|
234 | for (int i = 0; i < event.size(); ++i) |
---|
235 | if (event[i].isFinal()) { |
---|
236 | if (select > 2 && event[i].isNeutral() ) continue; |
---|
237 | if (select == 2 && !event[i].isVisible() ) continue; |
---|
238 | ++nStudy; |
---|
239 | |
---|
240 | // Store momenta. Use energy component for absolute momentum. |
---|
241 | Vec4 pNow = event[i].p(); |
---|
242 | pNow.e(pNow.pAbs()); |
---|
243 | pSum += pNow; |
---|
244 | pOrder.push_back(pNow); |
---|
245 | } |
---|
246 | |
---|
247 | // Very low multiplicities (0 or 1) not considered. |
---|
248 | if (nStudy < NSTUDYMIN) { |
---|
249 | if (nFew < TIMESTOPRINT) os << " PYTHIA Error in " << |
---|
250 | "Thrust::analyze: too few particles" << endl; |
---|
251 | ++nFew; |
---|
252 | return false; |
---|
253 | } |
---|
254 | |
---|
255 | // Try all combinations of reference vector orthogonal to two particles. |
---|
256 | for (int i1 = 0; i1 < nStudy - 1; ++i1) |
---|
257 | for (int i2 = i1 + 1; i2 < nStudy; ++i2) { |
---|
258 | nRef = cross3( pOrder[i1], pOrder[i2]); |
---|
259 | nRef /= nRef.pAbs(); |
---|
260 | pPart = 0.; |
---|
261 | |
---|
262 | // Add all momenta with sign; two choices for each reference particle. |
---|
263 | for (int i = 0; i < nStudy; ++i) if (i != i1 && i != i2) { |
---|
264 | if (dot3(pOrder[i], nRef) > 0.) pPart += pOrder[i]; |
---|
265 | else pPart -= pOrder[i]; |
---|
266 | } |
---|
267 | for (int j = 0; j < 4; ++j) { |
---|
268 | if (j == 0) pFull = pPart + pOrder[i1] + pOrder[i2]; |
---|
269 | else if (j == 1) pFull = pPart + pOrder[i1] - pOrder[i2]; |
---|
270 | else if (j == 2) pFull = pPart - pOrder[i1] + pOrder[i2]; |
---|
271 | else pFull = pPart - pOrder[i1] - pOrder[i2]; |
---|
272 | pFull.e(pFull.pAbs()); |
---|
273 | if (pFull.e() > pMax.e()) pMax = pFull; |
---|
274 | } |
---|
275 | } |
---|
276 | |
---|
277 | // Maximum gives thrust axis and value. |
---|
278 | eVal1 = pMax.e() / pSum.e(); |
---|
279 | eVec1 = pMax / pMax.e(); |
---|
280 | eVec1.e(0.); |
---|
281 | |
---|
282 | // Subtract momentum along thrust axis. |
---|
283 | double pAbsSum = 0.; |
---|
284 | for (int i = 0; i < nStudy; ++i) { |
---|
285 | pOrder[i] -= dot3( eVec1, pOrder[i]) * eVec1; |
---|
286 | pOrder[i].e(pOrder[i].pAbs()); |
---|
287 | pAbsSum += pOrder[i].e(); |
---|
288 | } |
---|
289 | |
---|
290 | // Simpleminded major and minor axes if too little transverse left. |
---|
291 | if (pAbsSum < MAJORMIN * pSum.e()) { |
---|
292 | if ( abs(eVec1.pz()) > 0.5) eVec2 = Vec4( 1., 0., 0., 0.); |
---|
293 | else eVec2 = Vec4( 0., 0., 1., 0.); |
---|
294 | eVec2 -= dot3( eVec1, eVec2) * eVec1; |
---|
295 | eVec2 /= eVec2.pAbs(); |
---|
296 | eVec3 = cross3( eVec1, eVec2); |
---|
297 | return true; |
---|
298 | } |
---|
299 | |
---|
300 | // Try all reference vectors orthogonal to one particles. |
---|
301 | pMax = 0.; |
---|
302 | for (int i1 = 0; i1 < nStudy; ++i1) { |
---|
303 | nRef = cross3( pOrder[i1], eVec1); |
---|
304 | nRef /= nRef.pAbs(); |
---|
305 | pPart = 0.; |
---|
306 | |
---|
307 | // Add all momenta with sign; two choices for each reference particle. |
---|
308 | for (int i = 0; i < nStudy; ++i) if (i != i1) { |
---|
309 | if (dot3(pOrder[i], nRef) > 0.) pPart += pOrder[i]; |
---|
310 | else pPart -= pOrder[i]; |
---|
311 | } |
---|
312 | pFull = pPart + pOrder[i1]; |
---|
313 | pFull.e(pFull.pAbs()); |
---|
314 | if (pFull.e() > pMax.e()) pMax = pFull; |
---|
315 | pFull = pPart - pOrder[i1]; |
---|
316 | pFull.e(pFull.pAbs()); |
---|
317 | if (pFull.e() > pMax.e()) pMax = pFull; |
---|
318 | } |
---|
319 | |
---|
320 | // Maximum gives major axis and value. |
---|
321 | eVal2 = pMax.e() / pSum.e(); |
---|
322 | eVec2 = pMax / pMax.e(); |
---|
323 | eVec2.e(0.); |
---|
324 | |
---|
325 | // Orthogonal direction gives minor axis, and from there value. |
---|
326 | eVec3 = cross3( eVec1, eVec2); |
---|
327 | pAbsSum = 0.; |
---|
328 | for (int i = 0; i < nStudy; ++i) |
---|
329 | pAbsSum += abs( dot3(eVec3, pOrder[i]) ); |
---|
330 | eVal3 = pAbsSum / pSum.e(); |
---|
331 | |
---|
332 | // Done. |
---|
333 | return true; |
---|
334 | |
---|
335 | } |
---|
336 | |
---|
337 | //-------------------------------------------------------------------------- |
---|
338 | |
---|
339 | // Provide a listing of the info. |
---|
340 | |
---|
341 | void Thrust::list(ostream& os) const { |
---|
342 | |
---|
343 | // Header. |
---|
344 | os << "\n -------- PYTHIA Thrust Listing ------------ \n" |
---|
345 | << "\n value e_x e_y e_z \n"; |
---|
346 | |
---|
347 | // The thrust, major and minor values and related event axes. |
---|
348 | os << setprecision(5); |
---|
349 | os << " Thr" << setw(11) << eVal1 << setw(11) << eVec1.px() |
---|
350 | << setw(10) << eVec1.py() << setw(10) << eVec1.pz() << "\n"; |
---|
351 | os << " Maj" << setw(11) << eVal2 << setw(11) << eVec2.px() |
---|
352 | << setw(10) << eVec2.py() << setw(10) << eVec2.pz() << "\n"; |
---|
353 | os << " Min" << setw(11) << eVal3 << setw(11) << eVec3.px() |
---|
354 | << setw(10) << eVec3.py() << setw(10) << eVec3.pz() << "\n"; |
---|
355 | |
---|
356 | // Listing finished. |
---|
357 | os << "\n -------- End PYTHIA Thrust Listing --------" << endl; |
---|
358 | |
---|
359 | } |
---|
360 | |
---|
361 | //========================================================================== |
---|
362 | |
---|
363 | // SingleClusterJet class. |
---|
364 | // Simple helper class to ClusterJet for a jet and its contents. |
---|
365 | |
---|
366 | //-------------------------------------------------------------------------- |
---|
367 | |
---|
368 | // Constants: could be changed here if desired, but normally should not. |
---|
369 | // These are of technical nature, as described for each. |
---|
370 | |
---|
371 | // Assign minimal pAbs to avoid division by zero. |
---|
372 | const double SingleClusterJet::PABSMIN = 1e-10; |
---|
373 | |
---|
374 | //-------------------------------------------------------------------------- |
---|
375 | |
---|
376 | // Distance measures between two SingleClusterJet objects. |
---|
377 | |
---|
378 | double dist2Fun(int measure, const SingleClusterJet& j1, |
---|
379 | const SingleClusterJet& j2) { |
---|
380 | |
---|
381 | // JADE distance. |
---|
382 | if (measure == 2) return 2. * j1.pJet.e() * j2.pJet.e() |
---|
383 | * (1. - dot3( j1.pJet, j2.pJet) / (j1.pAbs * j2.pAbs) ); |
---|
384 | |
---|
385 | // Durham distance. |
---|
386 | if (measure == 3) return 2. * pow2( min( j1.pJet.e(), j2.pJet.e() ) ) |
---|
387 | * (1. - dot3( j1.pJet, j2.pJet) / (j1.pAbs * j2.pAbs) ); |
---|
388 | |
---|
389 | // Lund distance; "default". |
---|
390 | return (j1.pAbs * j2.pAbs - dot3( j1.pJet, j2.pJet)) |
---|
391 | * 2. * j1.pAbs * j2.pAbs / pow2(j1.pAbs + j2.pAbs); |
---|
392 | |
---|
393 | } |
---|
394 | |
---|
395 | //========================================================================== |
---|
396 | |
---|
397 | // ClusterJet class. |
---|
398 | // This class performs a jet clustering according to different |
---|
399 | // distance measures: Lund, JADE or Durham. |
---|
400 | |
---|
401 | //-------------------------------------------------------------------------- |
---|
402 | |
---|
403 | // Constants: could be changed here if desired, but normally should not. |
---|
404 | // These are of technical nature, as described for each. |
---|
405 | |
---|
406 | // Maximum number of times that an error warning will be printed. |
---|
407 | const int ClusterJet::TIMESTOPRINT = 1; |
---|
408 | |
---|
409 | // Assume the pi+- mass for all particles, except the photon, in one option. |
---|
410 | const double ClusterJet::PIMASS = 0.13957; |
---|
411 | |
---|
412 | // Assign minimal pAbs to avoid division by zero. |
---|
413 | const double ClusterJet::PABSMIN = 1e-10; |
---|
414 | |
---|
415 | // Initial pT/m preclustering scale as fraction of clustering one. |
---|
416 | const double ClusterJet::PRECLUSTERFRAC = 0.1; |
---|
417 | |
---|
418 | // Step with which pT/m is reduced if preclustering gives too few jets. |
---|
419 | const double ClusterJet::PRECLUSTERSTEP = 0.8; |
---|
420 | |
---|
421 | //-------------------------------------------------------------------------- |
---|
422 | |
---|
423 | // Analyze event. |
---|
424 | |
---|
425 | bool ClusterJet::analyze(const Event& event, double yScaleIn, |
---|
426 | double pTscaleIn, int nJetMinIn, int nJetMaxIn, ostream& os) { |
---|
427 | |
---|
428 | // Input values. Initial values zero. |
---|
429 | yScale = yScaleIn; |
---|
430 | pTscale = pTscaleIn; |
---|
431 | nJetMin = nJetMinIn; |
---|
432 | nJetMax = nJetMaxIn; |
---|
433 | particles.resize(0); |
---|
434 | jets.resize(0); |
---|
435 | Vec4 pSum; |
---|
436 | distances.clear(); |
---|
437 | |
---|
438 | // Loop over desired particles in the event. |
---|
439 | for (int i = 0; i < event.size(); ++i) |
---|
440 | if (event[i].isFinal()) { |
---|
441 | if (select > 2 && event[i].isNeutral() ) continue; |
---|
442 | if (select == 2 && !event[i].isVisible() ) continue; |
---|
443 | |
---|
444 | // Store them, possibly with modified mass => new energy. |
---|
445 | Vec4 pTemp = event[i].p(); |
---|
446 | if (massSet == 0 || massSet == 1) { |
---|
447 | double mTemp = (massSet == 0 || event[i].id() == 22) |
---|
448 | ? 0. : PIMASS; |
---|
449 | double eTemp = sqrt(pTemp.pAbs2() + pow2(mTemp)); |
---|
450 | pTemp.e(eTemp); |
---|
451 | } |
---|
452 | particles.push_back( SingleClusterJet(pTemp, i) ); |
---|
453 | pSum += pTemp; |
---|
454 | } |
---|
455 | |
---|
456 | // Very low multiplicities not considered. |
---|
457 | nParticles = particles.size(); |
---|
458 | if (nParticles < nJetMin) { |
---|
459 | if (nFew < TIMESTOPRINT) os << " PYTHIA Error in " << |
---|
460 | "ClusterJet::analyze: too few particles" << endl; |
---|
461 | ++nFew; |
---|
462 | return false; |
---|
463 | } |
---|
464 | |
---|
465 | // Squared maximum distance in GeV^2 for joining. |
---|
466 | double p2Sum = pSum.m2Calc(); |
---|
467 | dist2Join = max( yScale * p2Sum, pow2(pTscale)); |
---|
468 | dist2BigMin = 2. * max( dist2Join, p2Sum); |
---|
469 | |
---|
470 | // Do preclustering if desired and possible. |
---|
471 | if (doPrecluster && nParticles > nJetMin + 2) { |
---|
472 | precluster(); |
---|
473 | if (doReassign) reassign(); |
---|
474 | } |
---|
475 | |
---|
476 | // If no preclustering: each particle is a starting jet. |
---|
477 | else for (int i = 0; i < nParticles; ++i) { |
---|
478 | jets.push_back( SingleClusterJet(particles[i]) ); |
---|
479 | particles[i].daughter = i; |
---|
480 | } |
---|
481 | |
---|
482 | // Begin iteration towards fewer jets. |
---|
483 | for ( ; ; ) { |
---|
484 | |
---|
485 | // Find the two closest jets. |
---|
486 | double dist2Min = dist2BigMin; |
---|
487 | int jMin = 0; |
---|
488 | int kMin = 0; |
---|
489 | for (int j = 0; j < int(jets.size()) - 1; ++j) |
---|
490 | for (int k = j + 1; k < int(jets.size()); ++k) { |
---|
491 | double dist2 = dist2Fun( measure, jets[j], jets[k]); |
---|
492 | if (dist2 < dist2Min) { |
---|
493 | dist2Min = dist2; |
---|
494 | jMin = j; |
---|
495 | kMin = k; |
---|
496 | } |
---|
497 | } |
---|
498 | |
---|
499 | // Stop if no pair below cut and not more jets than allowed. |
---|
500 | if ( dist2Min > dist2Join |
---|
501 | && (nJetMax < nJetMin || int(jets.size()) <= nJetMax) ) break; |
---|
502 | |
---|
503 | // Stop if reached minimum allowed number of jets. Else continue. |
---|
504 | if (int(jets.size()) <= nJetMin) break; |
---|
505 | |
---|
506 | // Join two closest jets. |
---|
507 | jets[jMin].pJet += jets[kMin].pJet; |
---|
508 | jets[jMin].pAbs = max( PABSMIN, jets[jMin].pJet.pAbs()); |
---|
509 | jets[jMin].multiplicity += jets[kMin].multiplicity; |
---|
510 | for (int i = 0; i < nParticles; ++i) |
---|
511 | if (particles[i].daughter == kMin) particles[i].daughter = jMin; |
---|
512 | |
---|
513 | // Save the last 5 distances. |
---|
514 | distances.push_front(dist2Min); |
---|
515 | if (distances.size() > 5) distances.pop_back(); |
---|
516 | |
---|
517 | // Move up last jet to empty slot to shrink list. |
---|
518 | jets[kMin] = jets.back(); |
---|
519 | jets.pop_back(); |
---|
520 | int iEnd = jets.size(); |
---|
521 | for (int i = 0; i < nParticles; ++i) |
---|
522 | if (particles[i].daughter == iEnd) particles[i].daughter = kMin; |
---|
523 | |
---|
524 | // Do reassignments of particles to nearest jet if desired. |
---|
525 | if (doReassign) reassign(); |
---|
526 | } |
---|
527 | |
---|
528 | // Order jets in decreasing energy. |
---|
529 | for (int j = 0; j < int(jets.size()) - 1; ++j) |
---|
530 | for (int k = int(jets.size()) - 1; k > j; --k) |
---|
531 | if (jets[k].pJet.e() > jets[k-1].pJet.e()) { |
---|
532 | swap( jets[k], jets[k-1]); |
---|
533 | for (int i = 0; i < nParticles; ++i) { |
---|
534 | if (particles[i].daughter == k) particles[i].daughter = k-1; |
---|
535 | else if (particles[i].daughter == k-1) particles[i].daughter = k; |
---|
536 | } |
---|
537 | } |
---|
538 | |
---|
539 | // Done. |
---|
540 | return true; |
---|
541 | } |
---|
542 | |
---|
543 | //-------------------------------------------------------------------------- |
---|
544 | |
---|
545 | // Precluster nearby particles to save computer time. |
---|
546 | |
---|
547 | void ClusterJet::precluster() { |
---|
548 | |
---|
549 | // Begin iteration over preclustering scale. |
---|
550 | distPre = PRECLUSTERFRAC * sqrt(dist2Join) / PRECLUSTERSTEP; |
---|
551 | for ( ; ;) { |
---|
552 | distPre *= PRECLUSTERSTEP; |
---|
553 | dist2Pre = pow2(distPre); |
---|
554 | for (int i = 0; i < nParticles; ++i) { |
---|
555 | particles[i].daughter = -1; |
---|
556 | particles[i].isAssigned = false; |
---|
557 | } |
---|
558 | |
---|
559 | // Sum up low-momentum region. Jet if enough momentum. |
---|
560 | Vec4 pCentral; |
---|
561 | int multCentral = 0; |
---|
562 | for (int i = 0; i < nParticles; ++i) |
---|
563 | if (particles[i].pAbs < 2. * distPre) { |
---|
564 | pCentral += particles[i].pJet; |
---|
565 | multCentral += particles[i].multiplicity; |
---|
566 | particles[i].isAssigned = true; |
---|
567 | } |
---|
568 | if (pCentral.pAbs() > 2. * distPre) { |
---|
569 | jets.push_back( SingleClusterJet(pCentral) ); |
---|
570 | jets.back().multiplicity = multCentral; |
---|
571 | for (int i = 0; i < nParticles; ++i) |
---|
572 | if (particles[i].isAssigned) particles[i].daughter = 0; |
---|
573 | } |
---|
574 | |
---|
575 | // Find fastest remaining particle until none left. |
---|
576 | for ( ; ;) { |
---|
577 | int iMax = -1; |
---|
578 | double pMax = 0.; |
---|
579 | for (int i = 0; i < nParticles; ++i) |
---|
580 | if ( !particles[i].isAssigned && particles[i].pAbs > pMax) { |
---|
581 | iMax = i; |
---|
582 | pMax = particles[i].pAbs; |
---|
583 | } |
---|
584 | if (iMax == -1) break; |
---|
585 | |
---|
586 | // Sum up precluster around it according to distance function. |
---|
587 | Vec4 pPre; |
---|
588 | int multPre = 0; |
---|
589 | int nRemain = 0; |
---|
590 | for (int i = 0; i < nParticles; ++i) |
---|
591 | if ( !particles[i].isAssigned) { |
---|
592 | double dist2 = dist2Fun( measure, particles[iMax], |
---|
593 | particles[i]); |
---|
594 | if (dist2 < dist2Pre) { |
---|
595 | pPre += particles[i].pJet; |
---|
596 | ++multPre; |
---|
597 | particles[i].isAssigned = true; |
---|
598 | particles[i].daughter = jets.size(); |
---|
599 | } else ++nRemain; |
---|
600 | } |
---|
601 | jets.push_back( SingleClusterJet(pPre) ); |
---|
602 | jets.back().multiplicity = multPre; |
---|
603 | |
---|
604 | // Decide whether sensible starting configuration or iterate. |
---|
605 | if (int(jets.size()) + nRemain < nJetMin) break; |
---|
606 | } |
---|
607 | if (int(jets.size()) >= nJetMin) break; |
---|
608 | } |
---|
609 | |
---|
610 | } |
---|
611 | |
---|
612 | //-------------------------------------------------------------------------- |
---|
613 | |
---|
614 | // Reassign particles to nearest jet to correct misclustering. |
---|
615 | |
---|
616 | void ClusterJet::reassign() { |
---|
617 | |
---|
618 | // Reset clustered momenta. |
---|
619 | for (int j = 0; j < int(jets.size()); ++j) { |
---|
620 | jets[j].pTemp = 0.; |
---|
621 | jets[j].multiplicity = 0; |
---|
622 | } |
---|
623 | |
---|
624 | // Loop through particles to find closest jet. |
---|
625 | for (int i = 0; i < nParticles; ++i) { |
---|
626 | particles[i].daughter = -1; |
---|
627 | double dist2Min = dist2BigMin; |
---|
628 | int jMin = 0; |
---|
629 | for (int j = 0; j < int(jets.size()); ++j) { |
---|
630 | double dist2 = dist2Fun( measure, particles[i], jets[j]); |
---|
631 | if (dist2 < dist2Min) { |
---|
632 | dist2Min = dist2; |
---|
633 | jMin = j; |
---|
634 | } |
---|
635 | } |
---|
636 | jets[jMin].pTemp += particles[i].pJet; |
---|
637 | ++jets[jMin].multiplicity; |
---|
638 | particles[i].daughter = jMin; |
---|
639 | } |
---|
640 | |
---|
641 | // Replace old by new jet momenta. |
---|
642 | for (int j = 0; j < int(jets.size()); ++j) { |
---|
643 | jets[j].pJet = jets[j].pTemp; |
---|
644 | jets[j].pAbs = max( PABSMIN, jets[j].pJet.pAbs()); |
---|
645 | } |
---|
646 | |
---|
647 | // Check that no empty clusters after reassignments. |
---|
648 | for ( ; ; ) { |
---|
649 | |
---|
650 | // If no empty jets then done. |
---|
651 | int jEmpty = -1; |
---|
652 | for (int j = 0; j < int(jets.size()); ++j) |
---|
653 | if (jets[j].multiplicity == 0) jEmpty = j; |
---|
654 | if (jEmpty == -1) return; |
---|
655 | |
---|
656 | // Find particle assigned to jet with largest distance to it. |
---|
657 | int iSplit = -1; |
---|
658 | double dist2Max = 0.; |
---|
659 | for (int i = 0; i < nParticles; ++i) { |
---|
660 | int j = particles[i].daughter; |
---|
661 | double dist2 = dist2Fun( measure, particles[i], jets[j]); |
---|
662 | if (dist2 > dist2Max) { |
---|
663 | iSplit = i; |
---|
664 | dist2Max = dist2; |
---|
665 | } |
---|
666 | } |
---|
667 | |
---|
668 | // Let this particle form new jet and subtract off from existing. |
---|
669 | int jSplit = particles[iSplit].daughter; |
---|
670 | jets[jEmpty] = SingleClusterJet( particles[iSplit].pJet ); |
---|
671 | jets[jSplit].pJet -= particles[iSplit].pJet; |
---|
672 | jets[jSplit].pAbs = max( PABSMIN,jets[jSplit].pJet.pAbs()); |
---|
673 | particles[iSplit].daughter = jEmpty; |
---|
674 | --jets[jSplit].multiplicity; |
---|
675 | } |
---|
676 | |
---|
677 | } |
---|
678 | |
---|
679 | //-------------------------------------------------------------------------- |
---|
680 | |
---|
681 | // Provide a listing of the info. |
---|
682 | |
---|
683 | void ClusterJet::list(ostream& os) const { |
---|
684 | |
---|
685 | // Header. |
---|
686 | string method = (measure == 1) ? "Lund pT" |
---|
687 | : ( (measure == 2) ? "JADE m" : "Durham kT" ) ; |
---|
688 | os << "\n -------- PYTHIA ClusterJet Listing, " << setw(9) << method |
---|
689 | << " =" << fixed << setprecision(3) << setw(7) << sqrt(dist2Join) |
---|
690 | << " GeV --- \n \n no mult p_x p_y p_z " |
---|
691 | << " e m \n"; |
---|
692 | |
---|
693 | // The jets. |
---|
694 | for (int i = 0; i < int(jets.size()); ++i) { |
---|
695 | os << setw(4) << i << setw(6) << jets[i].multiplicity << setw(11) |
---|
696 | << jets[i].pJet.px() << setw(11) << jets[i].pJet.py() |
---|
697 | << setw(11) << jets[i].pJet.pz() << setw(11) |
---|
698 | << jets[i].pJet.e() << setw(11) << jets[i].pJet.mCalc() |
---|
699 | << "\n"; |
---|
700 | } |
---|
701 | |
---|
702 | // Listing finished. |
---|
703 | os << "\n -------- End PYTHIA ClusterJet Listing ---------------" |
---|
704 | << "--------" << endl; |
---|
705 | } |
---|
706 | |
---|
707 | //========================================================================== |
---|
708 | |
---|
709 | // CellJet class. |
---|
710 | // This class performs a cone jet search in (eta, phi, E_T) space. |
---|
711 | |
---|
712 | //-------------------------------------------------------------------------- |
---|
713 | |
---|
714 | // Constants: could be changed here if desired, but normally should not. |
---|
715 | // These are of technical nature, as described for each. |
---|
716 | |
---|
717 | // Minimum number of particles to perform study. |
---|
718 | const int CellJet::TIMESTOPRINT = 1; |
---|
719 | |
---|
720 | //-------------------------------------------------------------------------- |
---|
721 | |
---|
722 | // Analyze event. |
---|
723 | |
---|
724 | bool CellJet::analyze(const Event& event, double eTjetMinIn, |
---|
725 | double coneRadiusIn, double eTseedIn, ostream& ) { |
---|
726 | |
---|
727 | // Input values. Initial values zero. |
---|
728 | eTjetMin = eTjetMinIn; |
---|
729 | coneRadius = coneRadiusIn; |
---|
730 | eTseed = eTseedIn; |
---|
731 | jets.resize(0); |
---|
732 | vector<SingleCell> cells; |
---|
733 | |
---|
734 | // Loop over desired particles in the event. |
---|
735 | for (int i = 0; i < event.size(); ++i) |
---|
736 | if (event[i].isFinal()) { |
---|
737 | if (select > 2 && event[i].isNeutral() ) continue; |
---|
738 | if (select == 2 && !event[i].isVisible() ) continue; |
---|
739 | |
---|
740 | // Find particle position in (eta, phi, pT) space. |
---|
741 | double etaNow = event[i].eta(); |
---|
742 | if (abs(etaNow) > etaMax) continue; |
---|
743 | double phiNow = event[i].phi(); |
---|
744 | double pTnow = event[i].pT(); |
---|
745 | int iEtaNow = max(1, min( nEta, 1 + int(nEta * 0.5 |
---|
746 | * (1. + etaNow / etaMax) ) ) ); |
---|
747 | int iPhiNow = max(1, min( nPhi, 1 + int(nPhi * 0.5 |
---|
748 | * (1. + phiNow / M_PI) ) ) ); |
---|
749 | int iCell = nPhi * iEtaNow + iPhiNow; |
---|
750 | |
---|
751 | // Add pT to cell already hit or book a new cell. |
---|
752 | bool found = false; |
---|
753 | for (int j = 0; j < int(cells.size()); ++j) { |
---|
754 | if (iCell == cells[j].iCell) { |
---|
755 | found = true; |
---|
756 | ++cells[j].multiplicity; |
---|
757 | cells[j].eTcell += pTnow; |
---|
758 | continue; |
---|
759 | } |
---|
760 | } |
---|
761 | if (!found) { |
---|
762 | double etaCell = (etaMax / nEta) * (2 * iEtaNow - 1 - nEta); |
---|
763 | double phiCell = (M_PI / nPhi) * (2 * iPhiNow - 1 - nPhi); |
---|
764 | cells.push_back( SingleCell( iCell, etaCell, phiCell, pTnow, 1) ); |
---|
765 | } |
---|
766 | } |
---|
767 | |
---|
768 | // Smear true bin content by calorimeter resolution. |
---|
769 | if (smear > 0 && rndmPtr != 0) |
---|
770 | for (int j = 0; j < int(cells.size()); ++j) { |
---|
771 | double eTeConv = (smear < 2) ? 1. : cosh( cells[j].etaCell ); |
---|
772 | double eBef = cells[j].eTcell * eTeConv; |
---|
773 | double eAft = 0.; |
---|
774 | do eAft = eBef + resolution * sqrt(eBef) * rndmPtr->gauss(); |
---|
775 | while (eAft < 0 || eAft > upperCut * eBef); |
---|
776 | cells[j].eTcell = eAft / eTeConv; |
---|
777 | } |
---|
778 | |
---|
779 | // Remove cells below threshold for seed or for use at all. |
---|
780 | for (int j = 0; j < int(cells.size()); ++j) { |
---|
781 | if (cells[j].eTcell < eTseed) cells[j].canBeSeed = false; |
---|
782 | if (cells[j].eTcell < threshold) cells[j].isUsed = true; |
---|
783 | } |
---|
784 | |
---|
785 | // Find seed cell: the one with highest pT of not yet probed ones. |
---|
786 | for ( ; ; ) { |
---|
787 | int jMax = 0; |
---|
788 | double eTmax = 0.; |
---|
789 | for (int j = 0; j < int(cells.size()); ++j) |
---|
790 | if (cells[j].canBeSeed && cells[j].eTcell > eTmax) { |
---|
791 | jMax = j; |
---|
792 | eTmax = cells[j].eTcell; |
---|
793 | } |
---|
794 | |
---|
795 | // If too small cell eT then done, else start new trial jet. |
---|
796 | if (eTmax < eTseed) break; |
---|
797 | double etaCenterNow = cells[jMax].etaCell; |
---|
798 | double phiCenterNow = cells[jMax].phiCell; |
---|
799 | double eTjetNow = 0.; |
---|
800 | |
---|
801 | // Sum up unused cells within required distance of seed. |
---|
802 | for (int j = 0; j < int(cells.size()); ++j) { |
---|
803 | if (cells[j].isUsed) continue; |
---|
804 | double dEta = abs( cells[j].etaCell - etaCenterNow ); |
---|
805 | if (dEta > coneRadius) continue; |
---|
806 | double dPhi = abs( cells[j].phiCell - phiCenterNow ); |
---|
807 | if (dPhi > M_PI) dPhi = 2. * M_PI - dPhi; |
---|
808 | if (dPhi > coneRadius) continue; |
---|
809 | if (pow2(dEta) + pow2(dPhi) > pow2(coneRadius)) continue; |
---|
810 | cells[j].isAssigned = true; |
---|
811 | eTjetNow += cells[j].eTcell; |
---|
812 | } |
---|
813 | |
---|
814 | // Reject cluster below minimum ET. |
---|
815 | if (eTjetNow < eTjetMin) { |
---|
816 | cells[jMax].canBeSeed = false; |
---|
817 | for (int j = 0; j < int(cells.size()); ++j) |
---|
818 | cells[j].isAssigned = false; |
---|
819 | |
---|
820 | // Else find new jet properties. |
---|
821 | } else { |
---|
822 | double etaWeightedNow = 0.; |
---|
823 | double phiWeightedNow = 0.; |
---|
824 | int multiplicityNow = 0; |
---|
825 | Vec4 pMassiveNow; |
---|
826 | for (int j = 0; j < int(cells.size()); ++j) |
---|
827 | if (cells[j].isAssigned) { |
---|
828 | cells[j].canBeSeed = false; |
---|
829 | cells[j].isUsed = true; |
---|
830 | cells[j].isAssigned = false; |
---|
831 | etaWeightedNow += cells[j].eTcell * cells[j].etaCell; |
---|
832 | double phiCell = cells[j].phiCell; |
---|
833 | if (abs(phiCell - phiCenterNow) > M_PI) |
---|
834 | phiCell += (phiCenterNow > 0.) ? 2. * M_PI : -2. * M_PI; |
---|
835 | phiWeightedNow += cells[j].eTcell * phiCell; |
---|
836 | multiplicityNow += cells[j].multiplicity; |
---|
837 | pMassiveNow += cells[j].eTcell * Vec4( |
---|
838 | cos(cells[j].phiCell), sin(cells[j].phiCell), |
---|
839 | sinh(cells[j].etaCell), cosh(cells[j].etaCell) ); |
---|
840 | } |
---|
841 | etaWeightedNow /= eTjetNow; |
---|
842 | phiWeightedNow /= eTjetNow; |
---|
843 | |
---|
844 | // Bookkeep new jet, in decreasing ET order. |
---|
845 | jets.push_back( SingleCellJet( eTjetNow, etaCenterNow, phiCenterNow, |
---|
846 | etaWeightedNow, phiWeightedNow, multiplicityNow, pMassiveNow) ); |
---|
847 | for (int i = int(jets.size()) - 1; i > 0; --i) { |
---|
848 | if (jets[i-1].eTjet > jets[i].eTjet) break; |
---|
849 | swap( jets[i-1], jets[i]); |
---|
850 | } |
---|
851 | } |
---|
852 | } |
---|
853 | |
---|
854 | // Done. |
---|
855 | return true; |
---|
856 | } |
---|
857 | |
---|
858 | //-------------------------------------------------------------------------- |
---|
859 | |
---|
860 | // Provide a listing of the info. |
---|
861 | |
---|
862 | void CellJet::list(ostream& os) const { |
---|
863 | |
---|
864 | // Header. |
---|
865 | os << "\n -------- PYTHIA CellJet Listing, eTjetMin = " |
---|
866 | << fixed << setprecision(3) << setw(8) << eTjetMin |
---|
867 | << ", coneRadius = " << setw(5) << coneRadius |
---|
868 | << " ------------------------------ \n \n no " |
---|
869 | << " eTjet etaCtr phiCtr etaWt phiWt mult p_x" |
---|
870 | << " p_y p_z e m \n"; |
---|
871 | |
---|
872 | // The jets. |
---|
873 | for (int i = 0; i < int(jets.size()); ++i) { |
---|
874 | os << setw(4) << i << setw(10) << jets[i].eTjet << setw(8) |
---|
875 | << jets[i].etaCenter << setw(8) << jets[i].phiCenter << setw(8) |
---|
876 | << jets[i].etaWeighted << setw(8) << jets[i].phiWeighted |
---|
877 | << setw(5) << jets[i].multiplicity << setw(11) |
---|
878 | << jets[i].pMassive.px() << setw(11) << jets[i].pMassive.py() |
---|
879 | << setw(11) << jets[i].pMassive.pz() << setw(11) |
---|
880 | << jets[i].pMassive.e() << setw(11) |
---|
881 | << jets[i].pMassive.mCalc() << "\n"; |
---|
882 | } |
---|
883 | |
---|
884 | // Listing finished. |
---|
885 | os << "\n -------- End PYTHIA CellJet Listing ------------------" |
---|
886 | << "-------------------------------------------------" |
---|
887 | << endl; |
---|
888 | } |
---|
889 | |
---|
890 | //========================================================================== |
---|
891 | |
---|
892 | // SlowJet class. |
---|
893 | // This class performs clustering in (y, phi, pT) space. |
---|
894 | |
---|
895 | //-------------------------------------------------------------------------- |
---|
896 | |
---|
897 | // Constants: could be changed here if desired, but normally should not. |
---|
898 | // These are of technical nature, as described for each. |
---|
899 | |
---|
900 | // Minimum number of particles to perform study. |
---|
901 | const int SlowJet::TIMESTOPRINT = 1; |
---|
902 | |
---|
903 | // Assume the pi+- mass for all particles, except the photon, in one option. |
---|
904 | const double SlowJet::PIMASS = 0.13957; |
---|
905 | |
---|
906 | // Small number to avoid division by zero. |
---|
907 | const double SlowJet::TINY = 1e-20; |
---|
908 | |
---|
909 | //-------------------------------------------------------------------------- |
---|
910 | |
---|
911 | // Set up list of particles to analyze, and initial distances. |
---|
912 | |
---|
913 | bool SlowJet::setup(const Event& event) { |
---|
914 | |
---|
915 | // Initial values zero. |
---|
916 | clusters.resize(0); |
---|
917 | jets.resize(0); |
---|
918 | jtSize = 0; |
---|
919 | |
---|
920 | // Loop over final particles in the event. |
---|
921 | Vec4 pTemp; |
---|
922 | double mTemp, pT2Temp, mTTemp, yTemp, phiTemp; |
---|
923 | for (int i = 0; i < event.size(); ++i) |
---|
924 | if (event[i].isFinal()) { |
---|
925 | |
---|
926 | // Always apply selection options for visible or charged particles. |
---|
927 | if (chargedOnly && event[i].isNeutral() ) continue; |
---|
928 | else if (visibleOnly && !event[i].isVisible() ) continue; |
---|
929 | |
---|
930 | // Normally use built-in selection machinery. |
---|
931 | if (noHook) { |
---|
932 | |
---|
933 | // Pseudorapidity cut to describe detector range. |
---|
934 | if (cutInEta && abs(event[i].eta()) > etaMax) continue; |
---|
935 | |
---|
936 | // Optionally modify mass and energy. |
---|
937 | pTemp = event[i].p(); |
---|
938 | mTemp = event[i].m(); |
---|
939 | if (modifyMass) { |
---|
940 | mTemp = (massSet == 0 || event[i].id() == 22) ? 0. : PIMASS; |
---|
941 | pTemp.e( sqrt(pTemp.pAbs2() + mTemp*mTemp) ); |
---|
942 | } |
---|
943 | |
---|
944 | // Alternatively pass info to SlowJetHook for decision. |
---|
945 | // User can also modify pTemp and mTemp. |
---|
946 | } else { |
---|
947 | pTemp = event[i].p(); |
---|
948 | mTemp = event[i].m(); |
---|
949 | if ( !sjHookPtr->include( i, event, pTemp, mTemp) ) continue; |
---|
950 | } |
---|
951 | |
---|
952 | // Store particle momentum, including some derived quantities. |
---|
953 | pT2Temp = max( TINY*TINY, pTemp.pT2()); |
---|
954 | mTTemp = sqrt( mTemp*mTemp + pT2Temp); |
---|
955 | yTemp = (pTemp.pz() > 0) |
---|
956 | ? log( max( TINY, pTemp.e() + pTemp.pz() ) / mTTemp ) |
---|
957 | : log( mTTemp / max( TINY, pTemp.e() - pTemp.pz() ) ); |
---|
958 | phiTemp = pTemp.phi(); |
---|
959 | clusters.push_back( SingleSlowJet(pTemp, pT2Temp, yTemp, phiTemp, i) ); |
---|
960 | } |
---|
961 | |
---|
962 | // Resize arrays to store distances between clusters. |
---|
963 | origSize = clusters.size(); |
---|
964 | clSize = origSize; |
---|
965 | clLast = clSize - 1; |
---|
966 | diB.resize(clSize); |
---|
967 | dij.resize(clSize * (clSize - 1) / 2); |
---|
968 | |
---|
969 | // Loop through particles and find distance to beams. |
---|
970 | for (int i = 0; i < clSize; ++i) { |
---|
971 | if (isAnti) diB[i] = 1. / clusters[i].pT2; |
---|
972 | else if (isKT) diB[i] = clusters[i].pT2; |
---|
973 | else diB[i] = 1.; |
---|
974 | |
---|
975 | // Loop through pairs and find relative distance. |
---|
976 | for (int j = 0; j < i; ++j) { |
---|
977 | dPhi = abs( clusters[i].phi - clusters[j].phi ); |
---|
978 | if (dPhi > M_PI) dPhi = 2. * M_PI - dPhi; |
---|
979 | dijTemp = (pow2( clusters[i].y - clusters[j].y) + dPhi*dPhi) / R2; |
---|
980 | if (isAnti) dijTemp /= max(clusters[i].pT2, clusters[j].pT2); |
---|
981 | else if (isKT) dijTemp *= min(clusters[i].pT2, clusters[j].pT2); |
---|
982 | dij[i*(i-1)/2 + j] = dijTemp; |
---|
983 | |
---|
984 | // End of original-particle loops. |
---|
985 | } |
---|
986 | } |
---|
987 | |
---|
988 | // Find first particle pair to join. |
---|
989 | findNext(); |
---|
990 | |
---|
991 | // Done. |
---|
992 | return true; |
---|
993 | |
---|
994 | } |
---|
995 | |
---|
996 | //-------------------------------------------------------------------------- |
---|
997 | |
---|
998 | // Do one recombination step, possibly giving a jet. |
---|
999 | |
---|
1000 | bool SlowJet::doStep() { |
---|
1001 | |
---|
1002 | // Fail if no possibility to take a step. |
---|
1003 | if (clSize == 0) return false; |
---|
1004 | |
---|
1005 | // When distance to beam is smallest the cluster is promoted to jet. |
---|
1006 | if (jMin == -1) { |
---|
1007 | |
---|
1008 | // Store new jet if its pT is above pTMin. |
---|
1009 | if (clusters[iMin].pT2 > pT2jetMin) { |
---|
1010 | jets.push_back( SingleSlowJet(clusters[iMin]) ); |
---|
1011 | ++jtSize; |
---|
1012 | |
---|
1013 | // Order jets in decreasing pT. |
---|
1014 | for (int i = jtSize - 1; i > 0; --i) { |
---|
1015 | if (jets[i].pT2 < jets[i-1].pT2) break; |
---|
1016 | swap( jets[i], jets[i-1]); |
---|
1017 | } |
---|
1018 | } |
---|
1019 | } |
---|
1020 | |
---|
1021 | // When distance between two clusters is smallest they are joined. |
---|
1022 | else { |
---|
1023 | |
---|
1024 | // Add iMin cluster to jMin. |
---|
1025 | clusters[jMin].p += clusters[iMin].p; |
---|
1026 | clusters[jMin].pT2 = max( TINY*TINY, clusters[jMin].p.pT2()); |
---|
1027 | double mTTemp = sqrt(clusters[jMin].p.m2Calc() + clusters[jMin].pT2); |
---|
1028 | clusters[jMin].y = (clusters[jMin].p.pz() > 0) |
---|
1029 | ? log( max( TINY, clusters[jMin].p.e() + clusters[jMin].p.pz() ) |
---|
1030 | / mTTemp ) : log( mTTemp |
---|
1031 | / max( TINY, clusters[jMin].p.e() - clusters[jMin].p.pz() ) ); |
---|
1032 | clusters[jMin].phi = clusters[jMin].p.phi(); |
---|
1033 | clusters[jMin].mult += clusters[iMin].mult; |
---|
1034 | clusters[jMin].idx.insert(clusters[iMin].idx.begin(), |
---|
1035 | clusters[iMin].idx.end()); |
---|
1036 | |
---|
1037 | // Update distances for and to new jMin. |
---|
1038 | if (isAnti) diB[jMin] = 1. / clusters[jMin].pT2; |
---|
1039 | else if (isKT) diB[jMin] = clusters[jMin].pT2; |
---|
1040 | else diB[jMin] = 1.; |
---|
1041 | for (int i = 0; i < clSize; ++i) if (i != jMin && i != iMin) { |
---|
1042 | dPhi = abs( clusters[i].phi - clusters[jMin].phi ); |
---|
1043 | if (dPhi > M_PI) dPhi = 2. * M_PI - dPhi; |
---|
1044 | dijTemp = (pow2( clusters[i].y - clusters[jMin].y) + dPhi*dPhi) / R2; |
---|
1045 | if (isAnti) dijTemp /= max(clusters[i].pT2, clusters[jMin].pT2); |
---|
1046 | else if (isKT) dijTemp *= min(clusters[i].pT2, clusters[jMin].pT2); |
---|
1047 | if (i < jMin) dij[jMin*(jMin-1)/2 + i] = dijTemp; |
---|
1048 | else dij[i*(i-1)/2 + jMin] = dijTemp; |
---|
1049 | } |
---|
1050 | } |
---|
1051 | |
---|
1052 | // Move up last cluster and distances to vacated position iMin. |
---|
1053 | if (iMin < clLast) { |
---|
1054 | clusters[iMin] = clusters[clLast]; |
---|
1055 | diB[iMin] = diB[clLast]; |
---|
1056 | for (int j = 0; j < iMin; ++j) |
---|
1057 | dij[iMin*(iMin-1)/2 + j] = dij[clLast*(clLast-1)/2 + j]; |
---|
1058 | for (int j = iMin + 1; j < clLast; ++j) |
---|
1059 | dij[j*(j-1)/2 + iMin] = dij[clLast*(clLast-1)/2 + j]; |
---|
1060 | } |
---|
1061 | |
---|
1062 | // Shrink cluster list by one. |
---|
1063 | clusters.pop_back(); |
---|
1064 | --clSize; |
---|
1065 | --clLast; |
---|
1066 | |
---|
1067 | // Find next cluster pair to join. |
---|
1068 | findNext(); |
---|
1069 | |
---|
1070 | // Done. |
---|
1071 | return true; |
---|
1072 | |
---|
1073 | } |
---|
1074 | |
---|
1075 | //-------------------------------------------------------------------------- |
---|
1076 | |
---|
1077 | // Provide a listing of the info. |
---|
1078 | |
---|
1079 | void SlowJet::list(bool listAll, ostream& os) const { |
---|
1080 | |
---|
1081 | // Header. |
---|
1082 | os << "\n -------- PYTHIA SlowJet Listing, p = " << setw(2) |
---|
1083 | << power << ", R = " << fixed << setprecision(3) << setw(5) << R |
---|
1084 | << ", pTjetMin =" << setw(8) << pTjetMin << ", etaMax = " << setw(6) |
---|
1085 | << etaMax << " --- \n \n no pTjet y phi mult " |
---|
1086 | << " p_x p_y p_z e m \n"; |
---|
1087 | |
---|
1088 | // The jets. |
---|
1089 | for (int i = 0; i < jtSize; ++i) { |
---|
1090 | os << setw(4) << i << setw(11) << sqrt(jets[i].pT2) << setw(9) |
---|
1091 | << jets[i].y << setw(9) << jets[i].phi << setw(6) |
---|
1092 | << jets[i].mult << setw(11) << jets[i].p.px() << setw(11) |
---|
1093 | << jets[i].p.py() << setw(11) << jets[i].p.pz() << setw(11) |
---|
1094 | << jets[i].p.e() << setw(11) << jets[i].p.mCalc() << "\n"; |
---|
1095 | } |
---|
1096 | |
---|
1097 | // Optionally list also clusters not yet jets. |
---|
1098 | if (listAll && clSize > 0) { |
---|
1099 | os << " -------- Below this line follows remaining clusters," |
---|
1100 | << " still pT-unordered -------------------\n"; |
---|
1101 | for (int i = 0; i < clSize; ++i) { |
---|
1102 | os << setw(4) << i + jtSize << setw(11) << sqrt(clusters[i].pT2) |
---|
1103 | << setw(9) << clusters[i].y << setw(9) << clusters[i].phi |
---|
1104 | << setw(6) << clusters[i].mult << setw(11) << clusters[i].p.px() |
---|
1105 | << setw(11) << clusters[i].p.py() << setw(11) << clusters[i].p.pz() |
---|
1106 | << setw(11) << clusters[i].p.e() << setw(11) |
---|
1107 | << clusters[i].p.mCalc() << "\n"; |
---|
1108 | } |
---|
1109 | } |
---|
1110 | |
---|
1111 | // Listing finished. |
---|
1112 | os << "\n -------- End PYTHIA SlowJet Listing ------------------" |
---|
1113 | << "-------------------------------------" << endl; |
---|
1114 | |
---|
1115 | } |
---|
1116 | |
---|
1117 | //-------------------------------------------------------------------------- |
---|
1118 | |
---|
1119 | // Find next cluster pair to join. |
---|
1120 | |
---|
1121 | void SlowJet::findNext() { |
---|
1122 | |
---|
1123 | // Find smallest of diB, dij. |
---|
1124 | if (clSize > 0) { |
---|
1125 | iMin = 0; |
---|
1126 | jMin = -1; |
---|
1127 | dMin = diB[0]; |
---|
1128 | for (int i = 1; i < clSize; ++i) { |
---|
1129 | if (diB[i] < dMin) { |
---|
1130 | iMin = i; |
---|
1131 | jMin = -1; |
---|
1132 | dMin = diB[i]; |
---|
1133 | } |
---|
1134 | for (int j = 0; j < i; ++j) { |
---|
1135 | if (dij[i*(i-1)/2 + j] < dMin) { |
---|
1136 | iMin = i; |
---|
1137 | jMin = j; |
---|
1138 | dMin = dij[i*(i-1)/2 + j]; |
---|
1139 | } |
---|
1140 | } |
---|
1141 | } |
---|
1142 | |
---|
1143 | // If no clusters left then instead default values. |
---|
1144 | } else { |
---|
1145 | iMin = -1; |
---|
1146 | jMin = -1; |
---|
1147 | dMin = 0.; |
---|
1148 | } |
---|
1149 | |
---|
1150 | } |
---|
1151 | |
---|
1152 | //========================================================================== |
---|
1153 | |
---|
1154 | } // end namespace Pythia8 |
---|