StRoot  1
Analysis.cc
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
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
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
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
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
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;
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
865  os << "\n -------- PYTHIA CellJet Listing, eTjetMin = "
866  << fixed << setprecision(3) << setw(8) << eTjetMin
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) );
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
1035  // Update distances for and to new jMin.
1036  if (isAnti) diB[jMin] = 1. / clusters[jMin].pT2;
1037  else if (isKT) diB[jMin] = clusters[jMin].pT2;
1038  else diB[jMin] = 1.;
1039  for (int i = 0; i < clSize; ++i) if (i != jMin && i != iMin) {
1040  dPhi = abs( clusters[i].phi - clusters[jMin].phi );
1041  if (dPhi > M_PI) dPhi = 2. * M_PI - dPhi;
1042  dijTemp = (pow2( clusters[i].y - clusters[jMin].y) + dPhi*dPhi) / R2;
1043  if (isAnti) dijTemp /= max(clusters[i].pT2, clusters[jMin].pT2);
1044  else if (isKT) dijTemp *= min(clusters[i].pT2, clusters[jMin].pT2);
1045  if (i < jMin) dij[jMin*(jMin-1)/2 + i] = dijTemp;
1046  else dij[i*(i-1)/2 + jMin] = dijTemp;
1047  }
1048  }
1049
1050  // Move up last cluster and distances to vacated position iMin.
1051  if (iMin < clLast) {
1052  clusters[iMin] = clusters[clLast];
1053  diB[iMin] = diB[clLast];
1054  for (int j = 0; j < iMin; ++j)
1055  dij[iMin*(iMin-1)/2 + j] = dij[clLast*(clLast-1)/2 + j];
1056  for (int j = iMin + 1; j < clLast; ++j)
1057  dij[j*(j-1)/2 + iMin] = dij[clLast*(clLast-1)/2 + j];
1058  }
1059
1060  // Shrink cluster list by one.
1061  clusters.pop_back();
1062  --clSize;
1063  --clLast;
1064
1065  // Find next cluster pair to join.
1066  findNext();
1067
1068  // Done.
1069  return true;
1070
1071 }
1072
1073 //--------------------------------------------------------------------------
1074
1075 // Provide a listing of the info.
1076
1077 void SlowJet::list(bool listAll, ostream& os) const {
1078
1080  os << "\n -------- PYTHIA SlowJet Listing, p = " << setw(2)
1081  << power << ", R = " << fixed << setprecision(3) << setw(5) << R
1082  << ", pTjetMin =" << setw(8) << pTjetMin << ", etaMax = " << setw(6)
1083  << etaMax << " --- \n \n no pTjet y phi mult "
1084  << " p_x p_y p_z e m \n";
1085
1086  // The jets.
1087  for (int i = 0; i < jtSize; ++i) {
1088  os << setw(4) << i << setw(11) << sqrt(jets[i].pT2) << setw(9)
1089  << jets[i].y << setw(9) << jets[i].phi << setw(6)
1090  << jets[i].mult << setw(11) << jets[i].p.px() << setw(11)
1091  << jets[i].p.py() << setw(11) << jets[i].p.pz() << setw(11)
1092  << jets[i].p.e() << setw(11) << jets[i].p.mCalc() << "\n";
1093  }
1094
1095  // Optionally list also clusters not yet jets.
1096  if (listAll && clSize > 0) {
1097  os << " -------- Below this line follows remaining clusters,"
1098  << " still pT-unordered -------------------\n";
1099  for (int i = 0; i < clSize; ++i) {
1100  os << setw(4) << i + jtSize << setw(11) << sqrt(clusters[i].pT2)
1101  << setw(9) << clusters[i].y << setw(9) << clusters[i].phi
1102  << setw(6) << clusters[i].mult << setw(11) << clusters[i].p.px()
1103  << setw(11) << clusters[i].p.py() << setw(11) << clusters[i].p.pz()
1104  << setw(11) << clusters[i].p.e() << setw(11)
1105  << clusters[i].p.mCalc() << "\n";
1106  }
1107  }
1108
1109  // Listing finished.
1110  os << "\n -------- End PYTHIA SlowJet Listing ------------------"
1111  << "-------------------------------------" << endl;
1112
1113 }
1114
1115 //--------------------------------------------------------------------------
1116
1117 // Find next cluster pair to join.
1118
1119 void SlowJet::findNext() {
1120
1121  // Find smallest of diB, dij.
1122  if (clSize > 0) {
1123  iMin = 0;
1124  jMin = -1;
1125  dMin = diB[0];
1126  for (int i = 1; i < clSize; ++i) {
1127  if (diB[i] < dMin) {
1128  iMin = i;
1129  jMin = -1;
1130  dMin = diB[i];
1131  }
1132  for (int j = 0; j < i; ++j) {
1133  if (dij[i*(i-1)/2 + j] < dMin) {
1134  iMin = i;
1135  jMin = j;
1136  dMin = dij[i*(i-1)/2 + j];
1137  }
1138  }
1139  }
1140
1141  // If no clusters left then instead default values.
1142  } else {
1143  iMin = -1;
1144  jMin = -1;
1145  dMin = 0.;
1146  }
1147
1148 }
1149
1150 //==========================================================================
1151
1152 } // end namespace Pythia8
Definition: AgUStep.h:26