PYTHIA8/pythia8130/src/Analysis.cxx

   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
   6 // Function definitions (not found in the header) for the
   7 // Sphericity, Thrust, ClusJet and CellJet 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. 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
 180 void Sphericity::list(ostream& os) {
 181
 182   // Header.
 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.
 214 const int    Thrust::NSTUDYMIN    = 2;
 215
 216 // Maximum number of times that an error warning will be printed.
 217 const int    Thrust::TIMESTOPRINT = 1;
 218
 219 // Major not too low or not possible to find major axis.
 220 const double Thrust::MAJORMIN     = 1e-10;
 221
 222 //*********
 223
 224 // Analyze event.
 225
 226 bool 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
 343 void Thrust::list(ostream& os) {
 344
 345   // Header.
 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.
 374 const double SingleClusterJet::PABSMIN  = 1e-10;
 375
 376 //*********
 377
 378 // Distance measures between two SingleClusterJet objects.
 379
 380 double dist2Fun(int measure, const SingleClusterJet& j1,
 381   const SingleClusterJet& j2) {
 382
 383   // JADE distance.
 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.
 409 const int    ClusterJet::TIMESTOPRINT   = 1;
 410
 411 // Assign minimal pAbs to avoid division by zero.
 412 const double ClusterJet::PABSMIN        = 1e-10;
 413
 414 // Initial pT/m preclustering scale as fraction of clustering one.
 415 const double ClusterJet::PRECLUSTERFRAC = 0.1;
 416
 417 // Step with which pT/m is reduced if preclustering gives too few jets.
 418 const double ClusterJet::PRECLUSTERSTEP = 0.8;
 419
 420 //*********
 421
 422 // Analyze event.
 423
 424 bool 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
 541 void 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
 610 void 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
 677 void ClusterJet::list(ostream& os) {
 678
 679   // Header.
 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.
 712 const int CellJet::TIMESTOPRINT = 1;
 713
 714 //*********
 715
 716 // Analyze event.
 717
 718 bool CellJet::analyze(const Event& event, double eTjetMinIn,
 719   double coneRadiusIn, double eTseedIn, ostream& ) {
 720
 721   // Input values. Initial values zero.
 722   eTjetMin   = eTjetMinIn;
 723   coneRadius = coneRadiusIn;
 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
 856 void CellJet::list(ostream& os) {
 857
 858   // Header.
 859   os << "\n --------  PYTHIA CellJet Listing, eTjetMin = "
 860      << fixed << setprecision(3) << setw(8) << eTjetMin
 861      << ", coneRadius = " << setw(5) << coneRadius
 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