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.
94 HepMC::GenEvent*
evt =
const_cast<HepMC::GenEvent*
>(*ievt);
100 ATH_MSG_DEBUG(
"Please use MC::Loops::findLoops for this event to obtain all particles and vertices in the loops");
103 for (
auto ip: *
evt) {
107 if (newpid ==
m_pidmap.end())
continue;
108 ip->set_pdg_id(newpid->second);
113 if (
evt->weights().empty()) {
114 ATH_MSG_DEBUG(
"Adding a unit weight to empty event weight vector");
115 evt->weights().push_back(1);
119 const HepMC::FourVector nullpos;
120 for (
auto iv:
evt->vertices()) {
121 if (iv->position() == nullpos) {
122 ATH_MSG_DEBUG(
"Setting representative event position vertex");
130 std::vector<HepMC::GenParticlePtr> beams_t;
132 if (
p->status() == 4) beams_t.push_back(
p);
134 if (beams_t.size() > 2) {
135 ATH_MSG_INFO(
"Invalid number of beam particles " << beams_t.size() <<
". Will try to fix.");
136 std::vector<HepMC::GenParticlePtr> bparttoremove;
137 for (
const auto& bpart : beams_t) {
138 if (bpart->id() == 0 && bpart->production_vertex()) bparttoremove.push_back(bpart);
140 for (
auto bpart: bparttoremove) {
141 bpart->production_vertex()->remove_particle_out(bpart);
146 std::vector<HepMC::GenParticlePtr> semi_disconnected, decay_loop_particles;
147 for (
auto ip :
evt->particles()) {
150 bool particle_to_fix =
false;
151 int abspid = std::abs(
ip->pdg_id());
152 auto vProd =
ip->production_vertex();
153 auto vEnd =
ip->end_vertex();
155 if ( (!vProd || vProd->id() == 0) && vEnd &&
ip->status() != 4) {
156 particle_to_fix =
true;
157 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without production vertex! HepMC status = " <<
ip->status());
160 if (vProd && !vEnd &&
ip->status() != 1) {
161 particle_to_fix =
true;
162 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without decay vertex! HepMC status = " <<
ip->status());
164 if (particle_to_fix) semi_disconnected.push_back(
ip);
166 if (abspid == 43 || abspid == 44 || abspid == 30353 || abspid == 30343) {
167 decay_loop_particles.push_back(
ip);
174 if ( semi_disconnected.size() == 4 || semi_disconnected.size() == 3 || semi_disconnected.size() == 2) {
177 HepMC::FourVector
sum(0,0,0,0);
178 for (
const auto&
part : semi_disconnected) {
179 if (!
part->production_vertex() || !
part->production_vertex()->id()) {
180 no_prov++;
sum +=
part->momentum();
182 if (!
part->end_vertex()) {
183 no_endv++;
sum -=
part->momentum();
186 ATH_MSG_INFO(
"Heuristics: found " << semi_disconnected.size() <<
" semi-disconnected particles. Momentum sum is " <<
sum);
188 if (no_endv && no_prov && ( no_endv + no_prov == semi_disconnected.size() )) {
189 if (std::abs(
sum.px()) < 1
e-2 && std::abs(
sum.py()) < 1
e-2 && std::abs(
sum.pz()) < 1
e-2 ) {
190 ATH_MSG_INFO(
"Try " << no_endv <<
"->" << no_prov <<
" splitting/merging.");
192 for (
auto part : semi_disconnected) {
193 if (!
part->production_vertex() ||
part->production_vertex()->id() == 0)
v->add_particle_out(
part);
195 for (
auto part : semi_disconnected) {
196 if (!
part->end_vertex())
v->add_particle_in(
part);
207 for (
auto part: decay_loop_particles) {
209 auto vend =
part->end_vertex();
210 auto vprod =
part->production_vertex();
211 if (!vprod || !vend)
continue;
212 bool loop_in_decay =
true;
215 auto sisters = vend->particles_in();
216 for (
auto sister: sisters) {
217 if (vprod != sister->production_vertex()) loop_in_decay =
false;
219 if (!loop_in_decay)
continue;
222 auto daughters = vend->particles_out();
223 for (
auto p : daughters) vprod->add_particle_out(
p);
224 for (
auto sister : sisters) {
225 vprod->remove_particle_out(sister);
226 vend->remove_particle_in(sister);
227 evt->remove_particle(sister);
229 evt->remove_vertex(vend);
234 std::vector<HepMC::GenParticlePtr> toremove;
235 for (
auto ip:
evt->particles()) {
240 bool bad_particle =
false;
256 int abs_pdg_id = std::abs(
ip->pdg_id());
257 bool is_decayed_weak_boson = ( abs_pdg_id == 23 || abs_pdg_id == 24 || abs_pdg_id == 25 ) &&
ip->end_vertex();
265 if (bad_particle) toremove.push_back(
ip);
269 if (toremove.empty())
continue;
270 ATH_MSG_DEBUG(
"Cleaning event record of " << toremove.size() <<
" bad particles");
272 const int num_particles_orig =
evt->particles().size();
273 for (
auto part: toremove)
evt->remove_particle(
part);
274 const int num_particles_filt =
evt->particles().size();
276 ATH_MSG_INFO(
"Particles filtered: " << num_particles_orig <<
" -> " << num_particles_filt);
280 if (
evt->weights().empty()) {
281 ATH_MSG_DEBUG(
"Adding a unit weight to empty event weight vector");
282 evt->weights().push_back(1);
286 if (
evt->signal_process_vertex() == NULL) {
287 const HepMC::FourVector nullpos;
288 for (HepMC::GenEvent::vertex_const_iterator iv =
evt->vertices_begin(); iv !=
evt->vertices_end(); ++iv) {
289 if ((*iv)->position() == nullpos) {
290 ATH_MSG_DEBUG(
"Setting representative event position vertex");
291 evt->set_signal_process_vertex(
const_cast<HepMC::GenVertex*
>(*iv));
299 std::vector<HepMC::GenParticlePtr> semi_disconnected, decay_loop_particles;
300 for (
auto ip : *
evt) {
303 bool particle_to_fix =
false;
304 int abspid = std::abs(
ip->pdg_id());
305 auto vProd =
ip->production_vertex();
306 auto vEnd =
ip->end_vertex();
308 if ( (!vProd || vProd->id() == 0) && vEnd &&
ip->status() != 4) {
309 particle_to_fix =
true;
310 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without production vertex! HepMC status = " <<
ip->status());
313 if (vProd && !vEnd &&
ip->status() != 1) {
314 particle_to_fix =
true;
315 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without decay vertex! HepMC status = " <<
ip->status());
317 if (particle_to_fix) semi_disconnected.push_back(
ip);
319 if (abspid == 43 || abspid == 44 || abspid == 30353 || abspid == 30343) {
320 decay_loop_particles.push_back(
ip);
326 if (semi_disconnected.size() == 4 ||semi_disconnected.size() == 3 || semi_disconnected.size() == 2) {
329 double vsum[4] = {0,0,0,0};
330 for (
auto part : semi_disconnected) {
331 if (!
part->production_vertex() ) {
333 vsum[0] +=
part->momentum().px();
334 vsum[1] +=
part->momentum().py();
335 vsum[2] +=
part->momentum().pz();
336 vsum[3] +=
part->momentum().e();
338 if (!
part->end_vertex()) {
340 vsum[0] -=
part->momentum().px();
341 vsum[1] -=
part->momentum().py();
342 vsum[2] -=
part->momentum().pz();
343 vsum[3] -=
part->momentum().e();
346 HepMC::FourVector
sum(vsum[0],vsum[1],vsum[2],vsum[3]);
347 ATH_MSG_INFO(
"Heuristics: found " << semi_disconnected.size() <<
" semi-disconnected particles. Momentum sum is " << vsum[0] <<
" " << vsum[1] <<
" " << vsum[2] <<
" " << vsum[3]);
349 if (no_endv && no_prov && ( no_endv + no_prov == semi_disconnected.size() )) {
350 if (std::abs(
sum.px()) < 1
e-2 && std::abs(
sum.py()) < 1
e-2 && std::abs(
sum.pz()) < 1
e-2 ) {
351 ATH_MSG_INFO(
"Try " << no_endv <<
"->" << no_prov <<
" splitting/merging.");
353 for (
auto part : semi_disconnected) {
354 if (!
part->production_vertex())
v->add_particle_out(
part);
356 for (
auto part : semi_disconnected) {
357 if (!
part->end_vertex())
v->add_particle_in(
part);
368 for (
auto part: decay_loop_particles) {
370 auto vend =
part->end_vertex();
371 auto vprod =
part->production_vertex();
372 if (!vend || !vprod)
continue;
373 bool loop_in_decay =
true;
375 std::vector<HepMC::GenParticlePtr> sisters;
377 for (
auto sister : sisters) {
378 if (vprod != sister->production_vertex()) loop_in_decay =
false;
380 if (!loop_in_decay)
continue;
382 std::vector<HepMC::GenParticlePtr> daughters;
384 for (
auto p : daughters) vprod->add_particle_out(
p);
385 for (
auto sister : sisters) {
386 vprod->remove_particle(sister);
387 vend->remove_particle(sister);
389 evt->remove_vertex(vend);
395 std::vector<HepMC::GenParticlePtr> toremove; toremove.reserve(10);
396 for (HepMC::GenEvent::particle_const_iterator
ip =
evt->particles_begin();
ip !=
evt->particles_end(); ++
ip) {
398 if (*
ip == NULL)
continue;
402 bool bad_particle =
false;
421 int abs_pdg_id = std::abs((*ip)->pdg_id());
422 bool is_decayed_weak_boson = ( abs_pdg_id == 23 || abs_pdg_id == 24 || abs_pdg_id == 25 ) && (*ip)->end_vertex();
431 if (bad_particle) toremove.push_back(*
ip);
435 if (toremove.empty())
continue;
436 ATH_MSG_DEBUG(
"Cleaning event record of " << toremove.size() <<
" bad particles");
439 const int num_particles_orig =
evt->particles_size();
440 int num_orphan_vtxs_orig = 0;
441 int num_noparent_vtxs_orig = 0;
442 int num_nochild_vtxs_orig = 0;
443 for (
auto v =
evt->vertices_begin();
v !=
evt->vertices_end(); ++
v) {
444 if ((*v)->particles_in_size()==0&&(*v)->particles_out_size()==0) num_orphan_vtxs_orig++;
445 if ((*v)->particles_in_size()==0) num_noparent_vtxs_orig++;
446 if ((*v)->particles_out_size()==0) num_nochild_vtxs_orig++;
449 int signal_vertex_bc =
evt->signal_process_vertex() ?
evt->signal_process_vertex()->barcode() : 0;
452 if (
evt->barcode_to_vertex (signal_vertex_bc) ==
nullptr) {
453 evt->set_signal_process_vertex (
nullptr);
457 const int num_particles_filt =
evt->particles_size();
458 int num_orphan_vtxs_filt = 0;
459 int num_noparent_vtxs_filt = 0;
460 int num_nochild_vtxs_filt = 0;
461 for (
auto v =
evt->vertices_begin();
v !=
evt->vertices_end(); ++
v) {
462 if ((*v)->particles_in_size()==0&&(*v)->particles_out_size()==0) num_orphan_vtxs_filt++;
463 if ((*v)->particles_in_size()==0) num_noparent_vtxs_filt++;
464 if ((*v)->particles_out_size()==0) num_nochild_vtxs_filt++;
468 ATH_MSG_INFO(
"Particles filtered: " << num_particles_orig <<
" -> " << num_particles_filt);
470 if (num_orphan_vtxs_filt != num_orphan_vtxs_orig)
471 ATH_MSG_WARNING(
"Change in orphaned vertices: " << num_orphan_vtxs_orig <<
" -> " << num_orphan_vtxs_filt);
472 if (num_noparent_vtxs_filt != num_noparent_vtxs_orig)
473 ATH_MSG_WARNING(
"Change in no-parent vertices: " << num_noparent_vtxs_orig <<
" -> " << num_noparent_vtxs_filt);
474 if (num_nochild_vtxs_filt != num_nochild_vtxs_orig)
475 ATH_MSG_WARNING(
"Change in no-parent vertices: " << num_nochild_vtxs_orig <<
" -> " << num_nochild_vtxs_filt);
479 return StatusCode::SUCCESS;
1.8.18