36 if (gp->parent_event() != ge)
return;
39 HepMC::GenVertex* vstart = gp->production_vertex();
40 HepMC::GenVertex* vend = gp->end_vertex();
43 if (vstart !=
nullptr) vstart->remove_particle(gp);
44 if (vend !=
nullptr) vend->remove_particle(gp);
52 if (vstart !=
nullptr && vend !=
nullptr && vend != vstart) {
53 bool is_only_link =
true;
54 for (
auto pchild=vstart->particles_out_const_begin();pchild!=vstart->particles_out_const_end();++pchild) {
55 if ((*pchild)->end_vertex() == vend) is_only_link =
false;
58 if (vend->position() != HepMC::FourVector())
59 vstart->set_position(vend->position());
60 while (vend->particles_out_size() > 0) {
61 vstart->add_particle_out(*vend->particles_out_const_begin());
63 while (vend->particles_in_size() > 0) {
64 vstart->add_particle_in(*vend->particles_in_const_begin());
74 std::vector<HepMC::GenVertex*> orphaned_vtxs;
75 for (HepMC::GenEvent::vertex_const_iterator vi = ge->vertices_begin(); vi != ge->vertices_end(); ++vi) {
76 if ((*vi)->particles_in_size() == 0 && (*vi)->particles_out_size() == 0) orphaned_vtxs.push_back(*vi);
78 for (HepMC::GenVertex* gv : orphaned_vtxs)
delete gv;
82 inline void reduce(HepMC::GenEvent* ge, std::vector<HepMC::GenParticlePtr> toremove) {
83 while (toremove.size()) {
84 auto gp = toremove.back();
95 HepMC::GenEvent*
evt =
const_cast<HepMC::GenEvent*
>(*ievt);
101 if (newpid ==
m_pidmap.end())
continue;
102 ip->set_pdg_id(newpid->second);
107 if (
evt->weights().empty()) {
108 ATH_MSG_DEBUG(
"Adding a unit weight to empty event weight vector");
109 evt->weights().push_back(1);
113 const HepMC::FourVector nullpos;
114 for (
auto iv:
evt->vertices()) {
115 if (iv->position() == nullpos) {
116 ATH_MSG_DEBUG(
"Setting representative event position vertex");
124 std::vector<HepMC::GenParticlePtr> beams_t;
126 if (
p->status() == 4) beams_t.push_back(
p);
128 if (beams_t.size() > 2) {
129 ATH_MSG_INFO(
"Invalid number of beam particles " << beams_t.size() <<
". Will try to fix.");
130 std::vector<HepMC::GenParticlePtr> bparttoremove;
131 for (
const auto& bpart : beams_t) {
132 if (bpart->id() == 0 && bpart->production_vertex()) bparttoremove.push_back(bpart);
134 for (
auto bpart: bparttoremove) {
135 bpart->production_vertex()->remove_particle_out(bpart);
140 std::vector<HepMC::GenParticlePtr> semi_disconnected, decay_loop_particles;
141 for (
auto ip :
evt->particles()) {
144 bool particle_to_fix =
false;
145 int abspid = std::abs(
ip->pdg_id());
146 auto vProd =
ip->production_vertex();
147 auto vEnd =
ip->end_vertex();
149 if ( (!vProd || vProd->id() == 0) && vEnd &&
ip->status() != 4) {
150 particle_to_fix =
true;
151 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without production vertex! HepMC status = " <<
ip->status());
154 if (vProd && !vEnd &&
ip->status() != 1) {
155 particle_to_fix =
true;
156 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without decay vertex! HepMC status = " <<
ip->status());
158 if (particle_to_fix) semi_disconnected.push_back(
ip);
160 if (abspid == 43 || abspid == 44 || abspid == 30353 || abspid == 30343) {
161 decay_loop_particles.push_back(
ip);
168 if ( semi_disconnected.size() == 4 || semi_disconnected.size() == 3 || semi_disconnected.size() == 2) {
171 HepMC::FourVector
sum(0,0,0,0);
172 for (
const auto&
part : semi_disconnected) {
173 if (!
part->production_vertex() || !
part->production_vertex()->id()) {
174 no_prov++;
sum +=
part->momentum();
176 if (!
part->end_vertex()) {
177 no_endv++;
sum -=
part->momentum();
180 ATH_MSG_INFO(
"Heuristics: found " << semi_disconnected.size() <<
" semi-disconnected particles. Momentum sum is " <<
sum);
182 if (no_endv && no_prov && ( no_endv + no_prov == semi_disconnected.size() )) {
183 if (std::abs(
sum.px()) < 1
e-2 && std::abs(
sum.py()) < 1
e-2 && std::abs(
sum.pz()) < 1
e-2 ) {
184 ATH_MSG_INFO(
"Try " << no_endv <<
"->" << no_prov <<
" splitting/merging.");
186 for (
auto part : semi_disconnected) {
187 if (!
part->production_vertex() ||
part->production_vertex()->id() == 0)
v->add_particle_out(
part);
189 for (
auto part : semi_disconnected) {
190 if (!
part->end_vertex())
v->add_particle_in(
part);
201 for (
auto part: decay_loop_particles) {
203 auto vend =
part->end_vertex();
204 auto vprod =
part->production_vertex();
205 if (!vprod || !vend)
continue;
206 bool loop_in_decay =
true;
209 auto sisters = vend->particles_in();
210 for (
auto sister: sisters) {
211 if (vprod != sister->production_vertex()) loop_in_decay =
false;
213 if (!loop_in_decay)
continue;
216 auto daughters = vend->particles_out();
217 for (
auto p : daughters) vprod->add_particle_out(
p);
218 for (
auto sister : sisters) {
219 vprod->remove_particle_out(sister);
220 vend->remove_particle_in(sister);
221 evt->remove_particle(sister);
223 evt->remove_vertex(vend);
228 std::vector<HepMC::GenParticlePtr> toremove;
229 for (
auto ip:
evt->particles()) {
234 bool bad_particle =
false;
250 int abs_pdg_id = std::abs(
ip->pdg_id());
251 bool is_decayed_weak_boson = ( abs_pdg_id == 23 || abs_pdg_id == 24 || abs_pdg_id == 25 ) &&
ip->end_vertex();
259 if (bad_particle) toremove.push_back(
ip);
263 if (toremove.empty())
continue;
264 ATH_MSG_DEBUG(
"Cleaning event record of " << toremove.size() <<
" bad particles");
266 const int num_particles_orig =
evt->particles().size();
267 for (
auto part: toremove)
evt->remove_particle(
part);
268 const int num_particles_filt =
evt->particles().size();
270 ATH_MSG_INFO(
"Particles filtered: " << num_particles_orig <<
" -> " << num_particles_filt);
274 if (
evt->weights().empty()) {
275 ATH_MSG_DEBUG(
"Adding a unit weight to empty event weight vector");
276 evt->weights().push_back(1);
280 if (
evt->signal_process_vertex() == NULL) {
281 const HepMC::FourVector nullpos;
282 for (HepMC::GenEvent::vertex_const_iterator iv =
evt->vertices_begin(); iv !=
evt->vertices_end(); ++iv) {
283 if ((*iv)->position() == nullpos) {
284 ATH_MSG_DEBUG(
"Setting representative event position vertex");
285 evt->set_signal_process_vertex(
const_cast<HepMC::GenVertex*
>(*iv));
293 std::vector<HepMC::GenParticlePtr> semi_disconnected, decay_loop_particles;
294 for (
auto ip : *
evt) {
297 bool particle_to_fix =
false;
298 int abspid = std::abs(
ip->pdg_id());
299 auto vProd =
ip->production_vertex();
300 auto vEnd =
ip->end_vertex();
302 if ( (!vProd || vProd->id() == 0) && vEnd &&
ip->status() != 4) {
303 particle_to_fix =
true;
304 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without production vertex! HepMC status = " <<
ip->status());
307 if (vProd && !vEnd &&
ip->status() != 1) {
308 particle_to_fix =
true;
309 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without decay vertex! HepMC status = " <<
ip->status());
311 if (particle_to_fix) semi_disconnected.push_back(
ip);
313 if (abspid == 43 || abspid == 44 || abspid == 30353 || abspid == 30343) {
314 decay_loop_particles.push_back(
ip);
320 if (semi_disconnected.size() == 4 ||semi_disconnected.size() == 3 || semi_disconnected.size() == 2) {
323 double vsum[4] = {0,0,0,0};
324 for (
auto part : semi_disconnected) {
325 if (!
part->production_vertex() ) {
327 vsum[0] +=
part->momentum().px();
328 vsum[1] +=
part->momentum().py();
329 vsum[2] +=
part->momentum().pz();
330 vsum[3] +=
part->momentum().e();
332 if (!
part->end_vertex()) {
334 vsum[0] -=
part->momentum().px();
335 vsum[1] -=
part->momentum().py();
336 vsum[2] -=
part->momentum().pz();
337 vsum[3] -=
part->momentum().e();
340 HepMC::FourVector
sum(vsum[0],vsum[1],vsum[2],vsum[3]);
341 ATH_MSG_INFO(
"Heuristics: found " << semi_disconnected.size() <<
" semi-disconnected particles. Momentum sum is " << vsum[0] <<
" " << vsum[1] <<
" " << vsum[2] <<
" " << vsum[3]);
343 if (no_endv && no_prov && ( no_endv + no_prov == semi_disconnected.size() )) {
344 if (std::abs(
sum.px()) < 1
e-2 && std::abs(
sum.py()) < 1
e-2 && std::abs(
sum.pz()) < 1
e-2 ) {
345 ATH_MSG_INFO(
"Try " << no_endv <<
"->" << no_prov <<
" splitting/merging.");
347 for (
auto part : semi_disconnected) {
348 if (!
part->production_vertex())
v->add_particle_out(
part);
350 for (
auto part : semi_disconnected) {
351 if (!
part->end_vertex())
v->add_particle_in(
part);
362 for (
auto part: decay_loop_particles) {
364 auto vend =
part->end_vertex();
365 auto vprod =
part->production_vertex();
366 if (!vend || !vprod)
continue;
367 bool loop_in_decay =
true;
369 std::vector<HepMC::GenParticlePtr> sisters;
371 for (
auto sister : sisters) {
372 if (vprod != sister->production_vertex()) loop_in_decay =
false;
374 if (!loop_in_decay)
continue;
376 std::vector<HepMC::GenParticlePtr> daughters;
378 for (
auto p : daughters) vprod->add_particle_out(
p);
379 for (
auto sister : sisters) {
380 vprod->remove_particle(sister);
381 vend->remove_particle(sister);
383 evt->remove_vertex(vend);
389 std::vector<HepMC::GenParticlePtr> toremove; toremove.reserve(10);
390 for (HepMC::GenEvent::particle_const_iterator
ip =
evt->particles_begin();
ip !=
evt->particles_end(); ++
ip) {
392 if (*
ip == NULL)
continue;
396 bool bad_particle =
false;
415 int abs_pdg_id = std::abs((*ip)->pdg_id());
416 bool is_decayed_weak_boson = ( abs_pdg_id == 23 || abs_pdg_id == 24 || abs_pdg_id == 25 ) && (*ip)->end_vertex();
425 if (bad_particle) toremove.push_back(*
ip);
429 if (toremove.empty())
continue;
430 ATH_MSG_DEBUG(
"Cleaning event record of " << toremove.size() <<
" bad particles");
433 const int num_particles_orig =
evt->particles_size();
434 int num_orphan_vtxs_orig = 0;
435 int num_noparent_vtxs_orig = 0;
436 int num_nochild_vtxs_orig = 0;
437 for (
auto v =
evt->vertices_begin();
v !=
evt->vertices_end(); ++
v) {
438 if ((*v)->particles_in_size()==0&&(*v)->particles_out_size()==0) num_orphan_vtxs_orig++;
439 if ((*v)->particles_in_size()==0) num_noparent_vtxs_orig++;
440 if ((*v)->particles_out_size()==0) num_nochild_vtxs_orig++;
443 int signal_vertex_bc =
evt->signal_process_vertex() ?
evt->signal_process_vertex()->barcode() : 0;
446 if (
evt->barcode_to_vertex (signal_vertex_bc) ==
nullptr) {
447 evt->set_signal_process_vertex (
nullptr);
451 const int num_particles_filt =
evt->particles_size();
452 int num_orphan_vtxs_filt = 0;
453 int num_noparent_vtxs_filt = 0;
454 int num_nochild_vtxs_filt = 0;
455 for (
auto v =
evt->vertices_begin();
v !=
evt->vertices_end(); ++
v) {
456 if ((*v)->particles_in_size()==0&&(*v)->particles_out_size()==0) num_orphan_vtxs_filt++;
457 if ((*v)->particles_in_size()==0) num_noparent_vtxs_filt++;
458 if ((*v)->particles_out_size()==0) num_nochild_vtxs_filt++;
462 ATH_MSG_INFO(
"Particles filtered: " << num_particles_orig <<
" -> " << num_particles_filt);
464 if (num_orphan_vtxs_filt != num_orphan_vtxs_orig)
465 ATH_MSG_WARNING(
"Change in orphaned vertices: " << num_orphan_vtxs_orig <<
" -> " << num_orphan_vtxs_filt);
466 if (num_noparent_vtxs_filt != num_noparent_vtxs_orig)
467 ATH_MSG_WARNING(
"Change in no-parent vertices: " << num_noparent_vtxs_orig <<
" -> " << num_noparent_vtxs_filt);
468 if (num_nochild_vtxs_filt != num_nochild_vtxs_orig)
469 ATH_MSG_WARNING(
"Change in no-parent vertices: " << num_nochild_vtxs_orig <<
" -> " << num_nochild_vtxs_filt);
473 return StatusCode::SUCCESS;
482 return StatusCode::SUCCESS;
491 return p->pdg_id() == 0;
496 if (!
p)
return false;
497 auto v=
p->production_vertex();
498 if (!
v)
return false;
500 for (
const auto& anc:
v->particles_in()) {
503 if (storage->find(
p->id()) != storage->end())
return false;
504 storage->insert(
p->id());
505 for (
const auto& anc:
v->particles_in()) {
506 if (
fromDecay(anc, storage))
return true;
509 for (
auto anc=
v->particles_in_const_begin(); anc !=
v->particles_in_const_end(); ++anc) {
512 if (storage->find(
p->barcode()) != storage->end())
return false;
513 storage->insert(
p->barcode());
514 for (
auto anc=
v->particles_in_const_begin(); anc !=
v->particles_in_const_end(); ++anc){
515 if (
fromDecay(*anc, storage))
return true;
523 auto storage = std::make_shared<std::set<int>>();
529 if (
p->production_vertex() ==
p->end_vertex() &&
p->end_vertex() != NULL)
return true;
533 for (
auto itrParent: *(
p->production_vertex())) {