CommitLineData
5ad4eb21 1// Analysis.cc is a part of the PYTHIA event generator.
2// Copyright (C) 2008 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 and CellJet classes.
8
9#include "Analysis.h"
10
11namespace 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.
24const int Sphericity::NSTUDYMIN = 2;
25
26// Maximum number of times that an error warning will be printed.
27const int Sphericity::TIMESTOPRINT = 1;
28
29// Assign mimimum squared momentum in weight to avoid division by zero.
30const double Sphericity::P2MIN = 1e-20;
31
32// Second eigenvalue not too low or not possible to find eigenvectors.
33const double Sphericity::EIGENVALUEMIN = 1e-10;
34
35//*********
36
37// Analyze event.
38
39bool 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. Random sign.
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 if (Rndm::flat() > 0.5) length = -length;
159
160 // Store eigenvectors.
161 if (iVal == 0) eVec1 = Vec4( eVec[1] / length,
162 eVec[2] / length, eVec[3] / length, 0.);
163 else eVec3 = Vec4( eVec[1] / length,
164 eVec[2] / length, eVec[3] / length, 0.);
165 }
166
167 // Middle eigenvector is orthogonal to the other two.
168 eVec2 = cross3( eVec1, eVec3);
169 if (Rndm::flat() > 0.5) eVec2 = -eVec2;
170
171 // Done.
172 return true;
173
174}
175
176//*********
177
178// Provide a listing of the info.
179
180void Sphericity::list(ostream& os) {
181
183 os << "\n -------- PYTHIA Sphericity Listing -------- \n";
184 if (powerInt !=2) os << " Nonstandard momentum power = "
185 << fixed << setprecision(3) << setw(6) << power << "\n";
186 os << "\n no lambda e_x e_y e_z \n";
187
188 // The three eigenvalues and eigenvectors.
189 os << setprecision(5);
190 os << " 1" << setw(11) << eVal1 << setw(11) << eVec1.px()
191 << setw(10) << eVec1.py() << setw(10) << eVec1.pz() << "\n";
192 os << " 2" << setw(11) << eVal2 << setw(11) << eVec2.px()
193 << setw(10) << eVec2.py() << setw(10) << eVec2.pz() << "\n";
194 os << " 3" << setw(11) << eVal3 << setw(11) << eVec3.px()
195 << setw(10) << eVec3.py() << setw(10) << eVec3.pz() << "\n";
196
197 // Listing finished.
198 os << "\n -------- End PYTHIA Sphericity Listing ----" << endl;
199
200}
201
202
203//**************************************************************************
204
205// Thrust class.
206// This class finds thrust-related properties of an event.
207
208//*********
209
210// Constants: could be changed here if desired, but normally should not.
211// These are of technical nature, as described for each.
212
213// Minimum number of particles to perform study.
214const int Thrust::NSTUDYMIN = 2;
215
216// Maximum number of times that an error warning will be printed.
217const int Thrust::TIMESTOPRINT = 1;
218
219// Major not too low or not possible to find major axis.
220const double Thrust::MAJORMIN = 1e-10;
221
222//*********
223
224// Analyze event.
225
226bool Thrust::analyze(const Event& event, ostream& os) {
227
228 // Initial values and counters zero.
229 eVal1 = eVal2 = eVal3 = 0.;
230 eVec1 = eVec2 = eVec3 = 0.;
231 int nStudy = 0;
232 vector<Vec4> pOrder;
233 Vec4 pSum, nRef, pPart, pFull, pMax;
234
235 // Loop over desired particles in the event.
236 for (int i = 0; i < event.size(); ++i)
237 if (event[i].isFinal()) {
238 if (select > 2 && event[i].isNeutral() ) continue;
239 if (select == 2 && !event[i].isVisible() ) continue;
240 ++nStudy;
241
242 // Store momenta. Use energy component for absolute momentum.
243 Vec4 pNow = event[i].p();
244 pNow.e(pNow.pAbs());
245 pSum += pNow;
246 pOrder.push_back(pNow);
247 }
248
249 // Very low multiplicities (0 or 1) not considered.
250 if (nStudy < NSTUDYMIN) {
251 if (nFew < TIMESTOPRINT) os << " PYTHIA Error in " <<
252 "Thrust::analyze: too few particles" << endl;
253 ++nFew;
254 return false;
255 }
256
257 // Try all combinations of reference vector orthogonal to two particles.
258 for (int i1 = 0; i1 < nStudy - 1; ++i1)
259 for (int i2 = i1 + 1; i2 < nStudy; ++i2) {
260 nRef = cross3( pOrder[i1], pOrder[i2]);
261 nRef /= nRef.pAbs();
262 pPart = 0.;
263
264 // Add all momenta with sign; two choices for each reference particle.
265 for (int i = 0; i < nStudy; ++i) if (i != i1 && i != i2) {
266 if (dot3(pOrder[i], nRef) > 0.) pPart += pOrder[i];
267 else pPart -= pOrder[i];
268 }
269 for (int j = 0; j < 4; ++j) {
270 if (j == 0) pFull = pPart + pOrder[i1] + pOrder[i2];
271 else if (j == 1) pFull = pPart + pOrder[i1] - pOrder[i2];
272 else if (j == 2) pFull = pPart - pOrder[i1] + pOrder[i2];
273 else pFull = pPart - pOrder[i1] - pOrder[i2];
274 pFull.e(pFull.pAbs());
275 if (pFull.e() > pMax.e()) pMax = pFull;
276 }
277 }
278
279 // Maximum gives thrust axis and value.
280 eVal1 = pMax.e() / pSum.e();
281 eVec1 = pMax / pMax.e();
282 eVec1.e(0.);
283
284 // Subtract momentum along thrust axis.
285 double pAbsSum = 0.;
286 for (int i = 0; i < nStudy; ++i) {
287 pOrder[i] -= dot3( eVec1, pOrder[i]) * eVec1;
288 pOrder[i].e(pOrder[i].pAbs());
289 pAbsSum += pOrder[i].e();
290 }
291
292 // Simpleminded major and minor axes if too little transverse left.
293 if (pAbsSum < MAJORMIN * pSum.e()) {
294 if ( abs(eVec1.pz()) > 0.5) eVec2 = Vec4( 1., 0., 0., 0.);
295 else eVec2 = Vec4( 0., 0., 1., 0.);
296 eVec2 -= dot3( eVec1, eVec2) * eVec1;
297 eVec2 /= eVec2.pAbs();
298 eVec3 = cross3( eVec1, eVec2);
299 return true;
300 }
301
302 // Try all reference vectors orthogonal to one particles.
303 pMax = 0.;
304 for (int i1 = 0; i1 < nStudy; ++i1) {
305 nRef = cross3( pOrder[i1], eVec1);
306 nRef /= nRef.pAbs();
307 pPart = 0.;
308
309 // Add all momenta with sign; two choices for each reference particle.
310 for (int i = 0; i < nStudy; ++i) if (i != i1) {
311 if (dot3(pOrder[i], nRef) > 0.) pPart += pOrder[i];
312 else pPart -= pOrder[i];
313 }
314 pFull = pPart + pOrder[i1];
315 pFull.e(pFull.pAbs());
316 if (pFull.e() > pMax.e()) pMax = pFull;
317 pFull = pPart - pOrder[i1];
318 pFull.e(pFull.pAbs());
319 if (pFull.e() > pMax.e()) pMax = pFull;
320 }
321
322 // Maximum gives major axis and value.
323 eVal2 = pMax.e() / pSum.e();
324 eVec2 = pMax / pMax.e();
325 eVec2.e(0.);
326
327 // Orthogonal direction gives minor axis, and from there value.
328 eVec3 = cross3( eVec1, eVec2);
329 pAbsSum = 0.;
330 for (int i = 0; i < nStudy; ++i)
331 pAbsSum += abs( dot3(eVec3, pOrder[i]) );
332 eVal3 = pAbsSum / pSum.e();
333
334 // Done.
335 return true;
336
337}
338
339//*********
340
341// Provide a listing of the info.
342
343void Thrust::list(ostream& os) {
344
346 os << "\n -------- PYTHIA Thrust Listing ------------ \n"
347 << "\n value e_x e_y e_z \n";
348
349 // The thrust, major and minor values and related event axes.
350 os << setprecision(5);
351 os << " Thr" << setw(11) << eVal1 << setw(11) << eVec1.px()
352 << setw(10) << eVec1.py() << setw(10) << eVec1.pz() << "\n";
353 os << " Maj" << setw(11) << eVal2 << setw(11) << eVec2.px()
354 << setw(10) << eVec2.py() << setw(10) << eVec2.pz() << "\n";
355 os << " Min" << setw(11) << eVal3 << setw(11) << eVec3.px()
356 << setw(10) << eVec3.py() << setw(10) << eVec3.pz() << "\n";
357
358 // Listing finished.
359 os << "\n -------- End PYTHIA Thrust Listing --------" << endl;
360
361}
362
363//**************************************************************************
364
365// SingleClusterJet class.
366// Simple helper class to ClusterJet for a jet and its contents.
367
368//*********
369
370// Constants: could be changed here if desired, but normally should not.
371// These are of technical nature, as described for each.
372
373// Assign minimal pAbs to avoid division by zero.
374const double SingleClusterJet::PABSMIN = 1e-10;
375
376//*********
377
378// Distance measures between two SingleClusterJet objects.
379
380double dist2Fun(int measure, const SingleClusterJet& j1,
381 const SingleClusterJet& j2) {
382
384 if (measure == 2) return 2. * j1.pJet.e() * j2.pJet.e()
385 * (1. - dot3( j1.pJet, j2.pJet) / (j1.pAbs * j2.pAbs) );
386
387 // Durham distance.
388 if (measure == 3) return 2. * pow2( min( j1.pJet.e(), j2.pJet.e() ) )
389 * (1. - dot3( j1.pJet, j2.pJet) / (j1.pAbs * j2.pAbs) );
390
391 // Lund distance; "default".
392 return (j1.pAbs * j2.pAbs - dot3( j1.pJet, j2.pJet))
393 * 2. * j1.pAbs * j2.pAbs / pow2(j1.pAbs + j2.pAbs);
394
395}
396
397//**************************************************************************
398
399// ClusterJet class.
400// This class performs a jet clustering according to different
401// distance measures: Lund, JADE or Durham.
402
403//*********
404
405// Constants: could be changed here if desired, but normally should not.
406// These are of technical nature, as described for each.
407
408// Maximum number of times that an error warning will be printed.
409const int ClusterJet::TIMESTOPRINT = 1;
410
411// Assign minimal pAbs to avoid division by zero.
412const double ClusterJet::PABSMIN = 1e-10;
413
414// Initial pT/m preclustering scale as fraction of clustering one.
415const double ClusterJet::PRECLUSTERFRAC = 0.1;
416
417// Step with which pT/m is reduced if preclustering gives too few jets.
418const double ClusterJet::PRECLUSTERSTEP = 0.8;
419
420//*********
421
422// Analyze event.
423
424bool ClusterJet::analyze(const Event& event, double yScaleIn,
425 double pTscaleIn, int nJetMinIn, int nJetMaxIn, ostream& os) {
426
427 // Input values. Initial values zero.
428 yScale = yScaleIn;
429 pTscale = pTscaleIn;
430 nJetMin = nJetMinIn;
431 nJetMax = nJetMaxIn;
432 particles.resize(0);
433 jets.resize(0);
434 Vec4 pSum;
435
436 // Loop over desired particles in the event.
437 for (int i = 0; i < event.size(); ++i)
438 if (event[i].isFinal()) {
439 if (select > 2 && event[i].isNeutral() ) continue;
440 if (select == 2 && !event[i].isVisible() ) continue;
441
442 // Store them, possibly with modified mass => new energy.
443 Vec4 pTemp = event[i].p();
444 if (massSet == 0 || massSet == 1) {
445 double mTemp = (massSet == 0 || event[i].id() == 22)
446 ? 0. : piMass;
447 double eTemp = sqrt(pTemp.pAbs2() + pow2(mTemp));
448 pTemp.e(eTemp);
449 }
450 particles.push_back( SingleClusterJet(pTemp, i) );
451 pSum += pTemp;
452 }
453
454 // Very low multiplicities not considered.
455 nParticles = particles.size();
456 if (nParticles < nJetMin) {
457 if (nFew < TIMESTOPRINT) os << " PYTHIA Error in " <<
458 "ClusterJet::analyze: too few particles" << endl;
459 ++nFew;
460 return false;
461 }
462
463 // Squared maximum distance in GeV^2 for joining.
464 double p2Sum = pSum.m2Calc();
465 dist2Join = max( yScale * p2Sum, pow2(pTscale));
466 dist2BigMin = 2. * max( dist2Join, p2Sum);
467
468 // Do preclustering if desired and possible.
469 if (doPrecluster && nParticles > nJetMin + 2) {
470 precluster();
471 if (doReassign) reassign();
472 }
473
474 // If no preclustering: each particle is a starting jet.
475 else for (int i = 0; i < nParticles; ++i) {
476 jets.push_back( SingleClusterJet(particles[i]) );
477 particles[i].daughter = i;
478 }
479
480 // Begin iteration towards fewer jets.
481 for ( ; ; ) {
482
483 // Find the two closest jets.
484 double dist2Min = dist2BigMin;
485 int jMin = 0;
486 int kMin = 0;
487 for (int j = 0; j < int(jets.size()) - 1; ++j)
488 for (int k = j + 1; k < int(jets.size()); ++k) {
489 double dist2 = dist2Fun( measure, jets[j], jets[k]);
490 if (dist2 < dist2Min) {
491 dist2Min = dist2;
492 jMin = j;
493 kMin = k;
494 }
495 }
496
497 // Stop if no pair below cut and not more jets than allowed.
498 if ( dist2Min > dist2Join
499 && (nJetMax < nJetMin || int(jets.size()) <= nJetMax) ) break;
500
501 // Join two closest jets.
502 jets[jMin].pJet += jets[kMin].pJet;
503 jets[jMin].pAbs = max( PABSMIN, jets[jMin].pJet.pAbs());
504 jets[jMin].multiplicity += jets[kMin].multiplicity;
505 for (int i = 0; i < nParticles; ++i)
506 if (particles[i].daughter == kMin) particles[i].daughter = jMin;
507
508 // Move up last jet to empty slot to shrink list.
509 jets[kMin] = jets.back();
510 jets.pop_back();
511 int iEnd = jets.size();
512 for (int i = 0; i < nParticles; ++i)
513 if (particles[i].daughter == iEnd) particles[i].daughter = kMin;
514
515 // Do reassignments of particles to nearest jet if desired.
516 if (doReassign) reassign();
517
518 // Stop if reached minimum allowed number of jets. Else continue.
519 if (int(jets.size()) <= nJetMin) break;
520 }
521
522 // Order jets in decreasing energy.
523 for (int j = 0; j < int(jets.size()) - 1; ++j)
524 for (int k = int(jets.size()) - 1; k > j; --k)
525 if (jets[k].pJet.e() > jets[k-1].pJet.e()) {
526 swap( jets[k], jets[k-1]);
527 for (int i = 0; i < nParticles; ++i) {
528 if (particles[i].daughter == k) particles[i].daughter = k-1;
529 else if (particles[i].daughter == k-1) particles[i].daughter = k;
530 }
531 }
532
533 // Done.
534 return true;
535}
536
537//*********
538
539// Precluster nearby particles to save computer time.
540
541void ClusterJet::precluster() {
542
543 // Begin iteration over preclustering scale.
544 distPre = PRECLUSTERFRAC * sqrt(dist2Join) / PRECLUSTERSTEP;
545 for ( ; ;) {
546 distPre *= PRECLUSTERSTEP;
547 dist2Pre = pow2(distPre);
548 for (int i = 0; i < nParticles; ++i) {
549 particles[i].daughter = -1;
550 particles[i].isAssigned = false;
551 }
552
553 // Sum up low-momentum region. Jet if enough momentum.
554 Vec4 pCentral;
555 int multCentral = 0;
556 for (int i = 0; i < nParticles; ++i)
557 if (particles[i].pAbs < 2. * distPre) {
558 pCentral += particles[i].pJet;
559 multCentral += particles[i].multiplicity;
560 particles[i].isAssigned = true;
561 }
562 if (pCentral.pAbs() > 2. * distPre) {
563 jets.push_back( SingleClusterJet(pCentral) );
564 jets.back().multiplicity = multCentral;
565 for (int i = 0; i < nParticles; ++i)
566 if (particles[i].isAssigned) particles[i].daughter = 0;
567 }
568
569 // Find fastest remaining particle until none left.
570 for ( ; ;) {
571 int iMax = -1;
572 double pMax = 0.;
573 for (int i = 0; i < nParticles; ++i)
574 if ( !particles[i].isAssigned && particles[i].pAbs > pMax) {
575 iMax = i;
576 pMax = particles[i].pAbs;
577 }
578 if (iMax == -1) break;
579
580 // Sum up precluster around it according to distance function.
581 Vec4 pPre;
582 int multPre = 0;
583 int nRemain = 0;
584 for (int i = 0; i < nParticles; ++i)
585 if ( !particles[i].isAssigned) {
586 double dist2 = dist2Fun( measure, particles[iMax],
587 particles[i]);
588 if (dist2 < dist2Pre) {
589 pPre += particles[i].pJet;
590 ++multPre;
591 particles[i].isAssigned = true;
592 particles[i].daughter = jets.size();
593 } else ++nRemain;
594 }
595 jets.push_back( SingleClusterJet(pPre) );
596 jets.back().multiplicity = multPre;
597
598 // Decide whether sensible starting configuration or iterate.
599 if (int(jets.size()) + nRemain < nJetMin) break;
600 }
601 if (int(jets.size()) >= nJetMin) break;
602 }
603
604}
605
606//*********
607
608// Reassign particles to nearest jet to correct misclustering.
609
610void ClusterJet::reassign() {
611
612 // Reset clustered momenta.
613 for (int j = 0; j < int(jets.size()); ++j) {
614 jets[j].pTemp = 0.;
615 jets[j].multiplicity = 0;
616 }
617
618 // Loop through particles to find closest jet.
619 for (int i = 0; i < nParticles; ++i) {
620 particles[i].daughter = -1;
621 double dist2Min = dist2BigMin;
622 int jMin = 0;
623 for (int j = 0; j < int(jets.size()); ++j) {
624 double dist2 = dist2Fun( measure, particles[i], jets[j]);
625 if (dist2 < dist2Min) {
626 dist2Min = dist2;
627 jMin = j;
628 }
629 }
630 jets[jMin].pTemp += particles[i].pJet;
631 ++jets[jMin].multiplicity;
632 particles[i].daughter = jMin;
633 }
634
635 // Replace old by new jet momenta.
636 for (int j = 0; j < int(jets.size()); ++j) {
637 jets[j].pJet = jets[j].pTemp;
638 jets[j].pAbs = max( PABSMIN, jets[j].pJet.pAbs());
639 }
640
641 // Check that no empty clusters after reassignments.
642 for ( ; ; ) {
643
644 // If no empty jets then done.
645 int jEmpty = -1;
646 for (int j = 0; j < int(jets.size()); ++j)
647 if (jets[j].multiplicity == 0) jEmpty = j;
648 if (jEmpty == -1) return;
649
650 // Find particle assigned to jet with largest distance to it.
651 int iSplit = -1;
652 double dist2Max = 0.;
653 for (int i = 0; i < nParticles; ++i) {
654 int j = particles[i].daughter;
655 double dist2 = dist2Fun( measure, particles[i], jets[j]);
656 if (dist2 > dist2Max) {
657 iSplit = i;
658 dist2Max = dist2;
659 }
660 }
661
662 // Let this particle form new jet and subtract off from existing.
663 int jSplit = particles[iSplit].daughter;
664 jets[jEmpty] = SingleClusterJet( particles[iSplit].pJet );
665 jets[jSplit].pJet -= particles[iSplit].pJet;
666 jets[jSplit].pAbs = max( PABSMIN,jets[jSplit].pJet.pAbs());
667 particles[iSplit].daughter = jEmpty;
668 --jets[jSplit].multiplicity;
669 }
670
671}
672
673//*********
674
675// Provide a listing of the info.
676
677void ClusterJet::list(ostream& os) {
678
680 string method = (measure == 1) ? "Lund pT"
681 : ( (measure == 2) ? "JADE m" : "Durham kT" ) ;
682 os << "\n -------- PYTHIA ClusterJet Listing, " << setw(9) << method
683 << " =" << fixed << setprecision(3) << setw(7) << sqrt(dist2Join)
684 << " GeV --- \n \n no mult p_x p_y p_z "
685 << " e m \n";
686
687 // The jets.
688 for (int i = 0; i < int(jets.size()); ++i) {
689 os << setw(4) << i << setw(6) << jets[i].multiplicity << setw(11)
690 << jets[i].pJet.px() << setw(11) << jets[i].pJet.py()
691 << setw(11) << jets[i].pJet.pz() << setw(11)
692 << jets[i].pJet.e() << setw(11) << jets[i].pJet.mCalc()
693 << "\n";
694 }
695
696 // Listing finished.
697 os << "\n -------- End PYTHIA ClusterJet Listing ---------------"
698 << "--------" << endl;
699}
700
701//**************************************************************************
702
703// CellJet class.
704// This class performs a cone jet search in (eta, phi, E_T) space.
705
706//*********
707
708// Constants: could be changed here if desired, but normally should not.
709// These are of technical nature, as described for each.
710
711// Minimum number of particles to perform study.
712const int CellJet::TIMESTOPRINT = 1;
713
714//*********
715
716// Analyze event.
717
718bool CellJet::analyze(const Event& event, double eTjetMinIn,
719 double coneRadiusIn, double eTseedIn, ostream& ) {
720
721 // Input values. Initial values zero.
722 eTjetMin = eTjetMinIn;
724 eTseed = eTseedIn;
725 jets.resize(0);
726 vector<SingleCell> cells;
727
728 // Loop over desired particles in the event.
729 for (int i = 0; i < event.size(); ++i)
730 if (event[i].isFinal()) {
731 if (select > 2 && event[i].isNeutral() ) continue;
732 if (select == 2 && !event[i].isVisible() ) continue;
733
734 // Find particle position in (eta, phi, pT) space.
735 double etaNow = event[i].eta();
736 if (abs(etaNow) > etaMax) continue;
737 double phiNow = event[i].phi();
738 double pTnow = event[i].pT();
739 int iEtaNow = max(1, min( nEta, 1 + int(nEta * 0.5
740 * (1. + etaNow / etaMax) ) ) );
741 int iPhiNow = max(1, min( nPhi, 1 + int(nPhi * 0.5
742 * (1. + phiNow / M_PI) ) ) );
743 int iCell = nPhi * iEtaNow + iPhiNow;
744
745 // Add pT to cell already hit or book a new cell.
746 bool found = false;
747 for (int j = 0; j < int(cells.size()); ++j) {
748 if (iCell == cells[j].iCell) {
749 found = true;
750 ++cells[j].multiplicity;
751 cells[j].eTcell += pTnow;
752 continue;
753 }
754 }
755 if (!found) {
756 double etaCell = (etaMax / nEta) * (2 * iEtaNow - 1 - nEta);
757 double phiCell = (M_PI / nPhi) * (2 * iPhiNow - 1 - nPhi);
758 cells.push_back( SingleCell( iCell, etaCell, phiCell, pTnow, 1) );
759 }
760 }
761
762 // Smear true bin content by calorimeter resolution.
763 if (smear > 0)
764 for (int j = 0; j < int(cells.size()); ++j) {
765 double eTeConv = (smear < 2) ? 1. : cosh( cells[j].etaCell );
766 double eBef = cells[j].eTcell * eTeConv;
767 double eAft = 0.;
768 do eAft = eBef + resolution * sqrt(eBef) * Rndm::gauss();
769 while (eAft < 0 || eAft > upperCut * eBef);
770 cells[j].eTcell = eAft / eTeConv;
771 }
772
773 // Remove cells below threshold for seed or for use at all.
774 for (int j = 0; j < int(cells.size()); ++j) {
775 if (cells[j].eTcell < eTseed) cells[j].canBeSeed = false;
776 if (cells[j].eTcell < threshold) cells[j].isUsed = true;
777 }
778
779 // Find seed cell: the one with highest pT of not yet probed ones.
780 for ( ; ; ) {
781 int jMax = 0;
782 double eTmax = 0.;
783 for (int j = 0; j < int(cells.size()); ++j)
784 if (cells[j].canBeSeed && cells[j].eTcell > eTmax) {
785 jMax = j;
786 eTmax = cells[j].eTcell;
787 }
788
789 // If too small cell eT then done, else start new trial jet.
790 if (eTmax < eTseed) break;
791 double etaCenterNow = cells[jMax].etaCell;
792 double phiCenterNow = cells[jMax].phiCell;
793 double eTjetNow = 0.;
794
795 // Sum up unused cells within required distance of seed.
796 for (int j = 0; j < int(cells.size()); ++j) {
797 if (cells[j].isUsed) continue;
798 double dEta = abs( cells[j].etaCell - etaCenterNow );
799 if (dEta > coneRadius) continue;
800 double dPhi = abs( cells[j].phiCell - phiCenterNow );
801 if (dPhi > M_PI) dPhi = 2. * M_PI - dPhi;
802 if (dPhi > coneRadius) continue;
803 if (pow2(dEta) + pow2(dPhi) > pow2(coneRadius)) continue;
804 cells[j].isAssigned = true;
805 eTjetNow += cells[j].eTcell;
806 }
807
808 // Reject cluster below minimum ET.
809 if (eTjetNow < eTjetMin) {
810 cells[jMax].canBeSeed = false;
811 for (int j = 0; j < int(cells.size()); ++j)
812 cells[j].isAssigned = false;
813
814 // Else find new jet properties.
815 } else {
816 double etaWeightedNow = 0.;
817 double phiWeightedNow = 0.;
818 int multiplicityNow = 0;
819 Vec4 pMassiveNow;
820 for (int j = 0; j < int(cells.size()); ++j)
821 if (cells[j].isAssigned) {
822 cells[j].canBeSeed = false;
823 cells[j].isUsed = true;
824 cells[j].isAssigned = false;
825 etaWeightedNow += cells[j].eTcell * cells[j].etaCell;
826 double phiCell = cells[j].phiCell;
827 if (abs(phiCell - phiCenterNow) > M_PI)
828 phiCell += (phiCenterNow > 0.) ? 2. * M_PI : -2. * M_PI;
829 phiWeightedNow += cells[j].eTcell * phiCell;
830 multiplicityNow += cells[j].multiplicity;
831 pMassiveNow += cells[j].eTcell * Vec4(
832 cos(cells[j].phiCell), sin(cells[j].phiCell),
833 sinh(cells[j].etaCell), cosh(cells[j].etaCell) );
834 }
835 etaWeightedNow /= eTjetNow;
836 phiWeightedNow /= eTjetNow;
837
838 // Bookkeep new jet, in decreasing ET order.
839 jets.push_back( SingleCellJet( eTjetNow, etaCenterNow, phiCenterNow,
840 etaWeightedNow, phiWeightedNow, multiplicityNow, pMassiveNow) );
841 for (int i = int(jets.size()) - 1; i > 0; --i) {
842 if (jets[i-1].eTjet > jets[i].eTjet) break;
843 swap( jets[i-1], jets[i]);
844 }
845 }
846 }
847
848 // Done.
849 return true;
850}
851
852//*********
853
854// Provide a listing of the info.
855
856void CellJet::list(ostream& os) {
857
859 os << "\n -------- PYTHIA CellJet Listing, eTjetMin = "
860 << fixed << setprecision(3) << setw(8) << eTjetMin
862 << " ------------------------------ \n \n no "
863 << " eTjet etaCtr phiCtr etaWt phiWt mult p_x"
864 << " p_y p_z e m \n";
865
866 // The jets.
867 for (int i = 0; i < int(jets.size()); ++i) {
868 os << setw(4) << i << setw(10) << jets[i].eTjet << setw(8)
869 << jets[i].etaCenter << setw(8) << jets[i].phiCenter << setw(8)
870 << jets[i].etaWeighted << setw(8) << jets[i].phiWeighted
871 << setw(5) << jets[i].multiplicity << setw(11)
872 << jets[i].pMassive.px() << setw(11) << jets[i].pMassive.py()
873 << setw(11) << jets[i].pMassive.pz() << setw(11)
874 << jets[i].pMassive.e() << setw(11)
875 << jets[i].pMassive.mCalc() << "\n";
876 }
877
878 // Listing finished.
879 os << "\n -------- End PYTHIA CellJet Listing ------------------"
880 << "-------------------------------------------------"
881 << endl;
882}
883
884//**************************************************************************
885
886} // end namespace Pythia8