StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Event.cc
1 // Event.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, 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 "Pythia8/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 // Set pointer to the particle data species of the particle.
29 
30 void Particle::setPDEPtr(ParticleDataEntry* pdePtrIn) {
31  pdePtr = pdePtrIn; if (pdePtrIn != 0 || evtPtr == 0) return;
32  pdePtr = (*evtPtr).particleDataPtr->particleDataEntryPtr( idSave);}
33 
34 //--------------------------------------------------------------------------
35 
36 // Find out if polarization is (close to) an integer.
37 
38 int Particle::intPol() const {
39 
40  double smallDbls[6] = { 0., 1., -1., 2., -2., 9.};
41  int smallInts[6] = { 0, 1, -1, 2, -2, 9 };
42  for (int iPol = 0; iPol < 6; ++ iPol)
43  if (abs(polSave - smallDbls[iPol]) < 1e-10) return smallInts[iPol];
44  return -9;
45 
46 }
47 
48 //--------------------------------------------------------------------------
49 
50 // Functions for rapidity and pseudorapidity.
51 
52 double Particle::y() const {
53  double temp = log( ( pSave.e() + abs(pSave.pz()) ) / max( TINY, mT() ) );
54  return (pSave.pz() > 0.) ? temp : -temp;
55 }
56 
57 double Particle::eta() const {
58  double temp = log( ( pSave.pAbs() + abs(pSave.pz()) ) / max( TINY, pT() ) );
59  return (pSave.pz() > 0.) ? temp : -temp;
60 }
61 
62 // Rapidity with minimal transverse mass, and after rotation/boost.
63 
64 double Particle::y(double mCut) const {
65  double mTmin = max( mCut, mT() );
66  double eMin = sqrt( pow2(mTmin) + pow2(pSave.pz()) );
67  double temp = log( ( eMin + abs(pSave.pz()) ) / mTmin );
68  return (pSave.pz() > 0.) ? temp : -temp;
69 }
70 
71 double Particle::y(double mCut, RotBstMatrix& M) const {
72  Vec4 pCopy = pSave;
73  pCopy.rotbst(M);
74  double mTmin = max( mCut, sqrt( m2() + pCopy.pT2()) );
75  double eMin = sqrt( pow2(mTmin) + pow2(pCopy.pz()) );
76  double temp = log( ( eMin + abs(pCopy.pz()) ) / mTmin );
77  return (pCopy.pz() > 0.) ? temp : -temp;
78 }
79 
80 //--------------------------------------------------------------------------
81 
82 // Method to find the index of the particle in the event record.
83 
84 int Particle::index() const { if (evtPtr == 0) return -1;
85  return (long(this) - long(&((*evtPtr)[0]))) / sizeof(Particle);
86 }
87 
88 //--------------------------------------------------------------------------
89 
90 // Trace the first and last copy of one and the same particle.
91 
92 int Particle::iTopCopy() const {
93 
94  if (evtPtr == 0) return -1;
95  int iUp = index();
96  while ( iUp > 0 && (*evtPtr)[iUp].mother2() == (*evtPtr)[iUp].mother1()
97  && (*evtPtr)[iUp].mother1() > 0) iUp = (*evtPtr)[iUp].mother1();
98  return iUp;
99 
100 }
101 
102 int Particle::iBotCopy() const {
103 
104  if (evtPtr == 0) return -1;
105  int iDn = index();
106  while ( iDn > 0 && (*evtPtr)[iDn].daughter2() == (*evtPtr)[iDn].daughter1()
107  && (*evtPtr)[iDn].daughter1() > 0) iDn = (*evtPtr)[iDn].daughter1();
108  return iDn;
109 
110 }
111 
112 //--------------------------------------------------------------------------
113 
114 // Trace the first and last copy of one and the same particle,
115 // also through shower branchings, making use of flavour matches.
116 // Stops tracing when this gives ambiguities.
117 
118 int Particle::iTopCopyId( bool simplify) const {
119 
120  // Check that particle belongs to event record. Initialize.
121  if (evtPtr == 0) return -1;
122  int iUp = index();
123 
124  // Simple solution when only first and last mother are studied.
125  if (simplify) for ( ; ; ) {
126  int mother1up = (*evtPtr)[iUp].mother1();
127  int id1up = (mother1up > 0) ? (*evtPtr)[mother1up].id() : 0;
128  int mother2up = (*evtPtr)[iUp].mother2();
129  int id2up = (mother2up > 0) ? (*evtPtr)[mother2up].id() : 0;
130  if (mother2up != mother1up && id2up == id1up) return iUp;
131  if (id1up != idSave && id2up != idSave) return iUp;
132  iUp = (id1up == idSave) ? mother1up : mother2up;
133  }
134 
135  // Else full solution where all mothers are studied.
136  for ( ; ; ) {
137  int iUpTmp = 0;
138  vector<int> mothersTmp = (*evtPtr)[iUp].motherList();
139  for (unsigned int i = 0; i < mothersTmp.size(); ++i)
140  if ( (*evtPtr)[mothersTmp[i]].id() == idSave) {
141  if (iUpTmp != 0) return iUp;
142  iUpTmp = mothersTmp[i];
143  }
144  if (iUpTmp == 0) return iUp;
145  iUp = iUpTmp;
146  }
147 
148 }
149 
150 int Particle::iBotCopyId( bool simplify) const {
151 
152  // Check that particle belongs to event record. Initialize.
153  if (evtPtr == 0) return -1;
154  int iDn = index();
155 
156  // Simple solution when only first and last daughter are studied.
157  if (simplify) for ( ; ; ) {
158  int daughter1dn = (*evtPtr)[iDn].daughter1();
159  int id1dn = (daughter1dn > 0) ? (*evtPtr)[daughter1dn].id() : 0;
160  int daughter2dn = (*evtPtr)[iDn].daughter2();
161  int id2dn = (daughter2dn > 0) ? (*evtPtr)[daughter2dn].id() : 0;
162  if (daughter2dn != daughter1dn && id2dn == id1dn) return iDn;
163  if (id1dn != idSave && id2dn != idSave) return iDn;
164  iDn = (id1dn == idSave) ? daughter1dn : daughter2dn;
165  }
166 
167  // Else full solution where all daughters are studied.
168  for ( ; ; ) {
169  int iDnTmp = 0;
170  vector<int> daughtersTmp = (*evtPtr)[iDn].daughterList();
171  for (unsigned int i = 0; i < daughtersTmp.size(); ++i)
172  if ( (*evtPtr)[daughtersTmp[i]].id() == idSave) {
173  if (iDnTmp != 0) return iDn;
174  iDnTmp = daughtersTmp[i];
175  }
176  if (iDnTmp == 0) return iDn;
177  iDn = iDnTmp;
178  }
179 
180 }
181 
182 //--------------------------------------------------------------------------
183 
184 // Find complete list of mothers.
185 
186 vector<int> Particle::motherList() const {
187 
188  // Vector of all the mothers; created empty. Done if no event pointer.
189  vector<int> motherVec;
190  if (evtPtr == 0) return motherVec;
191 
192  // Special cases in the beginning, where the meaning of zero is unclear.
193  int statusSaveAbs = abs(statusSave);
194  if (statusSaveAbs == 11 || statusSaveAbs == 12) ;
195  else if (mother1Save == 0 && mother2Save == 0) motherVec.push_back(0);
196 
197  // One mother or a carbon copy.
198  else if (mother2Save == 0 || mother2Save == mother1Save)
199  motherVec.push_back(mother1Save);
200 
201  // A range of mothers from string fragmentation.
202  else if ( (statusSaveAbs > 80 && statusSaveAbs < 90)
203  || (statusSaveAbs > 100 && statusSaveAbs < 107) )
204  for (int iRange = mother1Save; iRange <= mother2Save; ++iRange)
205  motherVec.push_back(iRange);
206 
207  // Two separate mothers.
208  else {
209  motherVec.push_back( min(mother1Save, mother2Save) );
210  motherVec.push_back( max(mother1Save, mother2Save) );
211  }
212 
213  // Done.
214  return motherVec;
215 
216 }
217 
218 //--------------------------------------------------------------------------
219 
220 // Find complete list of daughters.
221 
222 vector<int> Particle::daughterList() const {
223 
224  // Vector of all the daughters; created empty. Done if no event pointer.
225  vector<int> daughterVec;
226  if (evtPtr == 0) return daughterVec;
227 
228  // Simple cases: no or one daughter.
229  if (daughter1Save == 0 && daughter2Save == 0) ;
230  else if (daughter2Save == 0 || daughter2Save == daughter1Save)
231  daughterVec.push_back(daughter1Save);
232 
233  // A range of daughters.
234  else if (daughter2Save > daughter1Save)
235  for (int iRange = daughter1Save; iRange <= daughter2Save; ++iRange)
236  daughterVec.push_back(iRange);
237 
238  // Two separated daughters.
239  else {
240  daughterVec.push_back(daughter2Save);
241  daughterVec.push_back(daughter1Save);
242  }
243 
244  // Special case for two incoming beams: attach further
245  // initiators and remnants that have beam as mother.
246  if (abs(statusSave) == 12 || abs(statusSave) == 13) {
247  int i = index();
248  for (int iDau = i + 1; iDau < evtPtr->size(); ++iDau)
249  if ((*evtPtr)[iDau].mother1() == i) {
250  bool isIn = false;
251  for (int iIn = 0; iIn < int(daughterVec.size()); ++iIn)
252  if (iDau == daughterVec[iIn]) isIn = true;
253  if (!isIn) daughterVec.push_back(iDau);
254  }
255  }
256 
257  // Done.
258  return daughterVec;
259 
260 }
261 
262 //--------------------------------------------------------------------------
263 
264 // Find complete list of daughters recursively, i.e. including subsequent
265 // generations. Is intended specifically for resonance decays.
266 
267 vector<int> Particle::daughterListRecursive() const {
268 
269  // Vector of all the daughters; created empty. Done if no event pointer.
270  vector<int> daughterVec;
271  if (evtPtr == 0) return daughterVec;
272 
273  // Find first generation of daughters.
274  daughterVec = daughterList();
275 
276  // Recursively add daughters of unstable particles.
277  int size = daughterVec.size();
278  for (int iDau = 0; iDau < size; ++iDau) {
279  Particle& partNow = (*evtPtr)[daughterVec[iDau]];
280  if (!partNow.isFinal()) {
281  vector<int> grandDauVec = partNow.daughterList();
282  for (int i = 0; i < int(grandDauVec.size()); ++i)
283  daughterVec.push_back( grandDauVec[i] );
284  size += grandDauVec.size();
285  }
286  }
287 
288  // Done.
289  return daughterVec;
290 
291 }
292 
293 //--------------------------------------------------------------------------
294 
295 // Find complete list of sisters. Optionally traces up with iTopCopy
296 // and down with iBotCopy to give sisters at same level of evolution.
297 
298 vector<int> Particle::sisterList(bool traceTopBot) const {
299 
300  // Vector of all the sisters; created empty.
301  vector<int> sisterVec;
302  if (evtPtr == 0 || abs(statusSave) == 11) return sisterVec;
303 
304  // Find all daughters of the mother.
305  int iUp = (traceTopBot) ? iTopCopy() : index();
306  int iMother = (*evtPtr)[iUp].mother1();
307  vector<int> daughterVec = (*evtPtr)[iMother].daughterList();
308 
309  // Copy all daughters, excepting the input particle itself.
310  for (int j = 0; j < int(daughterVec.size()); ++j) {
311  int iDau = daughterVec[j];
312  if (iDau != iUp) {
313  int iDn = (traceTopBot) ? (*evtPtr)[iDau].iBotCopy() : iDau;
314  sisterVec.push_back( iDn);
315  }
316  }
317 
318  // Done.
319  return sisterVec;
320 
321 }
322 
323 //--------------------------------------------------------------------------
324 
325 // Check whether a given particle is an arbitrarily-steps-removed
326 // mother to another. For the parton -> hadron transition, only
327 // first-rank hadrons are associated with the respective end quark.
328 
329 bool Particle::isAncestor(int iAncestor) const {
330 
331  // Begin loop to trace upwards from the daughter.
332  if (evtPtr == 0) return false;
333  int iUp = index();
334  int sizeNow = (*evtPtr).size();
335  for ( ; ; ) {
336 
337  // If positive match then done.
338  if (iUp == iAncestor) return true;
339 
340  // If out of range then failed to find match.
341  if (iUp <= 0 || iUp > sizeNow) return false;
342 
343  // If unique mother then keep on moving up the chain.
344  int mother1up = (*evtPtr)[iUp].mother1();
345  int mother2up = (*evtPtr)[iUp].mother2();
346  if (mother2up == mother1up || mother2up == 0) {iUp = mother1up; continue;}
347 
348  // If many mothers, except hadronization, then fail tracing.
349  int statusUp = (*evtPtr)[iUp].statusAbs();
350  if (statusUp < 81 || statusUp > 86) return false;
351 
352  // For hadronization step, fail if not first rank, else move up.
353  if (statusUp == 82) {
354  iUp = (iUp + 1 < sizeNow && (*evtPtr)[iUp + 1].mother1() == mother1up)
355  ? mother1up : mother2up; continue;
356  }
357  if (statusUp == 83) {
358  if ((*evtPtr)[iUp - 1].mother1() == mother1up) return false;
359  iUp = mother1up; continue;
360  }
361  if (statusUp == 84) {
362  if (iUp + 1 < sizeNow && (*evtPtr)[iUp + 1].mother1() == mother1up)
363  return false;
364  iUp = mother1up; continue;
365  }
366 
367  // Fail for ministring -> one hadron and for junctions.
368  return false;
369 
370  }
371  // End of loop. Should never reach beyond here.
372  return false;
373 
374 }
375 
376 //--------------------------------------------------------------------------
377 
378 // Convert internal Pythia status codes to the HepMC status conventions.
379 
380 int Particle::statusHepMC() const {
381 
382  // Positive codes are final particles. Status -12 are beam particles.
383  if (statusSave > 0) return 1;
384  if (statusSave == -12) return 4;
385  if (evtPtr == 0) return 0;
386 
387  // Hadrons, muons, taus that decay normally are status 2.
388  if (isHadron() || abs(idSave) == 13 || abs(idSave) == 15) {
389  // Particle should not decay into itself (e.g. Bose-Einstein).
390  if ( (*evtPtr)[daughter1Save].id() != idSave) {
391  int statusDau = (*evtPtr)[daughter1Save].statusAbs();
392  if (statusDau > 90 && statusDau < 95) return 2;
393  }
394  }
395 
396  // Other acceptable negative codes as their positive counterpart.
397  if (statusSave <= -11 && statusSave >= -200) return -statusSave;
398 
399  // Unacceptable codes as 0.
400  return 0;
401 
402 }
403 
404 //--------------------------------------------------------------------------
405 
406 // Check if particle belonged to the final state at the Parton Level.
407 
408 bool Particle::isFinalPartonLevel() const {
409 
410  if (index() >= evtPtr->savedPartonLevelSize) return false;
411  if (statusSave > 0) return true;
412  if (daughter1Save >= evtPtr->savedPartonLevelSize) return true;
413  return false;
414 
415 }
416 
417 //--------------------------------------------------------------------------
418 
419 // Recursively remove the decay products of a particle, update it to be
420 // undecayed, and update all mother/daughter indices to be correct.
421 // Warning: assumes that decay chains are nicely ordered.
422 
423 bool Particle::undoDecay() {
424 
425  // Do not remove daughters of a parton, i.e. entry carrying colour.
426  if (evtPtr == 0) return false;
427  int i = index();
428  if (i < 0 || i >= int((*evtPtr).size())) return false;
429  if (colSave != 0 || acolSave != 0) return false;
430 
431  // Find range of daughters to remove.
432  int dau1 = daughter1Save;
433  if (dau1 == 0) return false;
434  int dau2 = daughter2Save;
435  if (dau2 == 0) dau2 = dau1;
436 
437  // Refuse if any of the daughters have other mothers.
438  for (int j = dau1; j <= dau2; ++j) if ((*evtPtr)[j].mother1() != i
439  || ((*evtPtr)[j].mother2() != i && (*evtPtr)[j].mother2() != 0) )
440  return false;
441 
442  // Initialize range arrays for daughters and granddaughters.
443  vector<int> dauBeg, dauEnd;
444  dauBeg.push_back( dau1);
445  dauEnd.push_back( dau2);
446 
447  // Begin recursive search through all decay chains.
448  int iRange = 0;
449  do {
450  for (int j = dauBeg[iRange]; j <= dauEnd[iRange]; ++j)
451  if ((*evtPtr)[j].status() < 0) {
452 
453  // Find new daughter range, if present.
454  dau1 = (*evtPtr)[j].daughter1();
455  if (dau1 == 0) return false;
456  dau2 = (*evtPtr)[j].daughter2();
457  if (dau2 == 0) dau2 = dau1;
458 
459  // Check if the range duplicates or contradicts existing ones.
460  bool isNew = true;
461  for (int iR = 0; iR < int(dauBeg.size()); ++iR) {
462  if (dau1 == dauBeg[iR] && dau2 == dauEnd[iR]) isNew = false;
463  else if (dau1 >= dauBeg[iR] && dau1 <= dauEnd[iR]) return false;
464  else if (dau2 >= dauBeg[iR] && dau2 <= dauEnd[iR]) return false;
465  }
466 
467  // Add new range where relevant. Keep ranges ordered.
468  if (isNew) {
469  dauBeg.push_back( dau1);
470  dauEnd.push_back( dau2);
471  for (int iR = int(dauBeg.size()) - 1; iR > 0; --iR) {
472  if (dauBeg[iR] < dauBeg[iR - 1]) {
473  swap( dauBeg[iR], dauBeg[iR - 1]);
474  swap( dauEnd[iR], dauEnd[iR - 1]);
475  } else break;
476  }
477  }
478 
479  // End of recursive search all decay chains.
480  }
481  } while (++iRange < int(dauBeg.size()));
482 
483  // Join adjacent ranges to reduce number of erase steps.
484  if (int(dauBeg.size()) > 1) {
485  int iRJ = 0;
486  do {
487  if (dauEnd[iRJ] + 1 == dauBeg[iRJ + 1]) {
488  for (int iRB = iRJ + 1; iRB < int(dauBeg.size()) - 1; ++iRB)
489  dauBeg[iRB] = dauBeg[iRB + 1];
490  for (int iRE = iRJ; iRE < int(dauEnd.size()) - 1; ++iRE)
491  dauEnd[iRE] = dauEnd[iRE + 1];
492  dauBeg.pop_back();
493  dauEnd.pop_back();
494  } else ++iRJ;
495  } while (iRJ < int(dauBeg.size()) - 1);
496  }
497 
498  // Iterate over relevant ranges, from bottom up.
499  for (int iR = int(dauBeg.size()) - 1; iR >= 0; --iR) {
500 
501  // Remove daughters in each range and update mother and daughter indices.
502  (*evtPtr).remove( dauBeg[iR], dauEnd[iR]);
503  }
504 
505  // Update mother that has been undecayed.
506  statusSave = abs(statusSave);
507  daughter1Save = 0;
508  daughter2Save = 0;
509 
510  // Done.
511  return true;
512 
513 }
514 
515 //--------------------------------------------------------------------------
516 
517 // Particle name, with status but imposed maximum length -> may truncate.
518 
519 string Particle::nameWithStatus(int maxLen) const {
520 
521  if (pdePtr == 0) return " ";
522  string temp = (statusSave > 0) ? pdePtr->name(idSave)
523  : "(" + pdePtr->name(idSave) + ")";
524  while (int(temp.length()) > maxLen) {
525  // Remove from end, excluding closing bracket and charge.
526  int iRem = temp.find_last_not_of(")+-0");
527  temp.erase(iRem, 1);
528  }
529  return temp;
530 }
531 
532 //--------------------------------------------------------------------------
533 
534 // Add offsets to mother and daughter pointers (must be non-negative).
535 
536 void Particle::offsetHistory( int minMother, int addMother, int minDaughter,
537  int addDaughter) {
538 
539  if (addMother < 0 || addDaughter < 0) return;
540  if ( mother1Save > minMother ) mother1Save += addMother;
541  if ( mother2Save > minMother ) mother2Save += addMother;
542  if (daughter1Save > minDaughter) daughter1Save += addDaughter;
543  if (daughter2Save > minDaughter) daughter2Save += addDaughter;
544 
545 }
546 
547 //--------------------------------------------------------------------------
548 
549 // Add offsets to colour and anticolour (must be positive).
550 
551 void Particle::offsetCol( int addCol) {
552 
553  if (addCol < 0) return;
554  if ( colSave > 0) colSave += addCol;
555  if (acolSave > 0) acolSave += addCol;
556 
557 }
558 
559 //--------------------------------------------------------------------------
560 
561 // Invariant mass of a pair and its square.
562 // (Not part of class proper, but tightly linked.)
563 
564 double m(const Particle& pp1, const Particle& pp2) {
565  double m2 = pow2(pp1.e() + pp2.e()) - pow2(pp1.px() + pp2.px())
566  - pow2(pp1.py() + pp2.py()) - pow2(pp1.pz() + pp2.pz());
567  return (m2 > 0. ? sqrt(m2) : 0.);
568 }
569 
570 double m2(const Particle& pp1, const Particle& pp2) {
571  double m2 = pow2(pp1.e() + pp2.e()) - pow2(pp1.px() + pp2.px())
572  - pow2(pp1.py() + pp2.py()) - pow2(pp1.pz() + pp2.pz());
573  return m2;
574 }
575 
576 //==========================================================================
577 
578 // Event class.
579 // This class holds info on the complete event record.
580 
581 //--------------------------------------------------------------------------
582 
583 // Constants: could be changed here if desired, but normally should not.
584 // These are of technical nature, as described for each.
585 
586 // Maxmimum number of mothers or daughter indices per line in listing.
587 const int Event::IPERLINE = 20;
588 
589 //--------------------------------------------------------------------------
590 
591 // Copy all information from one event record to another.
592 
593 Event& Event::operator=( const Event& oldEvent) {
594 
595  // Do not copy if same.
596  if (this != &oldEvent) {
597 
598  // Reset all current info in the event.
599  clear();
600 
601  // Copy particle data table; needed for individual particles.
602  particleDataPtr = oldEvent.particleDataPtr;
603 
604  // Copy all the particles one by one.
605  maxColTag = 100;
606  for (int i = 0; i < oldEvent.size(); ++i) append( oldEvent[i] );
607 
608  // Copy all the junctions one by one.
609  for (int i = 0; i < oldEvent.sizeJunction(); ++i)
610  appendJunction( oldEvent.getJunction(i) );
611 
612  // Copy all other values.
613  startColTag = oldEvent.startColTag;
614  maxColTag = oldEvent.maxColTag;
615  savedSize = oldEvent.savedSize;
616  savedJunctionSize = oldEvent.savedJunctionSize;
617  scaleSave = oldEvent.scaleSave;
618  scaleSecondSave = oldEvent.scaleSecondSave;
619  headerList = oldEvent.headerList;
620 
621  // Done.
622  }
623  return *this;
624 
625 }
626 
627 //--------------------------------------------------------------------------
628 
629 // Add a copy of an existing particle at the end of the event record;
630 // return index. Three cases, depending on sign of new status code:
631 // Positive: copy is viewed as daughter, status of original is negated.
632 // Negative: copy is viewed as mother, status of original is unchanged.
633 // Zero: the new is a perfect carbon copy (maybe to be changed later).
634 
635 int Event::copy(int iCopy, int newStatus) {
636 
637  // Protect against attempt to copy negative entries (e.g., junction codes)
638  // or entries beyond end of record.
639  if (iCopy < 0 || iCopy >= size()) return -1;
640 
641  // Simple carbon copy.
642  entry.push_back(entry[iCopy]);
643  int iNew = entry.size() - 1;
644 
645  // Set up to make new daughter of old.
646  if (newStatus > 0) {
647  entry[iCopy].daughters(iNew,iNew);
648  entry[iCopy].statusNeg();
649  entry[iNew].mothers(iCopy, iCopy);
650  entry[iNew].status(newStatus);
651 
652  // Set up to make new mother of old.
653  } else if (newStatus < 0) {
654  entry[iCopy].mothers(iNew,iNew);
655  entry[iNew].daughters(iCopy, iCopy);
656  entry[iNew].status(newStatus);
657  }
658 
659  // Done.
660  return iNew;
661 
662 }
663 
664 //--------------------------------------------------------------------------
665 
666 // Remove entries from iFirst to iLast, including endpoints, and fix history.
667 // (To the extent possible; history pointers in removed range are zeroed.)
668 
669 void Event::remove(int iFirst, int iLast, bool shiftHistory) {
670 
671  // Check that removal range is sensible.
672  if (iFirst < 0 || iLast >= int(entry.size()) || iLast < iFirst) return;
673  int nRem = iLast + 1 - iFirst;
674 
675  // Remove the entries.
676  entry.erase( entry.begin() + iFirst, entry.begin() + iLast + 1);
677 
678  // Loop over remaining particles; read out mothers and daughters.
679  if (shiftHistory) for (int i = 0; i < int(entry.size()); ++i) {
680  int iMot1 = entry[i].mother1();
681  int iMot2 = entry[i].mother2();
682  int iDau1 = entry[i].daughter1();
683  int iDau2 = entry[i].daughter2();
684 
685  // Shift mother and daughter indices according to removed number.
686  // Set zero if in removed range.
687  if (iMot1 > iLast) iMot1 -= nRem;
688  else if (iMot1 >= iFirst) iMot1 = 0;
689  if (iMot2 > iLast) iMot2 -= nRem;
690  else if (iMot2 >= iFirst) iMot2 = 0;
691  if (iDau1 > iLast) iDau1 -= nRem;
692  else if (iDau1 >= iFirst) iDau1 = 0;
693  if (iDau2 > iLast) iDau2 -= nRem;
694  else if (iDau2 >= iFirst) iDau2 = 0;
695 
696  // Set the new values.
697  entry[i].mothers(iMot1, iMot2);
698  entry[i].daughters(iDau1, iDau2);
699  }
700 
701 }
702 
703 //--------------------------------------------------------------------------
704 
705 // Print an event.
706 
707 void Event::list(bool showScaleAndVertex, bool showMothersAndDaughters,
708  int precision) const {
709 
710  // Header.
711  cout << "\n -------- PYTHIA Event Listing " << headerList << "----------"
712  << "-------------------------------------------------\n \n no "
713  << " id name status mothers daughters colou"
714  << "rs p_x p_y p_z e m \n";
715  if (showScaleAndVertex)
716  cout << " scale pol "
717  << " xProd yProd zProd tProd "
718  << " tau\n";
719 
720  // Precision. At high energy switch to scientific format for momenta.
721  int prec = max( 3, precision);
722  bool useFixed = (entry.empty() || entry[0].e() < 1e5);
723 
724  // Listing of complete event.
725  Vec4 pSum;
726  double chargeSum = 0.;
727  for (int i = 0; i < int(entry.size()); ++i) {
728  const Particle& pt = entry[i];
729 
730  // Basic line for a particle, always printed.
731  cout << setw(6) << i << setw(11) << pt.id() << " " << left
732  << setw(18) << pt.nameWithStatus(18) << right << setw(4)
733  << pt.status() << setw(6) << pt.mother1() << setw(6)
734  << pt.mother2() << setw(6) << pt.daughter1() << setw(6)
735  << pt.daughter2() << setw(6) << pt.col() << setw(6) << pt.acol()
736  << ( (useFixed) ? fixed : scientific ) << setprecision(prec)
737  << setw(8+prec) << pt.px() << setw(8+prec) << pt.py()
738  << setw(8+prec) << pt.pz() << setw(8+prec) << pt.e()
739  << setw(8+prec) << pt.m() << "\n";
740 
741  // Optional extra line for scale value, polarization and production vertex.
742  if (showScaleAndVertex)
743  cout << " " << setw(8+prec) << pt.scale()
744  << " " << fixed << setprecision(prec) << setw(8+prec) << pt.pol()
745  << " " << scientific << setprecision(prec)
746  << setw(8+prec) << pt.xProd() << setw(8+prec) << pt.yProd()
747  << setw(8+prec) << pt.zProd() << setw(8+prec) << pt.tProd()
748  << setw(8+prec) << pt.tau() << "\n";
749 
750  // Optional extra line, giving a complete list of mothers and daughters.
751  if (showMothersAndDaughters) {
752  int linefill = 2;
753  cout << " mothers:";
754  vector<int> allMothers = pt.motherList();
755  for (int j = 0; j < int(allMothers.size()); ++j) {
756  cout << " " << allMothers[j];
757  if (++linefill == IPERLINE) {
758  cout << "\n ";
759  linefill = 0;
760  }
761  }
762  cout << "; daughters:";
763  vector<int> allDaughters = pt.daughterList();
764  for (int j = 0; j < int(allDaughters.size()); ++j) {
765  cout << " " << allDaughters[j];
766  if (++linefill == IPERLINE) {
767  cout << "\n ";
768  linefill = 0;
769  }
770  }
771  if (linefill !=0) cout << "\n";
772  }
773 
774  // Extra blank separation line when each particle spans more than one line.
775  if (showScaleAndVertex || showMothersAndDaughters) cout << "\n";
776 
777  // Statistics on momentum and charge.
778  if (entry[i].status() > 0) {
779  pSum += entry[i].p();
780  chargeSum += entry[i].charge();
781  }
782  }
783 
784  // Line with sum charge, momentum, energy and invariant mass.
785  cout << fixed << setprecision(3) << " "
786  << "Charge sum:" << setw(7) << chargeSum << " Momentum sum:"
787  << ( (useFixed) ? fixed : scientific ) << setprecision(prec)
788  << setw(8+prec) << pSum.px() << setw(8+prec) << pSum.py()
789  << setw(8+prec) << pSum.pz() << setw(8+prec) << pSum.e()
790  << setw(8+prec) << pSum.mCalc() << "\n";
791 
792  // Listing finished.
793  cout << "\n -------- End PYTHIA Event Listing ----------------------------"
794  << "-------------------------------------------------------------------"
795  << endl;
796 }
797 
798 //--------------------------------------------------------------------------
799 
800 // Erase junction stored in specified slot and move up the ones under.
801 
802 void Event::eraseJunction(int i) {
803 
804  for (int j = i; j < int(junction.size()) - 1; ++j)
805  junction[j] = junction[j + 1];
806  junction.pop_back();
807 
808 }
809 
810 //--------------------------------------------------------------------------
811 
812 // Print the junctions in an event.
813 
814 void Event::listJunctions() const {
815 
816  // Header.
817  cout << "\n -------- PYTHIA Junction Listing "
818  << headerList.substr(0,30) << "\n \n no kind col0 col1 col2 "
819  << "endc0 endc1 endc2 stat0 stat1 stat2\n";
820 
821  // Loop through junctions in event and list them.
822  for (int i = 0; i < sizeJunction(); ++i)
823  cout << setw(6) << i << setw(6) << kindJunction(i) << setw(6)
824  << colJunction(i, 0) << setw(6) << colJunction(i, 1) << setw(6)
825  << colJunction(i, 2) << setw(6) << endColJunction(i, 0) << setw(6)
826  << endColJunction(i, 1) << setw(6) << endColJunction(i, 2) << setw(6)
827  << statusJunction(i, 0) << setw(6) << statusJunction(i, 1) << setw(6)
828  << statusJunction(i, 2) << "\n";
829 
830  // Alternative if no junctions. Listing finished.
831  if (sizeJunction() == 0) cout << " no junctions present \n";
832  cout << "\n -------- End PYTHIA Junction Listing --------------------"
833  << "------" << endl;
834 }
835 
836 //--------------------------------------------------------------------------
837 
838 // Operator overloading allows to append one event to an existing one.
839 
840 Event& Event::operator+=( const Event& addEvent) {
841 
842  // Find offsets. One less since won't copy line 0.
843  int offsetIdx = entry.size() - 1;
844  int offsetCol = maxColTag;
845 
846  // Add energy to zeroth line and calculate new invariant mass.
847  entry[0].p( entry[0].p() + addEvent[0].p() );
848  entry[0].m( entry[0].mCalc() );
849 
850  // Read out particles from line 1 (not 0) onwards.
851  Particle temp;
852  for (int i = 1; i < addEvent.size(); ++i) {
853  temp = addEvent[i];
854 
855  // Add offset to nonzero mother, daughter and colour indices.
856  if (temp.mother1() > 0) temp.mother1( temp.mother1() + offsetIdx );
857  if (temp.mother2() > 0) temp.mother2( temp.mother2() + offsetIdx );
858  if (temp.daughter1() > 0) temp.daughter1( temp.daughter1() + offsetIdx );
859  if (temp.daughter2() > 0) temp.daughter2( temp.daughter2() + offsetIdx );
860  if (temp.col() > 0) temp.col( temp.col() + offsetCol );
861  if (temp.acol() > 0) temp.acol( temp.acol() + offsetCol );
862 
863  // Append particle to summed event.
864  append( temp );
865  }
866 
867  // Read out junctions one by one.
868  Junction tempJ;
869  int begCol, endCol;
870  for (int i = 0; i < addEvent.sizeJunction(); ++i) {
871  tempJ = addEvent.getJunction(i);
872 
873  // Add colour offsets to all three legs.
874  for (int j = 0; j < 3; ++j) {
875  begCol = tempJ.col(j);
876  endCol = tempJ.endCol(j);
877  if (begCol > 0) begCol += offsetCol;
878  if (endCol > 0) endCol += offsetCol;
879  tempJ.cols( j, begCol, endCol);
880  }
881 
882  // Append junction to summed event.
883  appendJunction( tempJ );
884  }
885 
886  // Set header that indicates character as sum of events.
887  headerList = "(combination of several events) -------";
888 
889  // Done.
890  return *this;
891 
892 }
893 
894 //==========================================================================
895 
896 } // end namespace Pythia8
Definition: AgUStep.h:26