119 bool doDebug =
false;
120 bool doExtra =
false;
128 return StatusCode::FAILURE;
131 if (mcEvts->
empty()) {
133 return StatusCode::SUCCESS;
140 delete thinnedMcEvts;
141 thinnedMcEvts =
nullptr;
142 return StatusCode::FAILURE;
147 const HepMC::GenEvent* mcEvt = mcEvts->
front();
148 auto wtCont = mcEvt->weights();
149 if (!wtCont.empty()) {
153 int inEvent = mcEvt->event_number();
160 HepMC::GenEvent* thinEvt =
new HepMC::GenEvent(*mcEvt);
161 int nEvent = thinEvt->event_number();
165 std::cout <<
"========== BEGIN EVENT BEFORE THINNING ==========" << std::endl;
167 std::cout <<
"========== END EVENT BEFORE THINNING ==========" << std::endl;
187 typedef std::pair<HepMC::GenVertexPtr, HepMC::GenParticlePtr> vpPair;
188 std::vector<vpPair> removePV;
189 std::vector<vpPair> addinPV;
190 std::vector<vpPair> addoutPV;
191 std::vector<HepMC::GenVertexPtr> removeV;
194 std::list<HepMC::GenParticlePtr> deleteP;
195 std::list<HepMC::GenVertexPtr> deleteV;
207 std::vector<HepMC::GenVertexPtr> hadVertices;
210 for (
auto& hadv: thinEvt->vertices()) {
212 if (hadv->particles_in().size() < 2)
continue;
213 if (hadv->particles_out().size() < 1)
continue;
217 bool isHadVtx =
true;
218 bool isHadOut =
false;
219 for (
const auto& inp: hadv->particles_in() ) {
222 for (
const auto& vp: hadv->particles_out()) {
226 isHadVtx = isHadVtx && isHadOut;
227 if (isHadVtx) hadVertices.push_back(hadv);
228 if (doDebug && isHadVtx)
ATH_MSG_VERBOSE(
"Hadronization vertex " << hadv);
231 HepMC::GenEvent::vertex_iterator hadv = thinEvt->vertices_begin();
232 HepMC::GenEvent::vertex_iterator hadvB = thinEvt->vertices_begin();
233 HepMC::GenEvent::vertex_iterator hadvE = thinEvt->vertices_end();
235 for (; hadv != hadvE; ++hadv) {
236 if (!(*hadv))
continue;
237 if ((*hadv)->particles_in_size() < 2)
continue;
238 if ((*hadv)->particles_out_size() < 1)
continue;
243 bool isHadVtx =
true;
244 bool isHadOut =
false;
245 HepMC::GenVertex::particles_in_const_iterator inp = (*hadv)->particles_in_const_begin();
246 HepMC::GenVertex::particles_in_const_iterator inpE = (*hadv)->particles_in_const_end();
247 for (; inp != inpE; ++inp) {
250 HepMC::GenVertex::particles_out_const_iterator vp = (*hadv)->particles_out_const_begin();
251 HepMC::GenVertex::particles_out_const_iterator vpE = (*hadv)->particles_out_const_end();
252 for (; vp != vpE; ++vp) {
256 isHadVtx = isHadVtx && isHadOut;
257 if (isHadVtx) hadVertices.push_back(*hadv);
258 if (doDebug && isHadVtx)
ATH_MSG_VERBOSE(
"Hadronization vertex " << *hadv);
262 if (hadVertices.empty()) {
266 return StatusCode::SUCCESS;
275 for (
unsigned int iv = 0; iv < hadVertices.size(); ++iv) {
277 if (doDebug)
ATH_MSG_DEBUG(
"Removing partons from hadVertex " << ivtx);
278 for (
auto pin: ivtx->particles_in()) {
279 removePV.push_back(vpPair(ivtx, pin));
284 for (
unsigned int iv = 0; iv < hadVertices.size(); ++iv) {
286 auto descendants = HepMC::descendant_particles(ivtx);
287 for (
auto pout: descendants) {
289 if (vpar) removePV.push_back(vpPair(vpar,
pout));
291 if (vend) removePV.push_back(vpPair(vend,
pout));
292 deleteP.push_back(
pout);
296 for (
unsigned int iv = 0; iv < hadVertices.size(); ++iv) {
297 HepMC::GenVertex* ivtx = hadVertices[iv];
298 if (doDebug)
ATH_MSG_DEBUG(
"Removing partons from hadVertex " << ivtx);
299 HepMC::GenVertex::particles_in_const_iterator pin = ivtx->particles_in_const_begin();
300 HepMC::GenVertex::particles_in_const_iterator pinE = ivtx->particles_in_const_end();
301 for (; pin != pinE; ++pin) {
302 removePV.emplace_back(ivtx, *pin);
308 for (
unsigned int iv = 0; iv < hadVertices.size(); ++iv) {
309 HepMC::GenVertex* ivtx = hadVertices[iv];
310 HepMC::GenVertex::particle_iterator
pout = ivtx->particles_begin(HepMC::descendants);
311 HepMC::GenVertex::particle_iterator poutE = ivtx->particles_end(HepMC::descendants);
313 HepMC::GenVertex* vpar = (*pout)->production_vertex();
314 if (vpar) removePV.emplace_back(vpar, *
pout);
315 HepMC::GenVertex* vend = (*pout)->end_vertex();
316 if (vend) removePV.emplace_back(vend, *
pout);
317 deleteP.push_back(*
pout);
328 for (
auto hadv:thinEvt->vertices()) {
330 if (hadv->particles_in().size() == 0 && hadv->particles_out().size() == 0) {
331 removeV.push_back(hadv);
332 deleteV.push_back(hadv);
336 for (
auto pin: hadv->particles_in()) {
337 removePV.push_back(vpPair(hadv, pin));
340 for (
auto pout: hadv->particles_out()) {
341 removePV.push_back(vpPair(hadv,
pout));
344 removeV.push_back(hadv);
345 deleteV.push_back(hadv);
348 for (hadv = hadvB; hadv != hadvE; ++hadv) {
351 if ((*hadv)->particles_in_size() == 0 && (*hadv)->particles_out_size() == 0) {
352 removeV.push_back(*hadv);
353 deleteV.push_back(*hadv);
358 HepMC::GenVertex::particles_in_const_iterator pin = (*hadv)->particles_in_const_begin();
359 HepMC::GenVertex::particles_in_const_iterator pinE = (*hadv)->particles_in_const_end();
360 for (; pin != pinE; ++pin) {
361 removePV.emplace_back(*hadv, *pin);
364 HepMC::GenVertex::particles_out_const_iterator
pout = (*hadv)->particles_out_const_begin();
365 HepMC::GenVertex::particles_out_const_iterator poutE = (*hadv)->particles_out_const_end();
367 removePV.emplace_back(*hadv, *
pout);
370 removeV.push_back(*hadv);
371 deleteV.push_back(*hadv);
377 for (
unsigned int i = 0;
i < removePV.size(); ++
i) {
381 v->remove_particle_in(
p);
382 v->remove_particle_out(
p);
384 v->remove_particle(
p);
388 for (
unsigned int i = 0;
i < addoutPV.size(); ++
i) {
391 v->add_particle_out(
p);
394 for (
unsigned int iv = 1; iv < hadVertices.size(); ++iv) {
396 if (
v->particles_in().size() != 0 ||
v->particles_out().size() != 0) {
398 <<
" with in/out particles " <<
v->particles_in().size() <<
" " <<
v->particles_out().size());
400 thinEvt->remove_vertex(hadVertices[iv]);
405 for (
unsigned int iv = 1; iv < hadVertices.size(); ++iv) {
406 HepMC::GenVertex*
v = hadVertices[iv];
407 if (
v->particles_in_size() != 0 ||
v->particles_out_size() != 0) {
408 ATH_MSG_WARNING(
"Removing vertex " <<
HepMC::uniqueID(
v) <<
" for event " << nEvent <<
" with in/out particles " <<
v->particles_in_size() <<
" " <<
v->particles_out_size());
415 if (doDebug)
ATH_MSG_DEBUG(
"Deleting hadronization vertices " << deleteV.size());
418 for (dvItr = deleteV.begin(); dvItr != deleteV.end(); ++dvItr) {
420 if (*dvItr)
delete (*dvItr);
425 for (dpItr = deleteP.begin(); dpItr != deleteP.end(); ++dpItr) {
427 if (*dpItr)
delete (*dpItr);
435 if (doDebug && doExtra) {
436 std::cout <<
"========== BEGIN EVENT BEFORE CLUSTER ==========" << std::endl;
438 std::cout <<
"========== END EVENT BEFORE CLUSTER ==========" << std::endl;
453 using vpPair = std::pair<HepMC::GenVertexPtr, HepMC::GenParticlePtr>;
461 using pkPair = std::pair<HepMC::GenParticlePtr, HepMC::FourVector>;
462 std::vector<pkPair> changePK;
466 if (doDebug)
ATH_MSG_DEBUG(
"New parton pass " << inEvent <<
" " << thinEvt->particles_size() <<
" " << thinEvt->vertices_size());
476 for (
auto fp: thinEvt->particles() ) {
492 if (pvtx->particles_in().size() == 1 && pvtx->particles_out().size() == 1) {
495 if (!pp || pp->parent_event() ==
nullptr) {
502 if (!ppvtx || ppvtx->parent_event() ==
nullptr) {
509 removePV.push_back(vpPair(ppvtx, pp));
510 removePV.push_back(vpPair(pvtx, pp));
511 deleteP.push_back(pp);
512 removeV.push_back(pvtx);
513 deleteV.push_back(pvtx);
514 addoutPV.push_back(vpPair(ppvtx,
fp));
515 if (doDebug) {
ATH_MSG_DEBUG(
"1->1: ppvtx,pp,pvtx,fp,evtx " << ppvtx <<
" " << pp <<
" " << pvtx <<
" " <<
fp); }
524 if (pvtx->particles_in().size() == 2 && pvtx->particles_out().size() == 1) {
530 if (std::abs(pp1->momentum().perp()) < 1.e-3)
continue;
531 if (std::abs(pp2->momentum().perp()) < 1.e-3)
continue;
535 if (!ppvtx1 || ppvtx1->parent_event() ==
nullptr) {
540 if (!ppvtx2 || ppvtx2->parent_event() ==
nullptr) {
547 removePV.push_back(vpPair(pvtx,
fp));
548 removePV.push_back(vpPair(pvtx, pp1));
549 removePV.push_back(vpPair(pvtx, pp2));
550 deleteP.push_back(
fp);
551 removeV.push_back(pvtx);
552 deleteV.push_back(pvtx);
554 ATH_MSG_DEBUG(
"2->1: ppvtx1,pp1,ppvtx2,pp2,pvtx,fp " << ppvtx1 <<
" " << pp1
555 <<
" " << ppvtx2 <<
" " << pp2 <<
" " << pvtx <<
" "<<
fp);
564 if (pvtx->particles_in().size() == 1 && pvtx->particles_out().size() == 2) {
573 }
else if (
fp == pout2) {
579 ATH_MSG_WARNING(
"1->2: No match found for branching " <<
fp <<
" " << pvtx <<
" " << pout1 <<
" " << pout2);
582 if (
fp != pout1)
continue;
586 if (pout1->momentum().m() < -1.0 || pout1->momentum().perp() < 1.0)
continue;
587 if (pout2->momentum().m() < -1.0 || pout2->momentum().perp() < 1.0)
continue;
590 double m12 = p12.m();
592 if (std::abs(m12) > 10. + 1.0
e-5 * p12.e()) {
593 ATH_MSG_WARNING(
"Spacelike mass^2 for parton sum " << m12 <<
" " << pp <<
" " << pvtx <<
" " << pout1 <<
" " << pout2);
597 if (doDebug)
ATH_MSG_DEBUG(
"1->2: parton pair mass " << m12);
600 if (doDebug)
ATH_MSG_DEBUG(
"Keeping 1->2: parton mass " << m12);
605 if (!ppvtx || ppvtx->parent_event() ==
nullptr) {
613 if (doDebug)
ATH_MSG_DEBUG(
"Merging 1->2: mass " << p12.m());
614 changePK.push_back(pkPair(pp, p12));
615 removePV.push_back(vpPair(pvtx, pp));
616 removePV.push_back(vpPair(pvtx, pout1));
617 removePV.push_back(vpPair(pvtx, pout2));
618 deleteP.push_back(pout1);
619 deleteP.push_back(pout2);
620 removeV.push_back(pvtx);
621 deleteV.push_back(pvtx);
623 ATH_MSG_DEBUG(
"Merge 1->2: ppvtx,pp,pvtx,pout1,pout2,evtx " << ppvtx <<
" " << pp <<
" " << pvtx <<
" " << pout1 <<
" "
625 ATH_MSG_DEBUG(
"Merge 1->2: id " << pp->pdg_id() <<
" " << pout1->pdg_id() <<
" " << pout2->pdg_id());
632 if (pvtx->particles_in().size() == 1) {
635 if (std::abs(pp->pdg_id()) == MC::PROTON) iCase = -1;
640 if (doDebug)
ATH_MSG_DEBUG(
"Case not found " << pvtx <<
" " <<
fp <<
" " << pvtx->particles_in().size() <<
" " << pvtx->particles_out().size());
645 HepMC::GenEvent::particle_iterator finp = thinEvt->particles_begin();
646 HepMC::GenEvent::particle_iterator finpB = thinEvt->particles_begin();
647 HepMC::GenEvent::particle_iterator finpE = thinEvt->particles_end();
650 for (finp = finpB; finp != finpE; ++finp) {
658 HepMC::GenVertex* pvtx =
fp->production_vertex();
672 if (pvtx->particles_in_size() == 1 && pvtx->particles_out_size() == 1) {
674 HepMC::GenVertex::particles_in_const_iterator pitr = pvtx->particles_in_const_begin();
682 HepMC::GenVertex* ppvtx = pp->production_vertex();
691 removePV.emplace_back(ppvtx, pp);
692 removePV.emplace_back(pvtx, pp);
693 deleteP.push_back(pp);
694 removeV.push_back(pvtx);
695 deleteV.push_back(pvtx);
696 addoutPV.emplace_back(ppvtx,
fp);
709 if (pvtx->particles_in_size() == 2 && pvtx->particles_out_size() == 1) {
711 HepMC::GenVertex::particles_in_const_iterator pitr = pvtx->particles_in_const_begin();
718 if (std::abs(pp1->momentum().perp()) < 1.e-3)
continue;
719 if (std::abs(pp2->momentum().perp()) < 1.e-3)
continue;
721 HepMC::GenVertex* ppvtx1 = pp1->production_vertex();
722 HepMC::GenVertex* ppvtx2 = pp2->production_vertex();
737 removePV.emplace_back(pvtx,
fp);
738 removePV.emplace_back(pvtx, pp1);
739 removePV.emplace_back(pvtx, pp2);
740 deleteP.push_back(
fp);
741 removeV.push_back(pvtx);
742 deleteV.push_back(pvtx);
758 if (pvtx->particles_in_size() == 1 && pvtx->particles_out_size() == 2) {
759 HepMC::GenVertex::particles_out_const_iterator poutitr = pvtx->particles_out_const_begin();
770 }
else if (
fp == pout2) {
779 if (
fp != pout1)
continue;
781 HepMC::GenVertex::particles_in_const_iterator pitr = pvtx->particles_in_const_begin();
785 if (pout1->momentum().m() < -1.0 || pout1->momentum().perp() < 1.0)
continue;
786 if (pout2->momentum().m() < -1.0 || pout2->momentum().perp() < 1.0)
continue;
790 double m12 = p12.m();
792 if (fabs(m12) > 10. + 1.0
e-5 * p12.e()) {
797 if (doDebug)
ATH_MSG_DEBUG(
"1->2: parton pair mass " << m12);
800 if (doDebug)
ATH_MSG_DEBUG(
"Keeping 1->2: parton mass " << m12);
805 HepMC::GenVertex* ppvtx = pp->production_vertex();
815 if (doDebug)
ATH_MSG_DEBUG(
"Merging 1->2: mass " << p12.m());
817 changePK.emplace_back(pp, p12);
818 removePV.emplace_back(pvtx, pp);
819 removePV.emplace_back(pvtx, pout1);
820 removePV.emplace_back(pvtx, pout2);
822 deleteP.push_back(pout1);
823 deleteP.push_back(pout2);
824 removeV.push_back(pvtx);
825 deleteV.push_back(pvtx);
828 ATH_MSG_DEBUG(
"Merge 1->2: ppvtx,pp,pvtx,pout1,pout2,evtx " << ppvtx <<
" " << pp <<
" " << pvtx <<
" " << pout1 <<
" "
830 ATH_MSG_DEBUG(
"Merge 1->2: id " << pp->pdg_id() <<
" " << pout1->pdg_id() <<
" " << pout2->pdg_id());
839 if (pvtx->particles_in_size() == 1) {
841 HepMC::GenVertex::particles_in_const_iterator pitr = pvtx->particles_in_const_begin();
843 if (std::abs(pp->pdg_id()) == MC::PROTON) iCase = -1;
857 if (doDebug)
ATH_MSG_DEBUG(
"Actually removing particles " << removePV.size());
859 for (
unsigned int i = 0;
i < removePV.size(); ++
i) {
864 v->remove_particle_in(
p);
865 v->remove_particle_out(
p);
867 v->remove_particle(
p);
872 if (doDebug)
ATH_MSG_DEBUG(
"Actually add particles in/out " << addinPV.size() <<
" " << addoutPV.size());
873 for (
unsigned int i = 0;
i < addoutPV.size(); ++
i) {
877 v->add_particle_out(
p);
881 for (
unsigned int i = 0;
i < changePK.size(); ++
i) {
884 pp->set_momentum(changePK[
i].
second);
888 if (doDebug)
ATH_MSG_DEBUG(
"Actually remove vertices " << removeV.size());
889 for (
unsigned int i = 0;
i < removeV.size(); ++
i) {
891 int nv = thinEvt->vertices().size();
892 if (doDebug) {
ATH_MSG_VERBOSE(
"Removing vertex " << removeV[
i] <<
" " << nv <<
" " << thinEvt->vertices().size()); }
893 thinEvt->remove_vertex(removeV[
i]);
895 int nv = thinEvt->vertices_size();
896 if (thinEvt->remove_vertex(removeV[
i])) {
897 if (doDebug) {
ATH_MSG_VERBOSE(
"Removed vertex " << removeV[
i] <<
" " << nv <<
" " << thinEvt->vertices_size()); }
903 if (doDebug)
ATH_MSG_DEBUG(
"End while(moreP) pass " << moreP);
911 if (doDebug)
ATH_MSG_DEBUG(
"Deleting vertices " << deleteV.size());
914 for (dvItr = deleteV.begin(); dvItr != deleteV.end(); ++dvItr) {
916 if (*dvItr)
delete (*dvItr);
919 if (doDebug)
ATH_MSG_DEBUG(
"Deleting particles " << deleteP.size());
922 for (dpItr = deleteP.begin(); dpItr != deleteP.end(); ++dpItr) {
924 if (*dpItr)
delete (*dpItr);
932 if (doDebug && doExtra) {
933 std::cout <<
"========== BEGIN EVENT BEFORE SOFT ==========" << std::endl;
935 std::cout <<
"========== END EVENT BEFORE SOFT ==========" << std::endl;
940 std::list<HepMC::GenParticlePtr> pNotHad;
941 std::list<HepMC::GenParticlePtr> pHard;
943 std::vector<HepMC::GenParticlePtr> beams=thinEvt->beams();
944 for (
auto fp: thinEvt->particles()) {
948 double ep =
fp->momentum().e();
949 double pzp =
fp->momentum().pz();
950 double mtmax = (ep + pzp) * (ep - pzp);
951 auto ancestors = HepMC::ancestor_particles(
fp);
953 for (
auto gpar: ancestors) {
954 double e = gpar->momentum().e();
955 double pz = gpar->momentum().pz();
960 pNotHad.push_back(
fp);
961 int ida = std::abs(
fp->pdg_id());
962 bool keepid = (ida > 10 && ida < 20) || (ida > 1000000 && ida < 9000000);
965 pHard.insert(pHard.end(),ancestors.begin(),ancestors.end());
972 bool keepid2 = ( ida == 6) || (ida >= 11 && ida <= 16) || (ida >= 23 && ida <= 37) || (ida > 1000000 && ida < 9000000);
973 if (keepid2 &&
fp->end_vertex()) {
974 auto descendants = HepMC::descendant_particles(
fp);
975 pHard.insert(pHard.end(),descendants.begin(),descendants.end());
988 for (
auto ph: pHard) {
992 if (doDebug)
ATH_MSG_DEBUG(
"Hard GenParticles total " << nhard);
998 beams[0] = thinEvt->beam_particles().first;
999 beams[1] = thinEvt->beam_particles().second;
1001 HepMC::GenEvent::particle_iterator finp = thinEvt->particles_begin();
1002 HepMC::GenEvent::particle_iterator finpE = thinEvt->particles_end();
1004 for (; finp != finpE; ++finp) {
1006 HepMC::GenVertex* pvtx =
fp->production_vertex();
1007 if (!pvtx)
continue;
1009 double ep =
fp->momentum().e();
1010 double pzp =
fp->momentum().pz();
1011 double mtmax = (ep + pzp) * (ep - pzp);
1012 HepMC::GenVertex::particle_iterator gpar =
fp->production_vertex()->particles_begin(HepMC::ancestors);
1013 HepMC::GenVertex::particle_iterator gparB = gpar;
1014 HepMC::GenVertex::particle_iterator gparE =
fp->production_vertex()->particles_end(HepMC::ancestors);
1016 for (; gpar != gparE; ++gpar) {
1017 double e = (*gpar)->momentum().e();
1018 double pz = (*gpar)->momentum().pz();
1023 pNotHad.push_back(
fp);
1024 int ida = std::abs(
fp->pdg_id());
1025 bool keepid = (ida > 10 && ida < 20) || (ida > 1000000 && ida < 9000000);
1027 pHard.push_back(
fp);
1028 for (gpar = gparB; gpar != gparE; ++gpar)
1029 pHard.push_back(*gpar);
1036 bool keepid2 = ida == 6 || (ida >= 11 && ida <= 16) || (ida >= 23 && ida <= 37) || (ida > 1000000 && ida < 9000000);
1037 if (keepid2 &&
fp->end_vertex()) {
1038 HepMC::GenVertex::particle_iterator des =
fp->end_vertex()->particles_begin(HepMC::descendants);
1039 HepMC::GenVertex::particle_iterator desE =
fp->end_vertex()->particles_end(HepMC::descendants);
1040 for (; des != desE; ++des)
1041 pHard.push_back(*des);
1056 for (; hItr2 != hItr2E; ++hItr2) {
1060 if (doDebug)
ATH_MSG_DEBUG(
"Hard GenParticles total " << nhard);
1068 for (
auto p: pNotHad) {
1070 bool isHard =
false;
1071 for (
auto h: pHard) {
1077 if (doDebug)
ATH_MSG_DEBUG(
"Particle bc/isHard " <<
p <<
" " << isHard);
1078 if (isHard)
continue;
1080 if (pvtx) pvtx->remove_particle_out(
p);
1082 if (evtx) evtx->remove_particle_in(
p);
1093 for (; pItr != pItrE; ++pItr) {
1097 bool isHard =
false;
1098 for (hItr = hItrB; hItr != hItrE; ++hItr) {
1105 if (isHard)
continue;
1106 HepMC::GenVertex* pvtx =
p->production_vertex();
1107 if (pvtx) pvtx->remove_particle(
p);
1108 HepMC::GenVertex* evtx =
p->end_vertex();
1109 if (evtx) evtx->remove_particle(
p);
1122 for (
auto vtx: thinEvt->vertices()) {
1123 if (vtx->particles_in().size() != 0)
continue;
1124 if (vtx->particles_out().size() != 0)
continue;
1125 removeV.push_back(vtx);
1126 deleteV.push_back(vtx);
1128 if (doDebug)
ATH_MSG_DEBUG(
"Removing/deleting 0-particle vertices " << removeV.size() <<
" " << deleteV.size());
1129 for (
unsigned int i = 0;
i < removeV.size(); ++
i) {
1131 thinEvt->remove_vertex(removeV[
i]);
1134 HepMC::GenEvent::vertex_iterator vtx = thinEvt->vertices_begin();
1135 HepMC::GenEvent::vertex_iterator vtxE = thinEvt->vertices_end();
1136 for (; vtx != vtxE; ++vtx) {
1137 if ((*vtx)->particles_in_size() != 0)
continue;
1138 if ((*vtx)->particles_out_size() != 0)
continue;
1139 removeV.push_back(*vtx);
1140 deleteV.push_back(*vtx);
1143 if (doDebug)
ATH_MSG_DEBUG(
"Removing/deleting 0-particle vertices " << removeV.size() <<
" " << deleteV.size());
1144 for (
unsigned int i = 0;
i < removeV.size(); ++
i) {
1145 if (thinEvt->remove_vertex(removeV[
i])) {
1154 for (dvItr = deleteV.begin(); dvItr != deleteV.end(); ++dvItr) {
1156 if (*dvItr)
delete (*dvItr);
1164 if (doDebug && doExtra) {
1165 std::cout <<
"========== BEGIN EVENT BEFORE 1-BODY ==========" << std::endl;
1167 std::cout <<
"========== END EVENT BEFORE 1-BODY ==========" << std::endl;
1183 for (
auto v: thinEvt->vertices()) {
1184 if (
v->particles_in().size() != 1)
continue;
1185 if (
v->particles_out().size() != 1)
continue;
1186 pin = *(
v->particles_in().begin());
1187 pout = *(
v->particles_out().begin());
1188 if (pin->pdg_id() !=
pout->pdg_id())
continue;
1190 if (pin == beams[0] || pin == beams[1])
continue;
1192 if (!pvtx || pvtx->parent_event() ==
nullptr) {
1204 pvtx->remove_particle_in(pin);
1205 pvtx->remove_particle_out(pin);
1206 pvtx->add_particle_out(
pout);
1207 vtx11->remove_particle_in(pin);
1208 vtx11->remove_particle_out(pin);
1209 vtx11->remove_particle_in(
pout);
1210 vtx11->remove_particle_out(
pout);
1211 thinEvt->remove_vertex(vtx11);
1212 if (doDebug)
ATH_MSG_DEBUG(
"One-body new pvtx " << pvtx <<
" " << pvtx->particles_in().size() <<
" " << pvtx->particles_out().size());
1217 HepMC::GenVertex* vtx11;
1224 HepMC::GenEvent::vertex_iterator
v = thinEvt->vertices_begin();
1225 HepMC::GenEvent::vertex_iterator vE = thinEvt->vertices_end();
1228 for (;
v != vE; ++
v) {
1229 if ((*v)->particles_in_size() != 1)
continue;
1230 if ((*v)->particles_out_size() != 1)
continue;
1231 pin = *((*v)->particles_in_const_begin());
1232 pout = *((*v)->particles_out_const_begin());
1233 if (pin->pdg_id() !=
pout->pdg_id())
continue;
1235 if (pin == beams[0] || pin == beams[1])
continue;
1236 HepMC::GenVertex* pvtx = pin->production_vertex();
1249 HepMC::GenVertex* pvtx = pin->production_vertex();
1250 pvtx->remove_particle(pin);
1251 pvtx->add_particle_out(
pout);
1252 vtx11->remove_particle(pin);
1253 vtx11->remove_particle(
pout);
1254 thinEvt->remove_vertex(vtx11);
1257 if (doDebug)
ATH_MSG_DEBUG(
"One-body new pvtx " <<
HepMC::uniqueID(pvtx) <<
" " << pvtx->particles_in_size() <<
" " << pvtx->particles_out_size());
1278 for (
auto badv: thinEvt->vertices()) {
1279 if (!badv)
continue;
1280 if (badv->particles_in().size() != 1 || badv->particles_out().size() != 0)
continue;
1281 auto pitr = badv->particles_in().begin();
1283 if (pp->production_vertex())
continue;
1284 double pt = pp->momentum().perp();
1288 if (doDebug)
ATH_MSG_DEBUG(
"1->0: removing pp,badv,pt " << pp <<
" " << badv <<
" " <<
pt);
1289 removePV.push_back(vpPair(badv, pp));
1290 deleteP.push_back(pp);
1291 removeV.push_back(badv);
1292 deleteV.push_back(badv);
1296 for (
unsigned int i = 0;
i < removePV.size(); ++
i) {
1299 v->remove_particle_in(
p);
1300 v->remove_particle_out(
p);
1303 for (
unsigned int i = 0;
i < removeV.size(); ++
i) {
1304 thinEvt->remove_vertex(removeV[
i]);
1308 HepMC::GenEvent::vertex_iterator badv = thinEvt->vertices_begin();
1309 HepMC::GenEvent::vertex_iterator badvE = thinEvt->vertices_end();
1311 for (; badv != badvE; ++badv) {
1312 if (!(*badv))
continue;
1313 if ((*badv)->particles_in_size() != 1 || (*badv)->particles_out_size() != 0)
continue;
1314 HepMC::GenVertex::particles_in_const_iterator pitr = (*badv)->particles_in_const_begin();
1316 if (pp->production_vertex())
continue;
1317 double pt = pp->momentum().perp();
1321 if (doDebug)
ATH_MSG_DEBUG(
"1->0: removing pp,badv,pt " << pp <<
" " << *badv <<
" " <<
pt);
1322 removePV.emplace_back(*badv, pp);
1323 deleteP.push_back(pp);
1324 removeV.push_back(*badv);
1325 deleteV.push_back(*badv);
1330 for (
unsigned int i = 0;
i < removePV.size(); ++
i) {
1331 HepMC::GenVertex*
v = removePV[
i].first;
1333 v->remove_particle(
p);
1337 for (
unsigned int i = 0;
i < removeV.size(); ++
i) {
1338 if (!thinEvt->remove_vertex(removeV[
i])) {
ATH_MSG_WARNING(
"1->0: Failed to remove vertex " << removeV[
i]); }
1344 for (dvItr = deleteV.begin(); dvItr != deleteV.end(); ++dvItr) {
1345 if (*dvItr)
delete (*dvItr);
1350 for (dpItr = deleteP.begin(); dpItr != deleteP.end(); ++dpItr) {
1351 if (*dpItr)
delete (*dpItr);
1361 std::cout <<
"========== BEGIN EVENT AFTER THINNING ==========" << std::endl;
1363 std::cout <<
"========== END EVENT AFTER THINNING ==========" << std::endl;
1374 return StatusCode::SUCCESS;