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 {
99
100 HepMC::GenEvent*
evt =
const_cast<HepMC::GenEvent*
>(*ievt);
102
103 if (!
m_looper.loop_particles().empty() || !
m_looper.loop_vertices().empty()) {
106 ATH_MSG_DEBUG(
"Please use MC::Loops::findLoops for this event to obtain all particles and vertices in the loops");
107 }
108#ifdef HEPMC3
109 auto old_momentum =
evt->momentum_unit();
110 auto old_length =
evt->length_unit();
116#else
117 std::string old_momentum = (
evt->momentum_unit() == HepMC::Units::MomentumUnit::MEV ?
"MEV" :
"GEV" );
118 std::string old_length = (
evt->length_unit() == HepMC::Units::LengthUnit::MM ?
"MM" :
"CM" );
124#endif
125
127 for (auto ip: *evt) {
128
129 if (!ip) continue;
131 if (newpid ==
m_pidmap.end())
continue;
132 ip->set_pdg_id(newpid->second);
133
136 }
137 }
138#ifdef HEPMC3
139
140 if (
evt->weights().empty()) {
141 ATH_MSG_DEBUG(
"Adding a unit weight to empty event weight vector");
142 evt->weights().push_back(1);
143 }
144
146 const HepMC::FourVector nullpos;
147 for (
auto iv:
evt->vertices()) {
148 if (iv->position() == nullpos) {
149 ATH_MSG_DEBUG(
"Setting representative event position vertex");
151 break;
152 }
153 }
154 }
155
156
157 std::vector<HepMC::GenParticlePtr> beams_t;
159 if (
p->status() == 4) beams_t.push_back(p);
160 }
161 if (beams_t.size() > 2) {
162 ATH_MSG_INFO(
"Invalid number of beam particles " << beams_t.size() <<
". Will try to fix.");
163 std::vector<HepMC::GenParticlePtr> bparttoremove;
164 for (const auto& bpart : beams_t) {
165 if (bpart->id() == 0 && bpart->production_vertex()) bparttoremove.push_back(bpart);
166 }
167 for (auto bpart: bparttoremove) {
168 bpart->production_vertex()->remove_particle_out(bpart);
169 }
170 }
171
172
173
174
176 double units_problem = -1.;
177
178 if (beams_t.size() > 0 ){
180 if (
p->momentum().pz() > 1e9) units_problem =
p->momentum().pz();
181 }
182
183 } else {
185 if (p &&
p->momentum().pz() > 1e9) units_problem =
p->momentum().pz();
186 }
187 }
188 if (units_problem>0){
189 ATH_MSG_INFO(
"Apparent units problem; beam particles have z-momentum " << units_problem <<
" in MeV. Will divide by 1000.");
191 }
192 }
193
194
195 std::vector<HepMC::GenParticlePtr> semi_disconnected, decay_loop_particles;
196 for (
auto ip :
evt->particles()) {
197
198 if (!ip) continue;
199 bool particle_to_fix = false;
200 int abspid = std::abs(
ip->pdg_id());
201 auto vProd =
ip->production_vertex();
202 auto vEnd =
ip->end_vertex();
204 if ( (!vProd || vProd->id() == 0) && vEnd &&
ip->status() != 4) {
205 particle_to_fix = true;
206 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without production vertex! HepMC status = " <<
ip->status());
207 }
209 if (vProd && !vEnd &&
ip->status() != 1) {
210 particle_to_fix = true;
211 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without decay vertex! HepMC status = " <<
ip->status());
212 }
213 if (particle_to_fix) semi_disconnected.push_back(ip);
214
215 if (abspid == 43 || abspid == 44 || abspid == 30353 || abspid == 30343) {
216 decay_loop_particles.push_back(std::move(ip));
217 }
218 }
219
223
228 if ( !
m_ignoreSemiDisconnected && (semi_disconnected.size() == 4 || semi_disconnected.size() == 3 || semi_disconnected.size() == 2)) {
229 size_t no_endv = 0;
230 size_t no_prov = 0;
231 HepMC::FourVector
sum(0,0,0,0);
233 for (const auto& part : semi_disconnected) {
234 if (!
part->production_vertex() || !
part->production_vertex()->id()) {
235 no_prov++;
sum +=
part->momentum();
236 }
237 if (!
part->end_vertex()) {
238 no_endv++;
sum -=
part->momentum();
239 }
242 }
243 ATH_MSG_INFO(
"Heuristics: found " << semi_disconnected.size() <<
" semi-disconnected particles. Momentum sum is " << sum <<
" Standalone vertices " <<
standalone.size());
244 bool standalonevertex = (
standalone.size() == 1 && (*
standalone.begin())->particles_in().size() + (*
standalone.begin())->particles_out().size() == semi_disconnected.size());
246 if (! standalonevertex && no_endv && no_prov && ( no_endv + no_prov == semi_disconnected.size() )) {
247 if (std::abs(
sum.px()) < 1e-2 && std::abs(
sum.py()) < 1e-2 && std::abs(
sum.pz()) < 1e-2 ) {
248 ATH_MSG_INFO(
"Try " << no_endv <<
"->" << no_prov <<
" splitting/merging.");
250 for (auto part : semi_disconnected) {
251 if (!
part->production_vertex() ||
part->production_vertex()->id() == 0)
v->add_particle_out(std::move(part));
252 }
253 for (auto part : semi_disconnected) {
254 if (!
part->end_vertex())
v->add_particle_in(std::move(part));
255 }
256 evt->add_vertex(std::move(v));
257 }
258 }
259 }
260
265 for (auto part: decay_loop_particles) {
267 auto vend =
part->end_vertex();
268 auto vprod =
part->production_vertex();
269 if (!vprod || !vend) continue;
270 bool loop_in_decay = true;
273 auto sisters = vend->particles_in();
274 for (auto sister: sisters) {
275 if (vprod != sister->production_vertex()) loop_in_decay = false;
276 }
277 if (!loop_in_decay) continue;
278
280 auto daughters = vend->particles_out();
281 for (auto p : daughters) vprod->add_particle_out(std::move(p));
282 for (auto sister : sisters) {
283 vprod->remove_particle_out(sister);
284 vend->remove_particle_in(sister);
285 evt->remove_particle(std::move(sister));
286 }
287 evt->remove_vertex(std::move(vend));
288
289 }
290
291
292 std::vector<HepMC::GenParticlePtr> toremove;
293 for (
auto ip:
evt->particles()) {
294
295 if (!ip) continue;
297
298 bool bad_particle = false;
299
301 bad_particle = true;
305 }
306
308 bad_particle = true;
312 }
313
314 int abs_pdg_id = std::abs(
ip->pdg_id());
315 bool is_decayed_weak_boson = ( abs_pdg_id == 23 || abs_pdg_id == 24 || abs_pdg_id == 25 ) &&
ip->end_vertex();
317 bad_particle = true;
321 }
322
323 if (bad_particle) toremove.push_back(std::move(ip));
324 }
325
326
327 const int num_particles_orig =
evt->particles().size();
328
329
330 if (!toremove.empty()) {
331 ATH_MSG_DEBUG(
"Cleaning event record of " << toremove.size() <<
" bad particles");
332 for (
auto part: toremove)
evt->remove_particle(std::move(part));
333 }
334
336 int purged=0;
337 do {
338 purged=0;
339 const std::vector <HepMC::GenParticlePtr> allParticles=
evt->particles();
340 for(auto p : allParticles) {
342 if(
p->status() == 2 && !end_v) {
343 evt->remove_particle(std::move(p));
344 ++purged;
346 }
347 }
348 }
349 while (purged>0);
350 }
351
352 const int num_particles_filt =
evt->particles().size();
353
354 if(num_particles_orig!=num_particles_filt) {
355
356 ATH_MSG_INFO(
"Particles filtered: " << num_particles_orig <<
" -> " << num_particles_filt);
357 }
358 #else
359
360
361 if (
evt->weights().empty()) {
362 ATH_MSG_DEBUG(
"Adding a unit weight to empty event weight vector");
363 evt->weights().push_back(1);
364 }
365
366
367 if (
evt->signal_process_vertex() == NULL) {
368 const HepMC::FourVector nullpos;
369 for (HepMC::GenEvent::vertex_const_iterator iv =
evt->vertices_begin(); iv !=
evt->vertices_end(); ++iv) {
370 if ((*iv)->position() == nullpos) {
371 ATH_MSG_DEBUG(
"Setting representative event position vertex");
372 evt->set_signal_process_vertex(
const_cast<HepMC::GenVertex*
>(*iv));
373 break;
374 }
375 }
376 }
377
378
379
380 std::vector<HepMC::GenParticlePtr> semi_disconnected, decay_loop_particles;
381 for (auto ip : *evt) {
382
383 if (!ip) continue;
384 bool particle_to_fix = false;
385 int abspid = std::abs(
ip->pdg_id());
386 auto vProd =
ip->production_vertex();
387 auto vEnd =
ip->end_vertex();
389 if ( (!vProd || vProd->id() == 0) && vEnd &&
ip->status() != 4) {
390 particle_to_fix = true;
391 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without production vertex! HepMC status = " <<
ip->status());
392 }
394 if (vProd && !vEnd &&
ip->status() != 1) {
395 particle_to_fix = true;
396 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without decay vertex! HepMC status = " <<
ip->status());
397 }
398 if (particle_to_fix) semi_disconnected.push_back(ip);
399
400 if (abspid == 43 || abspid == 44 || abspid == 30353 || abspid == 30343) {
401 decay_loop_particles.push_back(ip);
402 }
403 }
404
407 if (semi_disconnected.size() == 4 ||semi_disconnected.size() == 3 || semi_disconnected.size() == 2) {
408 size_t no_endv = 0;
409 size_t no_prov = 0;
410 double vsum[4] = {0,0,0,0};
411 for (auto part : semi_disconnected) {
412 if (!
part->production_vertex() ) {
413 no_prov++;
414 vsum[0] +=
part->momentum().px();
415 vsum[1] +=
part->momentum().py();
416 vsum[2] +=
part->momentum().pz();
417 vsum[3] +=
part->momentum().e();
418 }
419 if (!
part->end_vertex()) {
420 no_endv++;
421 vsum[0] -=
part->momentum().px();
422 vsum[1] -=
part->momentum().py();
423 vsum[2] -=
part->momentum().pz();
424 vsum[3] -=
part->momentum().e();
425 }
426 }
427 HepMC::FourVector
sum(vsum[0],vsum[1],vsum[2],vsum[3]);
428 ATH_MSG_INFO(
"Heuristics: found " << semi_disconnected.size() <<
" semi-disconnected particles. Momentum sum is " << vsum[0] <<
" " << vsum[1] <<
" " << vsum[2] <<
" " << vsum[3]);
430 if (no_endv && no_prov && ( no_endv + no_prov == semi_disconnected.size() )) {
431 if (std::abs(
sum.px()) < 1e-2 && std::abs(
sum.py()) < 1e-2 && std::abs(
sum.pz()) < 1e-2 ) {
432 ATH_MSG_INFO(
"Try " << no_endv <<
"->" << no_prov <<
" splitting/merging.");
434 for (auto part : semi_disconnected) {
435 if (!
part->production_vertex())
v->add_particle_out(part);
436 }
437 for (auto part : semi_disconnected) {
438 if (!
part->end_vertex())
v->add_particle_in(part);
439 }
441 }
442 }
443 }
444
449 for (auto part: decay_loop_particles) {
451 auto vend =
part->end_vertex();
452 auto vprod =
part->production_vertex();
453 if (!vend || !vprod) continue;
454 bool loop_in_decay = true;
456 std::vector<HepMC::GenParticlePtr> sisters;
457 for (auto p = vend->particles_begin(HepMC::parents); p!= vend->particles_end(HepMC::parents); ++p) sisters.push_back(*p);
458 for (auto sister : sisters) {
459 if (vprod != sister->production_vertex()) loop_in_decay = false;
460 }
461 if (!loop_in_decay) continue;
462
463 std::vector<HepMC::GenParticlePtr> daughters;
464 for (auto p = vend->particles_begin(HepMC::children); p!= vend->particles_end(HepMC::children); ++p) daughters.push_back(*p);
465 for (auto p : daughters) vprod->add_particle_out(p);
466 for (auto sister : sisters) {
467 vprod->remove_particle(sister);
468 vend->remove_particle(sister);
469 }
470 evt->remove_vertex(vend);
471
472 }
473
474
475
476 std::vector<HepMC::GenParticlePtr> toremove; toremove.reserve(10);
477 for (HepMC::GenEvent::particle_const_iterator ip =
evt->particles_begin(); ip !=
evt->particles_end(); ++ip) {
478
479 if (*ip == NULL) continue;
481
482
483 bool bad_particle = false;
484
485
487 bad_particle = true;
490 if (
msgLvl( MSG::DEBUG ) ) (*ip)->print();
491 }
492
493
495 bad_particle = true;
498 if (
msgLvl( MSG::DEBUG ) ) (*ip)->print();
499 }
500
501
502 int abs_pdg_id = std::abs((*ip)->pdg_id());
503 bool is_decayed_weak_boson = ( abs_pdg_id == 23 || abs_pdg_id == 24 || abs_pdg_id == 25 ) && (*ip)->end_vertex();
505 bad_particle = true;
508 if (
msgLvl( MSG::DEBUG ) ) (*ip)->print();
509 }
510
511
512 if (bad_particle) toremove.push_back(*ip);
513 }
514
515
516 const int num_particles_orig =
evt->particles_size();
517 int num_orphan_vtxs_orig = 0;
518 int num_noparent_vtxs_orig = 0;
519 int num_nochild_vtxs_orig = 0;
520 for (
auto v =
evt->vertices_begin(); v !=
evt->vertices_end(); ++v) {
521 if ((*v)->particles_in_size()==0&&(*v)->particles_out_size()==0) num_orphan_vtxs_orig++;
522 if ((*v)->particles_in_size()==0) num_noparent_vtxs_orig++;
523 if ((*v)->particles_out_size()==0) num_nochild_vtxs_orig++;
524 }
525
526
527 if (!toremove.empty()) {
528 ATH_MSG_DEBUG(
"Cleaning event record of " << toremove.size() <<
" bad particles");
529
530 int signal_vertex_bc =
evt->signal_process_vertex() ?
evt->signal_process_vertex()->barcode() : 0;
531
533 if (
evt->barcode_to_vertex (signal_vertex_bc) ==
nullptr) {
534 evt->set_signal_process_vertex (
nullptr);
535 }
536 }
537
539 int purged=0;
540 do {
541 for (HepMC::GenParticle* p : *evt) {
543 if (
p->status() == 2 && !end_v) {
544 delete p->production_vertex()->remove_particle(p);
545 ++purged;
547 }
548 }
549 }
550 while (purged>0);
551 }
552
553
554 const int num_particles_filt =
evt->particles_size();
555 if(num_particles_orig!=num_particles_filt) {
556 int num_orphan_vtxs_filt = 0;
557 int num_noparent_vtxs_filt = 0;
558 int num_nochild_vtxs_filt = 0;
559 for (
auto v =
evt->vertices_begin(); v !=
evt->vertices_end(); ++v) {
560 if ((*v)->particles_in_size()==0&&(*v)->particles_out_size()==0) num_orphan_vtxs_filt++;
561 if ((*v)->particles_in_size()==0) num_noparent_vtxs_filt++;
562 if ((*v)->particles_out_size()==0) num_nochild_vtxs_filt++;
563 }
564
565
566 ATH_MSG_INFO(
"Particles filtered: " << num_particles_orig <<
" -> " << num_particles_filt);
567
568 if (num_orphan_vtxs_filt != num_orphan_vtxs_orig)
569 ATH_MSG_WARNING(
"Change in orphaned vertices: " << num_orphan_vtxs_orig <<
" -> " << num_orphan_vtxs_filt);
570 if (num_noparent_vtxs_filt != num_noparent_vtxs_orig)
571 ATH_MSG_WARNING(
"Change in no-parent vertices: " << num_noparent_vtxs_orig <<
" -> " << num_noparent_vtxs_filt);
572 if (num_nochild_vtxs_filt != num_nochild_vtxs_orig)
573 ATH_MSG_WARNING(
"Change in no-parent vertices: " << num_nochild_vtxs_orig <<
" -> " << num_nochild_vtxs_filt);
574 }
575#endif
576
577 }
578 return StatusCode::SUCCESS;
579}
#define ATH_MSG_WARNING(x)
static void reduce(HepMC::GenEvent *ge, HepMC::GenParticle *gp)
Remove an unwanted particle from the event, collapsing the graph structure consistently.
bool msgLvl(const MSG::Level lvl) const
DataModel_detail::const_iterator< DataVector > const_iterator
bool isPID0(const HepMC::ConstGenParticlePtr &p) const
bool isNonTransportableInDecayChain(const HepMC::ConstGenParticlePtr &p) const
std::map< int, int > m_replacedpid_counts
map of counters of replacements.
MC::Loops< HepMC::GenEvent, HepMC::ConstGenParticlePtr, HepMC::ConstGenVertexPtr > m_looper
member to detect loops
bool isSimpleLoop(const HepMC::ConstGenParticlePtr &p) const
void MeVToGeV(HepMC::GenEvent *evt)
Scale event energies/momenta by x 1/1000.
void line(std::ostream &os, const GenEvent &e)
void set_signal_process_vertex(GenEvent *e, T v)
GenVertexPtr newGenVertexPtr(const HepMC::FourVector &pos=HepMC::FourVector(0.0, 0.0, 0.0, 0.0), const int i=0)
GenParticle * GenParticlePtr
GenVertex * signal_process_vertex(const GenEvent *e)
const HepMC::GenVertex * ConstGenVertexPtr