]>
Commit | Line | Data |
---|---|---|
9419eeef | 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 |