7 #include "GaudiKernel/DataSvc.h"
8 #include "GaudiKernel/SystemOfUnits.h"
19 m_thistSvc(
"THistSvc",
name)
188 ATH_MSG_INFO(
"No decay vertex - is ignored for particles with status (list):" );
197 std::ifstream G4file;
198 G4file.open(fileLocation);
202 while(std::getline(G4file,
line)){
203 std::stringstream
ss(
line);
210 ATH_MSG_WARNING(
"Failed to open G4particle_acceptlist.txt, checking that all particles are known by Genat4 cannot be performed");
217 while(std::getline(G4file,
line)){
218 std::stringstream
ss(
line);
229 std::ifstream susyFile;
230 susyFile.open(
"susyParticlePdgid.txt");
232 if (!susyFile.fail()){
233 while(getline(susyFile,
line)){
234 std::stringstream ss1(
line);
241 ATH_MSG_WARNING(
"Failed to open susyParticlePdgid.txt, listing particles not present in PDTTable");
245 std::ifstream pdgFile;
248 if (!pdgFile.fail()){
250 while(std::getline(pdgFile,
line)){
251 std::stringstream
ss(
line);
258 ATH_MSG_INFO(
"extra accept list for PDG IDs not provided");
265 return StatusCode::SUCCESS;
272 bool filter_pass =
true;
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;
291 ATH_MSG_DEBUG(
"Please use MC::Loops::findLoops for this event to obtain all particles and vertices in the loops");
296 const auto xsec =
evt->cross_section();
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.");
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));
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(
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);
321 ATH_MSG_WARNING(
"Invalid number of beam particles " << beams_t.size() <<
" this generator interface should be fixed");
326 auto beams =
evt->beam_particles();
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");
335 ATH_MSG_WARNING(
"Beam particles have incorrectly set status -- this generator interface should be fixed");
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);
342 if(beams.first->pdg_id() == MC::OXYGEN && beams.second->pdg_id() == MC::OXYGEN){
345 if(beams.first->pdg_id() == MC::LEAD && beams.second->pdg_id() == MC::LEAD){
348 if(beams.first->pdg_id() == MC::OXYGEN && beams.second->pdg_id() == MC::PROTON){
351 if(beams.first->pdg_id() == MC::PROTON && beams.second->pdg_id() == MC::OXYGEN){
354 if(beams.first->pdg_id() == MC::LEAD && beams.second->pdg_id() == MC::PROTON){
357 if(beams.first->pdg_id() == MC::PROTON && beams.second->pdg_id() == MC::LEAD){
360 if(beams.first->pdg_id() == MC::OXYGEN && beams.second->pdg_id() == MC::HELIUM){
363 if(beams.first->pdg_id() == MC::HELIUM && beams.second->pdg_id() == MC::OXYGEN){
369 setFilterPassed(
false);
377 return StatusCode::FAILURE;
382 int vtxDisplacedstatuscode12CheckRateCnt=0;
383 int vtxDisplacedstatuscodenot12CheckRateCnt=0;
384 int vtxDisplacedMoreThan_1m_CheckRateCnt=0;
386 for (
const auto& vtx:
evt->vertices()) {
388 for (
auto vitr =
evt->vertices_begin(); vitr !=
evt->vertices_end(); ++vitr ) {
389 const HepMC::GenVertex* vtx = *vitr;
391 const HepMC::FourVector
pos = vtx->position();
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");
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;
424 for (
const auto&
part: vtx->particles_in()) {
426 for (
auto part_it = vtx->particles_in_const_begin(); part_it != vtx->particles_in_const_end(); ++part_it) {
427 auto part=(*part_it);
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()) {
438 for(
const auto& p_parents:
part->production_vertex()->particles_in()) {
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);
444 msg(MSG::WARNING) <<
"\t";
450 if (
part->status() == 1 ||
part->status() == 2){
451 vtxDisplacedstatuscode12CheckRateCnt += 1;
453 vtxDisplacedstatuscodenot12CheckRateCnt += 1;
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);
491 for (
auto pitr: *
evt) {
494 const HepMC::FourVector pmom = pitr->momentum();
495 const int pstatus = pitr->status();
496 const int ppdgid = pitr->pdg_id();
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");
513 if (ppdgid == 111 && !pitr->end_vertex() ) {
514 unDecPi0.push_back( pitr);
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!!");
537 std::vector<int>::size_type
count = 0;
546 ATH_MSG_WARNING(
"Stable particle not found in PDT, no lifetime check done");
553 const MC::DecodedPID decodedPID(ppdgid);
554 const int first_dig = decodedPID(0);
559 std::vector<int>::size_type
count =0;
566 nonG4_energy += pmom.e();
567 ATH_MSG_WARNING(
"Interacting particle not known by Geant4 with ID " << ppdgid);
582 unstNoEnd.push_back(pitr);
589 totalPx += pmom.px();
590 totalPy += pmom.py();
591 totalPz += pmom.pz();
594 negEnPart.push_back(pitr);
597 const double aener = std::abs(pmom.e());
599 tachyons.push_back(pitr);
609 auto vtx = pitr->end_vertex();
613 for (
auto desc: HepMC::descendant_particles(vtx)) {
615 for (
auto desc_it = vtx->particles_begin(HepMC::descendants); desc_it != vtx->particles_end(HepMC::descendants); ++desc_it) {
616 auto desc=(*desc_it);
618 if (std::abs(
desc->pdg_id()) ==
m_pdg) tau_child = 1;
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);
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();
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 > 1
e-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 < 1
e-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(
ip);
666 const double mass = pitr->generated_mass();
667 if (std::abs(
mass) > 1.0) {
677 ATH_MSG_WARNING(
"The energy of interacting particles not known by Geant4 is = " << nonG4_energy <<
" MeV");
685 double lostE = std::abs(totalE - cmenergy);
687 ATH_MSG_WARNING(
"ENERGY BALANCE FAILED : E-difference = " << lostE <<
" MeV");
689 ATH_MSG_WARNING(
"balance " << totalPx <<
" " << totalPy <<
" " << totalPz <<
" " << totalE);
703 ATH_MSG_WARNING(
"MOMENTUM BALANCE FAILED : SumPx = " << totalPx <<
" SumPy = " << totalPy <<
" SumPz = " << totalPz <<
" MeV");
717 if (!negEnPart.empty()) {
718 std::stringstream
ss;
719 ss <<
"NEGATIVE ENERGY PARTICLES FOUND :";
720 for (
const auto &
b: negEnPart){
732 if (!tachyons.empty()) {
733 std::stringstream
ss;
734 ss <<
"PARTICLES WITH |E| < |Pi| (i=x,y,z) FOUND :";
735 for (
auto b: tachyons){
747 if (!unstNoEnd.empty()) {
748 std::stringstream
ss;
749 ss <<
"Unstable particle with no decay vertex found: ";
750 for (
auto b: unstNoEnd){
762 if (!unDecPi0.empty()) {
763 std::stringstream
ss;
764 ss <<
"pi0 with no decay vertex found:";
765 for (
auto b: unDecPi0){
777 if (!undisplaceds.empty()) {
778 std::stringstream
ss;
779 ss <<
"Undisplaced decay vertices from displaced particle: ";
780 for (
auto b: undisplaceds){
795 setFilterPassed(
false);
806 return StatusCode::FAILURE;
809 return StatusCode::SUCCESS;
824 if (!
m_vtxNaNTest)
ATH_MSG_INFO(
" The check for NaN or inf in vtx. record is switched off, so is not included in the final TestHepMC efficiency ");
827 if (!
m_vtxDisplacedTest)
ATH_MSG_INFO(
" The check for displaced vertices is switched off, so is not included in the final TestHepMC efficiency ");
829 if (!
m_momNaNTest)
ATH_MSG_INFO(
" The check for NaN/inf in momentum record is switched off, so is not included in the final TestHepMC efficiency ");
839 if (!
m_momImbalanceTest)
ATH_MSG_INFO(
" The check for momentum imbalance is switched off, so is not included in the final TestHepMC efficiency ");
841 if (!
m_negativeEnergyTest)
ATH_MSG_INFO(
" The check for particles with negative energy is switched off, so is not included in the final TestHepMC efficiency ");
843 if (!
m_tachyonsTest)
ATH_MSG_INFO(
" The check for tachyons is switched off, so is not included in the final TestHepMC efficiency ");
846 if (!
m_unstableNoVtxTest)
ATH_MSG_INFO(
" The check for unstable part. without end vertex is switched off, so is not included in the final TestHepMC efficiency ");
848 if (!
m_pi0NoVtxTest)
ATH_MSG_INFO(
" The check for undecayed pi0's is switched off, so is not included in the final TestHepMC efficiency ");
852 if (!
m_lifeTimeTest)
ATH_MSG_INFO(
" The check for status 1 particles with too short lifetime is switched off, so is not included in the final TestHepMC efficiency ");
854 if (!
m_energyG4Test)
ATH_MSG_INFO(
" The check for energy not known by G4 is switched off, so is not included in the final TestHepMC efficiency ");
861 return StatusCode::FAILURE;
870 return StatusCode::FAILURE;
880 return StatusCode::FAILURE;
885 return StatusCode::SUCCESS;