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.
99 HepMC::GenEvent*
evt =
const_cast<HepMC::GenEvent*
>(*ievt);
105 ATH_MSG_DEBUG(
"Please use MC::Loops::findLoops for this event to obtain all particles and vertices in the loops");
108 auto old_momentum =
evt->momentum_unit();
109 auto old_length =
evt->length_unit();
116 std::string old_momentum = (
evt->momentum_unit() == HepMC::Units::MomentumUnit::MEV ?
"MEV" :
"GEV" );
125 for (
auto ip: *
evt) {
129 if (newpid ==
m_pidmap.end())
continue;
130 ip->set_pdg_id(newpid->second);
137 if (
evt->weights().empty()) {
138 ATH_MSG_DEBUG(
"Adding a unit weight to empty event weight vector");
139 evt->weights().push_back(1);
143 const HepMC::FourVector nullpos;
144 for (
auto iv:
evt->vertices()) {
145 if (iv->position() == nullpos) {
146 ATH_MSG_DEBUG(
"Setting representative event position vertex");
154 std::vector<HepMC::GenParticlePtr> beams_t;
156 if (
p->status() == 4) beams_t.push_back(
p);
158 if (beams_t.size() > 2) {
159 ATH_MSG_INFO(
"Invalid number of beam particles " << beams_t.size() <<
". Will try to fix.");
160 std::vector<HepMC::GenParticlePtr> bparttoremove;
161 for (
const auto& bpart : beams_t) {
162 if (bpart->id() == 0 && bpart->production_vertex()) bparttoremove.push_back(bpart);
164 for (
auto bpart: bparttoremove) {
165 bpart->production_vertex()->remove_particle_out(bpart);
170 std::vector<HepMC::GenParticlePtr> semi_disconnected, decay_loop_particles;
171 for (
auto ip :
evt->particles()) {
174 bool particle_to_fix =
false;
175 int abspid = std::abs(
ip->pdg_id());
176 auto vProd =
ip->production_vertex();
177 auto vEnd =
ip->end_vertex();
179 if ( (!vProd || vProd->id() == 0) && vEnd &&
ip->status() != 4) {
180 particle_to_fix =
true;
181 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without production vertex! HepMC status = " <<
ip->status());
184 if (vProd && !vEnd &&
ip->status() != 1) {
185 particle_to_fix =
true;
186 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without decay vertex! HepMC status = " <<
ip->status());
188 if (particle_to_fix) semi_disconnected.push_back(
ip);
190 if (abspid == 43 || abspid == 44 || abspid == 30353 || abspid == 30343) {
191 decay_loop_particles.push_back(std::move(
ip));
203 if ( !
m_ignoreSemiDisconnected && (semi_disconnected.size() == 4 || semi_disconnected.size() == 3 || semi_disconnected.size() == 2)) {
206 HepMC::FourVector
sum(0,0,0,0);
208 for (
const auto&
part : semi_disconnected) {
209 if (!
part->production_vertex() || !
part->production_vertex()->id()) {
210 no_prov++;
sum +=
part->momentum();
212 if (!
part->end_vertex()) {
213 no_endv++;
sum -=
part->momentum();
218 ATH_MSG_INFO(
"Heuristics: found " << semi_disconnected.size() <<
" semi-disconnected particles. Momentum sum is " <<
sum <<
" Standalone vertices " <<
standalone.size());
219 bool standalonevertex = (
standalone.size() == 1 && (*
standalone.begin())->particles_in().size() + (*
standalone.begin())->particles_out().size() == semi_disconnected.size());
221 if (! standalonevertex && no_endv && no_prov && ( no_endv + no_prov == semi_disconnected.size() )) {
222 if (std::abs(
sum.px()) < 1
e-2 && std::abs(
sum.py()) < 1
e-2 && std::abs(
sum.pz()) < 1
e-2 ) {
223 ATH_MSG_INFO(
"Try " << no_endv <<
"->" << no_prov <<
" splitting/merging.");
225 for (
auto part : semi_disconnected) {
226 if (!
part->production_vertex() ||
part->production_vertex()->id() == 0)
v->add_particle_out(std::move(
part));
228 for (
auto part : semi_disconnected) {
229 if (!
part->end_vertex())
v->add_particle_in(std::move(
part));
231 evt->add_vertex(std::move(
v));
240 for (
auto part: decay_loop_particles) {
242 auto vend =
part->end_vertex();
243 auto vprod =
part->production_vertex();
244 if (!vprod || !vend)
continue;
245 bool loop_in_decay =
true;
248 auto sisters = vend->particles_in();
249 for (
auto sister: sisters) {
250 if (vprod != sister->production_vertex()) loop_in_decay =
false;
252 if (!loop_in_decay)
continue;
255 auto daughters = vend->particles_out();
256 for (
auto p : daughters) vprod->add_particle_out(std::move(
p));
257 for (
auto sister : sisters) {
258 vprod->remove_particle_out(sister);
259 vend->remove_particle_in(sister);
260 evt->remove_particle(std::move(sister));
262 evt->remove_vertex(std::move(vend));
267 std::vector<HepMC::GenParticlePtr> toremove;
268 for (
auto ip:
evt->particles()) {
273 bool bad_particle =
false;
289 int abs_pdg_id = std::abs(
ip->pdg_id());
290 bool is_decayed_weak_boson = ( abs_pdg_id == 23 || abs_pdg_id == 24 || abs_pdg_id == 25 ) &&
ip->end_vertex();
298 if (bad_particle) toremove.push_back(
ip);
302 const int num_particles_orig =
evt->particles().size();
305 if (!toremove.empty()) {
306 ATH_MSG_DEBUG(
"Cleaning event record of " << toremove.size() <<
" bad particles");
307 for (
auto part: toremove)
evt->remove_particle(
part);
314 const std::vector <HepMC::GenParticlePtr> allParticles=
evt->particles();
315 for(
auto p : allParticles) {
317 if(
p->status() == 2 && !end_v) {
318 evt->remove_particle(std::move(
p));
327 const int num_particles_filt =
evt->particles().size();
329 if(num_particles_orig!=num_particles_filt) {
331 ATH_MSG_INFO(
"Particles filtered: " << num_particles_orig <<
" -> " << num_particles_filt);
336 if (
evt->weights().empty()) {
337 ATH_MSG_DEBUG(
"Adding a unit weight to empty event weight vector");
338 evt->weights().push_back(1);
342 if (
evt->signal_process_vertex() == NULL) {
343 const HepMC::FourVector nullpos;
344 for (HepMC::GenEvent::vertex_const_iterator iv =
evt->vertices_begin(); iv !=
evt->vertices_end(); ++iv) {
345 if ((*iv)->position() == nullpos) {
346 ATH_MSG_DEBUG(
"Setting representative event position vertex");
347 evt->set_signal_process_vertex(
const_cast<HepMC::GenVertex*
>(*iv));
355 std::vector<HepMC::GenParticlePtr> semi_disconnected, decay_loop_particles;
356 for (
auto ip : *
evt) {
359 bool particle_to_fix =
false;
360 int abspid = std::abs(
ip->pdg_id());
361 auto vProd =
ip->production_vertex();
362 auto vEnd =
ip->end_vertex();
364 if ( (!vProd || vProd->id() == 0) && vEnd &&
ip->status() != 4) {
365 particle_to_fix =
true;
366 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without production vertex! HepMC status = " <<
ip->status());
369 if (vProd && !vEnd &&
ip->status() != 1) {
370 particle_to_fix =
true;
371 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without decay vertex! HepMC status = " <<
ip->status());
373 if (particle_to_fix) semi_disconnected.push_back(
ip);
375 if (abspid == 43 || abspid == 44 || abspid == 30353 || abspid == 30343) {
376 decay_loop_particles.push_back(
ip);
382 if (semi_disconnected.size() == 4 ||semi_disconnected.size() == 3 || semi_disconnected.size() == 2) {
385 double vsum[4] = {0,0,0,0};
386 for (
auto part : semi_disconnected) {
387 if (!
part->production_vertex() ) {
389 vsum[0] +=
part->momentum().px();
390 vsum[1] +=
part->momentum().py();
391 vsum[2] +=
part->momentum().pz();
392 vsum[3] +=
part->momentum().e();
394 if (!
part->end_vertex()) {
396 vsum[0] -=
part->momentum().px();
397 vsum[1] -=
part->momentum().py();
398 vsum[2] -=
part->momentum().pz();
399 vsum[3] -=
part->momentum().e();
402 HepMC::FourVector
sum(vsum[0],vsum[1],vsum[2],vsum[3]);
403 ATH_MSG_INFO(
"Heuristics: found " << semi_disconnected.size() <<
" semi-disconnected particles. Momentum sum is " << vsum[0] <<
" " << vsum[1] <<
" " << vsum[2] <<
" " << vsum[3]);
405 if (no_endv && no_prov && ( no_endv + no_prov == semi_disconnected.size() )) {
406 if (std::abs(
sum.px()) < 1
e-2 && std::abs(
sum.py()) < 1
e-2 && std::abs(
sum.pz()) < 1
e-2 ) {
407 ATH_MSG_INFO(
"Try " << no_endv <<
"->" << no_prov <<
" splitting/merging.");
409 for (
auto part : semi_disconnected) {
410 if (!
part->production_vertex())
v->add_particle_out(
part);
412 for (
auto part : semi_disconnected) {
413 if (!
part->end_vertex())
v->add_particle_in(
part);
424 for (
auto part: decay_loop_particles) {
426 auto vend =
part->end_vertex();
427 auto vprod =
part->production_vertex();
428 if (!vend || !vprod)
continue;
429 bool loop_in_decay =
true;
431 std::vector<HepMC::GenParticlePtr> sisters;
433 for (
auto sister : sisters) {
434 if (vprod != sister->production_vertex()) loop_in_decay =
false;
436 if (!loop_in_decay)
continue;
438 std::vector<HepMC::GenParticlePtr> daughters;
440 for (
auto p : daughters) vprod->add_particle_out(
p);
441 for (
auto sister : sisters) {
442 vprod->remove_particle(sister);
443 vend->remove_particle(sister);
445 evt->remove_vertex(vend);
451 std::vector<HepMC::GenParticlePtr> toremove; toremove.reserve(10);
452 for (HepMC::GenEvent::particle_const_iterator
ip =
evt->particles_begin();
ip !=
evt->particles_end(); ++
ip) {
454 if (*
ip == NULL)
continue;
458 bool bad_particle =
false;
477 int abs_pdg_id = std::abs((*ip)->pdg_id());
478 bool is_decayed_weak_boson = ( abs_pdg_id == 23 || abs_pdg_id == 24 || abs_pdg_id == 25 ) && (*ip)->end_vertex();
487 if (bad_particle) toremove.push_back(*
ip);
491 const int num_particles_orig =
evt->particles_size();
492 int num_orphan_vtxs_orig = 0;
493 int num_noparent_vtxs_orig = 0;
494 int num_nochild_vtxs_orig = 0;
495 for (
auto v =
evt->vertices_begin();
v !=
evt->vertices_end(); ++
v) {
496 if ((*v)->particles_in_size()==0&&(*v)->particles_out_size()==0) num_orphan_vtxs_orig++;
497 if ((*v)->particles_in_size()==0) num_noparent_vtxs_orig++;
498 if ((*v)->particles_out_size()==0) num_nochild_vtxs_orig++;
502 if (!toremove.empty()) {
503 ATH_MSG_DEBUG(
"Cleaning event record of " << toremove.size() <<
" bad particles");
505 int signal_vertex_bc =
evt->signal_process_vertex() ?
evt->signal_process_vertex()->barcode() : 0;
508 if (
evt->barcode_to_vertex (signal_vertex_bc) ==
nullptr) {
509 evt->set_signal_process_vertex (
nullptr);
518 if (
p->status() == 2 && !end_v) {
519 delete p->production_vertex()->remove_particle(
p);
529 const int num_particles_filt =
evt->particles_size();
530 if(num_particles_orig!=num_particles_filt) {
531 int num_orphan_vtxs_filt = 0;
532 int num_noparent_vtxs_filt = 0;
533 int num_nochild_vtxs_filt = 0;
534 for (
auto v =
evt->vertices_begin();
v !=
evt->vertices_end(); ++
v) {
535 if ((*v)->particles_in_size()==0&&(*v)->particles_out_size()==0) num_orphan_vtxs_filt++;
536 if ((*v)->particles_in_size()==0) num_noparent_vtxs_filt++;
537 if ((*v)->particles_out_size()==0) num_nochild_vtxs_filt++;
541 ATH_MSG_INFO(
"Particles filtered: " << num_particles_orig <<
" -> " << num_particles_filt);
543 if (num_orphan_vtxs_filt != num_orphan_vtxs_orig)
544 ATH_MSG_WARNING(
"Change in orphaned vertices: " << num_orphan_vtxs_orig <<
" -> " << num_orphan_vtxs_filt);
545 if (num_noparent_vtxs_filt != num_noparent_vtxs_orig)
546 ATH_MSG_WARNING(
"Change in no-parent vertices: " << num_noparent_vtxs_orig <<
" -> " << num_noparent_vtxs_filt);
547 if (num_nochild_vtxs_filt != num_nochild_vtxs_orig)
548 ATH_MSG_WARNING(
"Change in no-parent vertices: " << num_nochild_vtxs_orig <<
" -> " << num_nochild_vtxs_filt);
552 return StatusCode::SUCCESS;
1.8.18