StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
EvtGen.h
1 // EvtGen.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2018 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 // Author: Philip Ilten.
6 
7 // This file contains an EvtGen interface. HepMC and EvtGen must be enabled.
8 
9 #ifndef Pythia8_EvtGen_H
10 #define Pythia8_EvtGen_H
11 
12 #include "Pythia8/Pythia.h"
13 #include "EvtGen/EvtGen.hh"
14 #include "EvtGenBase/EvtRandomEngine.hh"
15 #include "EvtGenBase/EvtParticle.hh"
16 #include "EvtGenBase/EvtParticleFactory.hh"
17 #include "EvtGenBase/EvtPDL.hh"
18 #include "EvtGenBase/EvtDecayTable.hh"
19 #include "EvtGenBase/EvtParticleDecayList.hh"
20 #include "EvtGenBase/EvtDecayBase.hh"
21 #include "EvtGenExternal/EvtExternalGenList.hh"
22 
23 using namespace Pythia8;
24 
25 //==========================================================================
26 
27 // A class to wrap the Pythia random number generator for use by EvtGen.
28 
29 class EvtGenRandom : public EvtRandomEngine {
30 
31 public:
32 
33  // Constructor.
34  EvtGenRandom(Rndm *rndmPtrIn) {rndmPtr = rndmPtrIn;}
35 
36  // Return a random number.
37  double random() {if (rndmPtr) return rndmPtr->flat(); else return -1.0;}
38 
39  // The random number pointer.
40  Rndm *rndmPtr;
41 
42 };
43 
44 //==========================================================================
45 
46 // A class to perform decays via the external EvtGen decay program,
47 // see http://evtgen.warwick.ac.uk/, the program manual provided with
48 // the EvtGen distribution, and D. J. Lange,
49 // Nucl. Instrum. Meth. A462, 152 (2001) for details.
50 
51 // EvtGen performs a series of decays from some initial particle
52 // decay, rather than just a single decay, and so EvtGen cannot be
53 // interfaced through the standard external DecayHandler class without
54 // considerable complication. Consequently, EvtGen is called on the
55 // complete event record after all steps of Pythia are completed.
56 
57 // Oftentimes a specific "signal" decay is needed to occur once in an
58 // event, and all other decays performed normally. This is possible
59 // via reading in a user decay file (with readDecayFile) and creating
60 // aliased particles with names ending with signalSuffix. By default,
61 // this is "_SIGNAL". When decay() is called, all particles in the
62 // Pythia event record that are of the same types as the signal
63 // particles are collected. One is selected at random and decayed via
64 // the channel(s) defined for that aliased signal particle. All other
65 // particles are decayed normally. The weight for the event is
66 // calculated and returned.
67 
68 // It is also possible to specify a status needed to consider a
69 // particle as a signal candidate. This can be done by modifying the
70 // signals map, e.g. if the tau- is a signal candidate, then
71 // EvtGenDecays.signals[15].status = 201
72 // will only only select as candidates any tau- with this status. This
73 // allows the event record to be changed before decays, so only
74 // certain particles are selected as possible signal candidates
75 // (e.g. passing kinematic requirements).
76 
77 // Please note that particles produced from a signal candidate decay
78 // are not searched for additional signal candidates. This means that
79 // if B0 and tau- have been designated as signal, then a tau- from a
80 // W- decay would be a signal candidate, while a tau- from a B0 decay
81 // would not. This restriction arises from the additional complexity
82 // of allowing recursive signal decays. The following statuses are
83 // used: 93 for particles decayed with EvtGen, 94 for particles
84 // oscillated with EvtGen, 95 for signal particles, and 96 for signal
85 // particles from an oscillation.
86 
87 class EvtGenDecays {
88 
89 public:
90 
91  // Constructor.
92  EvtGenDecays(Pythia *pythiaPtrIn, string decayFile, string particleDataFile,
93  EvtExternalGenList *extPtrIn = 0, EvtAbsRadCorr *fsrPtrIn = 0,
94  int mixing = 1, bool xml = false, bool limit = true,
95  bool extUse = true, bool fsrUse = true);
96 
97  // Destructor.
98  ~EvtGenDecays() {
99  if (evtgen) delete evtgen;
100  if (extOwner && extPtr) delete extPtr;
101  if (fsrOwner && fsrPtr) delete fsrPtr;
102  }
103 
104  // Perform all decays and return the event weight.
105  double decay();
106 
107  // Stop EvtGen decaying a particle.
108  void exclude(int id) {excIds.insert(id);}
109 
110  // Update the Pythia particle database from EvtGen.
111  void updatePythia();
112 
113  // Update the EvtGen particle database from Pythia.
114  void updateEvtGen();
115 
116  // Read an EvtGen user decay file.
117  void readDecayFile(string decayFile, bool xml = false) {
118  evtgen->readUDecay(decayFile.c_str(), xml); updateData();}
119 
120  // External model pointer and FSR model pointer.
121  bool extOwner, fsrOwner;
122  EvtExternalGenList *extPtr;
123  EvtAbsRadCorr *fsrPtr;
124  std::list<EvtDecayBase*> models;
125 
126  // Map of signal particle info.
127  struct Signal {int status; EvtId egId; vector<double> bfs; vector<int> map;
128  EvtParticleDecayList modes;};
129  map<int, Signal> signals;
130 
131  // The suffix indicating an EvtGen particle or alias is signal.
132  string signalSuffix;
133 
134 protected:
135 
136  // Update the particles to decay with EvtGen, and the selected signals.
137  void updateData(bool final = false);
138 
139  // Update the Pythia event record with an EvtGen decay tree.
140  void updateEvent(Particle *pyPro, EvtParticle *egPro,
141  vector<int> *pySigs = 0, vector<EvtParticle*> *egSigs = 0,
142  vector<double> *bfs = 0, double *wgt = 0);
143 
144  // Check if a particle can decay.
145  bool checkVertex(Particle *pyPro);
146 
147  // Check if a particle is signal.
148  bool checkSignal(Particle *pyPro);
149 
150  // Check if an EvtGen particle has oscillated.
151  bool checkOsc(EvtParticle *egPro);
152 
153  // Number of times to try a decay sampling (constant).
154  static const int NTRYDECAY = 1000;
155 
156  // The pointer to the associated Pythia object.
157  Pythia *pythiaPtr;
158 
159  // Random number wrapper for EvtGen.
160  EvtGenRandom rndm;
161 
162  // The EvtGen object.
163  EvtGen *evtgen;
164 
165  // Set of particle IDs to include and exclude decays with EvtGen.
166  set<int> incIds, excIds;
167 
168  // Flag whether the final particle update has been performed.
169  bool updated;
170 
171  // The selected signal iterator.
172  map<int, Signal>::iterator signal;
173 
174  // Parameters used to check if a particle should decay (as set via Pythia).
175  double tau0Max, tauMax, rMax, xyMax, zMax;
176  bool limitTau0, limitTau, limitRadius, limitCylinder, limitDecay;
177 
178 };
179 
180 //--------------------------------------------------------------------------
181 
182 // The constructor.
183 
184 // The EvtGenDecays object is associated with a single Pythia
185 // instance. This is to ensure a consistent random number generator
186 // across the two, as well as any updates to particle data, etc. Note
187 // that if multiple EvtGenDecays objects exist, that they will modify
188 // one anothers particle databases due to the design of EvtGen.
189 
190 // This constructor also sets all particles to be decayed by EvtGen as
191 // stable within Pythia. The parameters within Pythia used to check if
192 // a particle should be decayed, as described in the "Particle Decays"
193 // section of the Pythia manual, are set. Note that if the variable
194 // "limit" is set to "false", then no check will be made before
195 // decaying a particle with EvtGen.
196 
197 // The constructor is designed to have the exact same form as the
198 // EvtGen constructor except for these five differences.
199 // (1) The first variable is the pointer to the Pythia object.
200 // (2) The third last argument is a flag to limit decays based on the
201 // Pythia criteria (based on the particle decay vertex).
202 // (3) The second last argument is a flag if external models should be
203 // passed to EvtGen (default is true).
204 // (4) The last argument is a flag if an FSR model should be passed
205 // to EvtGen (default is true).
206 // (5) No random engine pointer is passed, as this is obtained from
207 // Pythia.
208 
209 // pythiaPtrIn: the pointer to the associated Pythia object.
210 // decayFile: the name of the decay file to pass to EvtGen.
211 // particleDataFile: the name of the particle data file to pass to EvtGen.
212 // extPtrIn: the optional EvtExternalGenList pointer, this must be
213 // be provided if fsrPtrIn is provided to avoid double
214 // initializations.
215 // fsrPtrIn: the EvtAbsRadCorr pointer to pass to EvtGen.
216 // mixing: the mixing type to pass to EvtGen.
217 // xml: flag to use XML files to pass to EvtGen.
218 // limit: flag to limit particle decays based on Pythia criteria.
219 // extUse: flag to use external models with EvtGen.
220 // fsrUse: flag to use radiative correction engine with EvtGen.
221 
222 EvtGenDecays::EvtGenDecays(Pythia *pythiaPtrIn, string decayFile,
223  string particleDataFile, EvtExternalGenList *extPtrIn,
224  EvtAbsRadCorr *fsrPtrIn, int mixing, bool xml, bool limit,
225  bool extUse, bool fsrUse) : extPtr(extPtrIn), fsrPtr(fsrPtrIn),
226  signalSuffix("_SIGNAL"), pythiaPtr(pythiaPtrIn), rndm(&pythiaPtr->rndm),
227  updated(false) {
228 
229  // Initialize EvtGen.
230  if (!extPtr && fsrPtr) {
231  if (pythiaPtr) pythiaPtr->info.errorMsg("Error in EvtGenDecays::"
232  "EvtGenDecays: extPtr is null but fsrPtr is provided");
233  return;
234  }
235  if (extPtr) extOwner = false;
236  else {extOwner = true; extPtr = new EvtExternalGenList();}
237  if (fsrPtr) fsrOwner = false;
238  else {fsrOwner = true; fsrPtr = extPtr->getPhotosModel();}
239  models = extPtr->getListOfModels();
240  evtgen = new EvtGen(decayFile.c_str(), particleDataFile.c_str(),
241  &rndm, fsrUse ? fsrPtr : 0, extUse ? &models : 0, mixing, xml);
242 
243  // Get the Pythia decay limits.
244  if (!pythiaPtr) return;
245  limitTau0 = pythiaPtr->settings.flag("ParticleDecays:limitTau0");
246  tau0Max = pythiaPtr->settings.parm("ParticleDecays:tau0Max");
247  limitTau = pythiaPtr->settings.flag("ParticleDecays:limitTau");
248  tauMax = pythiaPtr->settings.parm("ParticleDecays:tauMax");
249  limitRadius = pythiaPtr->settings.flag("ParticleDecays:limitRadius");
250  rMax = pythiaPtr->settings.parm("ParticleDecays:rMax");
251  limitCylinder = pythiaPtr->settings.flag("ParticleDecays:limitCylinder");
252  xyMax = pythiaPtr->settings.parm("ParticleDecays:xyMax");
253  zMax = pythiaPtr->settings.parm("ParticleDecays:zMax");
254  limitDecay = limit && (limitTau0 || limitTau ||
255  limitRadius || limitCylinder);
256 
257 }
258 
259 //--------------------------------------------------------------------------
260 
261 // Perform all decays and return the event weight.
262 
263 // All particles in the event record that can be decayed by EvtGen are
264 // decayed. If a particle is a signal particle, then this is stored in
265 // a vector of signal particles. A signal particle is only stored if
266 // its status is the same as the status provided in the signals map. A
267 // negative status in the signal map indicates that all statuses
268 // should be accepted. After all signal particles are identified, one
269 // is randomly chosen and decayed as signal. The remainder are decayed
270 // normally.
271 
272 // Forcing a signal decay changes the weight of an event from unity,
273 // and so the relative event weight is returned, given the forced
274 // signal decay. A weight of 0 indicates no signal in the event, while
275 // a weight of -1 indicates something is wrong, e.g. either the Pythia
276 // or EvtGen pointers are not available or the number of tries has
277 // been exceeded. For the event weight to be valid, one should not
278 // change the absolute branching fractions in the signal and inclusive
279 // definitions, but rather just remove the unwanted decay channels
280 // from the signal decay definition.
281 
282 double EvtGenDecays::decay() {
283 
284  // Reset the signal and signal counters.
285  if (!pythiaPtr || !evtgen) return -1;
286  if (!updated) updateData(true);
287 
288  // Loop over all particles in the Pythia event.
289  Event &event = pythiaPtr->event;
290  vector<int> pySigs; vector<EvtParticle*> egSigs, egPrts;
291  vector<double> bfs; double wgt(1.);
292  for (int iPro = 0; iPro < event.size(); ++iPro) {
293 
294  // Check particle is final and can be decayed by EvtGen.
295  Particle *pyPro = &event[iPro];
296  if (!pyPro->isFinal()) continue;
297  if (incIds.find(pyPro->id()) == incIds.end()) continue;
298  if (pyPro->status() == 93 || pyPro->status() == 94) continue;
299 
300  // Decay the progenitor with EvtGen.
301  EvtParticle *egPro = EvtParticleFactory::particleFactory
302  (EvtPDL::evtIdFromStdHep(pyPro->id()),
303  EvtVector4R(pyPro->e(), pyPro->px(), pyPro->py(), pyPro->pz()));
304  egPrts.push_back(egPro);
305  egPro->setDiagonalSpinDensity();
306  evtgen->generateDecay(egPro);
307  pyPro->tau(egPro->getLifetime());
308  if (!checkVertex(pyPro)) continue;
309 
310  // Add oscillations to event record.
311  if (checkOsc(egPro))
312  updateEvent(pyPro, egPro, &pySigs, &egSigs, &bfs, &wgt);
313 
314  // Undo decay if signal (duplicate to stop oscillations).
315  else if (checkSignal(pyPro)) {
316  pySigs.push_back(pyPro->index());
317  egSigs.push_back(egPro);
318  bfs.push_back(signal->second.bfs[0]);
319  wgt *= 1 - bfs.back();
320  egPro->deleteDaughters();
321  EvtParticle *egDau = EvtParticleFactory::particleFactory
322  (EvtPDL::evtIdFromStdHep(pyPro->id()),
323  EvtVector4R(pyPro->e(), pyPro->px(), pyPro->py(), pyPro->pz()));
324  egDau->addDaug(egPro);
325  egDau->setDiagonalSpinDensity();
326 
327  // If not signal, add to event record.
328  } else updateEvent(pyPro, egPro, &pySigs, &egSigs, &bfs, &wgt);
329  }
330  if (pySigs.size() == 0) {
331  for (int iPrt = 0; iPrt < (int)egPrts.size(); ++iPrt)
332  egPrts[iPrt]->deleteTree();
333  return 0;
334  }
335 
336  // Determine the decays of the signal particles (signal or background).
337  vector<int> modes; int force(-1), n(0);
338  for (int iTry = 1; iTry <= NTRYDECAY; ++iTry) {
339  modes.clear(); force = pythiaPtr->rndm.pick(bfs); n = 0;
340  for (int iSig = 0; iSig < (int)pySigs.size(); ++iSig) {
341  if (iSig == force) modes.push_back(0);
342  else modes.push_back(pythiaPtr->rndm.flat() > bfs[iSig]);
343  if (modes.back() == 0) ++n;
344  }
345  if (pythiaPtr->rndm.flat() <= 1.0 / n) break;
346  if (iTry == NTRYDECAY) {
347  wgt = 2.;
348  pythiaPtr->info.errorMsg("Warning in EvtGenDecays::decay: maximum "
349  "number of signal decay attempts exceeded");
350  }
351  }
352 
353  // Decay the signal particles and mark forced decay.
354  for (int iSig = 0; iSig < (int)pySigs.size(); ++iSig) {
355  EvtParticle *egSig = egSigs[iSig];
356  Particle *pySig = &event[pySigs[iSig]];
357  signal = signals.find(pySig->id());
358  if (egSig->getNDaug() > 0) egSig = egSig->getDaug(0);
359  if (modes[iSig] == 0) egSig->setId(signal->second.egId);
360  else {
361  signal->second.modes.getDecayModel(egSig);
362  egSig->setChannel(signal->second.map[egSig->getChannel()]);
363  }
364  if (iSig == force)
365  pySig->status(pySig->status() == 92 || pySig->status() == 94 ? 96 : 95);
366  evtgen->generateDecay(egSig);
367  updateEvent(pySig, egSigs[iSig]);
368  }
369 
370  // Delete all EvtGen particles and return weight.
371  for (int iPrt = 0; iPrt < (int)egPrts.size(); ++iPrt)
372  egPrts[iPrt]->deleteTree();
373  return 1. - wgt;
374 
375 }
376 
377 //--------------------------------------------------------------------------
378 
379 // Update the Pythia particle database from EvtGen.
380 
381 // Note that only the particle spin type, charge type, nominal mass,
382 // width, minimum mass, maximum mass, and nominal lifetime are
383 // set. The name string is not set.
384 
385 void EvtGenDecays::updatePythia() {
386  if (!pythiaPtr || !evtgen) return;
387  for (int entry = 0; entry < (int)EvtPDL::entries(); ++entry) {
388  EvtId egId = EvtPDL::getEntry(entry);
389  int pyId = EvtPDL::getStdHep(egId);
390  pythiaPtr->particleData.spinType (pyId, EvtPDL::getSpinType(egId));
391  pythiaPtr->particleData.chargeType(pyId, EvtPDL::chg3(egId));
392  pythiaPtr->particleData.m0 (pyId, EvtPDL::getMass(egId));
393  pythiaPtr->particleData.mWidth (pyId, EvtPDL::getWidth(egId));
394  pythiaPtr->particleData.mMin (pyId, EvtPDL::getMinMass(egId));
395  pythiaPtr->particleData.mMax (pyId, EvtPDL::getMaxMass(egId));
396  pythiaPtr->particleData.tau0 (pyId, EvtPDL::getctau(egId));
397  }
398 }
399 
400 //--------------------------------------------------------------------------
401 
402 // Update the EvtGen particle database from Pythia.
403 
404 // The particle update database between Pythia and EvtGen is not
405 // symmetric. Particularly, it is not possible to update the spin
406 // type, charge, or nominal lifetime in the EvtGen particle database.
407 
408 void EvtGenDecays::updateEvtGen() {
409  if (!pythiaPtr || !evtgen) return;
410  int pyId = pythiaPtr->particleData.nextId(1);
411  while (pyId != 0) {
412  EvtId egId = EvtPDL::evtIdFromStdHep(pyId);
413  EvtPDL::reSetMass (egId, pythiaPtr->particleData.m0(pyId));
414  EvtPDL::reSetWidth (egId, pythiaPtr->particleData.mWidth(pyId));
415  EvtPDL::reSetMassMin(egId, pythiaPtr->particleData.mMin(pyId));
416  EvtPDL::reSetMassMax(egId, pythiaPtr->particleData.mMax(pyId));
417  pyId = pythiaPtr->particleData.nextId(pyId);
418  }
419 }
420 
421 //--------------------------------------------------------------------------
422 
423 // Update the particles to decay with EvtGen, and the selected signals.
424 
425 // If final is false, then only signals are initialized in the signal
426 // map. Any particle or alias that ends with signalSuffix is taken as
427 // a signal particle. If final is true all particle entries in EvtGen
428 // are checked to see if they should be set stable in Pythia. If an
429 // EvtGen particle has no decay modes, then Pythia is still allowed to
430 // decay the particle. Additionally, the signal decay channels are
431 // turned off for the non-aliased signal particle.
432 
433 void EvtGenDecays::updateData(bool final) {
434 
435  // Loop over the EvtGen entries.
436  if (!pythiaPtr) return;
437  EvtDecayTable *egTable = EvtDecayTable::getInstance();
438  if (!egTable) return;
439  for (int iEntry = 0; iEntry < (int)EvtPDL::entries(); ++iEntry) {
440  EvtId egId = EvtPDL::getEntry(iEntry);
441  int pyId = EvtPDL::getStdHep(egId);
442  if (egTable->getNModes(egId) == 0) continue;
443  if (excIds.find(pyId) != excIds.end()) continue;
444 
445  // Stop Pythia from decaying the particle and include in decay set.
446  if (final) {
447  incIds.insert(pyId);
448  pythiaPtr->particleData.mayDecay(pyId, false);
449  }
450 
451  // Check for signal.
452  string egName = EvtPDL::name(egId);
453  if (egName.size() <= signalSuffix.size() || egName.substr
454  (egName.size() - signalSuffix.size()) != signalSuffix) continue;
455  signal = signals.find(pyId);
456  if (signal == signals.end()) {
457  signals[pyId].status = -1;
458  signal = signals.find(pyId);
459  }
460  signal->second.egId = egId;
461  signal->second.bfs = vector<double>(2, 0);
462  if (!final) continue;
463 
464  // Get the signal and background decay modes.
465  vector<EvtParticleDecayList> egList = egTable->getDecayTable();
466  int sigIdx = egId.getAlias();
467  int bkgIdx = EvtPDL::evtIdFromStdHep(pyId).getAlias();
468  if (sigIdx > (int)egList.size() || bkgIdx > (int)egList.size()) continue;
469  EvtParticleDecayList sigModes = egList[sigIdx];
470  EvtParticleDecayList bkgModes = egList[bkgIdx];
471  EvtParticleDecayList allModes = egList[bkgIdx];
472  double sum(0);
473 
474  // Sum signal branching fractions.
475  for (int iMode = 0; iMode < sigModes.getNMode(); ++iMode) {
476  EvtDecayBase *mode = sigModes.getDecayModel(iMode);
477  if (!mode) continue;
478  signal->second.bfs[0] += mode->getBranchingFraction();
479  sum += mode->getBranchingFraction();
480  }
481 
482  // Sum remaining background branching fractions.
483  for (int iMode = 0; iMode < allModes.getNMode(); ++iMode) {
484  EvtDecayBase *mode = allModes.getDecayModel(iMode);
485  if (!mode) continue;
486  bool match(false);
487  for (int jMode = 0; jMode < sigModes.getNMode(); ++jMode)
488  if (mode->matchingDecay(*(sigModes.getDecayModel(jMode)))) {
489  match = true; break;}
490  if (match) bkgModes.removeMode(mode);
491  else {
492  signal->second.map.push_back(iMode);
493  signal->second.bfs[1] += mode->getBranchingFraction();
494  sum += mode->getBranchingFraction();
495  }
496  }
497  signal->second.modes = bkgModes;
498  for (int iBf = 0; iBf < (int)signal->second.bfs.size(); ++iBf)
499  signal->second.bfs[iBf] /= sum;
500  }
501  if (final) updated = true;
502 
503 }
504 
505 //--------------------------------------------------------------------------
506 
507 // Update the Pythia event record with an EvtGen decay tree.
508 
509 // The production vertex of each particle (which can also be obtained
510 // in EvtGen via EvtParticle::get4Pos()) is set by the decay vertex of
511 // its mother, which in turn is calculated from the mother's
512 // lifetime. The status code 93 is used to indicate an external decay,
513 // while the status code 94 is used to indicate an oscillated external
514 // decay. If the progenitor has a single daughter with the same ID,
515 // this daughter is used as the progenitor. This is used to prevent
516 // double oscillations.
517 
518 // If the arguments after egPro are no NULL and a particle in the
519 // decay tree is a signal particle, the decay for this particle is
520 // removed and the particle is stored as a signal candidate in the
521 // pySigs and egSigs vectors, to be decayed later. However, if any of
522 // these arguments is NULL then the entire tree is written.
523 
524 void EvtGenDecays::updateEvent(Particle *pyPro, EvtParticle *egPro,
525  vector<int> *pySigs, vector<EvtParticle*> *egSigs,
526  vector<double> *bfs, double *wgt) {
527 
528  // Set up the mother vector.
529  if (!pythiaPtr) return;
530  EvtParticle* egMom = egPro;
531  if (egPro->getNDaug() == 1 && egPro->getPDGId() ==
532  egPro->getDaug(0)->getPDGId()) egMom = egPro->getDaug(0);
533  Event &event = pythiaPtr->event;
534  vector< pair<EvtParticle*, int> >
535  moms(1, pair<EvtParticle*, int>(egMom, pyPro->index()));
536 
537  // Loop over the mothers.
538  while (moms.size() != 0) {
539 
540  // Check if particle can decay.
541  egMom = moms.back().first;
542  int iMom = moms.back().second;
543  Particle* pyMom = &event[iMom];
544  moms.pop_back();
545  if (!checkVertex(pyMom)) continue;
546  bool osc(checkOsc(egMom));
547 
548  // Set the children of the mother.
549  pyMom->daughters(event.size(), event.size() + egMom->getNDaug() - 1);
550  pyMom->statusNeg();
551  Vec4 vProd = pyMom->vDec();
552  for (int iDtr = 0 ; iDtr < (int)egMom->getNDaug(); ++iDtr) {
553  EvtParticle *egDtr = egMom->getDaug(iDtr);
554  int id = egDtr->getPDGId();
555  EvtVector4R p = egDtr->getP4Lab();
556  int idx = event.append(id, 93, iMom, 0, 0, 0, 0, 0, p.get(1),
557  p.get(2), p.get(3), p.get(0), egDtr->mass());
558  Particle *pyDtr = &event.back();
559  pyDtr->vProd(vProd);
560  pyDtr->tau(egDtr->getLifetime());
561  if (pySigs && egSigs && bfs && wgt && checkSignal(pyDtr)) {
562  pySigs->push_back(pyDtr->index());
563  egSigs->push_back(egDtr);
564  bfs->push_back(signal->second.bfs[0]);
565  (*wgt) *= 1 - bfs->back();
566  egDtr->deleteDaughters();
567  }
568  if (osc) pyDtr->status(94);
569  if (egDtr->getNDaug() > 0)
570  moms.push_back(pair<EvtParticle*, int>(egDtr, idx));
571  }
572  }
573 
574 }
575 
576 //--------------------------------------------------------------------------
577 
578 // Check if a particle can decay.
579 
580 // Modified slightly from ParticleDecays::checkVertex.
581 
582 bool EvtGenDecays::checkVertex(Particle *pyPro) {
583  if (!limitDecay) return true;
584  if (limitTau0 && pyPro->tau0() > tau0Max) return false;
585  if (limitTau && pyPro->tau() > tauMax) return false;
586  if (limitRadius && pow2(pyPro->xDec()) + pow2(pyPro->yDec())
587  + pow2(pyPro->zDec()) > pow2(rMax)) return false;
588  if (limitCylinder && (pow2(pyPro->xDec()) + pow2(pyPro->yDec())
589  > pow2(xyMax) || abs(pyPro->zDec()) > zMax) ) return false;
590  return true;
591 }
592 
593 //--------------------------------------------------------------------------
594 
595 // Check if a particle is signal.
596 
597 bool EvtGenDecays::checkSignal(Particle *pyPro) {
598  signal = signals.find(pyPro->id());
599  if (signal != signals.end() && (signal->second.status < 0 ||
600  signal->second.status == pyPro->status())) return true;
601  else return false;
602 }
603 
604 //--------------------------------------------------------------------------
605 
606 // Check if an EvtGen particle has oscillated.
607 
608 // The criteria defined here for oscillation is a single daughter but
609 // with a different ID from the mother.
610 
611 bool EvtGenDecays::checkOsc(EvtParticle *egPro) {
612  if (egPro && egPro->getNDaug() == 1 &&
613  egPro->getPDGId() != egPro->getDaug(0)->getPDGId()) return true;
614  else return false;
615 }
616 
617 //==========================================================================
618 
619 #endif // end Pythia8_EvtGen_H
void setChannel(int i)
double mass() const
EvtParticle * getDaug(int i)
int getChannel() const
EvtVector4R getP4Lab() const
size_t getNDaug() const
Definition: EvtId.hh:27
Definition: EvtGen.hh:46
int getPDGId() const
void addDaug(EvtParticle *node)
double getLifetime()
void setDiagonalSpinDensity()