PYTHIA8/pythia8145/src/Event.cxx

   1 // Event.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2010 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 // Particle and Event classes, and some related global functions.
   8
   9 #include "Event.h"
  10
  11 namespace Pythia8 {
  12
  13 //==========================================================================
  14
  15 // Particle class.
  16 // This class holds info on a particle in general.
  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 // Small number to avoid division by zero.
  24 const double Particle::TINY = 1e-20;
  25
  26 //--------------------------------------------------------------------------
  27
  28 // Functions for rapidity and pseudorapidity.
  29
  30 double Particle::y() const {
  31   double temp = log( ( pSave.e() + abs(pSave.pz()) ) / max( TINY, mT() ) );
  32   return (pSave.pz() > 0) ? temp : -temp;
  33 }
  34
  35 double Particle::eta() const {
  36   double temp = log( ( pSave.pAbs() + abs(pSave.pz()) ) / max( TINY, pT() ) );
  37   return (pSave.pz() > 0) ? temp : -temp;
  38 }
  39
  40 //--------------------------------------------------------------------------
  41
  42 // Particle name, with status but imposed maximum length -> may truncate.
  43
  44 string Particle::nameWithStatus(int maxLen) const {
  45
  46   if (pdePtr == 0) return " ";
  47   string temp = (statusSave > 0) ? pdePtr->name(idSave)
  48     : "(" + pdePtr->name(idSave) + ")";
  49   while (int(temp.length()) > maxLen) {
  50     // Remove from end, excluding closing bracket and charge.
  51     int iRem = temp.find_last_not_of(")+-0");
  52     temp.erase(iRem, 1);
  53   }
  54   return temp;
  55 }
  56
  57 //--------------------------------------------------------------------------
  58
  59 // Add offsets to mother and daughter pointers (must be non-negative).
  60
  61 void Particle::offsetHistory( int minMother, int addMother, int minDaughter,
  62   int addDaughter) {
  63
  64   if (addMother < 0 || addDaughter < 0) return;
  65   if (  mother1Save > minMother  )   mother1Save += addMother;
  66   if (  mother2Save > minMother  )   mother2Save += addMother;
  67   if (daughter1Save > minDaughter) daughter1Save += addDaughter;
  68   if (daughter2Save > minDaughter) daughter2Save += addDaughter;
  69
  70 }
  71
  72 //--------------------------------------------------------------------------
  73
  74 // Add offsets to colour and anticolour (must be positive).
  75
  76 void Particle::offsetCol( int addCol) {
  77
  78   if (addCol < 0) return;
  79   if ( colSave > 0)  colSave += addCol;
  80   if (acolSave > 0) acolSave += addCol;
  81
  82 }
  83
  84 //--------------------------------------------------------------------------
  85
  86 // Invariant mass of a pair and its square.
  87 // (Not part of class proper, but tightly linked.)
  88
  89 double m(const Particle& pp1, const Particle& pp2) {
  90   double m2 = pow2(pp1.e() + pp2.e()) - pow2(pp1.px() + pp2.px())
  91      - pow2(pp1.py() + pp2.py()) - pow2(pp1.pz() + pp2.pz());
  92   return (m2 > 0. ? sqrt(m2) : 0.);
  93 }
  94
  95 double m2(const Particle& pp1, const Particle& pp2) {
  96   double m2 = pow2(pp1.e() + pp2.e()) - pow2(pp1.px() + pp2.px())
  97      - pow2(pp1.py() + pp2.py()) - pow2(pp1.pz() + pp2.pz());
  98   return m2;
  99 }
 100
 101 //==========================================================================
 102
 103 // Event class.
 104 // This class holds info on the complete event record.
 105
 106 //--------------------------------------------------------------------------
 107
 108 // Constants: could be changed here if desired, but normally should not.
 109 // These are of technical nature, as described for each.
 110
 111 // Maxmimum number of mothers or daughter indices per line in listing.
 112 const int Event::IPERLINE = 20;
 113
 114 //--------------------------------------------------------------------------
 115
 116 // Copy all information from one event record to another.
 117
 118 Event& Event::operator=( const Event& oldEvent) {
 119
 120   // Do not copy if same.
 121   if (this != &oldEvent) {
 122
 123     // Reset all current info in the event.
 124     clear();
 125
 126     // Copy all the particles one by one.
 127     for (int i = 0; i < oldEvent.size(); ++i) append( oldEvent[i] );
 128
 129     // Copy all the junctions one by one.
 130     for (int i = 0; i < oldEvent.sizeJunction(); ++i)
 131       appendJunction( oldEvent.getJunction(i) );
 132
 133     // Copy all other values.
 134     startColTag         = oldEvent.startColTag;
 135     maxColTag           = oldEvent.maxColTag;
 136     savedSize           = oldEvent.savedSize;
 137     savedJunctionSize   = oldEvent.savedJunctionSize;
 138     scaleSave           = oldEvent.scaleSave;
 139     scaleSecondSave     = oldEvent.scaleSecondSave;
 140     headerList          = oldEvent.headerList;
 141     particleDataPtr     = oldEvent.particleDataPtr;
 142
 143   // Done.
 144   }
 145   return *this;
 146
 147 }
 148
 149 //--------------------------------------------------------------------------
 150
 151 // Add a copy of an existing particle at the end of the event record;
 152 // return index. Three cases, depending on sign of new status code:
 153 // Positive: copy is viewed as daughter, status of original is negated.
 154 // Negative: copy is viewed as mother, status of original is unchanged.
 155 // Zero: the new is a perfect carbon copy (maybe to be changed later).
 156
 157 int Event::copy(int iCopy, int newStatus) {
 158
 159   // Simple carbon copy.
 160   entry.push_back(entry[iCopy]);
 161   int iNew = entry.size() - 1;
 162
 163   // Set up to make new daughter of old.
 164   if (newStatus > 0) {
 165     entry[iCopy].daughters(iNew,iNew);
 166     entry[iCopy].statusNeg();
 167     entry[iNew].mothers(iCopy, iCopy);
 168     entry[iNew].status(newStatus);
 169
 170   // Set up to make new mother of old.
 171   } else if (newStatus < 0) {
 172     entry[iCopy].mothers(iNew,iNew);
 173     entry[iNew].daughters(iCopy, iCopy);
 174     entry[iNew].status(newStatus);
 175   }
 176
 177   // Done.
 178   return iNew;
 179
 180 }
 181
 182 //--------------------------------------------------------------------------
 183
 184 // Print an event - special cases that rely on the general method.
 185 // Not inline to make them directly callable in (some) debuggers.
 186
 187 void Event::list() const {
 188   list(false, false, cout);
 189 }
 190
 191 void Event::list(ostream& os) const {
 192   list(false, false, os);
 193 }
 194
 195 void Event::list(bool showScaleAndVertex, bool showMothersAndDaughters)
 196 const {
 197   list(showScaleAndVertex, showMothersAndDaughters, cout);
 198 }
 199
 200 //--------------------------------------------------------------------------
 201
 202 // Print an event.
 203
 204 void Event::list(bool showScaleAndVertex, bool showMothersAndDaughters,
 205   ostream& os) const {
 206
 207   // Header.
 208   os << "\n --------  PYTHIA Event Listing  " << headerList << "----------"
 209      << "------------------------------------------------- \n \n    no    "
 210      << "    id   name            status     mothers   daughters     colou"
 211      << "rs      p_x        p_y        p_z         e          m \n";
 212   if (showScaleAndVertex)
 213     os << "                                    scale                      "
 214        << "                   xProd      yProd      zProd      tProd      "
 215        << " tau\n";
 216
 217   // At high energy switch to scientific format for momenta.
 218   bool useFixed = (entry[0].e() < 1e5);
 219
 220   // Listing of complete event.
 221   Vec4 pSum;
 222   double chargeSum = 0.;
 223   for (int i = 0; i < int(entry.size()); ++i) {
 224     const Particle& pt = entry[i];
 225
 226     // Basic line for a particle, always printed.
 227     os << setw(6) << i << setw(10) << pt.id() << "   " << left
 228        << setw(18) << pt.nameWithStatus(18) << right << setw(4)
 229        << pt.status() << setw(6) << pt.mother1() << setw(6)
 230        << pt.mother2() << setw(6) << pt.daughter1() << setw(6)
 231        << pt.daughter2() << setw(6) << pt.col() << setw(6) << pt.acol()
 232        << ( (useFixed) ? fixed : scientific ) << setprecision(3)
 233        << setw(11) << pt.px() << setw(11) << pt.py() << setw(11)
 234        << pt.pz() << setw(11) << pt.e() << setw(11) << pt.m() << "\n";
 235
 236     // Optional extra line for scale value and production vertex.
 237     if (showScaleAndVertex)
 238       os << "                              " << setw(11) << pt.scale()
 239          << "                                    " << scientific
 240          << setprecision(3) << setw(11) << pt.xProd() << setw(11)
 241          << pt.yProd() << setw(11) << pt.zProd() << setw(11)
 242          << pt.tProd() << setw(11) << pt.tau() << "\n";
 243
 244     // Optional extra line, giving a complete list of mothers and daughters.
 245     if (showMothersAndDaughters) {
 246       int linefill = 2;
 247       os << "                mothers:";
 248       vector<int> allMothers = motherList(i);
 249       for (int j = 0; j < int(allMothers.size()); ++j) {
 250         os << " " <<  allMothers[j];
 251         if (++linefill == IPERLINE) {os << "\n                "; linefill = 0;}
 252       }
 253       os << ";   daughters:";
 254       vector<int> allDaughters = daughterList(i);
 255       for (int j = 0; j < int(allDaughters.size()); ++j) {
 256         os << " " <<  allDaughters[j];
 257         if (++linefill == IPERLINE) {os << "\n                "; linefill = 0;}
 258       }
 259       if (linefill !=0) os << "\n";
 260     }
 261
 262     // Statistics on momentum and charge.
 263     if (entry[i].status() > 0) {
 264       pSum += entry[i].p();
 265       chargeSum += entry[i].charge();
 266     }
 267   }
 268
 269   // Line with sum charge, momentum, energy and invariant mass.
 270   os << fixed << setprecision(3) << "                                   "
 271      << "Charge sum:" << setw(7) << chargeSum << "           Momentum sum:"
 272      << ( (useFixed) ? fixed : scientific ) << setprecision(3)
 273      << setw(11) << pSum.px() << setw(11) << pSum.py() << setw(11)
 274      << pSum.pz() << setw(11) << pSum.e() << setw(11) << pSum.mCalc()
 275      << "\n";
 276
 277   // Listing finished.
 278   os << "\n --------  End PYTHIA Event Listing  ----------------------------"
 279      << "-------------------------------------------------------------------"
 280      << endl;
 281 }
 282
 283 //--------------------------------------------------------------------------
 284
 285 // Find complete list of mothers.
 286
 287 vector<int> Event::motherList(int i) const {
 288
 289   // Vector of all the mothers; created empty.
 290   vector<int> mothers;
 291
 292   // Read out the two official mother indices and status code.
 293   int mother1   = entry[i].mother1();
 294   int mother2   = entry[i].mother2();
 295   int statusAbs = entry[i].statusAbs();
 296
 297   // Special cases in the beginning, where the meaning of zero is unclear.
 298   if  (statusAbs == 11 || statusAbs == 12) ;
 299   else if (mother1 == 0 && mother2 == 0) mothers.push_back(0);
 300
 301   // One mother or a carbon copy
 302   else if (mother2 == 0 || mother2 == mother1) mothers.push_back(mother1);
 303
 304   // A range of mothers from string fragmentation.
 305   else if ( statusAbs > 80 &&  statusAbs < 90)
 306     for (int iRange = mother1; iRange <= mother2; ++iRange)
 307       mothers.push_back(iRange);
 308
 309   // Two separate mothers.
 310   else {
 311     mothers.push_back( min(mother1, mother2) );
 312     mothers.push_back( max(mother1, mother2) );
 313   }
 314
 315   // Done.
 316   return mothers;
 317
 318 }
 319
 320 //--------------------------------------------------------------------------
 321
 322 // Find complete list of daughters.
 323
 324 vector<int> Event::daughterList(int i) const {
 325
 326   // Vector of all the daughters; created empty.
 327   vector<int> daughters;
 328
 329   // Read out the two official daughter indices.
 330   int daughter1 = entry[i].daughter1();
 331   int daughter2 = entry[i].daughter2();
 332
 333   // Simple cases: no or one daughter.
 334   if (daughter1 == 0 && daughter2 == 0) ;
 335   else if (daughter2 == 0 || daughter2 == daughter1)
 336     daughters.push_back(daughter1);
 337
 338   // A range of daughters.
 339   else if (daughter2 > daughter1)
 340     for (int iRange = daughter1; iRange <= daughter2; ++iRange)
 341       daughters.push_back(iRange);
 342
 343   // Two separated daughters.
 344   else {
 345     daughters.push_back(daughter2);
 346     daughters.push_back(daughter1);
 347   }
 348
 349   // Special case for two incoming beams: attach further
 350   // initiators and remnants that have beam as mother.
 351   if (entry[i].statusAbs() == 12 || entry[i].statusAbs() == 13)
 352   for (int iDau = 3; iDau < size(); ++iDau)
 353   if (entry[iDau].mother1() == i) {
 354     bool isIn = false;
 355     for (int iIn = 0; iIn < int(daughters.size()); ++iIn)
 356       if (iDau == daughters[iIn]) isIn = true;
 357     if (!isIn) daughters.push_back(iDau);
 358   }
 359
 360   // Done.
 361   return daughters;
 362
 363 }
 364
 365 //--------------------------------------------------------------------------
 366
 367 // Convert internal Pythia status codes to the HepMC status conventions.
 368
 369 int Event::statusHepMC(int i) const {
 370
 371   // Positive codes are final particles. Status -12 are beam particles.
 372   int statusNow     = entry[i].status();
 373   if (statusNow > 0)    return 1;
 374   if (statusNow == -12) return 4;
 375
 376   // Hadrons, muons, taus that decay normally are status 2.
 377   int idNow         = entry[i].id();
 378   if (entry[i].isHadron() || abs(idNow) == 13 || abs(idNow) == 15) {
 379     int iDau        = entry[i].daughter1();
 380     // Particle should not decay into itself  (e.g. Bose-Einstein).
 381     if ( entry[iDau].id() != idNow) {
 382       int statusDau = entry[ iDau ].statusAbs();
 383       if (statusDau > 90 && statusDau < 95) return 2;
 384     }
 385   }
 386
 387   // Other acceptable negative codes as their positive counterpart.
 388   if (statusNow <= -11 && statusNow >= -200) return -statusNow;
 389
 390   // Unacceptable codes as 0.
 391   return 0;
 392
 393 }
 394
 395 //--------------------------------------------------------------------------
 396
 397 // Trace the first and last copy of one and the same particle.
 398
 399 int Event::iTopCopy( int i) const {
 400
 401   int iUp = i;
 402   while ( iUp > 0 && entry[iUp].mother2() == entry[iUp].mother1()
 403     && entry[iUp].mother1() > 0) iUp = entry[iUp].mother1();
 404   return iUp;
 405
 406 }
 407
 408 int Event::iBotCopy( int i) const {
 409
 410   int iDn = i;
 411   while ( iDn > 0 && entry[iDn].daughter2() == entry[iDn].daughter1()
 412     && entry[iDn].daughter1() > 0) iDn = entry[iDn].daughter1();
 413   return iDn;
 414
 415 }
 416
 417 //--------------------------------------------------------------------------
 418
 419 // Trace the first and last copy of one and the same particle,
 420 // also through shower branchings, making use of flavour matches.
 421 // Stops tracing when this gives ambiguities.
 422
 423 int Event::iTopCopyId( int i) const {
 424
 425   int id = entry[i].id();
 426   int iUp = i;
 427   for ( ; ; ) {
 428     int mother1 = entry[iUp].mother1();
 429     int id1     = (mother1 > 0) ? entry[mother1].id() : 0;
 430     int mother2 = entry[iUp].mother2();
 431     int id2     = (mother2 > 0) ? entry[mother2].id() : 0;
 432     if (mother2 != mother1 && id2 == id1) break;
 433     if (id1 == id) {
 434       iUp = mother1;
 435       continue;
 436     }
 437     if (id2 == id) {
 438       iUp = mother2;
 439       continue;
 440     }
 441     break;
 442   }
 443   return iUp;
 444
 445 }
 446
 447 int Event::iBotCopyId( int i) const {
 448
 449   int id = entry[i].id();
 450   int iDn = i;
 451   for ( ; ; ) {
 452     int daughter1 = entry[iDn].daughter1();
 453     int id1       = (daughter1 > 0) ? entry[daughter1].id() : 0;
 454     int daughter2 = entry[iDn].daughter2();
 455     int id2       = (daughter2 > 0) ? entry[daughter2].id() : 0;
 456     if (daughter2 != daughter1 && id2 == id1) break;
 457     if (id1 == id) {
 458       iDn = daughter1;
 459       continue;
 460     }
 461     if (id2 == id) {
 462       iDn = daughter2;
 463       continue;
 464     }
 465     break;
 466   }
 467   return iDn;
 468
 469 }
 470
 471 //--------------------------------------------------------------------------
 472
 473 // Find complete list of sisters.
 474
 475 vector<int> Event::sisterList(int i) const {
 476
 477   // Vector of all the sisters; created empty.
 478   vector<int> sisters;
 479   if (entry[i].statusAbs() == 11) return sisters;
 480
 481   // Find mother and all its daughters.
 482   int iMother = entry[i].mother1();
 483   vector<int> daughters = daughterList(iMother);
 484
 485   // Copy all daughters, excepting the input particle itself.
 486   for (int j = 0; j < int(daughters.size()); ++j)
 487   if (daughters[j] != i) sisters.push_back( daughters[j] );
 488
 489   // Done.
 490   return sisters;
 491
 492 }
 493
 494 //--------------------------------------------------------------------------
 495
 496 // Find complete list of sisters. Traces up with iTopCopy and
 497 // down with iBotCopy to give sisters at same level of evolution.
 498 // Should this not give any final particles, the search is widened.
 499
 500 vector<int> Event::sisterListTopBot(int i, bool widenSearch) const {
 501
 502   // Vector of all the sisters; created empty.
 503   vector<int> sisters;
 504   if (entry[i].statusAbs() == 11) return sisters;
 505
 506   // Trace up to first copy of current particle.
 507   int iUp = iTopCopy(i);
 508
 509   // Find mother and all its daughters.
 510   int iMother = entry[iUp].mother1();
 511   vector<int> daughters = daughterList(iMother);
 512
 513   // Trace all daughters down, excepting the input particle itself.
 514   for (int jD = 0; jD < int(daughters.size()); ++jD)
 515   if (daughters[jD] != iUp)
 516     sisters.push_back( iBotCopy( daughters[jD] ) );
 517
 518   // Prune any non-final particles from list.
 519   int jP = 0;
 520   while (jP < int(sisters.size())) {
 521     if (entry[sisters[jP]].status() > 0) ++jP;
 522     else {
 523       sisters[jP] = sisters.back();
 524       sisters.pop_back();
 525     }
 526   }
 527
 528   // If empty list then restore immediate daughters.
 529   if (sisters.size() == 0 && widenSearch) {
 530     for (int jR = 0; jR < int(daughters.size()); ++jR)
 531     if (daughters[jR] != iUp)
 532       sisters.push_back( iBotCopy( daughters[jR] ) );
 533
 534     // Then trace all daughters, not only bottom copy.
 535     for (int jT = 0; jT < int(sisters.size()); ++jT) {
 536       daughters = daughterList( sisters[jT] );
 537       for (int k = 0; k < int(daughters.size()); ++k)
 538         sisters.push_back( daughters[k] );
 539     }
 540
 541     // And then prune any non-final particles from list.
 542     int jN = 0;
 543     while (jN < int(sisters.size())) {
 544       if (entry[sisters[jN]].status() > 0) ++jN;
 545       else {
 546         sisters[jN] = sisters.back();
 547         sisters.pop_back();
 548       }
 549     }
 550   }
 551
 552   // Done.
 553   return sisters;
 554
 555 }
 556
 557 //--------------------------------------------------------------------------
 558
 559 // Check whether a given particle is an arbitrarily-steps-removed
 560 // mother to another. For the parton -> hadron transition, only
 561 // first-rank hadrons are associated with the respective end quark.
 562
 563 bool Event::isAncestor(int i, int iAncestor) const {
 564
 565   // Begin loop to trace upwards from the daughter.
 566   int iUp = i;
 567   for ( ; ; ) {
 568
 569     // If positive match then done.
 570     if (iUp == iAncestor) return true;
 571
 572     // If out of range then failed to find match.
 573     if (iUp <= 0 || iUp > size()) return false;
 574
 575     // If unique mother then keep on moving up the chain.
 576     int mother1 = entry[iUp].mother1();
 577     int mother2 = entry[iUp].mother2();
 578     if (mother2 == mother1 || mother2 == 0) {iUp = mother1; continue;}
 579
 580     // If many mothers, except hadronization, then fail tracing.
 581     int status = entry[iUp].statusAbs();
 582     if (status < 81 || status > 86) return false;
 583
 584     // For hadronization step, fail if not first rank, else move up.
 585     if (status == 82) {
 586       iUp = (iUp + 1 < size() && entry[iUp + 1].mother1() == mother1)
 587           ? mother1 : mother2; continue;
 588     }
 589     if (status == 83) {
 590       if (entry[iUp - 1].mother1() == mother1) return false;
 591       iUp = mother1; continue;
 592     }
 593     if (status == 84) {
 594       if (iUp + 1 < size() && entry[iUp + 1].mother1() == mother1)
 595         return false;
 596       iUp = mother1; continue;
 597     }
 598
 599     // Fail for ministring -> one hadron and for junctions.
 600     return false;
 601
 602   }
 603   // End of loop. Should never reach beyond here.
 604   return false;
 605
 606 }
 607
 608 //--------------------------------------------------------------------------
 609
 610 // Erase junction stored in specified slot and move up the ones under.
 611
 612 void Event::eraseJunction(int i) {
 613
 614   for (int j = i; j < int(junction.size()) - 1; ++j)
 615     junction[j] = junction[j + 1];
 616   junction.pop_back();
 617
 618 }
 619
 620 //--------------------------------------------------------------------------
 621
 622 // Print the junctions in an event.
 623
 624 void Event::listJunctions(ostream& os) const {
 625
 626   // Header.
 627   os << "\n --------  PYTHIA Junction Listing  "
 628      << headerList.substr(0,30) << "\n \n    no  kind  col0  col1  col2 "
 629      << "endc0 endc1 endc2 stat0 stat1 stat2\n";
 630
 631   // Loop through junctions in event and list them.
 632   for (int i = 0; i < sizeJunction(); ++i)
 633     os << setw(6) << i << setw(6) << kindJunction(i) << setw(6)
 634        << colJunction(i, 0) << setw(6) << colJunction(i, 1) << setw(6)
 635        << colJunction(i, 2) << setw(6) << endColJunction(i, 0) << setw(6)
 636        << endColJunction(i, 1) << setw(6) << endColJunction(i, 2) << setw(6)
 637        << statusJunction(i, 0) << setw(6) << statusJunction(i, 1) << setw(6)
 638        << statusJunction(i, 2) << "\n";
 639
 640   // Alternative if no junctions. Listing finished.
 641   if (sizeJunction() == 0) os << "    no junctions present \n";
 642   os << "\n --------  End PYTHIA Junction Listing  --------------------"
 643      << "------" << endl;
 644 }
 645
 646 //--------------------------------------------------------------------------
 647
 648 // Operator overloading allows to append one event to an existing one.
 649
 650 Event& Event::operator+=( const Event& addEvent) {
 651
 652   // Find offsets. One less since won't copy line 0.
 653   int offsetIdx = entry.size() - 1;
 654   int offsetCol = maxColTag;
 655
 656   // Add energy to zeroth line and calculate new invariant mass.
 657   entry[0].p( entry[0].p() + addEvent[0].p() );
 658   entry[0].m( entry[0].mCalc() );
 659
 660   // Read out particles from line 1 (not 0) onwards.
 661   Particle temp;
 662   for (int i = 1; i < addEvent.size(); ++i) {
 663     temp = addEvent[i];
 664
 665     // Add offset to nonzero mother, daughter and colour indices.
 666     if (temp.mother1() > 0) temp.mother1( temp.mother1() + offsetIdx );
 667     if (temp.mother2() > 0) temp.mother2( temp.mother2() + offsetIdx );
 668     if (temp.daughter1() > 0) temp.daughter1( temp.daughter1() + offsetIdx );
 669     if (temp.daughter2() > 0) temp.daughter2( temp.daughter2() + offsetIdx );
 670     if (temp.col() > 0) temp.col( temp.col() + offsetCol );
 671     if (temp.acol() > 0) temp.acol( temp.acol() + offsetCol );
 672
 673     // Append particle to summed event.
 674     append( temp );
 675   }
 676
 677   // Read out junctions one by one.
 678   Junction tempJ;
 679   int begCol, endCol;
 680   for (int i = 0; i < addEvent.sizeJunction(); ++i) {
 681     tempJ = addEvent.getJunction(i);
 682
 683     // Add colour offsets to all three legs.
 684     for (int  j = 0; j < 3; ++j) {
 685       begCol = tempJ.col(j);
 686       endCol = tempJ.endCol(j);
 687       if (begCol > 0) begCol += offsetCol;
 688       if (endCol > 0) endCol += offsetCol;
 689       tempJ.cols( j, begCol, endCol);
 690     }
 691
 692     // Append junction to summed event.
 693     appendJunction( tempJ );
 694   }
 695
 696   // Set header that indicates character as sum of events.
 697   headerList = "(combination of several events)  -------";
 698
 699   // Done.
 700   return *this;
 701
 702 }
 703
 704 //==========================================================================
 705
 706 } // end namespace Pythia8