Execute method.
270 {
271
272
273 bool filter_pass = true;
274
275
277 double totalPx = 0;
278 double totalPy = 0;
279 double totalPz = 0;
280 double totalE = 0;
281 double nonG4_energy = 0;
282 std::vector<HepMC::ConstGenParticlePtr> negEnPart;
283 std::vector<HepMC::ConstGenParticlePtr> tachyons;
284 std::vector<HepMC::ConstGenParticlePtr> unstNoEnd;
285 std::vector<HepMC::ConstGenParticlePtr> unDecPi0;
286 std::vector<HepMC::ConstGenParticlePtr> undisplaceds;
287
289 if (!
m_looper.loop_particles().empty() || !
m_looper.loop_vertices().empty()) {
292 ATH_MSG_DEBUG(
"Please use MC::Loops::findLoops for this event to obtain all particles and vertices in the loops");
294 }
295
296 const auto xsec =
evt->cross_section();
297 if (!xsec) {
298 ATH_MSG_WARNING(
"WATCH OUT: event is missing the generator cross-section!");
301 ATH_MSG_WARNING(
"-> Adding a dummy cross-section for debugging purposes.");
302
303 std::shared_ptr<HepMC3::GenCrossSection> dummy_xsec = std::make_shared<HepMC3::GenCrossSection>();
304 dummy_xsec->set_cross_section(1.0,0.0);
306 evt_nonconst->set_cross_section(std::move(dummy_xsec));
307 }
308 else {
310 }
311 }
312
313
314 std::vector<std::shared_ptr<const HepMC3::GenParticle>> beams_t;
315 for (
auto p :
evt->beams()) {
if (
p->status() == 4) beams_t.push_back(std::move(p)); }
316 std::pair<std::shared_ptr<const HepMC3::GenParticle>,std::shared_ptr<const HepMC3::GenParticle>> beams;
317 if (beams_t.size() == 2) {
318 beams.first=beams_t.at(0);
319 beams.second=beams_t.at(1);
320 } else {
321 ATH_MSG_WARNING(
"Invalid number of beam particles " << beams_t.size() <<
" this generator interface should be fixed");
323 for (const auto& part: beams_t) HepMC3::Print::line(part);
324 }
327 ATH_MSG_WARNING(
"Invalid beam particle pointers -- this generator interface should be fixed");
328 if (cmenergy < 0)
ATH_MSG_WARNING(
"Invalid expected beam energy: " << cmenergy <<
" MeV");
330 } else {
332 ATH_MSG_WARNING(
"Beam particles have incorrectly set status -- this generator interface should be fixed");
334 }
335 const double sumE = beams.first->momentum().e() + beams.second->momentum().e();
336 const double sumP = beams.first->momentum().pz() + beams.second->momentum().pz();
337 cmenergy = std::sqrt(sumE*sumE - sumP*sumP);
338
341 }
342 if(beams.first->pdg_id() ==
MC::LEAD && beams.second->pdg_id() ==
MC::LEAD){
344 }
347 }
350 }
353 }
356 }
359 }
362 }
363
365 ATH_MSG_FATAL(
"Beam particles have incorrect energy: " <<
m_cm_energy/Gaudi::Units::GeV <<
" GeV expected, vs. " << cmenergy/Gaudi::Units::GeV <<
" GeV found");
371 }
373
374 return StatusCode::FAILURE;
375 }
376 }
377
378
379 int vtxDisplacedstatuscode12CheckRateCnt=0;
380 int vtxDisplacedstatuscodenot12CheckRateCnt=0;
381 int vtxDisplacedMoreThan_1m_CheckRateCnt=0;
382 for (
const auto& vtx:
evt->vertices()) {
384
385
386 if ( std::isnan(
pos.x()) || std::isinf(
pos.x()) ||
387 std::isnan(
pos.y()) || std::isinf(
pos.y()) ||
388 std::isnan(
pos.z()) || std::isinf(
pos.z()) ) {
389 ATH_MSG_WARNING(
"NaN (Not A Number) or inf found in the event record vertex positions");
390
392 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
394 filter_pass = false;
395 }
396 }
397
398
399
400 const double dist_trans2 =
pos.x()*
pos.x() +
pos.y()*
pos.y();
401 const double dist2 = dist_trans2 +
pos.z()*
pos.z();
402 const double dist_trans = std::sqrt(dist_trans2);
403 const double dist = std::sqrt(dist2);
406 ++vtxDisplacedMoreThan_1m_CheckRateCnt;
407
409 filter_pass = false;
410 }
411 }
414 << "mm in transverse distance: " << dist_trans << "mm");
415
416 for (const auto& part: vtx->particles_in()) {
419 HepMC::Print::line(
msg( MSG::WARNING ).
stream(), part);
420 }
421 }
422
424 filter_pass = false;
425 }
426 }
427
430
431 for (const auto& part: vtx->particles_in()) {
434 HepMC::Print::line(
msg( MSG::WARNING ).
stream(),part);
435 }
436 ATH_MSG_WARNING(
"production vertex = " <<
part->production_vertex()->position().x() <<
", " <<
part->production_vertex()->position().y() <<
", " <<
part->production_vertex()->position().z());
437 ATH_MSG_WARNING(
"end vertex = " <<
part->end_vertex()->position().x() <<
", " <<
part->end_vertex()->position().y() <<
", " <<
part->end_vertex()->position().z());
439 if (
part->production_vertex()) {
440 for(
const auto& p_parents:
part->production_vertex()->particles_in()) {
442 msg(MSG::WARNING) <<
"\t";
443 HepMC::Print::line(
msg( MSG::WARNING ).
stream() , p_parents );
444 }
445 }
446 }
447
448 if (
part->status() == 1 ||
part->status() == 2){
449 vtxDisplacedstatuscode12CheckRateCnt += 1;
450 } else {
451 vtxDisplacedstatuscodenot12CheckRateCnt += 1;
452 }
453
456 if (
part->momentum().e()/Gaudi::Units::GeV < 10.) {
458 }
470 double endvx =
part->end_vertex()->position().x();
471 double endvy =
part->end_vertex()->position().y();
472 double endvz =
part->end_vertex()->position().z();
473 double prodvx =
part->production_vertex()->position().x();
474 double prodvy =
part->production_vertex()->position().y();
475 double prodvz =
part->production_vertex()->position().z();
476 double enddis = std::sqrt(endvx*endvx + endvy*endvy + endvz*endvz);
477 double proddis = std::sqrt(prodvx*prodvx + prodvy*prodvy + prodvz*prodvz);
480 }
481 }
482 }
483 }
487
488
489 for (auto pitr: *evt) {
490
491
493 const int pstatus = pitr->status();
494 const int ppdgid = pitr->pdg_id();
495
496 if ( std::isnan(pmom.px()) || std::isinf(pmom.px()) ||
497 std::isnan(pmom.py()) || std::isinf(pmom.py()) ||
498 std::isnan(pmom.pz()) || std::isinf(pmom.pz()) ||
499 std::isnan(pmom.e()) || std::isinf(pmom.e()) ) {
500 ATH_MSG_WARNING(
"NaN (Not A Number) or inf found in the event record momenta");
502
503 if (
m_dumpEvent) HepMC::Print::line(std::cout,pitr);
505 filter_pass = false;
506 }
507 }
508
509
511 if (ppdgid == 111 && !pitr->end_vertex() ) {
512 unDecPi0.push_back( pitr);
514 }
515 }
516
517
520 if (pd != NULL) {
521 double plifetime =
pd->lifetime()*1
e+12;
522 if (plifetime != 0 && plifetime <
m_min_tau) {
523 ATH_MSG_WARNING(
"Stable particle found with lifetime = " << plifetime <<
"~ns!!");
524 if (
m_dumpEvent) HepMC::Print::line(std::cout,pitr);
525
527
529 filter_pass = false;
530 }
531 }
532 }
533 else{
534 int susyPart = 0;
535 std::vector<int>::size_type
count = 0;
537
539 susyPart=1;
540 }
542 }
543 if (susyPart==0){
544 ATH_MSG_WARNING(
"Stable particle not found in PDT, no lifetime check done");
545 if (
m_dumpEvent) HepMC::Print::line(std::cout,pitr);
546 }
547 }
548 }
549
550
551 const MC::DecodedPID decodedPID(ppdgid);
552 const int first_dig = decodedPID(0);
553
555
556 int known_byG4 = 0;
557 std::vector<int>::size_type
count =0;
558
562 }
563 if(known_byG4==0){
564 nonG4_energy += pmom.e();
565 ATH_MSG_WARNING(
"Interacting particle not known by Geant4 with ID " << ppdgid);
566 }
567 }
568
569
574 filter_pass = false;
576 }
577 }
578
579
581 unstNoEnd.push_back(pitr);
583 }
584
585
586
588 totalPx += pmom.px();
589 totalPy += pmom.py();
590 totalPz += pmom.pz();
591 totalE += pmom.e();
592 if (pmom.e() < 0) {
593 negEnPart.push_back(pitr);
595 }
596 const double aener = std::abs(pmom.e());
598 tachyons.push_back(pitr);
600 }
601 }
602
603
605 int tau_child = 0;
608 auto vtx = pitr->end_vertex();
609 if (vtx) {
610 double p_energy = 0;
611 for (auto desc: HepMC::descendant_particles(vtx)) {
612 if (std::abs(
desc->pdg_id()) ==
m_pdg) tau_child = 1;
614 }
615 if (std::abs( p_energy - pmom.e()) >
m_energy_diff && !tau_child) {
617 <<
"Energy (original particle) > " <<
m_energy_diff <<
" MeV, "
618 <<
"Event #" <<
evt->event_number() <<
", "
619 << "The original particle = " << pitr);
621 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
622 }
623
625 const double tau_displacement = tau_decaypos.x()*tau_decaypos.x() + tau_decaypos.y()*tau_decaypos.y() + tau_decaypos.z()*tau_decaypos.z();
626
628 } else {
631 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
632 }
633 }
634
635
636 if (pitr->end_vertex()) {
637 auto decayvtx = pitr->end_vertex();
639 const double displacement = decaypos.x()*decaypos.x() + decaypos.y()*decaypos.y() + decaypos.z()*decaypos.z();
640 if (displacement > 1e-6) {
641 for (auto ip: *decayvtx) {
643 const double displacement2 = pos2.x()*pos2.x() + pos2.y()*pos2.y() + pos2.z()*pos2.z();
644 if (displacement2 < 1e-6) {
646 <<
" has undisplaced vertex (" <<
ip->production_vertex()
647 << " @ " << displacement2 << "mm) "
648 << " but parent vertex is displaced (" << decayvtx
649 << " @ " << displacement << "mm)");
650 undisplaceds.push_back(std::move(ip));
652 }
653 }
654 }
655 }
656
657
660 const double mass = pitr->generated_mass();
661 if (std::abs(mass) > 1.0) {
662 ATH_MSG_WARNING(
"Photon with non-zero mass found! Mass: " << mass <<
" MeV" << pitr);
664 }
665 }
666
667 }
668
669
671 ATH_MSG_WARNING(
"The energy of interacting particles not known by Geant4 is = " << nonG4_energy <<
" MeV");
673 filter_pass = false;
674 }
676 }
677
678
679 double lostE = std::abs(totalE - cmenergy);
681 ATH_MSG_WARNING(
"ENERGY BALANCE FAILED : E-difference = " << lostE <<
" MeV");
682
683 ATH_MSG_WARNING(
"balance " << totalPx <<
" " << totalPy <<
" " << totalPz <<
" " << totalE);
684
687 }
688 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
690 filter_pass = false;
691 }
693 }
694
695
697 ATH_MSG_WARNING(
"MOMENTUM BALANCE FAILED : SumPx = " << totalPx <<
" SumPy = " << totalPy <<
" SumPz = " << totalPz <<
" MeV");
702 }
703 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
705 filter_pass = false;
706 }
708 }
709
710
711 if (!negEnPart.empty()) {
712 std::stringstream
ss;
713 ss <<
"NEGATIVE ENERGY PARTICLES FOUND :";
714 for (const auto &b: negEnPart){
716 }
718 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
720 filter_pass = false;
721 }
723 }
724
725
726 if (!tachyons.empty()) {
727 std::stringstream
ss;
728 ss <<
"PARTICLES WITH |E| < |Pi| (i=x,y,z) FOUND :";
729 for (auto b: tachyons){
731 }
733 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
735 filter_pass = false;
736 }
738 }
739
740
741 if (!unstNoEnd.empty()) {
742 std::stringstream
ss;
743 ss <<
"Unstable particle with no decay vertex found: ";
744 for (auto b: unstNoEnd){
746 }
748 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
750 filter_pass = false;
751 }
753 }
754
755
756 if (!unDecPi0.empty()) {
757 std::stringstream
ss;
758 ss <<
"pi0 with no decay vertex found:";
759 for (auto b: unDecPi0){
761 }
763 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
765 filter_pass = false;
766 }
768 }
769
770
771 if (!undisplaceds.empty()) {
772 std::stringstream
ss{
"Undisplaced decay vertices from displaced particle: "};
774
776 }
778 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
780 filter_pass = false;
781 }
783 }
784
785 }
786
787
788 if (!filter_pass){
791 } else {
793 }
794
795
796
800 return StatusCode::FAILURE;
801 }
802
803 return StatusCode::SUCCESS;
804}
#define ATH_MSG_WARNING(x)
void setFilterPassed(bool state) const
const HepPDT::ParticleData * particleData(int pid) const
Access an element in the particle data table.
std::vector< int > m_uknownPDGID_tab
MC::Loops< HepMC::GenEvent, HepMC::ConstGenParticlePtr, HepMC::ConstGenVertexPtr > m_looper
member to detect loops
std::vector< int > m_G4pdgID_tab
std::vector< int > m_SusyPdgID_tab
int count(std::string s, const std::string ®x)
count how many occurances of a regx are in a string
HepMC3::FourVector FourVector
HepMC3::ConstGenParticlePtr ConstGenParticlePtr
bool valid_beam_particles(const GenEvent *e)
HepMC3::GenEvent GenEvent
int numberOfProtons(const T &p)
bool isPhoton(const T &p)
bool isStable(const T &p)
Identify if the particle is stable, i.e. has not decayed.
bool isSimInteracting(const T &p)
Identify if the particle could interact with the detector during the simulation, e....
bool isDecayed(const T &p)
Identify if the particle decayed.
bool isBeam(const T &p)
Identify if the particle is beam particle.
bool isValid(const T &p)
Av: we implement here an ATLAS-sepcific convention: all particles which are 99xxxxx are fine.
bool isNucleus(const T &p)
PDG rule 16 Nuclear codes are given as 10-digit numbers ±10LZZZAAAI.
double baryonNumber(const T &p)
1.16.1