AV: In case we have 3 particles, we try to add a vertex that correspond to 1->2 and 1->1 splitting. AV: we can try to do that for 4 particles as well.
Remove loops inside decay chains AV: Please note that this approach would distort some branching ratios. If some particle would have decay products with bad PDG ids, after the operation below the visible branching ratio of these decays would be zero.
Check that all particles coming into the decay vertex of bad particle cam from the same production vertex.
97 HepMC::GenEvent*
evt =
const_cast<HepMC::GenEvent*
>(*ievt);
103 ATH_MSG_DEBUG(
"Please use MC::Loops::findLoops for this event to obtain all particles and vertices in the loops");
106 for (
auto ip: *
evt) {
110 if (newpid ==
m_pidmap.end())
continue;
111 ip->set_pdg_id(newpid->second);
116 if (
evt->weights().empty()) {
117 ATH_MSG_DEBUG(
"Adding a unit weight to empty event weight vector");
118 evt->weights().push_back(1);
122 const HepMC::FourVector nullpos;
123 for (
auto iv:
evt->vertices()) {
124 if (iv->position() == nullpos) {
125 ATH_MSG_DEBUG(
"Setting representative event position vertex");
133 std::vector<HepMC::GenParticlePtr> beams_t;
135 if (
p->status() == 4) beams_t.push_back(
p);
137 if (beams_t.size() > 2) {
138 ATH_MSG_INFO(
"Invalid number of beam particles " << beams_t.size() <<
". Will try to fix.");
139 std::vector<HepMC::GenParticlePtr> bparttoremove;
140 for (
const auto& bpart : beams_t) {
141 if (bpart->id() == 0 && bpart->production_vertex()) bparttoremove.push_back(bpart);
143 for (
auto bpart: bparttoremove) {
144 bpart->production_vertex()->remove_particle_out(bpart);
149 std::vector<HepMC::GenParticlePtr> semi_disconnected, decay_loop_particles;
150 for (
auto ip :
evt->particles()) {
153 bool particle_to_fix =
false;
154 int abspid = std::abs(
ip->pdg_id());
155 auto vProd =
ip->production_vertex();
156 auto vEnd =
ip->end_vertex();
158 if ( (!vProd || vProd->id() == 0) && vEnd &&
ip->status() != 4) {
159 particle_to_fix =
true;
160 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without production vertex! HepMC status = " <<
ip->status());
163 if (vProd && !vEnd &&
ip->status() != 1) {
164 particle_to_fix =
true;
165 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without decay vertex! HepMC status = " <<
ip->status());
167 if (particle_to_fix) semi_disconnected.push_back(
ip);
169 if (abspid == 43 || abspid == 44 || abspid == 30353 || abspid == 30343) {
170 decay_loop_particles.push_back(std::move(
ip));
182 if ( !
m_ignoreSemiDisconnected && (semi_disconnected.size() == 4 || semi_disconnected.size() == 3 || semi_disconnected.size() == 2)) {
185 HepMC::FourVector
sum(0,0,0,0);
187 for (
const auto&
part : semi_disconnected) {
188 if (!
part->production_vertex() || !
part->production_vertex()->id()) {
189 no_prov++;
sum +=
part->momentum();
191 if (!
part->end_vertex()) {
192 no_endv++;
sum -=
part->momentum();
197 ATH_MSG_INFO(
"Heuristics: found " << semi_disconnected.size() <<
" semi-disconnected particles. Momentum sum is " <<
sum <<
" Standalone vertices " <<
standalone.size());
198 bool standalonevertex = (
standalone.size() == 1 && (*
standalone.begin())->particles_in().size() + (*
standalone.begin())->particles_out().size() == semi_disconnected.size());
200 if (! standalonevertex && no_endv && no_prov && ( no_endv + no_prov == semi_disconnected.size() )) {
201 if (std::abs(
sum.px()) < 1
e-2 && std::abs(
sum.py()) < 1
e-2 && std::abs(
sum.pz()) < 1
e-2 ) {
202 ATH_MSG_INFO(
"Try " << no_endv <<
"->" << no_prov <<
" splitting/merging.");
204 for (
auto part : semi_disconnected) {
205 if (!
part->production_vertex() ||
part->production_vertex()->id() == 0)
v->add_particle_out(std::move(
part));
207 for (
auto part : semi_disconnected) {
208 if (!
part->end_vertex())
v->add_particle_in(std::move(
part));
210 evt->add_vertex(std::move(
v));
219 for (
auto part: decay_loop_particles) {
221 auto vend =
part->end_vertex();
222 auto vprod =
part->production_vertex();
223 if (!vprod || !vend)
continue;
224 bool loop_in_decay =
true;
227 auto sisters = vend->particles_in();
228 for (
auto sister: sisters) {
229 if (vprod != sister->production_vertex()) loop_in_decay =
false;
231 if (!loop_in_decay)
continue;
234 auto daughters = vend->particles_out();
235 for (
auto p : daughters) vprod->add_particle_out(std::move(
p));
236 for (
auto sister : sisters) {
237 vprod->remove_particle_out(sister);
238 vend->remove_particle_in(sister);
239 evt->remove_particle(std::move(sister));
241 evt->remove_vertex(std::move(vend));
246 std::vector<HepMC::GenParticlePtr> toremove;
247 for (
auto ip:
evt->particles()) {
252 bool bad_particle =
false;
268 int abs_pdg_id = std::abs(
ip->pdg_id());
269 bool is_decayed_weak_boson = ( abs_pdg_id == 23 || abs_pdg_id == 24 || abs_pdg_id == 25 ) &&
ip->end_vertex();
277 if (bad_particle) toremove.push_back(
ip);
281 const int num_particles_orig =
evt->particles().size();
284 if (!toremove.empty()) {
285 ATH_MSG_DEBUG(
"Cleaning event record of " << toremove.size() <<
" bad particles");
286 for (
auto part: toremove)
evt->remove_particle(
part);
293 const std::vector <HepMC::GenParticlePtr> allParticles=
evt->particles();
294 for(
auto p : allParticles) {
296 if(
p->status() == 2 && !end_v) {
297 evt->remove_particle(std::move(
p));
306 const int num_particles_filt =
evt->particles().size();
308 if(num_particles_orig!=num_particles_filt) {
310 ATH_MSG_INFO(
"Particles filtered: " << num_particles_orig <<
" -> " << num_particles_filt);
315 if (
evt->weights().empty()) {
316 ATH_MSG_DEBUG(
"Adding a unit weight to empty event weight vector");
317 evt->weights().push_back(1);
321 if (
evt->signal_process_vertex() == NULL) {
322 const HepMC::FourVector nullpos;
323 for (HepMC::GenEvent::vertex_const_iterator iv =
evt->vertices_begin(); iv !=
evt->vertices_end(); ++iv) {
324 if ((*iv)->position() == nullpos) {
325 ATH_MSG_DEBUG(
"Setting representative event position vertex");
326 evt->set_signal_process_vertex(
const_cast<HepMC::GenVertex*
>(*iv));
334 std::vector<HepMC::GenParticlePtr> semi_disconnected, decay_loop_particles;
335 for (
auto ip : *
evt) {
338 bool particle_to_fix =
false;
339 int abspid = std::abs(
ip->pdg_id());
340 auto vProd =
ip->production_vertex();
341 auto vEnd =
ip->end_vertex();
343 if ( (!vProd || vProd->id() == 0) && vEnd &&
ip->status() != 4) {
344 particle_to_fix =
true;
345 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without production vertex! HepMC status = " <<
ip->status());
348 if (vProd && !vEnd &&
ip->status() != 1) {
349 particle_to_fix =
true;
350 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without decay vertex! HepMC status = " <<
ip->status());
352 if (particle_to_fix) semi_disconnected.push_back(
ip);
354 if (abspid == 43 || abspid == 44 || abspid == 30353 || abspid == 30343) {
355 decay_loop_particles.push_back(
ip);
361 if (semi_disconnected.size() == 4 ||semi_disconnected.size() == 3 || semi_disconnected.size() == 2) {
364 double vsum[4] = {0,0,0,0};
365 for (
auto part : semi_disconnected) {
366 if (!
part->production_vertex() ) {
368 vsum[0] +=
part->momentum().px();
369 vsum[1] +=
part->momentum().py();
370 vsum[2] +=
part->momentum().pz();
371 vsum[3] +=
part->momentum().e();
373 if (!
part->end_vertex()) {
375 vsum[0] -=
part->momentum().px();
376 vsum[1] -=
part->momentum().py();
377 vsum[2] -=
part->momentum().pz();
378 vsum[3] -=
part->momentum().e();
381 HepMC::FourVector
sum(vsum[0],vsum[1],vsum[2],vsum[3]);
382 ATH_MSG_INFO(
"Heuristics: found " << semi_disconnected.size() <<
" semi-disconnected particles. Momentum sum is " << vsum[0] <<
" " << vsum[1] <<
" " << vsum[2] <<
" " << vsum[3]);
384 if (no_endv && no_prov && ( no_endv + no_prov == semi_disconnected.size() )) {
385 if (std::abs(
sum.px()) < 1
e-2 && std::abs(
sum.py()) < 1
e-2 && std::abs(
sum.pz()) < 1
e-2 ) {
386 ATH_MSG_INFO(
"Try " << no_endv <<
"->" << no_prov <<
" splitting/merging.");
388 for (
auto part : semi_disconnected) {
389 if (!
part->production_vertex())
v->add_particle_out(
part);
391 for (
auto part : semi_disconnected) {
392 if (!
part->end_vertex())
v->add_particle_in(
part);
403 for (
auto part: decay_loop_particles) {
405 auto vend =
part->end_vertex();
406 auto vprod =
part->production_vertex();
407 if (!vend || !vprod)
continue;
408 bool loop_in_decay =
true;
410 std::vector<HepMC::GenParticlePtr> sisters;
412 for (
auto sister : sisters) {
413 if (vprod != sister->production_vertex()) loop_in_decay =
false;
415 if (!loop_in_decay)
continue;
417 std::vector<HepMC::GenParticlePtr> daughters;
419 for (
auto p : daughters) vprod->add_particle_out(
p);
420 for (
auto sister : sisters) {
421 vprod->remove_particle(sister);
422 vend->remove_particle(sister);
424 evt->remove_vertex(vend);
430 std::vector<HepMC::GenParticlePtr> toremove; toremove.reserve(10);
431 for (HepMC::GenEvent::particle_const_iterator
ip =
evt->particles_begin();
ip !=
evt->particles_end(); ++
ip) {
433 if (*
ip == NULL)
continue;
437 bool bad_particle =
false;
456 int abs_pdg_id = std::abs((*ip)->pdg_id());
457 bool is_decayed_weak_boson = ( abs_pdg_id == 23 || abs_pdg_id == 24 || abs_pdg_id == 25 ) && (*ip)->end_vertex();
466 if (bad_particle) toremove.push_back(*
ip);
470 const int num_particles_orig =
evt->particles_size();
471 int num_orphan_vtxs_orig = 0;
472 int num_noparent_vtxs_orig = 0;
473 int num_nochild_vtxs_orig = 0;
474 for (
auto v =
evt->vertices_begin();
v !=
evt->vertices_end(); ++
v) {
475 if ((*v)->particles_in_size()==0&&(*v)->particles_out_size()==0) num_orphan_vtxs_orig++;
476 if ((*v)->particles_in_size()==0) num_noparent_vtxs_orig++;
477 if ((*v)->particles_out_size()==0) num_nochild_vtxs_orig++;
481 if (!toremove.empty()) {
482 ATH_MSG_DEBUG(
"Cleaning event record of " << toremove.size() <<
" bad particles");
484 int signal_vertex_bc =
evt->signal_process_vertex() ?
evt->signal_process_vertex()->barcode() : 0;
487 if (
evt->barcode_to_vertex (signal_vertex_bc) ==
nullptr) {
488 evt->set_signal_process_vertex (
nullptr);
497 if (
p->status() == 2 && !end_v) {
498 delete p->production_vertex()->remove_particle(
p);
508 const int num_particles_filt =
evt->particles_size();
509 if(num_particles_orig!=num_particles_filt) {
510 int num_orphan_vtxs_filt = 0;
511 int num_noparent_vtxs_filt = 0;
512 int num_nochild_vtxs_filt = 0;
513 for (
auto v =
evt->vertices_begin();
v !=
evt->vertices_end(); ++
v) {
514 if ((*v)->particles_in_size()==0&&(*v)->particles_out_size()==0) num_orphan_vtxs_filt++;
515 if ((*v)->particles_in_size()==0) num_noparent_vtxs_filt++;
516 if ((*v)->particles_out_size()==0) num_nochild_vtxs_filt++;
520 ATH_MSG_INFO(
"Particles filtered: " << num_particles_orig <<
" -> " << num_particles_filt);
522 if (num_orphan_vtxs_filt != num_orphan_vtxs_orig)
523 ATH_MSG_WARNING(
"Change in orphaned vertices: " << num_orphan_vtxs_orig <<
" -> " << num_orphan_vtxs_filt);
524 if (num_noparent_vtxs_filt != num_noparent_vtxs_orig)
525 ATH_MSG_WARNING(
"Change in no-parent vertices: " << num_noparent_vtxs_orig <<
" -> " << num_noparent_vtxs_filt);
526 if (num_nochild_vtxs_filt != num_nochild_vtxs_orig)
527 ATH_MSG_WARNING(
"Change in no-parent vertices: " << num_nochild_vtxs_orig <<
" -> " << num_nochild_vtxs_filt);
531 return StatusCode::SUCCESS;
1.8.18