ParticleDecay/src/PythiaDecays.cc

0001 #include "FastSimulation/ParticlePropagator/interface/ParticlePropagator.h"
0002 #include "FastSimulation/Particle/interface/makeParticle.h"
0003 #include "FastSimulation/ParticleDecay/interface/PythiaDecays.h"
0004 #include "FWCore/ServiceRegistry/interface/RandomEngineSentry.h"
0005
0006 #include <Pythia8/Pythia.h>
0007
0008 #include <memory>
0009
0010 #include "Pythia8Plugins/HepMC2.h"
0011
0012 PythiaDecays::PythiaDecays() {
0013   // inspired by method Pythia8Hadronizer::residualDecay() in GeneratorInterface/Pythia8Interface/src/Py8GunBase.cc
0014   decayer = std::make_unique<Pythia8::Pythia>();
0015   p8RndmEngine = std::make_shared<gen::P8RndmEngine>();
0016   decayer->setRndmEnginePtr(p8RndmEngine);
0017   decayer->settings.flag("ProcessLevel:all", false);
0018   decayer->settings.flag("PartonLevel:FSRinResonances", false);
0019   decayer->settings.flag("ProcessLevel:resonanceDecays", false);
0020   decayer->init();
0021
0022   // forbid all decays
0023   // (decays are allowed selectively in the particleDaughters function)
0024   Pythia8::ParticleData& pdt = decayer->particleData;
0025   int pid = 1;
0026   while (pdt.nextId(pid) > pid) {
0027     pid = pdt.nextId(pid);
0028     pdt.mayDecay(pid, false);
0029   }
0030 }
0031
0032 PythiaDecays::~PythiaDecays() {}
0033
0034 const DaughterParticleList& PythiaDecays::particleDaughters(ParticlePropagator& particle,
0035                                                             CLHEP::HepRandomEngine* engine) {
0036   edm::RandomEngineSentry<gen::P8RndmEngine> sentry(p8RndmEngine.get(), engine);
0037
0038   theList.clear();
0039
0040   // inspired by method Pythia8Hadronizer::residualDecay() in GeneratorInterface/Pythia8Interface/src/Py8GunBase.cc
0041   int pid = particle.particle().pid();
0042   decayer->event.reset();
0043   Pythia8::Particle py8part(pid,
0044                             93,
0045                             0,
0046                             0,
0047                             0,
0048                             0,
0049                             0,
0050                             0,
0051                             particle.particle().momentum().x(),  // note: momentum().x() and Px() are the same
0052                             particle.particle().momentum().y(),
0053                             particle.particle().momentum().z(),
0054                             particle.particle().momentum().t(),
0055                             particle.particle().mass());
0056   py8part.vProd(particle.particle().X(), particle.particle().Y(), particle.particle().Z(), particle.particle().T());
0057   decayer->event.append(py8part);
0058
0059   int nentries_before = decayer->event.size();
0060   decayer->particleData.mayDecay(pid,
0061                                  true);  // switch on the decay of this and only this particle (avoid double decays)
0062   decayer->next();                       // do the decay
0063   decayer->particleData.mayDecay(pid, false);  // switch it off again
0064   int nentries_after = decayer->event.size();
0065   if (nentries_after <= nentries_before)
0066     return theList;
0067
0068   theList.reserve(nentries_after - nentries_before);
0069
0070   for (int ipart = nentries_before; ipart < nentries_after; ipart++) {
0071     Pythia8::Particle& py8daughter = decayer->event[ipart];
0072     theList
0073         .emplace_back(makeParticle(
0074             particle.particleDataTable(),
0075             py8daughter.id(),
0076             XYZTLorentzVector(py8daughter.px(), py8daughter.py(), py8daughter.pz(), py8daughter.e()),
0077             XYZTLorentzVector(py8daughter.xProd(), py8daughter.yProd(), py8daughter.zProd(), py8daughter.tProd())))
0078         .setMass(py8daughter.m());
0079   }
0080
0081   return theList;
0082 }