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#ifdef HEPMC3
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);
306 HepMC::GenEvent* evt_nonconst =
const_cast<HepMC::GenEvent*
>(
evt);
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 }
326#else
327 auto beams =
evt->beam_particles();
328#endif
331 ATH_MSG_WARNING(
"Invalid beam particle pointers -- this generator interface should be fixed");
332 if (cmenergy < 0)
ATH_MSG_WARNING(
"Invalid expected beam energy: " << cmenergy <<
" MeV");
334 } else {
336 ATH_MSG_WARNING(
"Beam particles have incorrectly set status -- this generator interface should be fixed");
338 }
339 const double sumE = beams.first->momentum().e() + beams.second->momentum().e();
340 const double sumP = beams.first->momentum().pz() + beams.second->momentum().pz();
341 cmenergy = std::sqrt(sumE*sumE - sumP*sumP);
342
345 }
346 if(beams.first->pdg_id() ==
MC::LEAD && beams.second->pdg_id() ==
MC::LEAD){
348 }
351 }
354 }
357 }
360 }
363 }
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");
370 setFilterPassed(false);
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#ifdef HEPMC3
387 for (
const auto& vtx:
evt->vertices()) {
388#else
389 for (
auto vitr =
evt->vertices_begin(); vitr !=
evt->vertices_end(); ++vitr ) {
390 const HepMC::GenVertex* vtx = *vitr;
391#endif
392 const HepMC::FourVector
pos = vtx->position();
393
394
395 if ( std::isnan(
pos.x()) || std::isinf(
pos.x()) ||
396 std::isnan(
pos.y()) || std::isinf(
pos.y()) ||
397 std::isnan(
pos.z()) || std::isinf(
pos.z()) ) {
398 ATH_MSG_WARNING(
"NaN (Not A Number) or inf found in the event record vertex positions");
399
403 filter_pass = false;
404 }
405 }
406
407
408
409 const double dist_trans2 =
pos.x()*
pos.x() +
pos.y()*
pos.y();
410 const double dist2 = dist_trans2 +
pos.z()*
pos.z();
411 const double dist_trans = std::sqrt(dist_trans2);
412 const double dist = std::sqrt(dist2);
415 ++vtxDisplacedMoreThan_1m_CheckRateCnt;
416
418 filter_pass = false;
419 }
420 }
423 << "mm in transverse distance: " << dist_trans << "mm");
424
425#ifdef HEPMC3
426 for (const auto& part: vtx->particles_in()) {
427#else
428 for (auto part_it = vtx->particles_in_const_begin(); part_it != vtx->particles_in_const_end(); ++part_it) {
429 auto part = (*part_it);
430#endif
434 }
435 }
436
438 filter_pass = false;
439 }
440 }
441
444
445#ifdef HEPMC3
446 for (const auto& part: vtx->particles_in()) {
447#else
448 for (auto part_it = vtx->particles_in_const_begin(); part_it != vtx->particles_in_const_end(); ++part_it) {
449 auto part=(*part_it);
450#endif
454 }
455 ATH_MSG_WARNING(
"production vertex = " <<
part->production_vertex()->position().x() <<
", " <<
part->production_vertex()->position().y() <<
", " <<
part->production_vertex()->position().z());
456 ATH_MSG_WARNING(
"end vertex = " <<
part->end_vertex()->position().x() <<
", " <<
part->end_vertex()->position().y() <<
", " <<
part->end_vertex()->position().z());
458 if (
part->production_vertex()) {
459#ifdef HEPMC3
460 for(
const auto& p_parents:
part->production_vertex()->particles_in()) {
461#else
462 for(
auto p_parents_it =
part->production_vertex()->particles_in_const_begin(); p_parents_it !=
part->production_vertex()->particles_in_const_end(); ++p_parents_it) {
463 auto p_parents=(*p_parents_it);
464#endif
466 msg(MSG::WARNING) <<
"\t";
468 }
469 }
470 }
471
472 if (
part->status() == 1 ||
part->status() == 2){
473 vtxDisplacedstatuscode12CheckRateCnt += 1;
474 } else {
475 vtxDisplacedstatuscodenot12CheckRateCnt += 1;
476 }
477
480 if (
part->momentum().e()/Gaudi::Units::GeV < 10.) {
482 }
494 double endvx =
part->end_vertex()->position().x();
495 double endvy =
part->end_vertex()->position().y();
496 double endvz =
part->end_vertex()->position().z();
497 double prodvx =
part->production_vertex()->position().x();
498 double prodvy =
part->production_vertex()->position().y();
499 double prodvz =
part->production_vertex()->position().z();
500 double enddis = std::sqrt(endvx*endvx + endvy*endvy + endvz*endvz);
501 double proddis = std::sqrt(prodvx*prodvx + prodvy*prodvy + prodvz*prodvz);
504 }
505 }
506 }
507 }
511
512
513 for (auto pitr: *evt) {
514
515
516 const HepMC::FourVector pmom = pitr->momentum();
517 const int pstatus = pitr->status();
518 const int ppdgid = pitr->pdg_id();
519
520 if ( std::isnan(pmom.px()) || std::isinf(pmom.px()) ||
521 std::isnan(pmom.py()) || std::isinf(pmom.py()) ||
522 std::isnan(pmom.pz()) || std::isinf(pmom.pz()) ||
523 std::isnan(pmom.e()) || std::isinf(pmom.e()) ) {
524 ATH_MSG_WARNING(
"NaN (Not A Number) or inf found in the event record momenta");
526
529 filter_pass = false;
530 }
531 }
532
533
535 if (ppdgid == 111 && !pitr->end_vertex() ) {
536 unDecPi0.push_back( pitr);
538 }
539 }
540
541
544 if (pd != NULL) {
545 double plifetime =
pd->lifetime()*1
e+12;
546 if (plifetime != 0 && plifetime <
m_min_tau) {
547 ATH_MSG_WARNING(
"Stable particle found with lifetime = " << plifetime <<
"~ns!!");
549
551
553 filter_pass = false;
554 }
555 }
556 }
557 else{
558 int susyPart = 0;
559 std::vector<int>::size_type
count = 0;
561
563 susyPart=1;
564 }
566 }
567 if (susyPart==0){
568 ATH_MSG_WARNING(
"Stable particle not found in PDT, no lifetime check done");
570 }
571 }
572 }
573
574
575 const MC::DecodedPID decodedPID(ppdgid);
576 const int first_dig = decodedPID(0);
577
579
580 int known_byG4 = 0;
581 std::vector<int>::size_type
count =0;
582
586 }
587 if(known_byG4==0){
588 nonG4_energy += pmom.e();
589 ATH_MSG_WARNING(
"Interacting particle not known by Geant4 with ID " << ppdgid);
590 }
591 }
592
593
598 filter_pass = false;
600 }
601 }
602
603
605 unstNoEnd.push_back(pitr);
607 }
608
609
610
612 totalPx += pmom.px();
613 totalPy += pmom.py();
614 totalPz += pmom.pz();
615 totalE += pmom.e();
616 if (pmom.e() < 0) {
617 negEnPart.push_back(pitr);
619 }
620 const double aener = std::abs(pmom.e());
622 tachyons.push_back(pitr);
624 }
625 }
626
627
629 int tau_child = 0;
632 auto vtx = pitr->end_vertex();
633 if (vtx) {
634 double p_energy = 0;
635#ifdef HEPMC3
636 for (auto desc: HepMC::descendant_particles(vtx)) {
637#else
638 for (auto desc_it = vtx->particles_begin(HepMC::descendants); desc_it != vtx->particles_end(HepMC::descendants); ++desc_it) {
639 auto desc=(*desc_it);
640#endif
641 if (std::abs(
desc->pdg_id()) ==
m_pdg) tau_child = 1;
643 }
644 if (std::abs( p_energy - pmom.e()) >
m_energy_diff && !tau_child) {
646 <<
"Energy (original particle) > " <<
m_energy_diff <<
" MeV, "
647 <<
"Event #" <<
evt->event_number() <<
", "
648 << "The original particle = " << pitr);
651 }
652
653 const HepMC::FourVector tau_decaypos = vtx->position();
654 const double tau_displacement = tau_decaypos.x()*tau_decaypos.x() + tau_decaypos.y()*tau_decaypos.y() + tau_decaypos.z()*tau_decaypos.z();
655
657 } else {
661 }
662 }
663
664
665 if (pitr->end_vertex()) {
666 auto decayvtx = pitr->end_vertex();
667 const HepMC::FourVector decaypos = decayvtx->position();
668 const double displacement = decaypos.x()*decaypos.x() + decaypos.y()*decaypos.y() + decaypos.z()*decaypos.z();
669 if (displacement > 1e-6) {
670 for (auto ip: *decayvtx) {
671 const HepMC::FourVector pos2 =
ip->production_vertex()->position();
672 const double displacement2 = pos2.x()*pos2.x() + pos2.y()*pos2.y() + pos2.z()*pos2.z();
673 if (displacement2 < 1e-6) {
675 <<
" has undisplaced vertex (" <<
ip->production_vertex()
676 << " @ " << displacement2 << "mm) "
677 << " but parent vertex is displaced (" << decayvtx
678 << " @ " << displacement << "mm)");
679 undisplaceds.push_back(std::move(ip));
681 }
682 }
683 }
684 }
685
686
689 const double mass = pitr->generated_mass();
690 if (std::abs(mass) > 1.0) {
691 ATH_MSG_WARNING(
"Photon with non-zero mass found! Mass: " << mass <<
" MeV" << pitr);
693 }
694 }
695
696 }
697
698
700 ATH_MSG_WARNING(
"The energy of interacting particles not known by Geant4 is = " << nonG4_energy <<
" MeV");
702 filter_pass = false;
703 }
705 }
706
707
708 double lostE = std::abs(totalE - cmenergy);
710 ATH_MSG_WARNING(
"ENERGY BALANCE FAILED : E-difference = " << lostE <<
" MeV");
711
712 ATH_MSG_WARNING(
"balance " << totalPx <<
" " << totalPy <<
" " << totalPz <<
" " << totalE);
713
716 }
719 filter_pass = false;
720 }
722 }
723
724
726 ATH_MSG_WARNING(
"MOMENTUM BALANCE FAILED : SumPx = " << totalPx <<
" SumPy = " << totalPy <<
" SumPz = " << totalPz <<
" MeV");
731 }
734 filter_pass = false;
735 }
737 }
738
739
740 if (!negEnPart.empty()) {
741 std::stringstream
ss;
742 ss <<
"NEGATIVE ENERGY PARTICLES FOUND :";
743 for (const auto &b: negEnPart){
745 }
749 filter_pass = false;
750 }
752 }
753
754
755 if (!tachyons.empty()) {
756 std::stringstream
ss;
757 ss <<
"PARTICLES WITH |E| < |Pi| (i=x,y,z) FOUND :";
758 for (auto b: tachyons){
760 }
764 filter_pass = false;
765 }
767 }
768
769
770 if (!unstNoEnd.empty()) {
771 std::stringstream
ss;
772 ss <<
"Unstable particle with no decay vertex found: ";
773 for (auto b: unstNoEnd){
775 }
779 filter_pass = false;
780 }
782 }
783
784
785 if (!unDecPi0.empty()) {
786 std::stringstream
ss;
787 ss <<
"pi0 with no decay vertex found:";
788 for (auto b: unDecPi0){
790 }
794 filter_pass = false;
795 }
797 }
798
799
800 if (!undisplaceds.empty()) {
801 std::stringstream
ss{
"Undisplaced decay vertices from displaced particle: "};
803
805 }
809 filter_pass = false;
810 }
812 }
813
814 }
815
816
817 if (!filter_pass){
818 setFilterPassed(false);
820 } else {
822 }
823
824
825
829 return StatusCode::FAILURE;
830 }
831
832 return StatusCode::SUCCESS;
833}
#define ATH_MSG_WARNING(x)
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
void line(std::ostream &os, const GenEvent &e)
void content(std::ostream &os, const GenEvent &e)
const GenParticle * ConstGenParticlePtr
bool valid_beam_particles(const GenEvent *e)
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.14.0