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