PYTHIA8/pythia8130/src/Event.cxx

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