122 bool doDebug =
false;
123 bool doExtra =
false;
131 return StatusCode::FAILURE;
134 if (mcEvts->
empty()) {
136 return StatusCode::SUCCESS;
143 delete thinnedMcEvts;
144 thinnedMcEvts =
nullptr;
145 return StatusCode::FAILURE;
150 const HepMC::GenEvent* mcEvt = mcEvts->
front();
151 auto wtCont = mcEvt->weights();
152 if (!wtCont.empty()) {
156 int inEvent = mcEvt->event_number();
163 HepMC::GenEvent* thinEvt =
new HepMC::GenEvent(*mcEvt);
164 int nEvent = thinEvt->event_number();
168 std::cout <<
"========== BEGIN EVENT BEFORE THINNING ==========" << std::endl;
170 std::cout <<
"========== END EVENT BEFORE THINNING ==========" << std::endl;
190 typedef std::pair<HepMC::GenVertexPtr, HepMC::GenParticlePtr> vpPair;
191 std::vector<vpPair> removePV;
192 std::vector<vpPair> addinPV;
193 std::vector<vpPair> addoutPV;
194 std::vector<HepMC::GenVertexPtr> removeV;
197 std::list<HepMC::GenParticlePtr> deleteP;
198 std::list<HepMC::GenVertexPtr> deleteV;
210 std::vector<HepMC::GenVertexPtr> hadVertices;
213 for (
auto& hadv: thinEvt->vertices()) {
215 if (hadv->particles_in().size() < 2)
continue;
216 if (hadv->particles_out().size() < 1)
continue;
220 bool isHadVtx =
true;
221 bool isHadOut =
false;
222 for (
const auto& inp: hadv->particles_in() ) {
225 for (
const auto& vp: hadv->particles_out()) {
229 isHadVtx = isHadVtx && isHadOut;
230 if (isHadVtx) hadVertices.push_back(hadv);
231 if (doDebug && isHadVtx)
ATH_MSG_VERBOSE(
"Hadronization vertex " << hadv);
234 HepMC::GenEvent::vertex_iterator hadv = thinEvt->vertices_begin();
235 HepMC::GenEvent::vertex_iterator hadvB = thinEvt->vertices_begin();
236 HepMC::GenEvent::vertex_iterator hadvE = thinEvt->vertices_end();
238 for (; hadv != hadvE; ++hadv) {
239 if (!(*hadv))
continue;
240 if ((*hadv)->particles_in_size() < 2)
continue;
241 if ((*hadv)->particles_out_size() < 1)
continue;
246 bool isHadVtx =
true;
247 bool isHadOut =
false;
248 HepMC::GenVertex::particles_in_const_iterator inp = (*hadv)->particles_in_const_begin();
249 HepMC::GenVertex::particles_in_const_iterator inpE = (*hadv)->particles_in_const_end();
250 for (; inp != inpE; ++inp) {
253 HepMC::GenVertex::particles_out_const_iterator vp = (*hadv)->particles_out_const_begin();
254 HepMC::GenVertex::particles_out_const_iterator vpE = (*hadv)->particles_out_const_end();
255 for (; vp != vpE; ++vp) {
259 isHadVtx = isHadVtx && isHadOut;
260 if (isHadVtx) hadVertices.push_back(*hadv);
261 if (doDebug && isHadVtx)
ATH_MSG_VERBOSE(
"Hadronization vertex " << *hadv);
265 if (hadVertices.empty()) {
269 return StatusCode::SUCCESS;
278 for (
unsigned int iv = 0; iv < hadVertices.size(); ++iv) {
280 if (doDebug)
ATH_MSG_DEBUG(
"Removing partons from hadVertex " << ivtx);
281 for (
auto pin: ivtx->particles_in()) {
282 removePV.push_back(vpPair(ivtx, pin));
287 for (
unsigned int iv = 0; iv < hadVertices.size(); ++iv) {
289 auto descendants = HepMC::descendant_particles(ivtx);
290 for (
auto pout: descendants) {
292 if (vpar) removePV.push_back(vpPair(vpar,
pout));
294 if (vend) removePV.push_back(vpPair(vend,
pout));
295 deleteP.push_back(
pout);
299 for (
unsigned int iv = 0; iv < hadVertices.size(); ++iv) {
300 HepMC::GenVertex* ivtx = hadVertices[iv];
301 if (doDebug)
ATH_MSG_DEBUG(
"Removing partons from hadVertex " << ivtx);
302 HepMC::GenVertex::particles_in_const_iterator pin = ivtx->particles_in_const_begin();
303 HepMC::GenVertex::particles_in_const_iterator pinE = ivtx->particles_in_const_end();
304 for (; pin != pinE; ++pin) {
305 removePV.emplace_back(ivtx, *pin);
311 for (
unsigned int iv = 0; iv < hadVertices.size(); ++iv) {
312 HepMC::GenVertex* ivtx = hadVertices[iv];
313 HepMC::GenVertex::particle_iterator
pout = ivtx->particles_begin(HepMC::descendants);
314 HepMC::GenVertex::particle_iterator poutE = ivtx->particles_end(HepMC::descendants);
316 HepMC::GenVertex* vpar = (*pout)->production_vertex();
317 if (vpar) removePV.emplace_back(vpar, *
pout);
318 HepMC::GenVertex* vend = (*pout)->end_vertex();
319 if (vend) removePV.emplace_back(vend, *
pout);
320 deleteP.push_back(*
pout);
331 for (
auto hadv:thinEvt->vertices()) {
333 if (hadv->particles_in().size() == 0 && hadv->particles_out().size() == 0) {
334 removeV.push_back(hadv);
335 deleteV.push_back(hadv);
339 for (
auto pin: hadv->particles_in()) {
340 removePV.push_back(vpPair(hadv, pin));
343 for (
auto pout: hadv->particles_out()) {
344 removePV.push_back(vpPair(hadv,
pout));
347 removeV.push_back(hadv);
348 deleteV.push_back(hadv);
351 for (hadv = hadvB; hadv != hadvE; ++hadv) {
354 if ((*hadv)->particles_in_size() == 0 && (*hadv)->particles_out_size() == 0) {
355 removeV.push_back(*hadv);
356 deleteV.push_back(*hadv);
361 HepMC::GenVertex::particles_in_const_iterator pin = (*hadv)->particles_in_const_begin();
362 HepMC::GenVertex::particles_in_const_iterator pinE = (*hadv)->particles_in_const_end();
363 for (; pin != pinE; ++pin) {
364 removePV.emplace_back(*hadv, *pin);
367 HepMC::GenVertex::particles_out_const_iterator
pout = (*hadv)->particles_out_const_begin();
368 HepMC::GenVertex::particles_out_const_iterator poutE = (*hadv)->particles_out_const_end();
370 removePV.emplace_back(*hadv, *
pout);
373 removeV.push_back(*hadv);
374 deleteV.push_back(*hadv);
380 for (
unsigned int i = 0;
i < removePV.size(); ++
i) {
384 v->remove_particle_in(
p);
385 v->remove_particle_out(
p);
387 v->remove_particle(
p);
391 for (
unsigned int i = 0;
i < addoutPV.size(); ++
i) {
394 v->add_particle_out(
p);
397 for (
unsigned int iv = 1; iv < hadVertices.size(); ++iv) {
399 if (
v->particles_in().size() != 0 ||
v->particles_out().size() != 0) {
401 <<
" with in/out particles " <<
v->particles_in().size() <<
" " <<
v->particles_out().size());
403 thinEvt->remove_vertex(hadVertices[iv]);
408 for (
unsigned int iv = 1; iv < hadVertices.size(); ++iv) {
409 HepMC::GenVertex*
v = hadVertices[iv];
410 if (
v->particles_in_size() != 0 ||
v->particles_out_size() != 0) {
411 ATH_MSG_WARNING(
"Removing vertex " <<
HepMC::uniqueID(
v) <<
" for event " << nEvent <<
" with in/out particles " <<
v->particles_in_size() <<
" " <<
v->particles_out_size());
418 if (doDebug)
ATH_MSG_DEBUG(
"Deleting hadronization vertices " << deleteV.size());
421 for (dvItr = deleteV.begin(); dvItr != deleteV.end(); ++dvItr) {
423 if (*dvItr)
delete (*dvItr);
428 for (dpItr = deleteP.begin(); dpItr != deleteP.end(); ++dpItr) {
430 if (*dpItr)
delete (*dpItr);
438 if (doDebug && doExtra) {
439 std::cout <<
"========== BEGIN EVENT BEFORE CLUSTER ==========" << std::endl;
441 std::cout <<
"========== END EVENT BEFORE CLUSTER ==========" << std::endl;
456 using vpPair = std::pair<HepMC::GenVertexPtr, HepMC::GenParticlePtr>;
464 using pkPair = std::pair<HepMC::GenParticlePtr, HepMC::FourVector>;
465 std::vector<pkPair> changePK;
469 if (doDebug)
ATH_MSG_DEBUG(
"New parton pass " << inEvent <<
" " << thinEvt->particles_size() <<
" " << thinEvt->vertices_size());
479 for (
auto fp: thinEvt->particles() ) {
495 if (pvtx->particles_in().size() == 1 && pvtx->particles_out().size() == 1) {
498 if (!pp || pp->parent_event() ==
nullptr) {
505 if (!ppvtx || ppvtx->parent_event() ==
nullptr) {
512 removePV.push_back(vpPair(ppvtx, pp));
513 removePV.push_back(vpPair(pvtx, pp));
514 deleteP.push_back(pp);
515 removeV.push_back(pvtx);
516 deleteV.push_back(pvtx);
517 addoutPV.push_back(vpPair(ppvtx,
fp));
518 if (doDebug) {
ATH_MSG_DEBUG(
"1->1: ppvtx,pp,pvtx,fp,evtx " << ppvtx <<
" " << pp <<
" " << pvtx <<
" " <<
fp); }
527 if (pvtx->particles_in().size() == 2 && pvtx->particles_out().size() == 1) {
533 if (std::abs(pp1->momentum().perp()) < 1.e-3)
continue;
534 if (std::abs(pp2->momentum().perp()) < 1.e-3)
continue;
538 if (!ppvtx1 || ppvtx1->parent_event() ==
nullptr) {
543 if (!ppvtx2 || ppvtx2->parent_event() ==
nullptr) {
550 removePV.push_back(vpPair(pvtx,
fp));
551 removePV.push_back(vpPair(pvtx, pp1));
552 removePV.push_back(vpPair(pvtx, pp2));
553 deleteP.push_back(
fp);
554 removeV.push_back(pvtx);
555 deleteV.push_back(pvtx);
557 ATH_MSG_DEBUG(
"2->1: ppvtx1,pp1,ppvtx2,pp2,pvtx,fp " << ppvtx1 <<
" " << pp1
558 <<
" " << ppvtx2 <<
" " << pp2 <<
" " << pvtx <<
" "<<
fp);
567 if (pvtx->particles_in().size() == 1 && pvtx->particles_out().size() == 2) {
576 }
else if (
fp == pout2) {
582 ATH_MSG_WARNING(
"1->2: No match found for branching " <<
fp <<
" " << pvtx <<
" " << pout1 <<
" " << pout2);
585 if (
fp != pout1)
continue;
589 if (pout1->momentum().m() < -1.0 || pout1->momentum().perp() < 1.0)
continue;
590 if (pout2->momentum().m() < -1.0 || pout2->momentum().perp() < 1.0)
continue;
593 double m12 = p12.m();
595 if (std::abs(m12) > 10. + 1.0
e-5 * p12.e()) {
596 ATH_MSG_WARNING(
"Spacelike mass^2 for parton sum " << m12 <<
" " << pp <<
" " << pvtx <<
" " << pout1 <<
" " << pout2);
600 if (doDebug)
ATH_MSG_DEBUG(
"1->2: parton pair mass " << m12);
603 if (doDebug)
ATH_MSG_DEBUG(
"Keeping 1->2: parton mass " << m12);
608 if (!ppvtx || ppvtx->parent_event() ==
nullptr) {
616 if (doDebug)
ATH_MSG_DEBUG(
"Merging 1->2: mass " << p12.m());
617 changePK.push_back(pkPair(pp, p12));
618 removePV.push_back(vpPair(pvtx, pp));
619 removePV.push_back(vpPair(pvtx, pout1));
620 removePV.push_back(vpPair(pvtx, pout2));
621 deleteP.push_back(pout1);
622 deleteP.push_back(pout2);
623 removeV.push_back(pvtx);
624 deleteV.push_back(pvtx);
626 ATH_MSG_DEBUG(
"Merge 1->2: ppvtx,pp,pvtx,pout1,pout2,evtx " << ppvtx <<
" " << pp <<
" " << pvtx <<
" " << pout1 <<
" "
628 ATH_MSG_DEBUG(
"Merge 1->2: id " << pp->pdg_id() <<
" " << pout1->pdg_id() <<
" " << pout2->pdg_id());
635 if (pvtx->particles_in().size() == 1) {
638 if (std::abs(pp->pdg_id()) == 2212) iCase = -1;
643 if (doDebug)
ATH_MSG_DEBUG(
"Case not found " << pvtx <<
" " <<
fp <<
" " << pvtx->particles_in().size() <<
" " << pvtx->particles_out().size());
648 HepMC::GenEvent::particle_iterator finp = thinEvt->particles_begin();
649 HepMC::GenEvent::particle_iterator finpB = thinEvt->particles_begin();
650 HepMC::GenEvent::particle_iterator finpE = thinEvt->particles_end();
653 for (finp = finpB; finp != finpE; ++finp) {
661 HepMC::GenVertex* pvtx =
fp->production_vertex();
675 if (pvtx->particles_in_size() == 1 && pvtx->particles_out_size() == 1) {
677 HepMC::GenVertex::particles_in_const_iterator pitr = pvtx->particles_in_const_begin();
685 HepMC::GenVertex* ppvtx = pp->production_vertex();
694 removePV.emplace_back(ppvtx, pp);
695 removePV.emplace_back(pvtx, pp);
696 deleteP.push_back(pp);
697 removeV.push_back(pvtx);
698 deleteV.push_back(pvtx);
699 addoutPV.emplace_back(ppvtx,
fp);
712 if (pvtx->particles_in_size() == 2 && pvtx->particles_out_size() == 1) {
714 HepMC::GenVertex::particles_in_const_iterator pitr = pvtx->particles_in_const_begin();
721 if (std::abs(pp1->momentum().perp()) < 1.e-3)
continue;
722 if (std::abs(pp2->momentum().perp()) < 1.e-3)
continue;
724 HepMC::GenVertex* ppvtx1 = pp1->production_vertex();
725 HepMC::GenVertex* ppvtx2 = pp2->production_vertex();
740 removePV.emplace_back(pvtx,
fp);
741 removePV.emplace_back(pvtx, pp1);
742 removePV.emplace_back(pvtx, pp2);
743 deleteP.push_back(
fp);
744 removeV.push_back(pvtx);
745 deleteV.push_back(pvtx);
761 if (pvtx->particles_in_size() == 1 && pvtx->particles_out_size() == 2) {
762 HepMC::GenVertex::particles_out_const_iterator poutitr = pvtx->particles_out_const_begin();
773 }
else if (
fp == pout2) {
782 if (
fp != pout1)
continue;
784 HepMC::GenVertex::particles_in_const_iterator pitr = pvtx->particles_in_const_begin();
788 if (pout1->momentum().m() < -1.0 || pout1->momentum().perp() < 1.0)
continue;
789 if (pout2->momentum().m() < -1.0 || pout2->momentum().perp() < 1.0)
continue;
793 double m12 = p12.m();
795 if (fabs(m12) > 10. + 1.0
e-5 * p12.e()) {
800 if (doDebug)
ATH_MSG_DEBUG(
"1->2: parton pair mass " << m12);
803 if (doDebug)
ATH_MSG_DEBUG(
"Keeping 1->2: parton mass " << m12);
808 HepMC::GenVertex* ppvtx = pp->production_vertex();
818 if (doDebug)
ATH_MSG_DEBUG(
"Merging 1->2: mass " << p12.m());
820 changePK.emplace_back(pp, p12);
821 removePV.emplace_back(pvtx, pp);
822 removePV.emplace_back(pvtx, pout1);
823 removePV.emplace_back(pvtx, pout2);
825 deleteP.push_back(pout1);
826 deleteP.push_back(pout2);
827 removeV.push_back(pvtx);
828 deleteV.push_back(pvtx);
831 ATH_MSG_DEBUG(
"Merge 1->2: ppvtx,pp,pvtx,pout1,pout2,evtx " << ppvtx <<
" " << pp <<
" " << pvtx <<
" " << pout1 <<
" "
833 ATH_MSG_DEBUG(
"Merge 1->2: id " << pp->pdg_id() <<
" " << pout1->pdg_id() <<
" " << pout2->pdg_id());
842 if (pvtx->particles_in_size() == 1) {
844 HepMC::GenVertex::particles_in_const_iterator pitr = pvtx->particles_in_const_begin();
846 if (std::abs(pp->pdg_id()) == 2212) iCase = -1;
860 if (doDebug)
ATH_MSG_DEBUG(
"Actually removing particles " << removePV.size());
862 for (
unsigned int i = 0;
i < removePV.size(); ++
i) {
867 v->remove_particle_in(
p);
868 v->remove_particle_out(
p);
870 v->remove_particle(
p);
875 if (doDebug)
ATH_MSG_DEBUG(
"Actually add particles in/out " << addinPV.size() <<
" " << addoutPV.size());
876 for (
unsigned int i = 0;
i < addoutPV.size(); ++
i) {
880 v->add_particle_out(
p);
884 for (
unsigned int i = 0;
i < changePK.size(); ++
i) {
887 pp->set_momentum(changePK[
i].
second);
891 if (doDebug)
ATH_MSG_DEBUG(
"Actually remove vertices " << removeV.size());
892 for (
unsigned int i = 0;
i < removeV.size(); ++
i) {
894 int nv = thinEvt->vertices().size();
895 if (doDebug) {
ATH_MSG_VERBOSE(
"Removing vertex " << removeV[
i] <<
" " << nv <<
" " << thinEvt->vertices().size()); }
896 thinEvt->remove_vertex(removeV[
i]);
898 int nv = thinEvt->vertices_size();
899 if (thinEvt->remove_vertex(removeV[
i])) {
900 if (doDebug) {
ATH_MSG_VERBOSE(
"Removed vertex " << removeV[
i] <<
" " << nv <<
" " << thinEvt->vertices_size()); }
906 if (doDebug)
ATH_MSG_DEBUG(
"End while(moreP) pass " << moreP);
914 if (doDebug)
ATH_MSG_DEBUG(
"Deleting vertices " << deleteV.size());
917 for (dvItr = deleteV.begin(); dvItr != deleteV.end(); ++dvItr) {
919 if (*dvItr)
delete (*dvItr);
922 if (doDebug)
ATH_MSG_DEBUG(
"Deleting particles " << deleteP.size());
925 for (dpItr = deleteP.begin(); dpItr != deleteP.end(); ++dpItr) {
927 if (*dpItr)
delete (*dpItr);
935 if (doDebug && doExtra) {
936 std::cout <<
"========== BEGIN EVENT BEFORE SOFT ==========" << std::endl;
938 std::cout <<
"========== END EVENT BEFORE SOFT ==========" << std::endl;
943 std::list<HepMC::GenParticlePtr> pNotHad;
944 std::list<HepMC::GenParticlePtr> pHard;
946 std::vector<HepMC::GenParticlePtr> beams=thinEvt->beams();
947 for (
auto fp: thinEvt->particles()) {
951 double ep =
fp->momentum().e();
952 double pzp =
fp->momentum().pz();
953 double mtmax = (ep + pzp) * (ep - pzp);
954 auto ancestors = HepMC::ancestor_particles(
fp);
956 for (
auto gpar: ancestors) {
957 double e = gpar->momentum().e();
958 double pz = gpar->momentum().pz();
963 pNotHad.push_back(
fp);
964 int ida = std::abs(
fp->pdg_id());
965 bool keepid = (ida > 10 && ida < 20) || (ida > 1000000 && ida < 9000000);
968 pHard.insert(pHard.end(),ancestors.begin(),ancestors.end());
975 bool keepid2 = ( ida == 6) || (ida >= 11 && ida <= 16) || (ida >= 23 && ida <= 37) || (ida > 1000000 && ida < 9000000);
976 if (keepid2 &&
fp->end_vertex()) {
977 auto descendants = HepMC::descendant_particles(
fp);
978 pHard.insert(pHard.end(),descendants.begin(),descendants.end());
991 for (
auto ph: pHard) {
995 if (doDebug)
ATH_MSG_DEBUG(
"Hard GenParticles total " << nhard);
1001 beams[0] = thinEvt->beam_particles().first;
1002 beams[1] = thinEvt->beam_particles().second;
1004 HepMC::GenEvent::particle_iterator finp = thinEvt->particles_begin();
1005 HepMC::GenEvent::particle_iterator finpE = thinEvt->particles_end();
1007 for (; finp != finpE; ++finp) {
1009 HepMC::GenVertex* pvtx =
fp->production_vertex();
1010 if (!pvtx)
continue;
1012 double ep =
fp->momentum().e();
1013 double pzp =
fp->momentum().pz();
1014 double mtmax = (ep + pzp) * (ep - pzp);
1015 HepMC::GenVertex::particle_iterator gpar =
fp->production_vertex()->particles_begin(HepMC::ancestors);
1016 HepMC::GenVertex::particle_iterator gparB = gpar;
1017 HepMC::GenVertex::particle_iterator gparE =
fp->production_vertex()->particles_end(HepMC::ancestors);
1019 for (; gpar != gparE; ++gpar) {
1020 double e = (*gpar)->momentum().e();
1021 double pz = (*gpar)->momentum().pz();
1026 pNotHad.push_back(
fp);
1027 int ida = std::abs(
fp->pdg_id());
1028 bool keepid = (ida > 10 && ida < 20) || (ida > 1000000 && ida < 9000000);
1030 pHard.push_back(
fp);
1031 for (gpar = gparB; gpar != gparE; ++gpar)
1032 pHard.push_back(*gpar);
1039 bool keepid2 = ida == 6 || (ida >= 11 && ida <= 16) || (ida >= 23 && ida <= 37) || (ida > 1000000 && ida < 9000000);
1040 if (keepid2 &&
fp->end_vertex()) {
1041 HepMC::GenVertex::particle_iterator
des =
fp->end_vertex()->particles_begin(HepMC::descendants);
1042 HepMC::GenVertex::particle_iterator desE =
fp->end_vertex()->particles_end(HepMC::descendants);
1043 for (;
des != desE; ++
des)
1044 pHard.push_back(*
des);
1059 for (; hItr2 != hItr2E; ++hItr2) {
1063 if (doDebug)
ATH_MSG_DEBUG(
"Hard GenParticles total " << nhard);
1071 for (
auto p: pNotHad) {
1073 bool isHard =
false;
1074 for (
auto h: pHard) {
1080 if (doDebug)
ATH_MSG_DEBUG(
"Particle bc/isHard " <<
p <<
" " << isHard);
1081 if (isHard)
continue;
1083 if (pvtx) pvtx->remove_particle_out(
p);
1085 if (evtx) evtx->remove_particle_in(
p);
1096 for (; pItr != pItrE; ++pItr) {
1100 bool isHard =
false;
1101 for (hItr = hItrB; hItr != hItrE; ++hItr) {
1108 if (isHard)
continue;
1109 HepMC::GenVertex* pvtx =
p->production_vertex();
1110 if (pvtx) pvtx->remove_particle(
p);
1111 HepMC::GenVertex* evtx =
p->end_vertex();
1112 if (evtx) evtx->remove_particle(
p);
1125 for (
auto vtx: thinEvt->vertices()) {
1126 if (vtx->particles_in().size() != 0)
continue;
1127 if (vtx->particles_out().size() != 0)
continue;
1128 removeV.push_back(vtx);
1129 deleteV.push_back(vtx);
1131 if (doDebug)
ATH_MSG_DEBUG(
"Removing/deleting 0-particle vertices " << removeV.size() <<
" " << deleteV.size());
1132 for (
unsigned int i = 0;
i < removeV.size(); ++
i) {
1134 thinEvt->remove_vertex(removeV[
i]);
1137 HepMC::GenEvent::vertex_iterator vtx = thinEvt->vertices_begin();
1138 HepMC::GenEvent::vertex_iterator vtxE = thinEvt->vertices_end();
1139 for (; vtx != vtxE; ++vtx) {
1140 if ((*vtx)->particles_in_size() != 0)
continue;
1141 if ((*vtx)->particles_out_size() != 0)
continue;
1142 removeV.push_back(*vtx);
1143 deleteV.push_back(*vtx);
1146 if (doDebug)
ATH_MSG_DEBUG(
"Removing/deleting 0-particle vertices " << removeV.size() <<
" " << deleteV.size());
1147 for (
unsigned int i = 0;
i < removeV.size(); ++
i) {
1148 if (thinEvt->remove_vertex(removeV[
i])) {
1157 for (dvItr = deleteV.begin(); dvItr != deleteV.end(); ++dvItr) {
1159 if (*dvItr)
delete (*dvItr);
1167 if (doDebug && doExtra) {
1168 std::cout <<
"========== BEGIN EVENT BEFORE 1-BODY ==========" << std::endl;
1170 std::cout <<
"========== END EVENT BEFORE 1-BODY ==========" << std::endl;
1186 for (
auto v: thinEvt->vertices()) {
1187 if (
v->particles_in().size() != 1)
continue;
1188 if (
v->particles_out().size() != 1)
continue;
1189 pin = *(
v->particles_in().begin());
1190 pout = *(
v->particles_out().begin());
1191 if (pin->pdg_id() !=
pout->pdg_id())
continue;
1193 if (pin == beams[0] || pin == beams[1])
continue;
1195 if (!pvtx || pvtx->parent_event() ==
nullptr) {
1207 pvtx->remove_particle_in(pin);
1208 pvtx->remove_particle_out(pin);
1209 pvtx->add_particle_out(
pout);
1210 vtx11->remove_particle_in(pin);
1211 vtx11->remove_particle_out(pin);
1212 vtx11->remove_particle_in(
pout);
1213 vtx11->remove_particle_out(
pout);
1214 thinEvt->remove_vertex(vtx11);
1215 if (doDebug)
ATH_MSG_DEBUG(
"One-body new pvtx " << pvtx <<
" " << pvtx->particles_in().size() <<
" " << pvtx->particles_out().size());
1220 HepMC::GenVertex* vtx11;
1227 HepMC::GenEvent::vertex_iterator
v = thinEvt->vertices_begin();
1228 HepMC::GenEvent::vertex_iterator vE = thinEvt->vertices_end();
1231 for (;
v != vE; ++
v) {
1232 if ((*v)->particles_in_size() != 1)
continue;
1233 if ((*v)->particles_out_size() != 1)
continue;
1234 pin = *((*v)->particles_in_const_begin());
1235 pout = *((*v)->particles_out_const_begin());
1236 if (pin->pdg_id() !=
pout->pdg_id())
continue;
1238 if (pin == beams[0] || pin == beams[1])
continue;
1239 HepMC::GenVertex* pvtx = pin->production_vertex();
1252 HepMC::GenVertex* pvtx = pin->production_vertex();
1253 pvtx->remove_particle(pin);
1254 pvtx->add_particle_out(
pout);
1255 vtx11->remove_particle(pin);
1256 vtx11->remove_particle(
pout);
1257 thinEvt->remove_vertex(vtx11);
1260 if (doDebug)
ATH_MSG_DEBUG(
"One-body new pvtx " <<
HepMC::uniqueID(pvtx) <<
" " << pvtx->particles_in_size() <<
" " << pvtx->particles_out_size());
1281 for (
auto badv: thinEvt->vertices()) {
1282 if (!badv)
continue;
1283 if (badv->particles_in().size() != 1 || badv->particles_out().size() != 0)
continue;
1284 auto pitr = badv->particles_in().begin();
1286 if (pp->production_vertex())
continue;
1287 double pt = pp->momentum().perp();
1291 if (doDebug)
ATH_MSG_DEBUG(
"1->0: removing pp,badv,pt " << pp <<
" " << badv <<
" " <<
pt);
1292 removePV.push_back(vpPair(badv, pp));
1293 deleteP.push_back(pp);
1294 removeV.push_back(badv);
1295 deleteV.push_back(badv);
1299 for (
unsigned int i = 0;
i < removePV.size(); ++
i) {
1302 v->remove_particle_in(
p);
1303 v->remove_particle_out(
p);
1306 for (
unsigned int i = 0;
i < removeV.size(); ++
i) {
1307 thinEvt->remove_vertex(removeV[
i]);
1311 HepMC::GenEvent::vertex_iterator badv = thinEvt->vertices_begin();
1312 HepMC::GenEvent::vertex_iterator badvE = thinEvt->vertices_end();
1314 for (; badv != badvE; ++badv) {
1315 if (!(*badv))
continue;
1316 if ((*badv)->particles_in_size() != 1 || (*badv)->particles_out_size() != 0)
continue;
1317 HepMC::GenVertex::particles_in_const_iterator pitr = (*badv)->particles_in_const_begin();
1319 if (pp->production_vertex())
continue;
1320 double pt = pp->momentum().perp();
1324 if (doDebug)
ATH_MSG_DEBUG(
"1->0: removing pp,badv,pt " << pp <<
" " << *badv <<
" " <<
pt);
1325 removePV.emplace_back(*badv, pp);
1326 deleteP.push_back(pp);
1327 removeV.push_back(*badv);
1328 deleteV.push_back(*badv);
1333 for (
unsigned int i = 0;
i < removePV.size(); ++
i) {
1334 HepMC::GenVertex*
v = removePV[
i].first;
1336 v->remove_particle(
p);
1340 for (
unsigned int i = 0;
i < removeV.size(); ++
i) {
1341 if (!thinEvt->remove_vertex(removeV[
i])) {
ATH_MSG_WARNING(
"1->0: Failed to remove vertex " << removeV[
i]); }
1347 for (dvItr = deleteV.begin(); dvItr != deleteV.end(); ++dvItr) {
1348 if (*dvItr)
delete (*dvItr);
1353 for (dpItr = deleteP.begin(); dpItr != deleteP.end(); ++dpItr) {
1354 if (*dpItr)
delete (*dpItr);
1364 std::cout <<
"========== BEGIN EVENT AFTER THINNING ==========" << std::endl;
1366 std::cout <<
"========== END EVENT AFTER THINNING ==========" << std::endl;
1377 return StatusCode::SUCCESS;