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