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