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.
96 HepMC::GenEvent*
evt =
const_cast<HepMC::GenEvent*
>(*ievt);
102 ATH_MSG_DEBUG(
"Please use MC::Loops::findLoops for this event to obtain all particles and vertices in the loops");
105 for (
auto ip: *
evt) {
109 if (newpid ==
m_pidmap.end())
continue;
110 ip->set_pdg_id(newpid->second);
115 if (
evt->weights().empty()) {
116 ATH_MSG_DEBUG(
"Adding a unit weight to empty event weight vector");
117 evt->weights().push_back(1);
121 const HepMC::FourVector nullpos;
122 for (
auto iv:
evt->vertices()) {
123 if (iv->position() == nullpos) {
124 ATH_MSG_DEBUG(
"Setting representative event position vertex");
132 std::vector<HepMC::GenParticlePtr> beams_t;
134 if (
p->status() == 4) beams_t.push_back(
p);
136 if (beams_t.size() > 2) {
137 ATH_MSG_INFO(
"Invalid number of beam particles " << beams_t.size() <<
". Will try to fix.");
138 std::vector<HepMC::GenParticlePtr> bparttoremove;
139 for (
const auto& bpart : beams_t) {
140 if (bpart->id() == 0 && bpart->production_vertex()) bparttoremove.push_back(bpart);
142 for (
auto bpart: bparttoremove) {
143 bpart->production_vertex()->remove_particle_out(bpart);
148 std::vector<HepMC::GenParticlePtr> semi_disconnected, decay_loop_particles;
149 for (
auto ip :
evt->particles()) {
152 bool particle_to_fix =
false;
153 int abspid = std::abs(
ip->pdg_id());
154 auto vProd =
ip->production_vertex();
155 auto vEnd =
ip->end_vertex();
157 if ( (!vProd || vProd->id() == 0) && vEnd &&
ip->status() != 4) {
158 particle_to_fix =
true;
159 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without production vertex! HepMC status = " <<
ip->status());
162 if (vProd && !vEnd &&
ip->status() != 1) {
163 particle_to_fix =
true;
164 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without decay vertex! HepMC status = " <<
ip->status());
166 if (particle_to_fix) semi_disconnected.push_back(
ip);
168 if (abspid == 43 || abspid == 44 || abspid == 30353 || abspid == 30343) {
169 decay_loop_particles.push_back(
ip);
176 if ( semi_disconnected.size() == 4 || semi_disconnected.size() == 3 || semi_disconnected.size() == 2) {
179 HepMC::FourVector
sum(0,0,0,0);
180 for (
const auto&
part : semi_disconnected) {
181 if (!
part->production_vertex() || !
part->production_vertex()->id()) {
182 no_prov++;
sum +=
part->momentum();
184 if (!
part->end_vertex()) {
185 no_endv++;
sum -=
part->momentum();
188 ATH_MSG_INFO(
"Heuristics: found " << semi_disconnected.size() <<
" semi-disconnected particles. Momentum sum is " <<
sum);
190 if (no_endv && no_prov && ( no_endv + no_prov == semi_disconnected.size() )) {
191 if (std::abs(
sum.px()) < 1
e-2 && std::abs(
sum.py()) < 1
e-2 && std::abs(
sum.pz()) < 1
e-2 ) {
192 ATH_MSG_INFO(
"Try " << no_endv <<
"->" << no_prov <<
" splitting/merging.");
194 for (
auto part : semi_disconnected) {
195 if (!
part->production_vertex() ||
part->production_vertex()->id() == 0)
v->add_particle_out(
part);
197 for (
auto part : semi_disconnected) {
198 if (!
part->end_vertex())
v->add_particle_in(
part);
209 for (
auto part: decay_loop_particles) {
211 auto vend =
part->end_vertex();
212 auto vprod =
part->production_vertex();
213 if (!vprod || !vend)
continue;
214 bool loop_in_decay =
true;
217 auto sisters = vend->particles_in();
218 for (
auto sister: sisters) {
219 if (vprod != sister->production_vertex()) loop_in_decay =
false;
221 if (!loop_in_decay)
continue;
224 auto daughters = vend->particles_out();
225 for (
auto p : daughters) vprod->add_particle_out(
p);
226 for (
auto sister : sisters) {
227 vprod->remove_particle_out(sister);
228 vend->remove_particle_in(sister);
229 evt->remove_particle(sister);
231 evt->remove_vertex(vend);
236 std::vector<HepMC::GenParticlePtr> toremove;
237 for (
auto ip:
evt->particles()) {
242 bool bad_particle =
false;
258 int abs_pdg_id = std::abs(
ip->pdg_id());
259 bool is_decayed_weak_boson = ( abs_pdg_id == 23 || abs_pdg_id == 24 || abs_pdg_id == 25 ) &&
ip->end_vertex();
267 if (bad_particle) toremove.push_back(
ip);
271 const int num_particles_orig =
evt->particles().size();
274 if (!toremove.empty()) {
275 ATH_MSG_DEBUG(
"Cleaning event record of " << toremove.size() <<
" bad particles");
276 for (
auto part: toremove)
evt->remove_particle(
part);
283 const std::vector <HepMC::GenParticlePtr> allParticles=
evt->particles();
284 for(
auto p : allParticles) {
286 if(
p->status() == 2 && !end_v) {
287 evt->remove_particle(
p);
296 const int num_particles_filt =
evt->particles().size();
298 if(num_particles_orig!=num_particles_filt) {
300 ATH_MSG_INFO(
"Particles filtered: " << num_particles_orig <<
" -> " << num_particles_filt);
305 if (
evt->weights().empty()) {
306 ATH_MSG_DEBUG(
"Adding a unit weight to empty event weight vector");
307 evt->weights().push_back(1);
311 if (
evt->signal_process_vertex() == NULL) {
312 const HepMC::FourVector nullpos;
313 for (HepMC::GenEvent::vertex_const_iterator iv =
evt->vertices_begin(); iv !=
evt->vertices_end(); ++iv) {
314 if ((*iv)->position() == nullpos) {
315 ATH_MSG_DEBUG(
"Setting representative event position vertex");
316 evt->set_signal_process_vertex(
const_cast<HepMC::GenVertex*
>(*iv));
324 std::vector<HepMC::GenParticlePtr> semi_disconnected, decay_loop_particles;
325 for (
auto ip : *
evt) {
328 bool particle_to_fix =
false;
329 int abspid = std::abs(
ip->pdg_id());
330 auto vProd =
ip->production_vertex();
331 auto vEnd =
ip->end_vertex();
333 if ( (!vProd || vProd->id() == 0) && vEnd &&
ip->status() != 4) {
334 particle_to_fix =
true;
335 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without production vertex! HepMC status = " <<
ip->status());
338 if (vProd && !vEnd &&
ip->status() != 1) {
339 particle_to_fix =
true;
340 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without decay vertex! HepMC status = " <<
ip->status());
342 if (particle_to_fix) semi_disconnected.push_back(
ip);
344 if (abspid == 43 || abspid == 44 || abspid == 30353 || abspid == 30343) {
345 decay_loop_particles.push_back(
ip);
351 if (semi_disconnected.size() == 4 ||semi_disconnected.size() == 3 || semi_disconnected.size() == 2) {
354 double vsum[4] = {0,0,0,0};
355 for (
auto part : semi_disconnected) {
356 if (!
part->production_vertex() ) {
358 vsum[0] +=
part->momentum().px();
359 vsum[1] +=
part->momentum().py();
360 vsum[2] +=
part->momentum().pz();
361 vsum[3] +=
part->momentum().e();
363 if (!
part->end_vertex()) {
365 vsum[0] -=
part->momentum().px();
366 vsum[1] -=
part->momentum().py();
367 vsum[2] -=
part->momentum().pz();
368 vsum[3] -=
part->momentum().e();
371 HepMC::FourVector
sum(vsum[0],vsum[1],vsum[2],vsum[3]);
372 ATH_MSG_INFO(
"Heuristics: found " << semi_disconnected.size() <<
" semi-disconnected particles. Momentum sum is " << vsum[0] <<
" " << vsum[1] <<
" " << vsum[2] <<
" " << vsum[3]);
374 if (no_endv && no_prov && ( no_endv + no_prov == semi_disconnected.size() )) {
375 if (std::abs(
sum.px()) < 1
e-2 && std::abs(
sum.py()) < 1
e-2 && std::abs(
sum.pz()) < 1
e-2 ) {
376 ATH_MSG_INFO(
"Try " << no_endv <<
"->" << no_prov <<
" splitting/merging.");
378 for (
auto part : semi_disconnected) {
379 if (!
part->production_vertex())
v->add_particle_out(
part);
381 for (
auto part : semi_disconnected) {
382 if (!
part->end_vertex())
v->add_particle_in(
part);
393 for (
auto part: decay_loop_particles) {
395 auto vend =
part->end_vertex();
396 auto vprod =
part->production_vertex();
397 if (!vend || !vprod)
continue;
398 bool loop_in_decay =
true;
400 std::vector<HepMC::GenParticlePtr> sisters;
402 for (
auto sister : sisters) {
403 if (vprod != sister->production_vertex()) loop_in_decay =
false;
405 if (!loop_in_decay)
continue;
407 std::vector<HepMC::GenParticlePtr> daughters;
409 for (
auto p : daughters) vprod->add_particle_out(
p);
410 for (
auto sister : sisters) {
411 vprod->remove_particle(sister);
412 vend->remove_particle(sister);
414 evt->remove_vertex(vend);
420 std::vector<HepMC::GenParticlePtr> toremove; toremove.reserve(10);
421 for (HepMC::GenEvent::particle_const_iterator
ip =
evt->particles_begin();
ip !=
evt->particles_end(); ++
ip) {
423 if (*
ip == NULL)
continue;
427 bool bad_particle =
false;
446 int abs_pdg_id = std::abs((*ip)->pdg_id());
447 bool is_decayed_weak_boson = ( abs_pdg_id == 23 || abs_pdg_id == 24 || abs_pdg_id == 25 ) && (*ip)->end_vertex();
456 if (bad_particle) toremove.push_back(*
ip);
460 const int num_particles_orig =
evt->particles_size();
461 int num_orphan_vtxs_orig = 0;
462 int num_noparent_vtxs_orig = 0;
463 int num_nochild_vtxs_orig = 0;
464 for (
auto v =
evt->vertices_begin();
v !=
evt->vertices_end(); ++
v) {
465 if ((*v)->particles_in_size()==0&&(*v)->particles_out_size()==0) num_orphan_vtxs_orig++;
466 if ((*v)->particles_in_size()==0) num_noparent_vtxs_orig++;
467 if ((*v)->particles_out_size()==0) num_nochild_vtxs_orig++;
471 if (!toremove.empty()) {
472 ATH_MSG_DEBUG(
"Cleaning event record of " << toremove.size() <<
" bad particles");
474 int signal_vertex_bc =
evt->signal_process_vertex() ?
evt->signal_process_vertex()->barcode() : 0;
477 if (
evt->barcode_to_vertex (signal_vertex_bc) ==
nullptr) {
478 evt->set_signal_process_vertex (
nullptr);
487 if (
p->status() == 2 && !end_v) {
488 delete p->production_vertex()->remove_particle(
p);
498 const int num_particles_filt =
evt->particles_size();
499 if(num_particles_orig!=num_particles_filt) {
500 int num_orphan_vtxs_filt = 0;
501 int num_noparent_vtxs_filt = 0;
502 int num_nochild_vtxs_filt = 0;
503 for (
auto v =
evt->vertices_begin();
v !=
evt->vertices_end(); ++
v) {
504 if ((*v)->particles_in_size()==0&&(*v)->particles_out_size()==0) num_orphan_vtxs_filt++;
505 if ((*v)->particles_in_size()==0) num_noparent_vtxs_filt++;
506 if ((*v)->particles_out_size()==0) num_nochild_vtxs_filt++;
510 ATH_MSG_INFO(
"Particles filtered: " << num_particles_orig <<
" -> " << num_particles_filt);
512 if (num_orphan_vtxs_filt != num_orphan_vtxs_orig)
513 ATH_MSG_WARNING(
"Change in orphaned vertices: " << num_orphan_vtxs_orig <<
" -> " << num_orphan_vtxs_filt);
514 if (num_noparent_vtxs_filt != num_noparent_vtxs_orig)
515 ATH_MSG_WARNING(
"Change in no-parent vertices: " << num_noparent_vtxs_orig <<
" -> " << num_noparent_vtxs_filt);
516 if (num_nochild_vtxs_filt != num_nochild_vtxs_orig)
517 ATH_MSG_WARNING(
"Change in no-parent vertices: " << num_nochild_vtxs_orig <<
" -> " << num_nochild_vtxs_filt);
521 return StatusCode::SUCCESS;
1.8.18