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
340 const int beam1PdgId = beams.first->pdg_id();
341 const int beam2PdgId = beams.second->pdg_id();
342
343
344
348 cmenergy /= std::sqrt(a1 * a2);
349 }
350 const bool nucleusProton =
352 const bool protonNucleus =
354
355 if (nucleusProton || protonNucleus) {
356 const int nucleusPdgId = nucleusProton ? beam1PdgId : beam2PdgId;
357 const auto& protonBeam = nucleusProton ? beams.second : beams.first;
358 const double direction = nucleusProton ? -1.0 : 1.0;
359
360 cmenergy =
361 direction * 2.0 * protonBeam->momentum().pz() *
362 std::sqrt(
365 );
366 }
367
369 ATH_MSG_FATAL(
"Beam particles have incorrect energy: " <<
m_cm_energy/Gaudi::Units::GeV <<
" GeV expected, vs. " << cmenergy/Gaudi::Units::GeV <<
" GeV found");
375 }
377
378 return StatusCode::FAILURE;
379 }
380 }
381
382
383 int vtxDisplacedstatuscode12CheckRateCnt=0;
384 int vtxDisplacedstatuscodenot12CheckRateCnt=0;
385 int vtxDisplacedMoreThan_1m_CheckRateCnt=0;
386 for (
const auto& vtx:
evt->vertices()) {
388
389
390 if ( std::isnan(
pos.x()) || std::isinf(
pos.x()) ||
391 std::isnan(
pos.y()) || std::isinf(
pos.y()) ||
392 std::isnan(
pos.z()) || std::isinf(
pos.z()) ) {
393 ATH_MSG_WARNING(
"NaN (Not A Number) or inf found in the event record vertex positions");
394
396 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
398 filter_pass = false;
399 }
400 }
401
402
403
404 const double dist_trans2 =
pos.x()*
pos.x() +
pos.y()*
pos.y();
405 const double dist2 = dist_trans2 +
pos.z()*
pos.z();
406 const double dist_trans = std::sqrt(dist_trans2);
407 const double dist = std::sqrt(dist2);
410 ++vtxDisplacedMoreThan_1m_CheckRateCnt;
411
413 filter_pass = false;
414 }
415 }
418 << "mm in transverse distance: " << dist_trans << "mm");
419
420 for (const auto& part: vtx->particles_in()) {
423 HepMC::Print::line(
msg( MSG::WARNING ).
stream(), part);
424 }
425 }
426
428 filter_pass = false;
429 }
430 }
431
434
435 for (const auto& part: vtx->particles_in()) {
438 HepMC::Print::line(
msg( MSG::WARNING ).
stream(),part);
439 }
440 ATH_MSG_WARNING(
"production vertex = " <<
part->production_vertex()->position().x() <<
", " <<
part->production_vertex()->position().y() <<
", " <<
part->production_vertex()->position().z());
441 ATH_MSG_WARNING(
"end vertex = " <<
part->end_vertex()->position().x() <<
", " <<
part->end_vertex()->position().y() <<
", " <<
part->end_vertex()->position().z());
443 if (
part->production_vertex()) {
444 for(
const auto& p_parents:
part->production_vertex()->particles_in()) {
446 msg(MSG::WARNING) <<
"\t";
447 HepMC::Print::line(
msg( MSG::WARNING ).
stream() , p_parents );
448 }
449 }
450 }
451
452 if (
part->status() == 1 ||
part->status() == 2){
453 vtxDisplacedstatuscode12CheckRateCnt += 1;
454 } else {
455 vtxDisplacedstatuscodenot12CheckRateCnt += 1;
456 }
457
460 if (
part->momentum().e()/Gaudi::Units::GeV < 10.) {
462 }
474 double endvx =
part->end_vertex()->position().x();
475 double endvy =
part->end_vertex()->position().y();
476 double endvz =
part->end_vertex()->position().z();
477 double prodvx =
part->production_vertex()->position().x();
478 double prodvy =
part->production_vertex()->position().y();
479 double prodvz =
part->production_vertex()->position().z();
480 double enddis = std::sqrt(endvx*endvx + endvy*endvy + endvz*endvz);
481 double proddis = std::sqrt(prodvx*prodvx + prodvy*prodvy + prodvz*prodvz);
484 }
485 }
486 }
487 }
491
492
493 for (auto pitr: *evt) {
494
495
497 const int pstatus = pitr->status();
498 const int ppdgid = pitr->pdg_id();
499
500 if ( std::isnan(pmom.px()) || std::isinf(pmom.px()) ||
501 std::isnan(pmom.py()) || std::isinf(pmom.py()) ||
502 std::isnan(pmom.pz()) || std::isinf(pmom.pz()) ||
503 std::isnan(pmom.e()) || std::isinf(pmom.e()) ) {
504 ATH_MSG_WARNING(
"NaN (Not A Number) or inf found in the event record momenta");
506
507 if (
m_dumpEvent) HepMC::Print::line(std::cout,pitr);
509 filter_pass = false;
510 }
511 }
512
513
515 if (ppdgid == 111 && !pitr->end_vertex() ) {
516 unDecPi0.push_back( pitr);
518 }
519 }
520
521
523 const auto plifetime =
m_gendata->particleLifetime(ppdgid);
524 if (plifetime) {
525 double lifetime = plifetime.value()*1
e+12;
526 if (lifetime != 0 && lifetime <
m_min_tau) {
527 ATH_MSG_WARNING(
"Stable particle found with lifetime = " << lifetime <<
"~ns!!");
528 if (
m_dumpEvent) HepMC::Print::line(std::cout,pitr);
529
531
533 filter_pass = false;
534 }
535 }
536 }
537 else{
538 int susyPart = 0;
539 std::vector<int>::size_type
count = 0;
541
543 susyPart=1;
544 }
546 }
547 if (susyPart==0){
548 ATH_MSG_WARNING(
"Stable particle not found in PDT, no lifetime check done");
549 if (
m_dumpEvent) HepMC::Print::line(std::cout,pitr);
550 }
551 }
552 }
553
554
555 const MC::DecodedPID decodedPID(ppdgid);
556 const int first_dig = decodedPID(0);
557
559
560 int known_byG4 = 0;
561 std::vector<int>::size_type
count =0;
562
566 }
567 if(known_byG4==0){
568 nonG4_energy += pmom.e();
569 ATH_MSG_WARNING(
"Interacting particle not known by Geant4 with ID " << ppdgid);
570 }
571 }
572
573
578 filter_pass = false;
580 }
581 }
582
583
585 unstNoEnd.push_back(pitr);
587 }
588
589
590
592 totalPx += pmom.px();
593 totalPy += pmom.py();
594 totalPz += pmom.pz();
595 totalE += pmom.e();
596 if (pmom.e() < 0) {
597 negEnPart.push_back(pitr);
599 }
600 const double aener = std::abs(pmom.e());
602 tachyons.push_back(pitr);
604 }
605 }
606
607
609 int tau_child = 0;
612 auto vtx = pitr->end_vertex();
613 if (vtx) {
614 double p_energy = 0;
615 for (auto desc: HepMC::descendant_particles(vtx)) {
616 if (std::abs(
desc->pdg_id()) ==
m_pdg) tau_child = 1;
618 }
619 if (std::abs( p_energy - pmom.e()) >
m_energy_diff && !tau_child) {
621 <<
"Energy (original particle) > " <<
m_energy_diff <<
" MeV, "
622 <<
"Event #" <<
evt->event_number() <<
", "
623 << "The original particle = " << pitr);
625 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
626 }
627
629 const double tau_displacement = tau_decaypos.x()*tau_decaypos.x() + tau_decaypos.y()*tau_decaypos.y() + tau_decaypos.z()*tau_decaypos.z();
630
632 } else {
635 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
636 }
637 }
638
639
640 if (pitr->end_vertex()) {
641 auto decayvtx = pitr->end_vertex();
643 const double displacement = decaypos.x()*decaypos.x() + decaypos.y()*decaypos.y() + decaypos.z()*decaypos.z();
644 if (displacement > 1e-6) {
645 for (auto ip: *decayvtx) {
647 const double displacement2 = pos2.x()*pos2.x() + pos2.y()*pos2.y() + pos2.z()*pos2.z();
648 if (displacement2 < 1e-6) {
650 <<
" has undisplaced vertex (" <<
ip->production_vertex()
651 << " @ " << displacement2 << "mm) "
652 << " but parent vertex is displaced (" << decayvtx
653 << " @ " << displacement << "mm)");
654 undisplaceds.push_back(std::move(ip));
656 }
657 }
658 }
659 }
660
661
664 const double mass = pitr->generated_mass();
665 if (std::abs(mass) > 1.0) {
666 ATH_MSG_WARNING(
"Photon with non-zero mass found! Mass: " << mass <<
" MeV" << pitr);
668 }
669 }
670
671 }
672
673
675 ATH_MSG_WARNING(
"The energy of interacting particles not known by Geant4 is = " << nonG4_energy <<
" MeV");
677 filter_pass = false;
678 }
680 }
681
682
683 double lostE = std::abs(totalE - cmenergy);
685 ATH_MSG_WARNING(
"ENERGY BALANCE FAILED : E-difference = " << lostE <<
" MeV");
686
687 ATH_MSG_WARNING(
"balance " << totalPx <<
" " << totalPy <<
" " << totalPz <<
" " << totalE);
688
691 }
692 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
694 filter_pass = false;
695 }
697 }
698
699
701 ATH_MSG_WARNING(
"MOMENTUM BALANCE FAILED : SumPx = " << totalPx <<
" SumPy = " << totalPy <<
" SumPz = " << totalPz <<
" MeV");
706 }
707 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
709 filter_pass = false;
710 }
712 }
713
714
715 if (!negEnPart.empty()) {
716 std::stringstream
ss;
717 ss <<
"NEGATIVE ENERGY PARTICLES FOUND :";
718 for (const auto &b: negEnPart){
720 }
722 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
724 filter_pass = false;
725 }
727 }
728
729
730 if (!tachyons.empty()) {
731 std::stringstream
ss;
732 ss <<
"PARTICLES WITH |E| < |Pi| (i=x,y,z) FOUND :";
733 for (auto b: tachyons){
735 }
737 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
739 filter_pass = false;
740 }
742 }
743
744
745 if (!unstNoEnd.empty()) {
746 std::stringstream
ss;
747 ss <<
"Unstable particle with no decay vertex found: ";
748 for (auto b: unstNoEnd){
750 }
752 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
754 filter_pass = false;
755 }
757 }
758
759
760 if (!unDecPi0.empty()) {
761 std::stringstream
ss;
762 ss <<
"pi0 with no decay vertex found:";
763 for (auto b: unDecPi0){
765 }
767 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
769 filter_pass = false;
770 }
772 }
773
774
775 if (!undisplaceds.empty()) {
776 std::stringstream
ss{
"Undisplaced decay vertices from displaced particle: "};
778
780 }
782 if (
m_dumpEvent) HepMC::Print::content(std::cout,*evt);
784 filter_pass = false;
785 }
787 }
788
789 }
790
791
792 if (!filter_pass){
795 } else {
797 }
798
799
800
804 return StatusCode::FAILURE;
805 }
806
807 return StatusCode::SUCCESS;
808}
#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)
const Amg::Vector3D & direction() const
Method to retrieve the direction at the Intersection.
1.17.0