Execute method.
271 {
272
273
274 bool filter_pass = true;
275
276
278 double totalPx = 0;
279 double totalPy = 0;
280 double totalPz = 0;
281 double totalE = 0;
282 double nonG4_energy = 0;
283 std::vector<HepMC::ConstGenParticlePtr> negEnPart;
284 std::vector<HepMC::ConstGenParticlePtr> tachyons;
285 std::vector<HepMC::ConstGenParticlePtr> unstNoEnd;
286 std::vector<HepMC::ConstGenParticlePtr> unDecPi0;
287 std::vector<HepMC::ConstGenParticlePtr> undisplaceds;
288
290 if (!
m_looper.loop_particles().empty() || !
m_looper.loop_vertices().empty()) {
293 ATH_MSG_DEBUG(
"Please use MC::Loops::findLoops for this event to obtain all particles and vertices in the loops");
295 }
296
297 const auto xsec =
evt->cross_section();
298 if (!xsec) {
299 ATH_MSG_WARNING(
"WATCH OUT: event is missing the generator cross-section!");
302 ATH_MSG_WARNING(
"-> Adding a dummy cross-section for debugging purposes.");
303
304 std::shared_ptr<HepMC3::GenCrossSection> dummy_xsec = std::make_shared<HepMC3::GenCrossSection>();
305 dummy_xsec->set_cross_section(1.0,0.0);
307 evt_nonconst->set_cross_section(std::move(dummy_xsec));
308 }
309 else {
311 }
312 }
313
314
315 std::vector<std::shared_ptr<const HepMC3::GenParticle>> beams_t;
316 for (
auto p :
evt->beams()) {
if (
p->status() == 4) beams_t.push_back(std::move(p)); }
317 std::pair<std::shared_ptr<const HepMC3::GenParticle>,std::shared_ptr<const HepMC3::GenParticle>> beams;
318 if (beams_t.size() == 2) {
319 beams.first=beams_t.at(0);
320 beams.second=beams_t.at(1);
321 } else {
322 ATH_MSG_WARNING(
"Invalid number of beam particles " << beams_t.size() <<
" this generator interface should be fixed");
324 for (const auto& part: beams_t) HepMC3::Print::line(part);
325 }
328 ATH_MSG_WARNING(
"Invalid beam particle pointers -- this generator interface should be fixed");
329 if (cmenergy < 0)
ATH_MSG_WARNING(
"Invalid expected beam energy: " << cmenergy <<
" MeV");
331 } else {
333 ATH_MSG_WARNING(
"Beam particles have incorrectly set status -- this generator interface should be fixed");
335 }
336 const double sumE = beams.first->momentum().e() + beams.second->momentum().e();
337 const double sumP = beams.first->momentum().pz() + beams.second->momentum().pz();
338 cmenergy = std::sqrt(sumE*sumE - sumP*sumP);
339
342 }
343 if(beams.first->pdg_id() ==
MC::LEAD && beams.second->pdg_id() ==
MC::LEAD){
345 }
348 }
351 }
354 }
357 }
360 }
363 }
364
366 ATH_MSG_FATAL(
"Beam particles have incorrect energy: " <<
m_cm_energy/Gaudi::Units::GeV <<
" GeV expected, vs. " << cmenergy/Gaudi::Units::GeV <<
" GeV found");
372 }
374
375 return StatusCode::FAILURE;
376 }
377 }
378
379
380 int vtxDisplacedstatuscode12CheckRateCnt=0;
381 int vtxDisplacedstatuscodenot12CheckRateCnt=0;
382 int vtxDisplacedMoreThan_1m_CheckRateCnt=0;
383 for (
const auto& vtx:
evt->vertices()) {
385
386
387 if ( std::isnan(
pos.x()) || std::isinf(
pos.x()) ||
388 std::isnan(
pos.y()) || std::isinf(
pos.y()) ||
389 std::isnan(
pos.z()) || std::isinf(
pos.z()) ) {
390 ATH_MSG_WARNING(
"NaN (Not A Number) or inf found in the event record vertex positions");
391
393 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
395 filter_pass = false;
396 }
397 }
398
399
400
401 const double dist_trans2 =
pos.x()*
pos.x() +
pos.y()*
pos.y();
402 const double dist2 = dist_trans2 +
pos.z()*
pos.z();
403 const double dist_trans = std::sqrt(dist_trans2);
404 const double dist = std::sqrt(dist2);
407 ++vtxDisplacedMoreThan_1m_CheckRateCnt;
408
410 filter_pass = false;
411 }
412 }
415 << "mm in transverse distance: " << dist_trans << "mm");
416
417 for (const auto& part: vtx->particles_in()) {
420 HepMC::Print::line(
msg( MSG::WARNING ).
stream(), part);
421 }
422 }
423
425 filter_pass = false;
426 }
427 }
428
431
432 for (const auto& part: vtx->particles_in()) {
435 HepMC::Print::line(
msg( MSG::WARNING ).
stream(),part);
436 }
437 ATH_MSG_WARNING(
"production vertex = " <<
part->production_vertex()->position().x() <<
", " <<
part->production_vertex()->position().y() <<
", " <<
part->production_vertex()->position().z());
438 ATH_MSG_WARNING(
"end vertex = " <<
part->end_vertex()->position().x() <<
", " <<
part->end_vertex()->position().y() <<
", " <<
part->end_vertex()->position().z());
440 if (
part->production_vertex()) {
441 for(
const auto& p_parents:
part->production_vertex()->particles_in()) {
443 msg(MSG::WARNING) <<
"\t";
444 HepMC::Print::line(
msg( MSG::WARNING ).
stream() , p_parents );
445 }
446 }
447 }
448
449 if (
part->status() == 1 ||
part->status() == 2){
450 vtxDisplacedstatuscode12CheckRateCnt += 1;
451 } else {
452 vtxDisplacedstatuscodenot12CheckRateCnt += 1;
453 }
454
457 if (
part->momentum().e()/Gaudi::Units::GeV < 10.) {
459 }
471 double endvx =
part->end_vertex()->position().x();
472 double endvy =
part->end_vertex()->position().y();
473 double endvz =
part->end_vertex()->position().z();
474 double prodvx =
part->production_vertex()->position().x();
475 double prodvy =
part->production_vertex()->position().y();
476 double prodvz =
part->production_vertex()->position().z();
477 double enddis = std::sqrt(endvx*endvx + endvy*endvy + endvz*endvz);
478 double proddis = std::sqrt(prodvx*prodvx + prodvy*prodvy + prodvz*prodvz);
481 }
482 }
483 }
484 }
488
489
490 for (auto pitr: *evt) {
491
492
494 const int pstatus = pitr->status();
495 const int ppdgid = pitr->pdg_id();
496
497 if ( std::isnan(pmom.px()) || std::isinf(pmom.px()) ||
498 std::isnan(pmom.py()) || std::isinf(pmom.py()) ||
499 std::isnan(pmom.pz()) || std::isinf(pmom.pz()) ||
500 std::isnan(pmom.e()) || std::isinf(pmom.e()) ) {
501 ATH_MSG_WARNING(
"NaN (Not A Number) or inf found in the event record momenta");
503
504 if (
m_dumpEvent) HepMC::Print::line(std::cout,pitr);
506 filter_pass = false;
507 }
508 }
509
510
512 if (ppdgid == 111 && !pitr->end_vertex() ) {
513 unDecPi0.push_back( pitr);
515 }
516 }
517
518
520 const auto plifetime =
m_gendata->particleLifetime(ppdgid);
521 if (plifetime) {
522 double lifetime = plifetime.value()*1
e+12;
523 if (lifetime != 0 && lifetime <
m_min_tau) {
524 ATH_MSG_WARNING(
"Stable particle found with lifetime = " << lifetime <<
"~ns!!");
525 if (
m_dumpEvent) HepMC::Print::line(std::cout,pitr);
526
528
530 filter_pass = false;
531 }
532 }
533 }
534 else{
535 int susyPart = 0;
536 std::vector<int>::size_type
count = 0;
538
540 susyPart=1;
541 }
543 }
544 if (susyPart==0){
545 ATH_MSG_WARNING(
"Stable particle not found in PDT, no lifetime check done");
546 if (
m_dumpEvent) HepMC::Print::line(std::cout,pitr);
547 }
548 }
549 }
550
551
552 const MC::DecodedPID decodedPID(ppdgid);
553 const int first_dig = decodedPID(0);
554
556
557 int known_byG4 = 0;
558 std::vector<int>::size_type
count =0;
559
563 }
564 if(known_byG4==0){
565 nonG4_energy += pmom.e();
566 ATH_MSG_WARNING(
"Interacting particle not known by Geant4 with ID " << ppdgid);
567 }
568 }
569
570
575 filter_pass = false;
577 }
578 }
579
580
582 unstNoEnd.push_back(pitr);
584 }
585
586
587
589 totalPx += pmom.px();
590 totalPy += pmom.py();
591 totalPz += pmom.pz();
592 totalE += pmom.e();
593 if (pmom.e() < 0) {
594 negEnPart.push_back(pitr);
596 }
597 const double aener = std::abs(pmom.e());
599 tachyons.push_back(pitr);
601 }
602 }
603
604
606 int tau_child = 0;
609 auto vtx = pitr->end_vertex();
610 if (vtx) {
611 double p_energy = 0;
612 for (auto desc: HepMC::descendant_particles(vtx)) {
613 if (std::abs(
desc->pdg_id()) ==
m_pdg) tau_child = 1;
615 }
616 if (std::abs( p_energy - pmom.e()) >
m_energy_diff && !tau_child) {
618 <<
"Energy (original particle) > " <<
m_energy_diff <<
" MeV, "
619 <<
"Event #" <<
evt->event_number() <<
", "
620 << "The original particle = " << pitr);
622 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
623 }
624
626 const double tau_displacement = tau_decaypos.x()*tau_decaypos.x() + tau_decaypos.y()*tau_decaypos.y() + tau_decaypos.z()*tau_decaypos.z();
627
629 } else {
632 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
633 }
634 }
635
636
637 if (pitr->end_vertex()) {
638 auto decayvtx = pitr->end_vertex();
640 const double displacement = decaypos.x()*decaypos.x() + decaypos.y()*decaypos.y() + decaypos.z()*decaypos.z();
641 if (displacement > 1e-6) {
642 for (auto ip: *decayvtx) {
644 const double displacement2 = pos2.x()*pos2.x() + pos2.y()*pos2.y() + pos2.z()*pos2.z();
645 if (displacement2 < 1e-6) {
647 <<
" has undisplaced vertex (" <<
ip->production_vertex()
648 << " @ " << displacement2 << "mm) "
649 << " but parent vertex is displaced (" << decayvtx
650 << " @ " << displacement << "mm)");
651 undisplaceds.push_back(std::move(ip));
653 }
654 }
655 }
656 }
657
658
661 const double mass = pitr->generated_mass();
662 if (std::abs(mass) > 1.0) {
663 ATH_MSG_WARNING(
"Photon with non-zero mass found! Mass: " << mass <<
" MeV" << pitr);
665 }
666 }
667
668 }
669
670
672 ATH_MSG_WARNING(
"The energy of interacting particles not known by Geant4 is = " << nonG4_energy <<
" MeV");
674 filter_pass = false;
675 }
677 }
678
679
680 double lostE = std::abs(totalE - cmenergy);
682 ATH_MSG_WARNING(
"ENERGY BALANCE FAILED : E-difference = " << lostE <<
" MeV");
683
684 ATH_MSG_WARNING(
"balance " << totalPx <<
" " << totalPy <<
" " << totalPz <<
" " << totalE);
685
688 }
689 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
691 filter_pass = false;
692 }
694 }
695
696
698 ATH_MSG_WARNING(
"MOMENTUM BALANCE FAILED : SumPx = " << totalPx <<
" SumPy = " << totalPy <<
" SumPz = " << totalPz <<
" MeV");
703 }
704 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
706 filter_pass = false;
707 }
709 }
710
711
712 if (!negEnPart.empty()) {
713 std::stringstream
ss;
714 ss <<
"NEGATIVE ENERGY PARTICLES FOUND :";
715 for (const auto &b: negEnPart){
717 }
719 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
721 filter_pass = false;
722 }
724 }
725
726
727 if (!tachyons.empty()) {
728 std::stringstream
ss;
729 ss <<
"PARTICLES WITH |E| < |Pi| (i=x,y,z) FOUND :";
730 for (auto b: tachyons){
732 }
734 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
736 filter_pass = false;
737 }
739 }
740
741
742 if (!unstNoEnd.empty()) {
743 std::stringstream
ss;
744 ss <<
"Unstable particle with no decay vertex found: ";
745 for (auto b: unstNoEnd){
747 }
749 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
751 filter_pass = false;
752 }
754 }
755
756
757 if (!unDecPi0.empty()) {
758 std::stringstream
ss;
759 ss <<
"pi0 with no decay vertex found:";
760 for (auto b: unDecPi0){
762 }
764 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
766 filter_pass = false;
767 }
769 }
770
771
772 if (!undisplaceds.empty()) {
773 std::stringstream
ss{
"Undisplaced decay vertices from displaced particle: "};
775
777 }
779 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
781 filter_pass = false;
782 }
784 }
785
786 }
787
788
789 if (!filter_pass){
792 } else {
794 }
795
796
797
801 return StatusCode::FAILURE;
802 }
803
804 return StatusCode::SUCCESS;
805}
#define ATH_MSG_WARNING(x)
virtual void setFilterPassed(bool state, const EventContext &ctx) const
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::shared_ptr< GenData > m_gendata
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