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