269 {
270
271
272 bool filter_pass = true;
273
274
276 double totalPx = 0;
277 double totalPy = 0;
278 double totalPz = 0;
279 double totalE = 0;
280 double nonG4_energy = 0;
281 std::vector<HepMC::ConstGenParticlePtr> negEnPart;
282 std::vector<HepMC::ConstGenParticlePtr> tachyons;
283 std::vector<HepMC::ConstGenParticlePtr> unstNoEnd;
284 std::vector<HepMC::ConstGenParticlePtr> unDecPi0;
285 std::vector<HepMC::ConstGenParticlePtr> undisplaceds;
286
288 if (!
m_looper.loop_particles().empty() || !
m_looper.loop_vertices().empty()) {
291 ATH_MSG_DEBUG(
"Please use MC::Loops::findLoops for this event to obtain all particles and vertices in the loops");
293 }
294
295#ifdef HEPMC3
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);
305 HepMC::GenEvent* evt_nonconst =
const_cast<HepMC::GenEvent*
>(
evt);
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 }
325#else
326 auto beams =
evt->beam_particles();
327#endif
330 ATH_MSG_WARNING(
"Invalid beam particle pointers -- this generator interface should be fixed");
331 if (cmenergy < 0)
ATH_MSG_WARNING(
"Invalid expected beam energy: " << cmenergy <<
" MeV");
333 } else {
335 ATH_MSG_WARNING(
"Beam particles have incorrectly set status -- this generator interface should be fixed");
337 }
338 const double sumE = beams.first->momentum().e() + beams.second->momentum().e();
339 const double sumP = beams.first->momentum().pz() + beams.second->momentum().pz();
340 cmenergy = std::sqrt(sumE*sumE - sumP*sumP);
341
344 }
345 if(beams.first->pdg_id() ==
MC::LEAD && beams.second->pdg_id() ==
MC::LEAD){
347 }
350 }
353 }
356 }
359 }
362 }
365 }
366
368 ATH_MSG_FATAL(
"Beam particles have incorrect energy: " <<
m_cm_energy/Gaudi::Units::GeV <<
" GeV expected, vs. " << cmenergy/Gaudi::Units::GeV <<
" GeV found");
369 setFilterPassed(false);
374 }
376
377 return StatusCode::FAILURE;
378 }
379 }
380
381
382 int vtxDisplacedstatuscode12CheckRateCnt=0;
383 int vtxDisplacedstatuscodenot12CheckRateCnt=0;
384 int vtxDisplacedMoreThan_1m_CheckRateCnt=0;
385#ifdef HEPMC3
386 for (
const auto& vtx:
evt->vertices()) {
387#else
388 for (
auto vitr =
evt->vertices_begin(); vitr !=
evt->vertices_end(); ++vitr ) {
389 const HepMC::GenVertex* vtx = *vitr;
390#endif
391 const HepMC::FourVector
pos = vtx->position();
392
393
394 if ( std::isnan(
pos.x()) || std::isinf(
pos.x()) ||
395 std::isnan(
pos.y()) || std::isinf(
pos.y()) ||
396 std::isnan(
pos.z()) || std::isinf(
pos.z()) ) {
397 ATH_MSG_WARNING(
"NaN (Not A Number) or inf found in the event record vertex positions");
398
402 filter_pass = false;
403 }
404 }
405
406
407
408 const double dist_trans2 =
pos.x()*
pos.x() +
pos.y()*
pos.y();
409 const double dist2 = dist_trans2 +
pos.z()*
pos.z();
410 const double dist_trans = std::sqrt(dist_trans2);
411 const double dist = std::sqrt(dist2);
414 ++vtxDisplacedMoreThan_1m_CheckRateCnt;
415
417 filter_pass = false;
418 }
419 }
422
423#ifdef HEPMC3
424 for (const auto& part: vtx->particles_in()) {
425#else
426 for (auto part_it = vtx->particles_in_const_begin(); part_it != vtx->particles_in_const_end(); ++part_it) {
427 auto part=(*part_it);
428#endif
432 }
433 ATH_MSG_WARNING(
"production vertex = " <<
part->production_vertex()->position().x() <<
", " <<
part->production_vertex()->position().y() <<
", " <<
part->production_vertex()->position().z());
434 ATH_MSG_WARNING(
"end vertex = " <<
part->end_vertex()->position().x() <<
", " <<
part->end_vertex()->position().y() <<
", " <<
part->end_vertex()->position().z());
436 if (
part->production_vertex()) {
437#ifdef HEPMC3
438 for(
const auto& p_parents:
part->production_vertex()->particles_in()) {
439#else
440 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) {
441 auto p_parents=(*p_parents_it);
442#endif
444 msg(MSG::WARNING) <<
"\t";
446 }
447 }
448 }
449
450 if (
part->status() == 1 ||
part->status() == 2){
451 vtxDisplacedstatuscode12CheckRateCnt += 1;
452 } else {
453 vtxDisplacedstatuscodenot12CheckRateCnt += 1;
454 }
455
458 if (
part->momentum().e()/Gaudi::Units::GeV < 10.) {
460 }
472 double endvx =
part->end_vertex()->position().x();
473 double endvy =
part->end_vertex()->position().y();
474 double endvz =
part->end_vertex()->position().z();
475 double prodvx =
part->production_vertex()->position().x();
476 double prodvy =
part->production_vertex()->position().y();
477 double prodvz =
part->production_vertex()->position().z();
478 double enddis = std::sqrt(endvx*endvx + endvy*endvy + endvz*endvz);
479 double proddis = std::sqrt(prodvx*prodvx + prodvy*prodvy + prodvz*prodvz);
482 }
483 }
484 }
485 }
489
490
491 for (auto pitr: *evt) {
492
493
494 const HepMC::FourVector pmom = pitr->momentum();
495 const int pstatus = pitr->status();
496 const int ppdgid = pitr->pdg_id();
497
498 if ( std::isnan(pmom.px()) || std::isinf(pmom.px()) ||
499 std::isnan(pmom.py()) || std::isinf(pmom.py()) ||
500 std::isnan(pmom.pz()) || std::isinf(pmom.pz()) ||
501 std::isnan(pmom.e()) || std::isinf(pmom.e()) ) {
502 ATH_MSG_WARNING(
"NaN (Not A Number) or inf found in the event record momenta");
504
507 filter_pass = false;
508 }
509 }
510
511
513 if (ppdgid == 111 && !pitr->end_vertex() ) {
514 unDecPi0.push_back( pitr);
516 }
517 }
518
519
522 if (pd != NULL) {
523 double plifetime =
pd->lifetime()*1
e+12;
524 if (plifetime != 0 && plifetime <
m_min_tau) {
525 ATH_MSG_WARNING(
"Stable particle found with lifetime = " << plifetime <<
"~ns!!");
527
529
531 filter_pass = false;
532 }
533 }
534 }
535 else{
536 int susyPart = 0;
537 std::vector<int>::size_type
count = 0;
539
541 susyPart=1;
542 }
544 }
545 if (susyPart==0){
546 ATH_MSG_WARNING(
"Stable particle not found in PDT, no lifetime check done");
548 }
549 }
550 }
551
552
553 const MC::DecodedPID decodedPID(ppdgid);
554 const int first_dig = decodedPID(0);
555
557
558 int known_byG4 = 0;
559 std::vector<int>::size_type
count =0;
560
564 }
565 if(known_byG4==0){
566 nonG4_energy += pmom.e();
567 ATH_MSG_WARNING(
"Interacting particle not known by Geant4 with ID " << ppdgid);
568 }
569 }
570
571
576 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#ifdef HEPMC3
613 for (auto desc: HepMC::descendant_particles(vtx)) {
614#else
615 for (auto desc_it = vtx->particles_begin(HepMC::descendants); desc_it != vtx->particles_end(HepMC::descendants); ++desc_it) {
616 auto desc=(*desc_it);
617#endif
618 if (std::abs(
desc->pdg_id()) ==
m_pdg) tau_child = 1;
620 }
621 if (std::abs( p_energy - pmom.e()) >
m_energy_diff && !tau_child) {
623 <<
"Energy (original particle) > " <<
m_energy_diff <<
" MeV, "
624 <<
"Event #" <<
evt->event_number() <<
", "
625 << "The original particle = " << pitr);
628 }
629
630 const HepMC::FourVector tau_decaypos = vtx->position();
631 const double tau_displacement = tau_decaypos.x()*tau_decaypos.x() + tau_decaypos.y()*tau_decaypos.y() + tau_decaypos.z()*tau_decaypos.z();
632
634 } else {
638 }
639 }
640
641
642 if (pitr->end_vertex()) {
643 auto decayvtx = pitr->end_vertex();
644 const HepMC::FourVector decaypos = decayvtx->position();
645 const double displacement = decaypos.x()*decaypos.x() + decaypos.y()*decaypos.y() + decaypos.z()*decaypos.z();
646 if (displacement > 1e-6) {
647 for (auto ip: *decayvtx) {
648 const HepMC::FourVector pos2 =
ip->production_vertex()->position();
649 const double displacement2 = pos2.x()*pos2.x() + pos2.y()*pos2.y() + pos2.z()*pos2.z();
650 if (displacement2 < 1e-6) {
652 <<
" has undisplaced vertex (" <<
ip->production_vertex()
653 << " @ " << displacement2 << "mm) "
654 << " but parent vertex is displaced (" << decayvtx
655 << " @ " << displacement << "mm)");
656 undisplaceds.push_back(std::move(ip));
658 }
659 }
660 }
661 }
662
663
666 const double mass = pitr->generated_mass();
667 if (std::abs(mass) > 1.0) {
668 ATH_MSG_WARNING(
"Photon with non-zero mass found! Mass: " << mass <<
" MeV" << pitr);
670 }
671 }
672
673 }
674
675
677 ATH_MSG_WARNING(
"The energy of interacting particles not known by Geant4 is = " << nonG4_energy <<
" MeV");
679 filter_pass = false;
680 }
682 }
683
684
685 double lostE = std::abs(totalE - cmenergy);
687 ATH_MSG_WARNING(
"ENERGY BALANCE FAILED : E-difference = " << lostE <<
" MeV");
688
689 ATH_MSG_WARNING(
"balance " << totalPx <<
" " << totalPy <<
" " << totalPz <<
" " << totalE);
690
693 }
696 filter_pass = false;
697 }
699 }
700
701
703 ATH_MSG_WARNING(
"MOMENTUM BALANCE FAILED : SumPx = " << totalPx <<
" SumPy = " << totalPy <<
" SumPz = " << totalPz <<
" MeV");
708 }
711 filter_pass = false;
712 }
714 }
715
716
717 if (!negEnPart.empty()) {
718 std::stringstream
ss;
719 ss <<
"NEGATIVE ENERGY PARTICLES FOUND :";
720 for (const auto &b: negEnPart){
722 }
726 filter_pass = false;
727 }
729 }
730
731
732 if (!tachyons.empty()) {
733 std::stringstream
ss;
734 ss <<
"PARTICLES WITH |E| < |Pi| (i=x,y,z) FOUND :";
735 for (auto b: tachyons){
737 }
741 filter_pass = false;
742 }
744 }
745
746
747 if (!unstNoEnd.empty()) {
748 std::stringstream
ss;
749 ss <<
"Unstable particle with no decay vertex found: ";
750 for (auto b: unstNoEnd){
752 }
756 filter_pass = false;
757 }
759 }
760
761
762 if (!unDecPi0.empty()) {
763 std::stringstream
ss;
764 ss <<
"pi0 with no decay vertex found:";
765 for (auto b: unDecPi0){
767 }
771 filter_pass = false;
772 }
774 }
775
776
777 if (!undisplaceds.empty()) {
778 std::stringstream
ss{
"Undisplaced decay vertices from displaced particle: "};
780
782 }
786 filter_pass = false;
787 }
789 }
790
791 }
792
793
794 if (!filter_pass){
795 setFilterPassed(false);
797 } else {
799 }
800
801
802
806 return StatusCode::FAILURE;
807 }
808
809 return StatusCode::SUCCESS;
810}
#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