283 const std::vector<int>& old2New,
285 std::vector<bool>&
done,
288 FirstLastRelative firstLast,
289 const std::string& what,
290 const std::string& ind)
293 auto findNew = [old2New](
size_t old) ->
int {
296 int newIdx = findNew(
index);
297 int hepmc =
event[
index].statusHepMC();
301 << newIdx <<
" (" << hepmc <<
") ";
303 LOG(
debug) << ind <<
" already done";
308 using IdList = std::pair<int, int>;
309 constexpr int invalid = 0xFFFFFFF;
310 IdList newRelatives = std::make_pair(invalid, -invalid);
313 auto addId = [](IdList& l,
size_t id) {
314 l.first = std::min(
int(
id), l.first);
315 l.second = std::max(
int(
id), l.second);
319 auto& particle =
event[
index];
320 auto relatives = getter(particle);
322 LOG(
debug) << ind <<
" Check " << what <<
"s ["
323 << std::setw(3) << firstLast(particle).first <<
","
324 << std::setw(3) << firstLast(particle).second <<
"] "
327 for (
auto relativeIdx : relatives) {
328 int newRelative = findNew(relativeIdx);
329 if (newRelative >= 0) {
334 << relativeIdx <<
" -> "
335 << newRelative <<
" to be kept" << std::endl;
336 addId(newRelatives, newRelative);
341 << relativeIdx <<
" not to be kept "
342 << (
done[relativeIdx] ?
"already done" :
"to be done")
346 auto& relative =
event[relativeIdx];
347 if (not
done[relativeIdx]) {
362 auto grandRelatives = firstLast(relative);
363 int grandRelative1 = grandRelatives.first;
364 int grandRelative2 = grandRelatives.second;
365 assert(grandRelative1 != invalid);
366 assert(grandRelative2 != -invalid);
367 if (grandRelative1 > 0) {
368 addId(newRelatives, grandRelative1);
370 if (grandRelative2 > 0) {
371 addId(newRelatives, grandRelative2);
375 << relativeIdx <<
" gave new relatives "
376 << grandRelative1 <<
" -> " << grandRelative2;
379 << (newRelatives.second - newRelatives.first + 1) <<
" new "
382 if (newRelatives.first != invalid) {
385 int newRelative1 = newRelatives.first;
386 int newRelative2 = newRelatives.second;
387 setter(particle, newRelative1, newRelative2);
388 LOG(
debug) << ind <<
" " << what <<
"s: "
389 << firstLast(particle).first <<
" ("
390 << newRelative1 <<
"),"
391 << firstLast(particle).second <<
" ("
392 << newRelative2 <<
")";
395 setter(particle, 0, 0);
404 std::vector<int> old2new(
event.size(), -1);
411 for (
size_t i = 1;
i <
event.size();
i++) {
412 auto& particle =
event[
i];
419 auto findNew = [old2new](
size_t old) ->
int {
424 auto getMothers = [](
const Pythia8::Particle& particle) {
return particle.motherList(); };
425 auto setMothers = [](Pythia8::Particle& particle,
int m1,
int m2) { particle.mothers(m1, m2); };
426 auto firstLastMothers = [](
const Pythia8::Particle& particle) {
return std::make_pair(particle.mother1(), particle.mother2()); };
428 std::vector<bool> motherDone(
event.size(),
false);
429 for (
size_t i = 1;
i <
event.size(); ++
i) {
441 auto getDaughters = [](
const Pythia8::Particle& particle) {
446 if (std::abs(particle.status()) == 13) {
447 int d1 = particle.daughter1();
448 int d2 = particle.daughter2();
449 if (d1 == 0 and d2 == 0) {
450 return std::vector<int>();
453 return std::vector<int>{d1};
456 std::vector<int> ret(d2-d1+1);
457 std::iota(ret.begin(), ret.end(), d1);
460 return std::vector<int>{d2,d1};
462 return particle.daughterList(); };
463 auto setDaughters = [](Pythia8::Particle& particle,
int d1,
int d2) { particle.daughters(d1, d2); };
464 auto firstLastDaughters = [](
const Pythia8::Particle& particle) {
return std::make_pair(particle.daughter1(), particle.daughter2()); };
466 std::vector<bool> daughterDone(
event.size(),
false);
467 for (
size_t i =
event.size() - 1;
i > 0; --
i) {
479 Pythia8::Event pruned;
480 pruned.init(
"Pruned event", &
mPythia.particleData);
483 for (
size_t i = 1;
i <
event.size();
i++) {
484 int newIdx = findNew(
i);
489 auto particle =
event[
i];
490 int realIdx = pruned.append(particle);
491 assert(realIdx == newIdx);
507 using IdList = std::pair<int, int>;
509 auto addId = [](IdList& l,
size_t id) {
510 l.first = std::min(
int(
id), l.first);
511 l.second = std::max(
int(
id), l.second);
513 constexpr int invalid = 0xFFFFFFF;
515 std::vector<bool> shareDone(pruned.size(),
false);
516 for (
size_t i = 1;
i < pruned.size();
i++) {
521 auto& particle = pruned[
i];
522 auto daughters = particle.daughterList();
523 IdList allDaughters = std::make_pair(invalid, -invalid);
524 IdList allMothers = std::make_pair(invalid, -invalid);
525 addId(allMothers,
i);
526 for (
auto daughterIdx : daughters) {
528 addId(allDaughters, daughterIdx);
529 auto& daughter = pruned[daughterIdx];
530 auto otherMothers = daughter.motherList();
531 for (
auto otherMotherIdx : otherMothers) {
536 addId(allMothers, otherMotherIdx);
539 auto& otherMother = pruned[otherMotherIdx];
540 int otherDaughter1 = otherMother.daughter1();
541 int otherDaughter2 = otherMother.daughter2();
542 if (otherDaughter1 > 0) {
543 addId(allDaughters, otherDaughter1);
545 if (otherDaughter2 > 0) {
546 addId(allDaughters, otherDaughter2);
557 int minDaughter = allDaughters.first;
558 int maxDaughter = allDaughters.second;
559 int minMother = allMothers.first;
560 int maxMother = allMothers.second;
561 if (minMother != invalid) {
563 for (
size_t motherIdx = minMother; motherIdx <= maxMother;
565 shareDone[motherIdx] =
true;
566 if (minDaughter == invalid) {
567 pruned[motherIdx].daughters(0, 0);
569 pruned[motherIdx].daughters(minDaughter, maxDaughter);
573 if (minDaughter != invalid) {
576 for (
size_t daughterIdx = minDaughter; daughterIdx <= maxDaughter;
578 if (minMother == invalid) {
579 pruned[daughterIdx].mothers(0, 0);
581 pruned[daughterIdx].mothers(minMother, maxMother);
587 LOG(info) <<
"Pythia event was pruned from " <<
event.size()
588 <<
" to " << pruned.size() <<
" particles";
627 std::function<bool(
const Pythia8::Particle&)> partonSelect = [](
const Pythia8::Particle&) {
return true; };
629 if (not includeParton) {
632 partonSelect = [](
const Pythia8::Particle& particle) {
633 switch (particle.statusHepMC()) {
647 auto finalSelect = [partonSelect,
this](
const Pythia8::Particle& p) {
return partonSelect(p) &&
mUserFilterFcn(p); };
652 auto nParticles =
event.size();
653 for (Int_t iparticle = 1; iparticle < nParticles; iparticle++) {
655 auto particle =
event[iparticle];
656 auto pdg = particle.id();
662 particle.mother1() - 1,
663 particle.mother2() - 1,
664 particle.daughter1() - 1,
665 particle.daughter2() - 1,
675 particle.statusHepMC() == 1);
689 eventHeader->
putInfo<std::string>(Key::generator,
"pythia8");
690 eventHeader->
putInfo<
int>(Key::generatorVersion, PYTHIA_VERSION_INTEGER);
691 eventHeader->
putInfo<std::string>(Key::processName,
mPythia.info.name());
697 eventHeader->
putInfo<
int>(Key::acceptedEvents, info.nAccepted());
698 eventHeader->
putInfo<
int>(Key::attemptedEvents, info.nTried());
701 eventHeader->
putInfo<
int>(Key::pdfParton1Id, info.id1pdf());
702 eventHeader->
putInfo<
int>(Key::pdfParton2Id, info.id2pdf());
703 eventHeader->
putInfo<
float>(Key::pdfX1, info.x1pdf());
704 eventHeader->
putInfo<
float>(Key::pdfX2, info.x2pdf());
705 eventHeader->
putInfo<
float>(Key::pdfScale, info.QFac());
706 eventHeader->
putInfo<
float>(Key::pdfXF1, info.pdf1());
707 eventHeader->
putInfo<
float>(Key::pdfXF2, info.pdf2());
710 eventHeader->
putInfo<
float>(Key::xSection, info.sigmaGen() * 1e9);
711 eventHeader->
putInfo<
float>(Key::xSectionError, info.sigmaErr() * 1e9);
714 eventHeader->
putInfo<
float>(Key::eventScale, info.QRen());
715 eventHeader->
putInfo<
int>(Key::mpi, info.nMPI());
719 auto xsecErr = info.weightContainerPtr->getTotalXsecErr();
720 for (
auto w : info.weightContainerPtr->getTotalXsec()) {
722 eventHeader->
putInfo<
float>(Key::weight + post, info.weightValueByIndex(iw));
723 eventHeader->
putInfo<
float>(Key::xSection + post,
w * 1e9);
724 eventHeader->
putInfo<
float>(Key::xSectionError + post, xsecErr[iw] * 1e9);
728#if PYTHIA_VERSION_INTEGER < 8300
729 auto hiinfo =
mPythia.info.hiinfo;
731 auto hiinfo =
mPythia.info.hiInfo;
736 eventHeader->SetB(hiinfo->b());
737 eventHeader->
putInfo<
double>(Key::impactParameter, hiinfo->b());
739#if PYTHIA_VERSION_INTEGER >= 8310
740 eventHeader->
putInfo<
double>(Key::planeAngle, hiinfo->phi());
742 auto bImp = hiinfo->b();
745 int nPartProtonProj, nPartNeutronProj, nPartProtonTarg, nPartNeutronTarg;
746 int nRemnProtonProj, nRemnNeutronProj, nRemnProtonTarg, nRemnNeutronTarg;
747 int nFreeNeutronProj, nFreeProtonProj, nFreeNeutronTarg, nFreeProtonTarg;
750 getNpart(nPartProtonProj, nPartNeutronProj, nPartProtonTarg, nPartNeutronTarg);
751 getNremn(nRemnProtonProj, nRemnNeutronProj, nRemnProtonTarg, nRemnNeutronTarg);
752 getNfreeSpec(nFreeNeutronProj, nFreeProtonProj, nFreeNeutronTarg, nFreeProtonTarg);
753 eventHeader->
putInfo<
int>(Key::nColl, nColl);
755 eventHeader->
putInfo<
int>(
"Npart", nPart);
756 eventHeader->
putInfo<
int>(
"Npart_proj_p", nPartProtonProj);
757 eventHeader->
putInfo<
int>(
"Npart_proj_n", nPartNeutronProj);
758 eventHeader->
putInfo<
int>(
"Npart_targ_p", nPartProtonTarg);
759 eventHeader->
putInfo<
int>(
"Npart_targ_n", nPartNeutronTarg);
760 eventHeader->
putInfo<
int>(
"Nremn_proj_p", nRemnProtonProj);
761 eventHeader->
putInfo<
int>(
"Nremn_proj_n", nRemnNeutronProj);
762 eventHeader->
putInfo<
int>(
"Nremn_targ_p", nRemnProtonTarg);
763 eventHeader->
putInfo<
int>(
"Nremn_targ_n", nRemnNeutronTarg);
764 eventHeader->
putInfo<
int>(
"Nfree_proj_n", nFreeNeutronProj);
765 eventHeader->
putInfo<
int>(
"Nfree_proj_p", nFreeProtonProj);
766 eventHeader->
putInfo<
int>(
"Nfree_targ_n", nFreeNeutronTarg);
767 eventHeader->
putInfo<
int>(
"Nfree_targ_p", nFreeProtonTarg);
775 eventHeader->
putInfo<
int>(Key::nPartProjectile,
776 hiinfo->nAbsProj() + hiinfo->nDiffProj());
777 eventHeader->
putInfo<
int>(Key::nPartTarget,
778 hiinfo->nAbsTarg() + hiinfo->nDiffTarg());
779 eventHeader->
putInfo<
int>(Key::nCollHard, hiinfo->nCollNDTot());
877#if PYTHIA_VERSION_INTEGER < 8300
878 auto hiinfo = info.hiinfo;
880 auto hiinfo = info.hiInfo;
888 nPart = (hiinfo->nAbsProj() + hiinfo->nDiffProj() +
889 hiinfo->nAbsTarg() + hiinfo->nDiffTarg());
890 if (not hiinfo->subCollisionsPtr()) {
894 int nProtonProj, nNeutronProj, nProtonTarg, nNeutronTarg;
895 getNpart(info, nProtonProj, nNeutronProj, nProtonTarg, nNeutronTarg);
896 nPart = nProtonProj + nNeutronProj + nProtonTarg + nNeutronTarg;
906#if PYTHIA_VERSION_INTEGER < 8300
907 auto hiinfo = info.hiinfo;
909 auto hiinfo = info.hiInfo;
912 nProtonProj = nNeutronProj = nProtonTarg = nNeutronTarg = 0;
917 nProtonProj = hiinfo->nAbsProj() + hiinfo->nDiffProj();
918 nProtonTarg = hiinfo->nAbsTarg() + hiinfo->nDiffTarg();
919 if (not hiinfo->subCollisionsPtr()) {
920#if PYTHIA_VERSION_INTEGER < 8310
921 LOG(fatal) <<
"No sub-collision pointer from Pythia";
927 std::vector<Pythia8::Nucleon*> projW;
928 std::vector<Pythia8::Nucleon*> targW;
931 auto scptr = hiinfo->subCollisionsPtr();
932 for (
auto sc : *scptr) {
935 auto pW = sc.proj->status() == Pythia8::Nucleon::ABS || sc.proj->status() == Pythia8::Nucleon::DIFF;
936 auto tW = sc.targ->status() == Pythia8::Nucleon::ABS || sc.targ->status() == Pythia8::Nucleon::DIFF;
939 if (pW && std::find(projW.begin(), projW.end(), sc.proj) == projW.end()) {
940 projW.push_back(sc.proj);
941 if (sc.proj->id() == 2212) {
943 }
else if (sc.proj->id() == 2112) {
949 if (tW && std::find(targW.begin(), targW.end(), sc.targ) == targW.end()) {
950 targW.push_back(sc.targ);
951 if (sc.targ->id() == 2212) {
953 }
else if (sc.targ->id() == 2112) {