CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
PythiaDecays.cc
Go to the documentation of this file.
5 
6 #include <Pythia8/Pythia.h>
7 
8 #include <memory>
9 
10 #include "Pythia8Plugins/HepMC2.h"
11 
13  // inspired by method Pythia8Hadronizer::residualDecay() in GeneratorInterface/Pythia8Interface/src/Py8GunBase.cc
14  decayer = std::make_unique<Pythia8::Pythia>();
15  p8RndmEngine = std::make_unique<gen::P8RndmEngine>();
16  decayer->setRndmEnginePtr(p8RndmEngine.get());
17  decayer->settings.flag("ProcessLevel:all", false);
18  decayer->settings.flag("PartonLevel:FSRinResonances", false);
19  decayer->settings.flag("ProcessLevel:resonanceDecays", false);
20  decayer->init();
21 
22  // forbid all decays
23  // (decays are allowed selectively in the particleDaughters function)
24  Pythia8::ParticleData& pdt = decayer->particleData;
25  int pid = 1;
26  while (pdt.nextId(pid) > pid) {
27  pid = pdt.nextId(pid);
28  pdt.mayDecay(pid, false);
29  }
30 }
31 
33 
35  CLHEP::HepRandomEngine* engine) {
37 
38  theList.clear();
39 
40  // inspired by method Pythia8Hadronizer::residualDecay() in GeneratorInterface/Pythia8Interface/src/Py8GunBase.cc
41  int pid = particle.particle().pid();
42  decayer->event.reset();
43  Pythia8::Particle py8part(pid,
44  93,
45  0,
46  0,
47  0,
48  0,
49  0,
50  0,
51  particle.particle().momentum().x(), // note: momentum().x() and Px() are the same
52  particle.particle().momentum().y(),
53  particle.particle().momentum().z(),
54  particle.particle().momentum().t(),
55  particle.particle().mass());
56  py8part.vProd(particle.particle().X(), particle.particle().Y(), particle.particle().Z(), particle.particle().T());
57  decayer->event.append(py8part);
58 
59  int nentries_before = decayer->event.size();
60  decayer->particleData.mayDecay(pid,
61  true); // switch on the decay of this and only this particle (avoid double decays)
62  decayer->next(); // do the decay
63  decayer->particleData.mayDecay(pid, false); // switch it off again
64  int nentries_after = decayer->event.size();
65  if (nentries_after <= nentries_before)
66  return theList;
67 
68  theList.reserve(nentries_after - nentries_before);
69 
70  for (int ipart = nentries_before; ipart < nentries_after; ipart++) {
71  Pythia8::Particle& py8daughter = decayer->event[ipart];
72  theList
73  .emplace_back(makeParticle(
74  particle.particleDataTable(),
75  py8daughter.id(),
76  XYZTLorentzVector(py8daughter.px(), py8daughter.py(), py8daughter.pz(), py8daughter.e()),
77  XYZTLorentzVector(py8daughter.xProd(), py8daughter.yProd(), py8daughter.zProd(), py8daughter.tProd())))
78  .setMass(py8daughter.m());
79  }
80 
81  return theList;
82 }
int pid() const
get the HEP particle ID number
Definition: RawParticle.h:277
double Z() const
z of vertex
Definition: RawParticle.h:288
const DaughterParticleList & particleDaughters(ParticlePropagator &particle, CLHEP::HepRandomEngine *)
Definition: PythiaDecays.cc:34
const XYZTLorentzVector & momentum() const
the momentum fourvector
Definition: RawParticle.h:321
const HepPDT::ParticleDataTable * particleDataTable() const
DaughterParticleList theList
Definition: PythiaDecays.h:35
RawParticle makeParticle(HepPDT::ParticleDataTable const *, int id, const math::XYZTLorentzVector &p)
Definition: makeParticle.cc:28
RawParticle const & particle() const
The particle being propagated.
HepPDT::ParticleData ParticleData
double T() const
vertex time
Definition: RawParticle.h:289
std::unique_ptr< Pythia8::Pythia > decayer
Definition: PythiaDecays.h:36
double mass() const
get the MEASURED mass
Definition: RawParticle.h:295
std::unique_ptr< gen::P8RndmEngine > p8RndmEngine
Definition: PythiaDecays.h:37
double Y() const
y of vertex
Definition: RawParticle.h:287
double X() const
x of vertex
Definition: RawParticle.h:286
math::XYZTLorentzVector XYZTLorentzVector
Definition: RawParticle.h:25
std::vector< RawParticle > DaughterParticleList
Definition: PythiaDecays.h:25