diff --git a/Simulation/Tools/McEventCollectionFilter/src/TruthClosureCheck.cxx b/Simulation/Tools/McEventCollectionFilter/src/TruthClosureCheck.cxx index 3d8d00a0b9ba398dd614d507449c45b105d0594e..4d7f27bf4e5b0e89ac55e50fb5fda39d10f2d8bb 100644 --- a/Simulation/Tools/McEventCollectionFilter/src/TruthClosureCheck.cxx +++ b/Simulation/Tools/McEventCollectionFilter/src/TruthClosureCheck.cxx @@ -1,11 +1,11 @@ /* - Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration */ #include "TruthClosureCheck.h" // -#include "HepMC/GenEvent.h" -#include "HepMC/GenVertex.h" +#include "AtlasHepMC/GenEvent.h" +#include "AtlasHepMC/GenVertex.h" #include "GeneratorObjects/McEventCollection.h" #include "StoreGate/ReadHandle.h" @@ -40,32 +40,28 @@ StatusCode TruthClosureCheck::initialize() { StatusCode TruthClosureCheck::sanityCheck(const HepMC::GenEvent& event) const { //Sanity check - HepMC::GenEvent::particle_const_iterator partIter = event.particles_begin(); - const HepMC::GenEvent::particle_const_iterator partEnd = event.particles_end(); bool resetProblem(false); - while (partIter!=partEnd) { - const HepMC::GenParticle& particle(**partIter); - if (particle.status()==1) { - if (!particle.production_vertex()) { + for (auto particle: event) { + if (particle->status() == 1) { + if (!particle->production_vertex()) { ATH_MSG_ERROR("Status 1 particle without a production vertex!! " << particle); resetProblem = true; } - if (!m_postSimulation && particle.end_vertex()) { + if (!m_postSimulation && particle->end_vertex()) { ATH_MSG_ERROR("Status 1 particle with an end vertex!! " << particle); resetProblem = true; } } - else if (particle.status()==2) { - if (!particle.production_vertex()) { + else if (particle->status() == 2) { + if (!particle->production_vertex()) { ATH_MSG_ERROR("Status 2 particle without a production vertex!! " << particle); resetProblem = true; } - if (!particle.end_vertex()) { + if (!particle->end_vertex()) { ATH_MSG_ERROR("Status 2 particle without an end vertex!! " << particle); resetProblem = true; } } - ++partIter; } if (resetProblem) { return StatusCode::FAILURE; @@ -74,6 +70,28 @@ StatusCode TruthClosureCheck::sanityCheck(const HepMC::GenEvent& event) const { return StatusCode::SUCCESS; } +#ifdef HEPMC3 +void TruthClosureCheck::printGenVertex(HepMC::ConstGenVertexPtr origVertex, + HepMC::ConstGenVertexPtr resetVertex) const +{ + ATH_MSG_INFO("----------------------------------"); + ATH_MSG_INFO("Original Vertex:"); + ATH_MSG_INFO( origVertex ); + ATH_MSG_INFO("Particles In:"); + for (auto originalPartIn: origVertex->particles_in()) ATH_MSG_INFO( originalPartIn ); + ATH_MSG_INFO("Particles Out:"); + for (auto originalPartOut: origVertex->particles_out()) ATH_MSG_INFO( originalPartOut ); + ATH_MSG_INFO("----------------------------------"); + ATH_MSG_INFO("Reset Vertex:"); + ATH_MSG_INFO( resetVertex ); + ATH_MSG_INFO("Particles In:"); + for (auto resetPartIn: resetVertex->particles_in()) ATH_MSG_INFO( resetPartIn ); + ATH_MSG_INFO("Particles Out:"); + for (auto resetPartOut: resetVertex->particles_out()) ATH_MSG_INFO( resetPartOut ); + ATH_MSG_INFO("----------------------------------"); + return; +} +#else void TruthClosureCheck::printGenVertex(const HepMC::GenVertex& origVertex, const HepMC::GenVertex& resetVertex) const @@ -115,7 +133,69 @@ void TruthClosureCheck::printGenVertex(const HepMC::GenVertex& origVertex, ATH_MSG_INFO("----------------------------------"); return; } +#endif +#ifdef HEPMC3 +StatusCode TruthClosureCheck::compareGenVertex(HepMC::ConstGenVertexPtr origVertex, + HepMC::ConstGenVertexPtr resetVertex) const +{ + if (!origVertex && !resetVertex) return StatusCode::SUCCESS; + if (!origVertex || !resetVertex) return StatusCode::FAILURE; + bool pass{true}; + + if (HepMC::barcode(origVertex) != HepMC::barcode(resetVertex)) { + ATH_MSG_ERROR ("vertex barcode differs! Original: "<<HepMC::barcode(origVertex)<<", Reset: "<<HepMC::barcode(resetVertex)); + pass = false; + } + if (origVertex->position() != resetVertex->position()) { + const HepMC::FourVector& oP = origVertex->position(); + const HepMC::FourVector& rP = resetVertex->position(); + ATH_MSG_ERROR("vertex position differs! Original: ("<<oP.x()<<","<<oP.y()<<","<<oP.z()<<","<<oP.t()<<"), Reset: ("<<rP.x()<<","<<rP.y()<<","<<rP.z()<<","<<rP.t()<<")"); + pass = false; + } + if (origVertex->particles_in().size() != resetVertex->particles_in().size() ) { + ATH_MSG_ERROR ("particles_in_size differs! Original: "<<origVertex->particles_in().size()<<", Reset: "<<resetVertex->particles_in().size()); + pass = false; + } + if (origVertex->particles_out().size() != resetVertex->particles_out().size() ) { + ATH_MSG_ERROR ("particles_out_size differs! Original: "<<origVertex->particles_out().size()<<", Reset: "<<resetVertex->particles_out().size()); + pass = false; + } + std::vector<HepMC::ConstGenParticlePtr> OriginalListOfParticlesIn = origVertex->particles_in(); + std::vector<HepMC::ConstGenParticlePtr> ResetListOfParticlesIn = resetVertex->particles_in(); + for ( std::vector<HepMC::ConstGenParticlePtr>::iterator originalPartInIter(OriginalListOfParticlesIn.begin()),resetPartInIter(ResetListOfParticlesIn.begin()); + originalPartInIter != OriginalListOfParticlesIn.end() && resetPartInIter != ResetListOfParticlesIn.end(); + ++originalPartInIter, ++resetPartInIter + ) { + if (compareGenParticle(*originalPartInIter,*resetPartInIter).isFailure()) { + ATH_MSG_ERROR ( "input particle properties differ!" ); + pass &= false; + } + ATH_MSG_VERBOSE ( "particles match!" ); + } + + + std::vector<HepMC::ConstGenParticlePtr> OriginalListOfParticlesOut = origVertex->particles_out(); + std::vector<HepMC::ConstGenParticlePtr> ResetListOfParticlesOut = resetVertex->particles_out(); + std::sort(OriginalListOfParticlesOut.begin(), OriginalListOfParticlesOut.end(), [](auto& a, auto& b) -> bool{return HepMC::barcode(a) > HepMC::barcode(b); }); + std::sort(ResetListOfParticlesOut.begin(), ResetListOfParticlesOut.end(), [](auto& a, auto& b) -> bool{return HepMC::barcode(a) > HepMC::barcode(b); }); + + for ( std::vector<HepMC::ConstGenParticlePtr>::iterator originalPartOutIter(OriginalListOfParticlesOut.begin()), resetPartOutIter(ResetListOfParticlesOut.begin()); + originalPartOutIter != OriginalListOfParticlesOut.end() && resetPartOutIter != ResetListOfParticlesOut.end(); + ++originalPartOutIter, ++resetPartOutIter + ) { + if (compareGenParticle(*originalPartOutIter,*resetPartOutIter).isFailure()) { + ATH_MSG_ERROR ( "output particle properties differ!" ); + pass &= false; + } + ATH_MSG_VERBOSE ( "particles match!" ); + } + if(!pass) { return StatusCode::FAILURE; } + return StatusCode::SUCCESS; +} + + +#else StatusCode TruthClosureCheck::compareGenVertex(const HepMC::GenVertex& origVertex, const HepMC::GenVertex& resetVertex) const { @@ -181,6 +261,7 @@ StatusCode TruthClosureCheck::compareGenVertex(const HepMC::GenVertex& origVerte if(!pass) { return StatusCode::FAILURE; } return StatusCode::SUCCESS; } +#endif StatusCode TruthClosureCheck::compareMomenta(const HepMC::FourVector& origMomenta, @@ -209,6 +290,37 @@ StatusCode TruthClosureCheck::compareMomenta(const HepMC::FourVector& origMoment return StatusCode::SUCCESS; } +#ifdef HEPMC3 +StatusCode TruthClosureCheck::compareGenParticle(HepMC::ConstGenParticlePtr origParticle, + HepMC::ConstGenParticlePtr resetParticle) const +{ + if (!origParticle && !resetParticle) return StatusCode::SUCCESS; + if (!origParticle || !resetParticle) return StatusCode::FAILURE; + + bool pass{true}; + if (HepMC::barcode(origParticle) != HepMC::barcode(resetParticle)) { + ATH_MSG_ERROR ("particle barcode differs! Original: "<<HepMC::barcode(origParticle)<<", Reset: "<<HepMC::barcode(resetParticle)); + pass &= false; + } + if (origParticle->status() != resetParticle->status()) { + ATH_MSG_ERROR ("particle status differs! Original: "<<origParticle->status()<<", Reset: "<<resetParticle->status()); + pass &= false; + } + if (origParticle->pdg_id() != resetParticle->pdg_id()) { + ATH_MSG_ERROR ("particle pdg_id differs! Original: "<<origParticle->pdg_id()<<", Reset: "<<resetParticle->pdg_id()); + pass &= false; + } + if (compareMomenta(origParticle->momentum(), resetParticle->momentum()).isFailure()) { + const HepMC::FourVector& oM = origParticle->momentum(); + const HepMC::FourVector& rM = resetParticle->momentum(); + ATH_MSG_DEBUG("particle momentum differs! Original: ("<<oM.x()<<","<<oM.y()<<","<<oM.z()<<","<<oM.t()<<"), Reset: ("<<rM.x()<<","<<rM.y()<<","<<rM.z()<<","<<rM.t()<<")"); + pass &= false; + } + if(!pass) { return StatusCode::FAILURE; } + return StatusCode::SUCCESS; +} + +#else StatusCode TruthClosureCheck::compareGenParticle(const HepMC::GenParticle& origParticle, const HepMC::GenParticle& resetParticle) const { @@ -234,6 +346,86 @@ StatusCode TruthClosureCheck::compareGenParticle(const HepMC::GenParticle& origP if(!pass) { return StatusCode::FAILURE; } return StatusCode::SUCCESS; } +#endif + +#ifdef HEPMC3 +//------------------------------------------------- +StatusCode TruthClosureCheck::execute() { +//------------------------------------------------- + + ATH_MSG_DEBUG( " execute..... " ); + SG::ReadHandle<McEventCollection> originalMcEventCollection(m_originalMcEventCollection); + if (!originalMcEventCollection.isValid()) { + ATH_MSG_ERROR("Could not find original McEventCollection called " << originalMcEventCollection.name() << " in store " << originalMcEventCollection.store() << "."); + return StatusCode::FAILURE; + } + const HepMC::GenEvent& originalEvent(**(originalMcEventCollection->begin())); + + if (sanityCheck(originalEvent).isFailure()) { + ATH_MSG_FATAL("Problems in original GenEvent - bailing out."); + return StatusCode::FAILURE; + } + + SG::ReadHandle<McEventCollection> resetMcEventCollection(m_resetMcEventCollection); + if (!resetMcEventCollection.isValid()) { + ATH_MSG_ERROR("Could not find reset McEventCollection called " << resetMcEventCollection.name() << " in store " << resetMcEventCollection.store() << "."); + return StatusCode::FAILURE; + } + const HepMC::GenEvent& resetEvent(**(resetMcEventCollection->begin())); + + if (sanityCheck(resetEvent).isFailure()) { + ATH_MSG_FATAL("Problems in reset GenEvent - bailing out."); + return StatusCode::FAILURE; + } + + if (originalEvent.event_number() != resetEvent.event_number() ) { + ATH_MSG_ERROR ("event_number differs! Original: "<<originalEvent.event_number()<<", Reset: "<<resetEvent.event_number()); + } + + if (HepMC::signal_process_id(originalEvent) != HepMC::signal_process_id(resetEvent) ) { + ATH_MSG_ERROR ("signal_process_id differs! Original: "<<HepMC::signal_process_id(originalEvent)<<", Reset: "<<HepMC::signal_process_id(resetEvent)); + } +///-> + std::vector<HepMC::ConstGenParticlePtr> OriginalBeams = originalEvent.beams(); + std::vector<HepMC::ConstGenParticlePtr> ResetBeams = resetEvent.beams(); + + for ( std::vector<HepMC::ConstGenParticlePtr>::iterator originalBeamIter(OriginalBeams.begin()), resetBeamIter(ResetBeams.begin()); + originalBeamIter != OriginalBeams.end() && resetBeamIter != ResetBeams.end(); + ++originalBeamIter, ++resetBeamIter + ) { + if (compareGenParticle(*originalBeamIter,*resetBeamIter).isFailure()) { + ATH_MSG_ERROR ( "Beam particle properties differ!" ); + } +///<= + if (originalEvent.particles().size() != resetEvent.particles().size() ) { + ATH_MSG_ERROR ("particles_size differs! Original: "<<originalEvent.particles().size()<<", Reset: "<<resetEvent.particles().size()); + } + + if (originalEvent.vertices().size() != resetEvent.vertices().size() ) { + ATH_MSG_ERROR ("vertices_size differs! Original: "<<originalEvent.vertices().size()<<", Reset: "<<resetEvent.vertices().size()); + } + } + + //loop over vertices in Background GenEvent + std::vector<HepMC::ConstGenVertexPtr> OriginalListOfVertices = originalEvent.vertices(); + std::vector<HepMC::ConstGenVertexPtr> ResetListOfVertices = resetEvent.vertices(); + + for( std::vector<HepMC::ConstGenVertexPtr>::iterator origVertexIter(OriginalListOfVertices.begin()),resetVertexIter(ResetListOfVertices.begin()); + origVertexIter != OriginalListOfVertices.end() && resetVertexIter != ResetListOfVertices.end(); + ++origVertexIter, ++resetVertexIter + ) { + if (compareGenVertex(*origVertexIter,*resetVertexIter).isFailure()) { + ATH_MSG_ERROR ( "vertices differ!" ); + printGenVertex(*origVertexIter,*resetVertexIter); + } + } + ATH_MSG_INFO( "Completed Truth reset closure check ..... " ); + + return StatusCode::SUCCESS; + +} + +#else //------------------------------------------------- @@ -318,3 +510,4 @@ StatusCode TruthClosureCheck::execute() { return StatusCode::SUCCESS; } +#endif diff --git a/Simulation/Tools/McEventCollectionFilter/src/TruthClosureCheck.h b/Simulation/Tools/McEventCollectionFilter/src/TruthClosureCheck.h index ed51fe64686a888afd6a59482cfe9c80c229b1f6..653dcdf56eab0bf3bd3ea0ec63002d67088fbb1d 100644 --- a/Simulation/Tools/McEventCollectionFilter/src/TruthClosureCheck.h +++ b/Simulation/Tools/McEventCollectionFilter/src/TruthClosureCheck.h @@ -22,14 +22,23 @@ public: private: StatusCode sanityCheck(const HepMC::GenEvent& event) const; +#ifdef HEPMC3 + StatusCode compareGenVertex(HepMC::ConstGenVertexPtr origVertex, + HepMC::ConstGenVertexPtr resetVertex) const; + StatusCode compareGenParticle(HepMC::ConstGenParticlePtr origParticle, + HepMC::ConstGenParticlePtr resetParticle) const; + void printGenVertex(HepMC::ConstGenVertexPtr origVertex, + HepMC::ConstGenVertexPtr resetVertex) const; +#else StatusCode compareGenVertex(const HepMC::GenVertex& origVertex, const HepMC::GenVertex& resetVertex) const; StatusCode compareGenParticle(const HepMC::GenParticle& origParticle, const HepMC::GenParticle& resetParticle) const; - StatusCode compareMomenta(const HepMC::FourVector& origMomenta, - const HepMC::FourVector& resetMomenta) const; void printGenVertex(const HepMC::GenVertex& origVertex, const HepMC::GenVertex& resetVertex) const; +#endif + StatusCode compareMomenta(const HepMC::FourVector& origMomenta, + const HepMC::FourVector& resetMomenta) const; SG::ReadHandleKey<McEventCollection> m_originalMcEventCollection; SG::ReadHandleKey<McEventCollection> m_resetMcEventCollection;