20 #include <generators/utilities/GeneratorConst.h>
21 #include "Pythia8/Pythia.h"
22 #include "EvtGen/EvtGen.hh"
23 #include "EvtGenBase/EvtRandomEngine.hh"
24 #include "EvtGenBase/EvtParticle.hh"
25 #include "EvtGenBase/EvtParticleFactory.hh"
26 #include "EvtGenBase/EvtPDL.hh"
27 #include "EvtGenBase/EvtDecayTable.hh"
28 #include "EvtGenBase/EvtParticleDecayList.hh"
29 #include "EvtGenBase/EvtDecayBase.hh"
30 #include "EvtGenExternal/EvtExternalGenList.hh"
99 EvtGenDecays(Pythia8::Pythia* pythiaPtrIn, EvtGen* evtGen,
bool limit =
true);
102 EvtGenDecays(Pythia8::Pythia* pythiaPtrIn, std::string decayFile, std::string particleDataFile,
103 EvtExternalGenList* extPtrIn = 0, EvtAbsRadCorr* fsrPtrIn = 0,
104 int mixing = 1,
bool xml =
false,
bool limit =
true,
105 bool extUse =
true,
bool fsrUse =
true);
168 void updateEvent(Pythia8::Particle* pyPro, EvtParticle* egPro,
169 std::vector<int>* pySigs = 0, std::vector<EvtParticle*>* egSigs = 0,
170 std::vector<double>* bfs = 0,
double* wgt = 0);
262 extOwner(false), fsrOwner(false), extPtr(nullptr), fsrPtr(nullptr),
263 signalSuffix(
"_SIGNAL"), pythiaPtr(pythiaPtrIn), rndm(nullptr),
264 evtgen(evtGen), updated(false)
270 std::string particleDataFile, EvtExternalGenList* extPtrIn,
271 EvtAbsRadCorr* fsrPtrIn,
int mixing,
bool xml,
bool limit,
272 bool extUse,
bool fsrUse) : extPtr(extPtrIn), fsrPtr(fsrPtrIn),
273 signalSuffix(
"_SIGNAL"), pythiaPtr(pythiaPtrIn), rndm(&pythiaPtr->rndm),
280 "EvtGenDecays: extPtr is null but fsrPtr is provided");
288 evtgen =
new EvtGen(decayFile.c_str(), particleDataFile.c_str(),
342 Pythia8::Event&
event =
pythiaPtr->event;
343 std::vector<int> pySigs; std::vector<EvtParticle*> egSigs, egPrts;
344 std::vector<double> bfs;
double wgt(1.);
345 for (
int iPro = 0; iPro <
event.size(); ++iPro) {
348 Pythia8::Particle* pyPro = &
event[iPro];
349 if (!pyPro->isFinal())
continue;
351 if (pyPro->status() == 93 || pyPro->status() == 94)
continue;
354 EvtParticle* egPro = EvtParticleFactory::particleFactory
355 (EvtPDL::evtIdFromStdHep(pyPro->id()),
356 EvtVector4R(pyPro->e(), pyPro->px(), pyPro->py(), pyPro->pz()));
357 egPrts.push_back(egPro);
358 egPro->setDiagonalSpinDensity();
359 evtgen->generateDecay(egPro);
360 pyPro->tau(egPro->getLifetime());
365 updateEvent(pyPro, egPro, &pySigs, &egSigs, &bfs, &wgt);
369 pySigs.push_back(pyPro->index());
370 egSigs.push_back(egPro);
371 bfs.push_back(
signal->second.bfs[0]);
372 wgt *= 1 - bfs.back();
373 egPro->deleteDaughters();
374 EvtParticle* egDau = EvtParticleFactory::particleFactory
375 (EvtPDL::evtIdFromStdHep(pyPro->id()),
376 EvtVector4R(pyPro->e(), pyPro->px(), pyPro->py(), pyPro->pz()));
377 egDau->addDaug(egPro);
378 egDau->setDiagonalSpinDensity();
381 }
else updateEvent(pyPro, egPro, &pySigs, &egSigs, &bfs, &wgt);
383 if (pySigs.size() == 0) {
384 for (
int iPrt = 0; iPrt < (int)egPrts.size(); ++iPrt)
385 egPrts[iPrt]->deleteTree();
390 std::vector<int> modes;
int force(-1);
391 for (
int iTry = 1; iTry <=
NTRYDECAY; ++iTry) {
394 modes.clear(); force =
pythiaPtr->rndm.pick(bfs); n = 0;
395 for (
int iSig = 0; iSig < (int)pySigs.size(); ++iSig) {
396 if (iSig == force) modes.push_back(0);
397 else modes.push_back(
pythiaPtr->rndm.flat() > bfs[iSig]);
398 if (modes.back() == 0) ++n;
400 if (
pythiaPtr->rndm.flat() <= 1.0 / n)
break;
403 pythiaPtr->info.errorMsg(
"Warning in EvtGenDecays::decay: maximum "
404 "number of signal decay attempts exceeded");
409 for (
int iSig = 0; iSig < (int)pySigs.size(); ++iSig) {
410 EvtParticle* egSig = egSigs[iSig];
411 Pythia8::Particle* pySig = &
event[pySigs[iSig]];
413 if (egSig->getNDaug() > 0) egSig = egSig->getDaug(0);
414 if (modes[iSig] == 0) egSig->setId(
signal->second.egId);
416 signal->second.modes.getDecayModel(egSig);
417 egSig->setChannel(
signal->second.map[egSig->getChannel()]);
420 pySig->status(pySig->status() == 92 || pySig->status() == 94 ? 96 : 95);
421 evtgen->generateDecay(egSig);
426 for (
int iPrt = 0; iPrt < (int)egPrts.size(); ++iPrt)
427 egPrts[iPrt]->deleteTree();
443 for (
int entry = 0; entry < (int)EvtPDL::entries(); ++entry) {
444 EvtId egId = EvtPDL::getEntry(entry);
445 int pyId = EvtPDL::getStdHep(egId);
446 pythiaPtr->particleData.spinType(pyId, EvtPDL::getSpinType(egId));
447 pythiaPtr->particleData.chargeType(pyId, EvtPDL::chg3(egId));
448 pythiaPtr->particleData.m0(pyId, EvtPDL::getMass(egId));
449 pythiaPtr->particleData.mWidth(pyId, EvtPDL::getWidth(egId));
450 pythiaPtr->particleData.mMin(pyId, EvtPDL::getMinMass(egId));
451 pythiaPtr->particleData.mMax(pyId, EvtPDL::getMaxMass(egId));
452 pythiaPtr->particleData.tau0(pyId, EvtPDL::getctau(egId));
467 int pyId =
pythiaPtr->particleData.nextId(0);
469 EvtId egId = EvtPDL::evtIdFromStdHep(pyId);
470 EvtPDL::reSetMass(egId,
pythiaPtr->particleData.m0(pyId));
471 EvtPDL::reSetWidth(egId,
pythiaPtr->particleData.mWidth(pyId));
472 EvtPDL::reSetMassMin(egId,
pythiaPtr->particleData.mMin(pyId));
473 EvtPDL::reSetMassMax(egId,
pythiaPtr->particleData.mMax(pyId));
474 pyId =
pythiaPtr->particleData.nextId(pyId);
495 EvtDecayTable* egTable = EvtDecayTable::getInstance();
496 if (!egTable)
return;
497 for (
int iEntry = 0; iEntry < (int)EvtPDL::entries(); ++iEntry) {
498 EvtId egId = EvtPDL::getEntry(iEntry);
499 int pyId = EvtPDL::getStdHep(egId);
500 if (egTable->getNModes(egId) == 0)
continue;
506 pythiaPtr->particleData.mayDecay(pyId,
false);
510 std::string egName = EvtPDL::name(egId);
511 if (egName.size() <=
signalSuffix.size() || egName.substr
518 signal->second.egId = egId;
519 signal->second.bfs = std::vector<double>(2, 0);
520 if (!
final)
continue;
523 std::vector<EvtParticleDecayList> egList = egTable->getDecayTable();
524 int sigIdx = egId.getAlias();
525 int bkgIdx = EvtPDL::evtIdFromStdHep(pyId).getAlias();
526 if (sigIdx > (
int)egList.size() || bkgIdx > (
int)egList.size())
continue;
527 EvtParticleDecayList sigModes = egList[sigIdx];
528 EvtParticleDecayList bkgModes = egList[bkgIdx];
529 EvtParticleDecayList allModes = egList[bkgIdx];
533 for (
int iMode = 0; iMode < sigModes.getNMode(); ++iMode) {
534 EvtDecayBase* mode = sigModes.getDecayModel(iMode);
536 signal->second.bfs[0] += mode->getBranchingFraction();
537 sum += mode->getBranchingFraction();
541 for (
int iMode = 0; iMode < allModes.getNMode(); ++iMode) {
542 EvtDecayBase* mode = allModes.getDecayModel(iMode);
545 for (
int jMode = 0; jMode < sigModes.getNMode(); ++jMode)
546 if (mode->matchingDecay(*(sigModes.getDecayModel(jMode)))) {
549 if (match) bkgModes.removeMode(mode);
551 signal->second.map.push_back(iMode);
552 signal->second.bfs[1] += mode->getBranchingFraction();
553 sum += mode->getBranchingFraction();
556 signal->second.modes = bkgModes;
557 for (
int iBf = 0; iBf < (int)
signal->second.bfs.size(); ++iBf)
558 signal->second.bfs[iBf] /= sum;
584 std::vector<int>* pySigs, std::vector<EvtParticle*>* egSigs,
585 std::vector<double>* bfs,
double* wgt)
590 EvtParticle* egMom = egPro;
591 if (egPro->getNDaug() == 1 && egPro->getPDGId() ==
592 egPro->getDaug(0)->getPDGId()) egMom = egPro->getDaug(0);
593 Pythia8::Event&
event =
pythiaPtr->event;
594 std::vector< std::pair<EvtParticle*, int> >
595 moms(1, std::pair<EvtParticle*, int>(egMom, pyPro->index()));
598 while (moms.size() != 0) {
601 egMom = moms.back().first;
602 int iMom = moms.back().second;
603 Pythia8::Particle* pyMom = &
event[iMom];
609 pyMom->daughters(event.size(), event.size() + egMom->getNDaug() - 1);
611 Pythia8::Vec4 vProd = pyMom->vDec();
612 for (
int iDtr = 0 ; iDtr < (int)egMom->getNDaug(); ++iDtr) {
613 EvtParticle* egDtr = egMom->getDaug(iDtr);
614 int id = egDtr->getPDGId();
615 EvtVector4R p = egDtr->getP4Lab();
616 int idx =
event.append(
id, 93, iMom, 0, 0, 0, 0, 0, p.get(1),
617 p.get(2), p.get(3), p.get(0), egDtr->mass());
618 Pythia8::Particle* pyDtr = &
event.back();
620 pyDtr->tau(egDtr->getLifetime());
621 if (pySigs && egSigs && bfs && wgt &&
checkSignal(pyDtr)) {
622 pySigs->push_back(pyDtr->index());
623 egSigs->push_back(egDtr);
624 bfs->push_back(
signal->second.bfs[0]);
625 (*wgt) *= 1 - bfs->back();
626 egDtr->deleteDaughters();
628 if (osc) pyDtr->status(94);
629 if (egDtr->getAttribute(
"FSR"))
630 pyDtr->status(Belle2::GeneratorConst::FSR_STATUS_CODE);
631 if (egDtr->getNDaug() > 0)
632 moms.push_back(std::pair<EvtParticle*, int>(egDtr, idx));
646 using ::Pythia8::pow2;
651 if (
limitRadius && pow2(pyPro->xDec()) + pow2(pyPro->yDec())
652 + pow2(pyPro->zDec()) > pow2(
rMax))
return false;
653 if (
limitCylinder && (pow2(pyPro->xDec()) + pow2(pyPro->yDec())
654 > pow2(
xyMax) || abs(pyPro->zDec()) >
zMax))
return false;
666 signal->second.status == pyPro->status()))
return true;
679 if (egPro && egPro->getNDaug() == 1 &&
680 egPro->getPDGId() != egPro->getDaug(0)->getPDGId())
return true;