271 bool filter_pass =
true;
279 double nonG4_energy = 0;
280 std::vector<HepMC::ConstGenParticlePtr> negEnPart;
281 std::vector<HepMC::ConstGenParticlePtr> tachyons;
282 std::vector<HepMC::ConstGenParticlePtr> unstNoEnd;
283 std::vector<HepMC::ConstGenParticlePtr> unDecPi0;
284 std::vector<HepMC::ConstGenParticlePtr> undisplaceds;
290 ATH_MSG_DEBUG(
"Please use MC::Loops::findLoops for this event to obtain all particles and vertices in the loops");
295 const auto xsec =
evt->cross_section();
297 ATH_MSG_WARNING(
"WATCH OUT: event is missing the generator cross-section!");
300 ATH_MSG_WARNING(
"-> Adding a dummy cross-section for debugging purposes.");
302 std::shared_ptr<HepMC3::GenCrossSection> dummy_xsec = std::make_shared<HepMC3::GenCrossSection>();
303 dummy_xsec->set_cross_section(1.0,0.0);
304 HepMC::GenEvent* evt_nonconst =
const_cast<HepMC::GenEvent*
>(
evt);
305 evt_nonconst->set_cross_section(dummy_xsec);
313 std::vector<std::shared_ptr<const HepMC3::GenParticle>> beams_t;
314 for (
auto p :
evt->beams()) {
if (
p->status() == 4) beams_t.push_back(
p); }
315 std::pair<std::shared_ptr<const HepMC3::GenParticle>,std::shared_ptr<const HepMC3::GenParticle>> beams;
316 if (beams_t.size() == 2) {
317 beams.first=beams_t.at(0);
318 beams.second=beams_t.at(1);
320 ATH_MSG_WARNING(
"Invalid number of beam particles " << beams_t.size() <<
" this generator interface should be fixed");
325 auto beams =
evt->beam_particles();
329 ATH_MSG_WARNING(
"Invalid beam particle pointers -- this generator interface should be fixed");
330 if (cmenergy < 0)
ATH_MSG_WARNING(
"Invalid expected beam energy: " << cmenergy <<
" MeV");
333 if (beams.first->status() != 4 || beams.second->status() != 4) {
334 ATH_MSG_WARNING(
"Beam particles have incorrectly set status -- this generator interface should be fixed");
337 const double sumE = beams.first->momentum().e() + beams.second->momentum().e();
338 const double sumP = beams.first->momentum().pz() + beams.second->momentum().pz();
339 cmenergy = std::sqrt(sumE*sumE - sumP*sumP);
342 ATH_MSG_FATAL(
"Beam particles have incorrect energy: " <<
m_cm_energy/1000. <<
" GeV expected, vs. " << cmenergy/1000. <<
" GeV found");
343 setFilterPassed(
false);
351 return StatusCode::FAILURE;
356 int vtxDisplacedstatuscode12CheckRateCnt=0;
357 int vtxDisplacedstatuscodenot12CheckRateCnt=0;
358 int vtxDisplacedMoreThan_1m_CheckRateCnt=0;
360 for (
const auto& vtx:
evt->vertices()) {
362 for (
auto vitr =
evt->vertices_begin(); vitr !=
evt->vertices_end(); ++vitr ) {
363 const HepMC::GenVertex* vtx = *vitr;
365 const HepMC::FourVector
pos = vtx->position();
368 if ( std::isnan(
pos.x()) || std::isinf(
pos.x()) ||
369 std::isnan(
pos.y()) || std::isinf(
pos.y()) ||
370 std::isnan(
pos.z()) || std::isinf(
pos.z()) ) {
371 ATH_MSG_WARNING(
"NaN (Not A Number) or inf found in the event record vertex positions");
382 const double dist_trans2 =
pos.x()*
pos.x() +
pos.y()*
pos.y();
383 const double dist2 = dist_trans2 +
pos.z()*
pos.z();
384 const double dist_trans = std::sqrt(dist_trans2);
385 const double dist = std::sqrt(dist2);
388 ++vtxDisplacedMoreThan_1m_CheckRateCnt;
398 for (
const auto&
part: vtx->particles_in()) {
400 for (
auto part_it = vtx->particles_in_const_begin(); part_it != vtx->particles_in_const_end(); ++part_it) {
401 auto part=(*part_it);
407 ATH_MSG_WARNING(
"production vertex = " <<
part->production_vertex()->position().x() <<
", " <<
part->production_vertex()->position().y() <<
", " <<
part->production_vertex()->position().z());
408 ATH_MSG_WARNING(
"end vertex = " <<
part->end_vertex()->position().x() <<
", " <<
part->end_vertex()->position().y() <<
", " <<
part->end_vertex()->position().z());
410 if (
part->production_vertex()) {
412 for(
const auto& p_parents:
part->production_vertex()->particles_in()) {
414 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) {
415 auto p_parents=(*p_parents_it);
418 msg(MSG::WARNING) <<
"\t";
424 if (
part->status() == 1 ||
part->status() == 2){
425 vtxDisplacedstatuscode12CheckRateCnt += 1;
427 vtxDisplacedstatuscodenot12CheckRateCnt += 1;
432 if (
part->momentum().e()*1
e-3 < 10.) {
446 double endvx =
part->end_vertex()->position().x();
447 double endvy =
part->end_vertex()->position().y();
448 double endvz =
part->end_vertex()->position().z();
449 double prodvx =
part->production_vertex()->position().x();
450 double prodvy =
part->production_vertex()->position().y();
451 double prodvz =
part->production_vertex()->position().z();
452 double enddis = std::sqrt(endvx*endvx + endvy*endvy + endvz*endvz);
453 double proddis = std::sqrt(prodvx*prodvx + prodvy*prodvy + prodvz*prodvz);
465 for (
auto pitr: *
evt) {
468 const HepMC::FourVector pmom = pitr->momentum();
469 const int pstatus = pitr->status();
470 const int ppdgid = pitr->pdg_id();
472 if ( std::isnan(pmom.px()) || std::isinf(pmom.px()) ||
473 std::isnan(pmom.py()) || std::isinf(pmom.py()) ||
474 std::isnan(pmom.pz()) || std::isinf(pmom.pz()) ||
475 std::isnan(pmom.e()) || std::isinf(pmom.e()) ) {
476 ATH_MSG_WARNING(
"NaN (Not A Number) or inf found in the event record momenta");
486 if (pstatus == 1 || pstatus == 2) {
487 if (ppdgid == 111 && !pitr->end_vertex() ) {
488 unDecPi0.push_back( pitr);
497 double plifetime =
pd->lifetime()*1
e+12;
498 if (plifetime != 0 && plifetime <
m_min_tau) {
499 ATH_MSG_WARNING(
"Stable particle found with lifetime = " << plifetime <<
"~ns!!");
511 std::vector<int>::size_type
count = 0;
520 ATH_MSG_WARNING(
"Stable particle not found in PDT, no lifetime check done");
528 int first_dig = ppdgid;
529 while (first_dig > 9) first_dig /= 10;
531 if ((pstatus == 1 ) && (!pitr->end_vertex()) && (
MC::isSimInteracting(pitr)) && (!
pid.isNucleus()) && (first_dig != 9) ) {
534 std::vector<int>::size_type
count =0;
541 nonG4_energy += pmom.e();
542 ATH_MSG_WARNING(
"Interacting particle not known by Geant4 with ID " << ppdgid);
556 if (!pitr->end_vertex() && pstatus == 2) {
557 unstNoEnd.push_back(pitr);
563 if ( pstatus == 1 && !pitr->end_vertex() ) {
564 totalPx += pmom.px();
565 totalPy += pmom.py();
566 totalPz += pmom.pz();
569 negEnPart.push_back(pitr);
572 const double aener = std::abs(pmom.e());
574 tachyons.push_back(pitr);
582 if (std::abs(ppdgid) ==
m_pdg && (pstatus == 1 || pstatus == 2)) {
584 auto vtx = pitr->end_vertex();
588 for (
auto desc: HepMC::descendant_particles(vtx)) {
590 for (
auto desc_it = vtx->particles_begin(HepMC::descendants); desc_it != vtx->particles_end(HepMC::descendants); ++desc_it) {
591 auto desc=(*desc_it);
593 if (std::abs(
desc->pdg_id()) ==
m_pdg) tau_child = 1;
596 if (std::abs( p_energy - pmom.e()) >
m_energy_diff && !tau_child) {
598 <<
"Energy (original particle) > " <<
m_energy_diff <<
" MeV, "
599 <<
"Event #" <<
evt->event_number() <<
", "
600 <<
"The original particle = " << pitr);
605 const HepMC::FourVector tau_decaypos = vtx->position();
606 const double tau_displacement = tau_decaypos.x()*tau_decaypos.x() + tau_decaypos.y()*tau_decaypos.y() + tau_decaypos.z()*tau_decaypos.z();
617 if (pitr->end_vertex()) {
618 auto decayvtx = pitr->end_vertex();
619 const HepMC::FourVector decaypos = decayvtx->position();
620 const double displacement = decaypos.x()*decaypos.x() + decaypos.y()*decaypos.y() + decaypos.z()*decaypos.z();
621 if (displacement > 1
e-6) {
622 for (
auto ip: *decayvtx) {
623 const HepMC::FourVector pos2 =
ip->production_vertex()->position();
624 const double displacement2 = pos2.x()*pos2.x() + pos2.y()*pos2.y() + pos2.z()*pos2.z();
625 if (displacement2 < 1
e-6) {
627 <<
" has undisplaced vertex (" <<
ip->production_vertex()
628 <<
" @ " << displacement2 <<
"mm) "
629 <<
" but parent vertex is displaced (" << decayvtx
630 <<
" @ " << displacement <<
"mm)");
631 undisplaceds.push_back(
ip);
640 if (ppdgid == 22 && pstatus == 1) {
641 const double mass = pitr->generated_mass();
642 if (std::abs(
mass) > 1.0) {
652 ATH_MSG_WARNING(
"The energy of interacting particles not known by Geant4 is = " << nonG4_energy <<
" MeV");
660 double lostE = std::abs(totalE - cmenergy);
662 ATH_MSG_WARNING(
"ENERGY BALANCE FAILED : E-difference = " << lostE <<
" MeV");
664 ATH_MSG_WARNING(
"balance " << totalPx <<
" " << totalPy <<
" " << totalPz <<
" " << totalE);
678 ATH_MSG_WARNING(
"MOMENTUM BALANCE FAILED : SumPx = " << totalPx <<
" SumPy = " << totalPy <<
" SumPz = " << totalPz <<
" MeV");
692 if (!negEnPart.empty()) {
693 std::stringstream
ss;
694 ss <<
"NEGATIVE ENERGY PARTICLES FOUND :";
695 for (
auto b: negEnPart){
707 if (!tachyons.empty()) {
708 std::stringstream
ss;
709 ss <<
"PARTICLES WITH |E| < |Pi| (i=x,y,z) FOUND :";
710 for (
auto b: tachyons){
722 if (!unstNoEnd.empty()) {
723 std::stringstream
ss;
724 ss <<
"Unstable particle with no decay vertex found: ";
725 for (
auto b: unstNoEnd){
737 if (!unDecPi0.empty()) {
738 std::stringstream
ss;
739 ss <<
"pi0 with no decay vertex found:";
740 for (
auto b: unDecPi0){
752 if (!undisplaceds.empty()) {
753 std::stringstream
ss;
754 ss <<
"Undisplaced decay vertices from displaced particle: ";
755 for (
auto b: undisplaceds){
770 setFilterPassed(
false);
781 return StatusCode::FAILURE;
784 return StatusCode::SUCCESS;
1.8.18