diff --git a/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/CMakeLists.txt b/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/CMakeLists.txt index 6c6df917f1bf1fb03cb6f04062c573bc4a409175..4815a6961e86a46e2ec5ae30cb5b22e2445c9946 100644 --- a/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/CMakeLists.txt +++ b/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/CMakeLists.txt @@ -11,7 +11,7 @@ atlas_add_library( PhotonVertexSelectionLib PhotonVertexSelection/*.h Root/*.cxx PUBLIC_HEADERS PhotonVertexSelection INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} - LINK_LIBRARIES ${ROOT_LIBRARIES} AsgDataHandlesLib AsgTools xAODEgamma xAODEventInfo xAODMetaData + LINK_LIBRARIES ${ROOT_LIBRARIES} AsgDataHandlesLib AsgTools xAODEgamma xAODEventInfo xAODMetaData AnaAlgorithmLib PATCoreLib xAODTracking PRIVATE_LINK_LIBRARIES PathResolver egammaUtils) diff --git a/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/PhotonVertexSelection/BuildVertexPointingAlg.h b/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/PhotonVertexSelection/BuildVertexPointingAlg.h new file mode 100644 index 0000000000000000000000000000000000000000..c2744433f22c3d4b5a39c1aa7024815fcb5f21a2 --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/PhotonVertexSelection/BuildVertexPointingAlg.h @@ -0,0 +1,89 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef PHOTONVERTEXSELECTION_BUILDVERTEXPOINTINGALG_H +#define PHOTONVERTEXSELECTION_BUILDVERTEXPOINTINGALG_H + +#include <AnaAlgorithm/AnaAlgorithm.h> +#include <PATCore/IAsgSelectionTool.h> + +#include <xAODTracking/VertexContainerFwd.h> +#include <xAODEgamma/EgammaContainerFwd.h> +#include <xAODEventInfo/EventInfo.h> + +#include "AsgDataHandles/ReadHandleKey.h" +#include <AsgDataHandles/WriteHandleKey.h> +#include <AsgDataHandles/WriteDecorHandleKey.h> +#include <AsgDataHandles/ReadDecorHandleKey.h> +#include <AsgTools/AsgTool.h> +#include <AsgTools/PropertyWrapper.h> + +#include <string> + + +/** + * @brief An algorithm to build a vertex from photon pointing + * + * This algorithm is needed by @c VertexSelectionToolBDT. + * + */ +class BuildVertexPointingAlg : public EL::AnaAlgorithm { + public: + BuildVertexPointingAlg(const std::string& name, + ISvcLocator* svcLoc = nullptr); + virtual StatusCode initialize(); + virtual StatusCode execute(); + + private: + SG::ReadHandleKey<xAOD::EgammaContainer> m_photonContainerKey{ + this, "PhotonContainerKey", "Photons", + "The input photons (or electrons) container"}; + SG::ReadHandleKey<xAOD::EventInfo> m_eventInfoKey{this, "EventInfoKey", + "EventInfo", ""}; + SG::WriteHandleKey<xAOD::VertexContainer> m_pointingVertexContainerKey{ + this, "PointingVertexContainerKey", "PhotonPointingVertices", ""}; + + // additional information written to the vertex as decorations + SG::WriteDecorHandleKey<xAOD::VertexContainer> m_nphotons_good_key{ + this, "NGoodPhotonsKey", m_pointingVertexContainerKey, "nphotons_good", ""}; + SG::WriteDecorHandleKey<xAOD::VertexContainer> m_z_common_nphotons_key{ + this, "ZCommonPhotonsKey", m_pointingVertexContainerKey, "z_common_nphotons", ""}; + SG::WriteDecorHandleKey<xAOD::VertexContainer> m_photons_px_key { + this, "PhotonsPxKey", m_pointingVertexContainerKey, "photons_px", ""}; + SG::WriteDecorHandleKey<xAOD::VertexContainer> m_photons_py_key { + this, "PhotonsPyKey", m_pointingVertexContainerKey, "photons_py", ""}; + SG::WriteDecorHandleKey<xAOD::VertexContainer> m_photons_pz_key { + this, "PhotonsPzKey", m_pointingVertexContainerKey, "photons_pz", ""}; + SG::WriteDecorHandleKey<xAOD::VertexContainer> m_photon_links_key { + this, "PhotonLinksKey", m_pointingVertexContainerKey, "PhotonLinks", ""}; + + // decoration read from Photon by xAOD::PVHelpers::getZCommonAndError function + // the user should not change them + SG::ReadDecorHandleKey<xAOD::EgammaContainer> m_photon_zvertex_key{ + this, "PhotonZVertexKey", m_photonContainerKey, "zvertex"}; + SG::ReadDecorHandleKey<xAOD::EgammaContainer> m_photon_errz_key{ + this, "PhotonErrZKey", m_photonContainerKey, "errz"}; + SG::ReadDecorHandleKey<xAOD::EgammaContainer> m_photon_HPV_zvertex_key{ + this, "PhotonHPVZVertexKey", m_photonContainerKey, "HPV_zvertex"}; + SG::ReadDecorHandleKey<xAOD::EgammaContainer> m_photon_HPV_errz_key{ + this, "PhotonHPVErrZKey", m_photonContainerKey, "HPV_errz"}; + + ToolHandle<IAsgSelectionTool> m_selectionTool{ + this, "selectionTool", "CP::AsgPtEtaSelectionTool/PhotonSelectionTool", + "the selection tool we apply to photons to compute pointing variables"}; + + ToolHandle<IAsgSelectionTool> m_goodPhotonSelectionTool{ + this, "goodPhotonSelectionTool", "", + "the selection tool to count good photons"}; + + Gaudi::Property<float> m_convPtCut {this, "convPtCut", 2E3, "the conversion pt cut"}; + Gaudi::Property<int> m_nphotons_to_use {this, "nphotonsToUse", -1, "the number of photons to use for the pointing computation"}; + + void selectPhotons(const xAOD::EgammaContainer& original_photons, + ConstDataVector<DataVector<xAOD::Egamma>>& photons_selected) const; + + +}; + +#endif // PHOTONVERTEXSELECTION_BUILDVERTEXPOINTINGALG_H \ No newline at end of file diff --git a/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/PhotonVertexSelection/DecoratePhotonPointingAlg.h b/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/PhotonVertexSelection/DecoratePhotonPointingAlg.h new file mode 100644 index 0000000000000000000000000000000000000000..b462aaa05f614ceab3c3507c8eae47d1d1e1e32a --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/PhotonVertexSelection/DecoratePhotonPointingAlg.h @@ -0,0 +1,38 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef PHOTONVERTEXSELECTION_DECORATEPHOTONPOINTINGALG_H +#define PHOTONVERTEXSELECTION_DECORATEPHOTONPOINTINGALG_H + +#include <AnaAlgorithm/AnaAlgorithm.h> +#include <AsgDataHandles/ReadHandleKey.h> +#include <AsgTools/ToolHandle.h> +#include <PhotonVertexSelection/IPhotonPointingTool.h> +#include <xAODEgamma/EgammaContainerFwd.h> + +/** + * @brief An algorithm to decorate photons (also electrons) with pointing variables. + * + * This algorithm uses the PhotonPointingTool to decorate photons with pointing + * variables. + * + */ +class DecoratePhotonPointingAlg : public EL::AnaAlgorithm { + public: + DecoratePhotonPointingAlg(const std::string& name, + ISvcLocator* svcLoc = nullptr); + + virtual StatusCode initialize() override; + StatusCode execute() override; + + private: + SG::ReadHandleKey<xAOD::EgammaContainer> m_photonContainerKey{ + this, "PhotonContainerKey", "Photons", "The input Photons container (it can be also Electrons)"}; + + ToolHandle<CP::IPhotonPointingTool> m_pointingTool{ + this, "PhotonPointingTool", + "CP::PhotonPointingTool/PhotonPointingTool", ""}; +}; + +#endif // PHOTONVERTEXSELECTION_DECORATEPHOTONPOINTINGALG_H \ No newline at end of file diff --git a/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/PhotonVertexSelection/PhotonVertexSelectionDict.h b/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/PhotonVertexSelection/PhotonVertexSelectionDict.h index d0cdc15e935dd345a384d79c67eb0cc9748718b5..eae3ba6449c886afe402a357d188cfc6162a4461 100644 --- a/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/PhotonVertexSelection/PhotonVertexSelectionDict.h +++ b/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/PhotonVertexSelection/PhotonVertexSelectionDict.h @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration */ #ifndef PhotonVertexSelection_PhotonVertexSelectionDict_H @@ -10,6 +10,9 @@ #include "PhotonVertexSelection/IPhotonPointingTool.h" #include "PhotonVertexSelection/PhotonPointingTool.h" +#include "PhotonVertexSelection/DecoratePhotonPointingAlg.h" +// #include "PhotonVertexSelection/BuildVertexPointingAlg.h" + #include "PhotonVertexSelection/PhotonVertexHelpers.h" diff --git a/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/PhotonVertexSelection/selection.xml b/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/PhotonVertexSelection/selection.xml index 1cc3db3d6691dd608e9c374ccc318e8e6cd5843d..518084bf3a0dc87b2e3f3c3ca0e7eb5045bdec5a 100644 --- a/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/PhotonVertexSelection/selection.xml +++ b/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/PhotonVertexSelection/selection.xml @@ -5,7 +5,8 @@ <class name="CP::IPhotonPointingTool" /> <class name="CP::PhotonPointingTool" /> - + <class name="DecoratePhotonPointingAlg" /> + <function pattern="xAOD::PVHelpers::*" /> <class name="std::pair<const xAOD::Vertex*, float>" /> diff --git a/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/Root/DecoratePhotonPointingAlg.cxx b/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/Root/DecoratePhotonPointingAlg.cxx new file mode 100644 index 0000000000000000000000000000000000000000..4fa271d27ef3d948974be0c248de6a75ecb7c951 --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/Root/DecoratePhotonPointingAlg.cxx @@ -0,0 +1,31 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ + +#include "PhotonVertexSelection/DecoratePhotonPointingAlg.h" + +#include <xAODEgamma/EgammaContainer.h> + +#include "AsgDataHandles/ReadHandle.h" + +DecoratePhotonPointingAlg::DecoratePhotonPointingAlg(const std::string& name, + ISvcLocator* svcLoc) + : EL::AnaAlgorithm(name, svcLoc) {} + +StatusCode DecoratePhotonPointingAlg::initialize() { + ATH_CHECK(m_photonContainerKey.initialize()); + ATH_CHECK(m_pointingTool.retrieve()); + return StatusCode::SUCCESS; +} + +StatusCode DecoratePhotonPointingAlg::execute() { + SG::ReadHandle<xAOD::EgammaContainer> egammas{m_photonContainerKey}; + if (!egammas.isValid()) { + ATH_MSG_ERROR("Invalid egamma container from " << m_photonContainerKey.key()); + return StatusCode::FAILURE; + } + + ATH_CHECK(m_pointingTool->updatePointingAuxdata(*egammas)); + + return StatusCode::SUCCESS; +} \ No newline at end of file diff --git a/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/python/PhotonVertexSelectionConfig.py b/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/python/PhotonVertexSelectionConfig.py index 90078bc087c030b6b8141c319e5a6a66eaa323f1..682bc93ad2c3b91d0973493e2d8ca462436ad29c 100644 --- a/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/python/PhotonVertexSelectionConfig.py +++ b/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/python/PhotonVertexSelectionConfig.py @@ -1,16 +1,65 @@ -# Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration +# Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator from AthenaConfiguration.ComponentFactory import CompFactory + def PhotonPointingToolCfg(flags, name="PhotonPointingTool", **kwargs): acc = ComponentAccumulator() kwargs.setdefault("isSimulation", flags.Input.isMC) acc.setPrivateTools(CompFactory.CP.PhotonPointingTool(name, **kwargs)) return acc -def PhotonVertexSelectionToolCfg( - flags, name="PhotonVertexSelectionTool", **kwargs): + +def PhotonVertexSelectionToolCfg(flags, name="PhotonVertexSelectionTool", **kwargs): acc = ComponentAccumulator() - acc.setPrivateTools( CompFactory.CP.PhotonVertexSelectionTool( **kwargs ) ) + acc.setPrivateTools(CompFactory.CP.PhotonVertexSelectionTool(**kwargs)) return acc + + +def DecoratePhotonPointingAlgCfg(flags, name="DecoratePhotonPointingAlg", **kwargs): + acc = ComponentAccumulator() + + kwargs.setdefault("PhotonContainerKey", "Photons") + + if not kwargs.get("PhotonPointingTool", None): + toolAcc = PhotonPointingToolCfg(flags, ContainerName=kwargs['PhotonContainerKey']) + tool = toolAcc.popPrivateTools() + acc.merge(toolAcc) + kwargs["PhotonPointingTool"] = tool + + alg = CompFactory.DecoratePhotonPointingAlg("DecoratePhotonPointingAlg", **kwargs) + acc.addEventAlgo(alg) + return acc + + +def BuildVertexPointingAlgCfg(flags, name="BuildVertexPointingAlg", **kwargs): + from AthenaCommon.SystemOfUnits import GeV + + acc = ComponentAccumulator() + kwargs.setdefault( + "selectionTool", + CompFactory.CP.AsgPtEtaSelectionTool( + "CutPhotonPtEta", + minPt=20 * GeV, + maxEta=2.37, + etaGapLow=1.37, + etaGapHigh=1.52, + useClusterEta=True, + ), + ) + kwargs.setdefault( + "goodPhotonSelectionTool", + CompFactory.CP.AsgFlagSelectionTool( + "GoodPhotonSelectionTool", + selectionFlags=[ + "Tight" + ], # TODO: why not DFCommonPhotonsIsEMTight (not working)? + ), + ) + kwargs.setdefault("nphotonsToUse", 2) + kwargs.setdefault("PhotonContainerKey", "Photons") + alg = CompFactory.BuildVertexPointingAlg("BuildVertexPointingAlg", **kwargs) + acc.addEventAlgo(alg) + return acc + diff --git a/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/src/BuildVertexPointingAlg.cxx b/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/src/BuildVertexPointingAlg.cxx new file mode 100644 index 0000000000000000000000000000000000000000..c06bc846fbf1c0524b3b051f330c011feb23be05 --- /dev/null +++ b/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/src/BuildVertexPointingAlg.cxx @@ -0,0 +1,182 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ + +#include "PhotonVertexSelection/BuildVertexPointingAlg.h" + +#include <AsgDataHandles/ReadHandle.h> +#include <AsgDataHandles/WriteDecorHandle.h> +#include <AsgDataHandles/WriteHandle.h> +#include <AthContainers/ConstDataVector.h> +#include <EventPrimitives/EventPrimitives.h> +#include <TLorentzVector.h> +#include <xAODEgamma/EgammaContainer.h> +#include <xAODTracking/VertexAuxContainer.h> +#include <xAODTracking/VertexContainer.h> + +#include <numeric> +#include <utility> + +#include "PhotonVertexSelection/PhotonVertexHelpers.h" + +BuildVertexPointingAlg::BuildVertexPointingAlg(const std::string& name, + ISvcLocator* svcLoc) + : EL::AnaAlgorithm(name, svcLoc) {} + +StatusCode BuildVertexPointingAlg::initialize() { + // read + ATH_CHECK(m_photonContainerKey.initialize()); + ATH_CHECK(m_eventInfoKey.initialize()); + + ATH_CHECK(m_photon_zvertex_key.initialize()); + ATH_CHECK(m_photon_errz_key.initialize()); + ATH_CHECK(m_photon_HPV_zvertex_key.initialize()); + ATH_CHECK(m_photon_HPV_errz_key.initialize()); + + // write + ATH_CHECK(m_pointingVertexContainerKey.initialize()); + + ATH_CHECK(m_nphotons_good_key.initialize()); + ATH_CHECK(m_z_common_nphotons_key.initialize()); + ATH_CHECK(m_photons_px_key.initialize()); + ATH_CHECK(m_photons_py_key.initialize()); + ATH_CHECK(m_photons_pz_key.initialize()); + ATH_CHECK(m_photon_links_key.initialize()); + + // tool + ATH_CHECK(m_selectionTool.retrieve()); + + ATH_MSG_INFO("Use at maximum " << m_nphotons_to_use + << " photons for pointing"); + ATH_MSG_INFO("Selecting photons for pointing with " + << m_selectionTool << " with " + << m_selectionTool->getAcceptInfo().getNCuts() << "cuts:"); + for (unsigned int icut = 0; + icut < m_selectionTool->getAcceptInfo().getNCuts(); ++icut) { + ATH_MSG_INFO( + "cut " << icut << ": " + << m_selectionTool->getAcceptInfo().getCutDescription(icut)); + } + return StatusCode::SUCCESS; +} + +StatusCode BuildVertexPointingAlg::execute() { + // the code has been written with photons in mind, but it is important + // the code to work also with electrons since electrons are used as a + // proxy of photons in performance studies + SG::ReadHandle<xAOD::EgammaContainer> original_photons{m_photonContainerKey}; + ATH_CHECK(original_photons.isValid()); + + ATH_MSG_DEBUG("number of photon in the container (" + << m_photonContainerKey << "): " << original_photons->size()); + + // these photons will be used to compute pointing variables + ConstDataVector<DataVector<xAOD::Egamma>> photons_selected(SG::VIEW_ELEMENTS); + selectPhotons(*original_photons, photons_selected); + + const xAOD::EventInfo* ei = SG::get(m_eventInfoKey); + const std::pair<float, float> z_common_and_error = + xAOD::PVHelpers::getZCommonAndError(ei, photons_selected.asDataVector(), + m_convPtCut); + + // create new vertex container + auto pointingVertices = std::make_unique<xAOD::VertexContainer>(); + auto pointingVerticesAux = std::make_unique<xAOD::VertexAuxContainer>(); + pointingVertices->setStore(pointingVerticesAux.get()); + + // create new vertex + xAOD::Vertex* vtx_pointing = + (*pointingVertices).push_back(std::make_unique<xAOD::Vertex>()); + + vtx_pointing->setVertexType(xAOD::VxType::PriVtx); + + vtx_pointing->setZ(z_common_and_error.first); + vtx_pointing->setX(0.0); + vtx_pointing->setY(0.0); + AmgSymMatrix(3) cov; + cov.setZero(); + cov(2, 2) = z_common_and_error.second * z_common_and_error.second; + vtx_pointing->setCovariancePosition(cov); + + auto writeHandle = SG::makeHandle(m_pointingVertexContainerKey); + ANA_CHECK(writeHandle.record(std::move(pointingVertices), + std::move(pointingVerticesAux))); + + // below only decorations to the vertex + + const unsigned int n_good_photons = + std::count_if(photons_selected.begin(), photons_selected.end(), + [this](const xAOD::Egamma* photon) { + return m_goodPhotonSelectionTool->accept(photon); + }); + + std::vector<ElementLink<xAOD::EgammaContainer>> ph_links; + ph_links.reserve(photons_selected.size()); + for (const xAOD::Egamma* photon : photons_selected) { + ph_links.emplace_back(*original_photons, photon->index()); + } + + auto dec_nphotons = SG::makeHandle<unsigned int>(m_z_common_nphotons_key); + auto dec_ph_links = SG::makeHandle<decltype(ph_links)>(m_photon_links_key); + auto dec_nphotons_good = SG::makeHandle<unsigned int>(m_nphotons_good_key); + + dec_nphotons(*vtx_pointing) = photons_selected.size(); + dec_ph_links(*vtx_pointing) = ph_links; + dec_nphotons_good(*vtx_pointing) = n_good_photons; + + // compute common vertex from the pointing + const TLorentzVector p4_photons = std::accumulate( + photons_selected.begin(), photons_selected.end(), + TLorentzVector(0, 0, 0, 0), + [](TLorentzVector sum, const xAOD::Egamma* photon) { + TLorentzVector v_ph; + const xAOD::CaloCluster* cluster = photon->caloCluster(); + v_ph.SetPtEtaPhiM(photon->e() / cosh(cluster->etaBE(2)), + cluster->etaBE(2), cluster->phiBE(2), 0.0); + sum += v_ph; + return sum; + }); + + auto dec_photons_px = SG::makeHandle<float>(m_photons_px_key); + auto dec_photons_py = SG::makeHandle<float>(m_photons_py_key); + auto dec_photons_pz = SG::makeHandle<float>(m_photons_pz_key); + + dec_photons_px(*vtx_pointing) = p4_photons.X(); + dec_photons_py(*vtx_pointing) = p4_photons.Y(); + dec_photons_pz(*vtx_pointing) = p4_photons.Z(); + + return StatusCode::SUCCESS; +} + +void BuildVertexPointingAlg::selectPhotons( + const xAOD::EgammaContainer& original_photons, + ConstDataVector<DataVector<xAOD::Egamma>>& photons_selected) const { + for (const auto& photon : original_photons) { + if (m_selectionTool->accept(photon)) { + photons_selected.push_back(photon); + ATH_MSG_DEBUG("photon selected for pointing: pT = " << photon->pt() + << " GeV"); + } + } + ATH_MSG_DEBUG( + "number of photons selected for pointing: " << photons_selected.size()); + + // limit the number of photons to use + if (m_nphotons_to_use >= 0 and + photons_selected.size() > static_cast<std::size_t>(m_nphotons_to_use)) { + std::sort(photons_selected.begin(), photons_selected.end(), + [](const xAOD::IParticle* a, const xAOD::IParticle* b) { + return a->pt() > b->pt(); // sort from the highest pT + }); + + photons_selected.resize(m_nphotons_to_use); + } + ATH_MSG_DEBUG("number of photons selected for pointing after limiting: " + << photons_selected.size()); + if (ATH_UNLIKELY(this->msgLvl(MSG::DEBUG))) { + for (const xAOD::Egamma* photon : photons_selected) { + ATH_MSG_DEBUG("photon selected for pointing: pT = " << photon->pt() + << " GeV"); + } + } +} diff --git a/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/src/components/PhotonVertexSelection_entries.cxx b/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/src/components/PhotonVertexSelection_entries.cxx index 35aa49809177a1a057c2c4064f05c1bb7dd2c3fb..c3e4f5c27aaa7d122404e5af2a7512632f755810 100644 --- a/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/src/components/PhotonVertexSelection_entries.cxx +++ b/PhysicsAnalysis/ElectronPhotonID/PhotonVertexSelection/src/components/PhotonVertexSelection_entries.cxx @@ -1,8 +1,12 @@ #include "PhotonVertexSelection/PhotonVertexSelectionTool.h" #include "PhotonVertexSelection/PhotonPointingTool.h" #include "../PhotonVertexSelectionAlg.h" +#include "PhotonVertexSelection/DecoratePhotonPointingAlg.h" +#include "PhotonVertexSelection/BuildVertexPointingAlg.h" DECLARE_COMPONENT( CP::PhotonVertexSelectionTool ) DECLARE_COMPONENT( CP::PhotonPointingTool ) DECLARE_COMPONENT( CP::PhotonVertexSelectionAlg ) +DECLARE_COMPONENT( DecoratePhotonPointingAlg ) +DECLARE_COMPONENT( BuildVertexPointingAlg ) diff --git a/Tracking/TrkConfig/python/TrkVertexWeightCalculatorsConfig.py b/Tracking/TrkConfig/python/TrkVertexWeightCalculatorsConfig.py index c790a3df3110936c3695a0aecbeaa334ad225187..39ca6ba093b8894a87fe218165b619eaaabd0909 100644 --- a/Tracking/TrkConfig/python/TrkVertexWeightCalculatorsConfig.py +++ b/Tracking/TrkConfig/python/TrkVertexWeightCalculatorsConfig.py @@ -1,8 +1,9 @@ -# Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration +# Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration # Configuration of TrkVertexWeightCalculators package from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator from AthenaConfiguration.ComponentFactory import CompFactory + def SumPt2VertexWeightCalculatorCfg(flags, name="SumPt2VertexWeightCalculator", **kwargs): acc = ComponentAccumulator() @@ -11,6 +12,7 @@ def SumPt2VertexWeightCalculatorCfg(flags, name="SumPt2VertexWeightCalculator", CompFactory.Trk.SumPtVertexWeightCalculator(name, **kwargs)) return acc + def SumPtVertexWeightCalculatorCfg(flags, name="SumPtVertexWeightCalculator", **kwargs): acc = ComponentAccumulator() @@ -18,3 +20,55 @@ def SumPtVertexWeightCalculatorCfg(flags, name="SumPtVertexWeightCalculator", acc.setPrivateTools( CompFactory.Trk.SumPtVertexWeightCalculator(name, **kwargs)) return acc + + +def BDTVertexWeightCalculatorCfg(flags, **kwargs): + """ + Configure the BDTVertexWeightCalculator. Note: this tool needs to be run after the + DecoratePhotonPointingAlg and BuildVertexPointingAlg. + Use BDTVertexWeightCalculatorSeqCfg to have the full sequence. + """ + acc = ComponentAccumulator() + kwargs.setdefault( + "BDTFile", "PhotonVertexSelection/BDT/2023-02-28/global_ggHW_phcount_BDT.root" + ) + kwargs.setdefault("BDTName", "lgbm") + kwargs.setdefault("PointingVertexContainerKey", "PhotonPointingVertices") + tool = CompFactory.BDTVertexWeightCalculator( + "BDTVertexWeightCalculator", **kwargs + ) + acc.setPrivateTools(tool) + return acc + + +def BDTVertexWeightCalculatorSeqCfg(flags, container='Photons', **kwargs): + """ + Configure BDTVertexWeightCalculator and the algorithms that are needed to run before. + Optional parameters are passed only to the tool. + """ + acc = ComponentAccumulator() + + from PhotonVertexSelection.PhotonVertexSelectionConfig import ( + DecoratePhotonPointingAlgCfg, + BuildVertexPointingAlgCfg, + ) + + # this algorithm decorates the photons with the pointing information + acc.merge(DecoratePhotonPointingAlgCfg(flags, "DecoratePhotonPointingAlg", PhotonContainerKey=container)) + + # this algorithm creates the vertex from photon pointing + acc.merge( + BuildVertexPointingAlgCfg( + flags, + "BuildVertexPointingAlg", + PhotonContainerKey=container, + PointingVertexContainerKey=kwargs.get( + "PointingVertexContainerKey", "PhotonPointingVertices" + ), + ) + ) + + accTool = BDTVertexWeightCalculatorCfg(flags, **kwargs) + tool = acc.popToolsAndMerge(accTool) + acc.setPrivateTools(tool) + return acc \ No newline at end of file diff --git a/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/CMakeLists.txt b/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/CMakeLists.txt index 2650c449a89ad7a9a6bd75faabe4e643a880daff..374f3be583a2cffe4a3bc2b9aaeb2326e5a852ca 100644 --- a/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/CMakeLists.txt +++ b/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/CMakeLists.txt @@ -1,10 +1,10 @@ -# Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration +# Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration # Declare the package name: atlas_subdir( TrkVertexWeightCalculators ) # External dependencies: -find_package( ROOT COMPONENTS Core MathCore Hist Matrix ) +find_package( ROOT COMPONENTS Tree Core MathCore Hist Matrix ) # Component(s) in the package: atlas_add_component( TrkVertexWeightCalculators @@ -12,7 +12,17 @@ atlas_add_component( TrkVertexWeightCalculators src/components/*.cxx INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps GeoPrimitives xAODTracking GaudiKernel - TrkVertexFitterInterfaces GeneratorObjects TrkParameters VxVertex TrkNeuralNetworkUtilsLib) + TrkVertexFitterInterfaces GeneratorObjects VxVertex TrkNeuralNetworkUtilsLib + AsgTools MVAUtils AnaAlgorithmLib PhotonVertexSelectionLib PATCoreLib + PRIVATE_LINK_LIBRARIES PathResolver PhotonVertexSelectionLib TrkParameters) -# Install files from the package: +# Install files from the package atlas_install_runtime( share/*.root ) +atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} ) +atlas_install_scripts( scripts/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} ) + +# Add test +atlas_add_test( TrkVertexWeightCalculatorsBDT_test + SCRIPT python -m TrkVertexWeightCalculators.TrkVertexWeightCalculatorsConfig + POST_EXEC_SCRIPT noerror.sh ) + diff --git a/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/TrkVertexWeightCalculators/BDTVertexWeightCalculator.h b/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/TrkVertexWeightCalculators/BDTVertexWeightCalculator.h new file mode 100644 index 0000000000000000000000000000000000000000..0198571f3a10b7bb4b1bf5befa62c90d742fdc35 --- /dev/null +++ b/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/TrkVertexWeightCalculators/BDTVertexWeightCalculator.h @@ -0,0 +1,72 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TRKVERTEXWEIGHTCALCULATORS_BDTVERTEXWEIGHTCALCULATOR_H +#define TRKVERTEXWEIGHTCALCULATORS_BDTVERTEXWEIGHTCALCULATOR_H + +#include <AsgDataHandles/ReadDecorHandleKey.h> +#include <AsgDataHandles/ReadHandleKey.h> +#include <AsgTools/AsgTool.h> +#include <TrkVertexFitterInterfaces/IVertexWeightCalculator.h> +#include <xAODTracking/VertexContainer.h> + +#include <memory> +#include <string> + +#include "AsgTools/PropertyWrapper.h" + +namespace MVAUtils { +class BDT; +} + +/** + * @brief A tool to select the primary vertex associated to the hard-scatter + * using a BDT + * + * This tool uses a BDT to select the hard-scatter vertex from a collection of + * primary vertices. + * It reads the pointing information from the pointing vertex created by + * @c BuildVertexPointingAlg. + * + */ +class BDTVertexWeightCalculator final + : public asg::AsgTool, + virtual public Trk::IVertexWeightCalculator { + ASG_TOOL_CLASS(BDTVertexWeightCalculator, IVertexWeightCalculator) + public: + BDTVertexWeightCalculator(const std::string& name); + virtual ~BDTVertexWeightCalculator() override; + virtual StatusCode initialize() override; + /** + * @brief Estimate the compatibility of the vertex with a hard scatter vertex, + * with respect to pileup vertices + */ + virtual double estimateSignalCompatibility( + const xAOD::Vertex& vertex) const override; + const xAOD::Vertex* getVertex(const xAOD::VertexContainer& vertices) const; + + private: + SG::ReadHandleKey<xAOD::VertexContainer> m_pointingVertexContainerKey{ + this, "PointingVertexContainerKey", "PhotonPointingVertices", + "The container with the vertex build with (photon) pointing"}; + SG::ReadDecorHandleKey<xAOD::VertexContainer> m_nphotons_good_key{ + this, "NGoodPhotonsKey", m_pointingVertexContainerKey, "nphotons_good"}; + SG::ReadDecorHandleKey<xAOD::VertexContainer> m_photons_px_key{ + this, "PhotonsPxKey", m_pointingVertexContainerKey, "photons_px"}; + SG::ReadDecorHandleKey<xAOD::VertexContainer> m_photons_py_key{ + this, "PhotonsPyKey", m_pointingVertexContainerKey, "photons_py"}; + SG::ReadDecorHandleKey<xAOD::VertexContainer> m_photons_pz_key{ + this, "PhotonsPzKey", m_pointingVertexContainerKey, "photons_pz"}; + Gaudi::Property<std::string> m_bdt_file{this, "BDTFile", "", + "BDT ROOT filename"}; + Gaudi::Property<std::string> m_bdt_name{this, "BDTName", "BDT", + "BDT tree name"}; + + std::unique_ptr<MVAUtils::BDT> m_bdt; + + std::vector<float> get_input_values(const xAOD::Vertex& vertex) const; + StatusCode initialize_BDT(); +}; + +#endif // TRKVERTEXWEIGHTCALCULATORS_BDTVERTEXWEIGHTCALCULATOR_H \ No newline at end of file diff --git a/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/TrkVertexWeightCalculators/DecorateVertexScoreAlg.h b/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/TrkVertexWeightCalculators/DecorateVertexScoreAlg.h new file mode 100644 index 0000000000000000000000000000000000000000..6d98a280bee59a08446b5b539620bbcffd95e344 --- /dev/null +++ b/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/TrkVertexWeightCalculators/DecorateVertexScoreAlg.h @@ -0,0 +1,46 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef DECORATEVERTEXSCOREALG_H +#define DECORATEVERTEXSCOREALG_H + +#include <AnaAlgorithm/AnaAlgorithm.h> +#include <AsgDataHandles/ReadHandleKey.h> +#include <AsgDataHandles/WriteDecorHandleKey.h> +#include <AsgTools/ToolHandle.h> +#include <TrkVertexFitterInterfaces/IVertexWeightCalculator.h> +#include <xAODTracking/VertexContainerFwd.h> + +#include <string> + +/** + * @brief An algorithm to decorate vertices with a score + * + * This algorithm uses the VertexWeightCalculator to decorate vertices with a + * score. The score should be related to the classification of the vertex as + * coming from the hard-scatter interaction or pileup. + * + */ +class DecorateVertexScoreAlg : public EL::AnaAlgorithm { + public: + DecorateVertexScoreAlg(const std::string& name, + ISvcLocator* svcLoc = nullptr); + + virtual StatusCode initialize() override; + StatusCode execute() override; + + private: + SG::ReadHandleKey<xAOD::VertexContainer> m_vertexContainer_key{ + this, "VertexContainer", "PrimaryVertices", "The input Vertex container"}; + + SG::WriteDecorHandleKey<xAOD::VertexContainer> m_vertexScoreDecor_key{ + this, "VertexScoreDecor", m_vertexContainer_key, "score", + "The output decoration for the score vertices"}; + + ToolHandle<Trk::IVertexWeightCalculator> m_vertexWeightCalculator{ + this, "VertexWeightCalculator", "BDTVertexWeightCalculator", + "The tool to compute the score"}; +}; + +#endif // DECORATEVERTEXSCOREALG_H \ No newline at end of file diff --git a/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/TrkVertexWeightCalculators/TrkVertexWeightCalculatorsDict.h b/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/TrkVertexWeightCalculators/TrkVertexWeightCalculatorsDict.h new file mode 100644 index 0000000000000000000000000000000000000000..0277e03d0fcb03f7fcfa6ff504fdf55ec18586f5 --- /dev/null +++ b/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/TrkVertexWeightCalculators/TrkVertexWeightCalculatorsDict.h @@ -0,0 +1,10 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef __TRKVERTEXWEIGHTCALCULATORS_DICT__ +#define __TRKVERTEXWEIGHTCALCULATORS_DICT__ +#include "TrkVertexWeightCalculators/BDTVertexWeightCalculator.h" +#include "TrkVertexWeightCalculators/DecorateVertexScoreAlg.h" +#endif + diff --git a/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/TrkVertexWeightCalculators/selection.xml b/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/TrkVertexWeightCalculators/selection.xml new file mode 100644 index 0000000000000000000000000000000000000000..2e383158444fdf98303aa152604d158c929070fd --- /dev/null +++ b/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/TrkVertexWeightCalculators/selection.xml @@ -0,0 +1,4 @@ + <lcgdict> + <class name="BDTVertexWeightCalculator" /> + <class name="DecorateVertexScoreAlg" /> + </lcgdict> diff --git a/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/doc/packagedoc.h b/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/doc/packagedoc.h new file mode 100644 index 0000000000000000000000000000000000000000..997fc3721a6cd8b6789f52cf5f97614a6742b3ad --- /dev/null +++ b/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/doc/packagedoc.h @@ -0,0 +1,17 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ + +/** + +@page TrkVertexWeightCalculators_page TrkVertexWeightCalculators Package + +This package provides tools to set a weight to reconstructed primary vertices, in order to select the one +coming from the hard-scatter interaction. + +- SumPtVertexWeightCalculator: return sumPT / sumPT2 of the tracks associated to the vertex. +- TrueVertexDistanceWeightCalculator: return 1 / deltaZ, where deltaZ is the distance between the reconstructed vertex and the true vertex. +- BDTVertexWeightCalculator: BDT-based, developed by Davide Riva during his bachelor thesis. @author Ruggero Turra <ruggero.turra@cern.ch> + +Most of the configuration is in Tracking/TrkConfig/python/TrkVertexWeightCalculatorsConfig.py +**/ \ No newline at end of file diff --git a/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/python/TrkVertexWeightCalculatorsConfig.py b/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/python/TrkVertexWeightCalculatorsConfig.py new file mode 100644 index 0000000000000000000000000000000000000000..4e3b9b55141528d3d020024e4030bad9f9ec4ac7 --- /dev/null +++ b/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/python/TrkVertexWeightCalculatorsConfig.py @@ -0,0 +1,150 @@ +# Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration + +from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator +from AthenaConfiguration.ComponentFactory import CompFactory + + +def DecorateVertexScoreAlgCfg(flags, name="DecorateVertexScoreAlg", **kwargs): + """ + This algorithm decorates all the vertices with the score computed by a tool. + """ + acc = ComponentAccumulator() + if not kwargs.get("VertexWeightCalculator", None): + from TrkConfig.TrkVertexWeightCalculatorsConfig import BDTVertexWeightCalculatorSeqCfg + toolBDTAcc = BDTVertexWeightCalculatorSeqCfg(flags) + tool = toolBDTAcc.popPrivateTools() + acc.merge(toolBDTAcc) + kwargs["VertexWeightCalculator"] = tool + + kwargs.setdefault("VertexScoreDecor", "score") + alg = CompFactory.DecorateVertexScoreAlg(name, **kwargs) + acc.addEventAlgo(alg) + return acc + + +def TrkVertexWeightCalculatorDebugCfg(flags, **kwargs): + """ + This is a test configuration for the TrkVertexWeightCalculator. It is not meant to be run in production. + It produces a ROOT file with a tree containing relevant information to check the performance of the tool. + """ + from AthenaCommon.Constants import DEBUG + from AthenaCommon.Logging import logging + + mlog = logging.getLogger("TrkVertexWeightCalculatorBDTDebugCfg") + mlog.warning( + "This is a test algorithm, it is not meant to be run in production." + ) + + acc = ComponentAccumulator() + acc.merge(DecorateVertexScoreAlgCfg(flags, **kwargs)) + mlog.info( + "Setting the output level of BuildVertexPointingAlg, DecorateVertexScoreAlg, DecorateVertexScoreAlg/VertexSelectionTool to DEBUG" + ) + acc.getEventAlgo("BuildVertexPointingAlg").OutputLevel = DEBUG + acc.getEventAlgo("DecorateVertexScoreAlg").OutputLevel = DEBUG + acc.getEventAlgo( + "DecorateVertexScoreAlg" + ).VertexWeightCalculator.OutputLevel = DEBUG + + from TrkConfig.TrkVertexWeightCalculatorsConfig import ( + SumPt2VertexWeightCalculatorCfg, + ) + + tool_pt2 = acc.popToolsAndMerge(SumPt2VertexWeightCalculatorCfg(flags)) + acc.merge( + DecorateVertexScoreAlgCfg( + flags, + "DecorateVertexScoreAlgSumPt2", + VertexWeightCalculator=tool_pt2, + VertexScoreDecor="score_sumpt2", + ) + ) + acc.getEventAlgo("DecorateVertexScoreAlgSumPt2").OutputLevel = DEBUG + + tool = CompFactory.Trk.TrueVertexDistanceWeightCalculator() + acc.merge( + DecorateVertexScoreAlgCfg( + flags, + "DecorateVertexScoreAlgTrueVertexDistance", + VertexWeightCalculator=tool, + VertexScoreDecor="score_true_vertex_distance", + ) + ) + acc.getEventAlgo("DecorateVertexScoreAlgTrueVertexDistance").OutputLevel = DEBUG + + sysService = CompFactory.CP.SystematicsSvc("SystematicsSvc", sigmaRecommended=0) + acc.addService(sysService) + + histoSvc = CompFactory.THistSvc( + Output=[ + f"ANALYSIS DATAFILE='{flags.Output.HISTFileName}' TYPE='ROOT' OPT='RECREATE'" + ] + ) + acc.addService(histoSvc) + acc.setAppProperty("HistogramPersistency", "ROOT") + + acc.addEventAlgo( + CompFactory.CP.TreeMakerAlg("TreeMaker", TreeName=flags.Output.TreeName) + ) + branches = [ + "EventInfo.runNumber -> runNumber", + "EventInfo.eventNumber -> eventNumber", + "EventInfo.actualInteractionsPerCrossing -> actualInteractionsPerCrossing", + "EventInfo.averageInteractionsPerCrossing -> averageInteractionsPerCrossing", + "PrimaryVertices.x -> vtx_x", + "PrimaryVertices.y -> vtx_y", + "PrimaryVertices.z -> vtx_z", + "PrimaryVertices.score -> vtx_score", + "PhotonPointingVertices.z -> z_common", + "PhotonPointingVertices.nphotons_good -> nphotons_good", + "PrimaryVertices.score_sumpt2 -> vtx_score_sumpt2", + "PrimaryVertices.score_true_vertex_distance -> vtx_score_true_vertex_distance", + ] + acc.addEventAlgo( + CompFactory.CP.AsgxAODNTupleMakerAlg( + "NTupleMaker", TreeName=flags.Output.TreeName, Branches=branches + ) + ) + acc.addEventAlgo( + CompFactory.CP.TreeFillerAlg("TreeFiller", TreeName=flags.Output.TreeName) + ) + + return acc + + +def TrkVertexWeightCalculatorBDTDebugRunCfg(): + from AthenaConfiguration.AllConfigFlags import initConfigFlags + + flags = initConfigFlags() + flags.Exec.MaxEvents = 100 + from AthenaConfiguration.TestDefaults import defaultTestFiles + + flags.Input.Files = defaultTestFiles.AOD_RUN3_MC + flags.Output.HISTFileName = "test_tree.root" + flags.addFlag("Output.TreeName", "tree") + + flags.fillFromArgs() + flags.lock() + flags.dump() + + from AthenaConfiguration.MainServicesConfig import MainServicesCfg + + acc = MainServicesCfg(flags) + + from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg + + acc.merge(PoolReadCfg(flags)) + + acc.merge(TrkVertexWeightCalculatorDebugCfg(flags)) + acc.printConfig(withDetails=True, summariseProps=True) + acc.store(open("TrkVertexWeightCalculatorBDTConfig.pkl", "wb")) + return acc + + +if __name__ == "__main__": + acc = TrkVertexWeightCalculatorBDTDebugRunCfg() + status = acc.run() + + import sys + + sys.exit(not status.isSuccess()) diff --git a/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/python/__init__.py b/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/python/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/scripts/TrkVertexWeightCalculatorsExample.py b/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/scripts/TrkVertexWeightCalculatorsExample.py new file mode 100755 index 0000000000000000000000000000000000000000..55c894959de286944e55051da94f41a53dcb1b19 --- /dev/null +++ b/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/scripts/TrkVertexWeightCalculatorsExample.py @@ -0,0 +1,20 @@ +#!/usr/bin/env python +# Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration + +""" +This script is equivalent to python -m TrkVertexWeightCalculators.TrkVertexWeightCalculatorsConfig +but it is an installed executable script. This simplify the grid submission +since there is no easy way to run CA with pathena / prun. + +pathena --trf "TrkVertexWeighCalculatorsExample.py --filesInput=%IN Output.HISTFileName=%OUT.root --evtMax=-1" --inDS mc21_13p6TeV:mc21_13p6TeV.601521.PhPy8EG_PDF4LHC21_ggH_MINLO_gammagamma.merge.AOD.e8472_e8455_s3873_s3874_r13829_r13831 --outDS user.turra.mc21_13p6TeV.601521.PhPy8EG_PDF4LHC21_ggH_MINLO_gammagamma.dumpvertex.NTUP.e8472_s3873_r13829_v9 +""" + +import sys + +from TrkVertexWeightCalculators.TrkVertexWeightCalculatorsConfig import TrkVertexWeightCalculatorBDTDebugRunCfg + +acc = TrkVertexWeightCalculatorBDTDebugRunCfg() +status = acc.run() + +sys.exit(not status.isSuccess()) + diff --git a/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/src/BDTVertexWeightCalculator.cxx b/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/src/BDTVertexWeightCalculator.cxx new file mode 100644 index 0000000000000000000000000000000000000000..cf0628fa948ee66ce27aa0badbc107827f03e544 --- /dev/null +++ b/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/src/BDTVertexWeightCalculator.cxx @@ -0,0 +1,147 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ +#include "TrkVertexWeightCalculators/BDTVertexWeightCalculator.h" + +#include "MVAUtils/BDT.h" +#include "PathResolver/PathResolver.h" +#include "PhotonVertexSelection/PhotonVertexHelpers.h" +#include "TFile.h" +#include "TTree.h" +#include "TVector3.h" +#include "xAODTracking/Vertex.h" +#include "xAODTracking/VertexContainer.h" +#include <AsgDataHandles/ReadHandle.h> +#include <AsgDataHandles/ReadDecorHandle.h> + +#include <cmath> +#include <limits> + +BDTVertexWeightCalculator::BDTVertexWeightCalculator(const std::string& name) + : asg::AsgTool(name) { } + +BDTVertexWeightCalculator::~BDTVertexWeightCalculator() = default; + +StatusCode BDTVertexWeightCalculator::initialize() { + ATH_MSG_DEBUG("Initializing " << name() << "..."); + + ATH_CHECK(m_pointingVertexContainerKey.initialize()); + ATH_CHECK(m_nphotons_good_key.initialize()); + ATH_CHECK(m_photons_px_key.initialize()); + ATH_CHECK(m_photons_py_key.initialize()); + ATH_CHECK(m_photons_pz_key.initialize()); + + ATH_CHECK(initialize_BDT()); + + return StatusCode::SUCCESS; +} + +StatusCode BDTVertexWeightCalculator::initialize_BDT() { + const std::string fullPathToFile = PathResolverFindCalibFile(m_bdt_file); + if (fullPathToFile.empty()) { + ATH_MSG_ERROR("Can not find BDT file: '" << m_bdt_file + << "'. Please check if the " + "file exists."); + return StatusCode::FAILURE; + } + + std::unique_ptr<TFile> rootFile(TFile::Open(fullPathToFile.c_str(), "READ")); + if (!rootFile) { + ATH_MSG_ERROR("Can not open BDT file root file: '" + << fullPathToFile + << "'. Please check if the file " + "exists and is not corrupted."); + + return StatusCode::FAILURE; + } + std::unique_ptr<TTree> tree((TTree*)rootFile->Get(m_bdt_name.value().c_str())); + if (!tree) { + ATH_MSG_ERROR("Can not find BDT tree '" << m_bdt_name << "' in file: '" + << fullPathToFile << "'."); + } + + ATH_MSG_INFO("Loading BDT '" << m_bdt_name << "' from file: '" + << fullPathToFile << "'."); + m_bdt = std::make_unique<MVAUtils::BDT>(tree.get()); + + return StatusCode::SUCCESS; +} + +std::vector<float> BDTVertexWeightCalculator::get_input_values( + const xAOD::Vertex& vertex) const { + SG::ReadHandle<xAOD::VertexContainer> pointingVertexContainer{ + m_pointingVertexContainerKey}; + if (pointingVertexContainer->size() != 1) { + ATH_MSG_ERROR("Pointing vertex container has size " + << pointingVertexContainer->size() + << " instead of 1"); + return {}; + } + const xAOD::Vertex* pointingVertex = pointingVertexContainer->front(); + + const float zcommon = pointingVertex->z(); + const float zcommon_error = + std::sqrt(pointingVertex->covariancePosition()(2, 2)); + + auto acc_nph_good = SG::makeHandle<unsigned int>(m_nphotons_good_key); + const float nph_good = acc_nph_good(*pointingVertex); + const float vtx_sumpt = xAOD::PVHelpers::getVertexSumPt(&vertex, 1); + const float vtx_sumpt2 =xAOD::PVHelpers::getVertexSumPt(&vertex, 2); + const float vtx_z = vertex.z(); + + const float deltaz_sig_photons = std::abs(vtx_z - zcommon) / zcommon_error; + + auto acc_photons_px = SG::makeHandle<float>(m_photons_px_key); + auto acc_photons_py = SG::makeHandle<float>(m_photons_py_key); + auto acc_photons_pz = SG::makeHandle<float>(m_photons_pz_key); + TVector3 photon_sump(acc_photons_px(*pointingVertex), + acc_photons_py(*pointingVertex), + acc_photons_pz(*pointingVertex)); + TLorentzVector vertex_p4 = xAOD::PVHelpers::getVertexMomentum(&vertex, true); + const float deltaphi_sig_photons = vertex_p4.Vect().DeltaPhi(photon_sump); + + ATH_MSG_DEBUG("Input variables: nph_good: " + << nph_good << ", vtx_sumpt: " << vtx_sumpt + << ", vtx_sumpt2: " << vtx_sumpt2 + << ", deltaz_sig_photons: " << deltaz_sig_photons + << ", deltaphi_sig_photons: " << deltaphi_sig_photons); + + return {nph_good, vtx_sumpt, vtx_sumpt2, deltaz_sig_photons, + deltaphi_sig_photons}; +} + +double BDTVertexWeightCalculator::estimateSignalCompatibility( + const xAOD::Vertex& vertex) const { + if (vertex.vertexType() != xAOD::VxType::VertexType::PriVtx and + vertex.vertexType() != xAOD::VxType::VertexType::PileUp) { + return 0.; + } + const std::vector<float> input_values = get_input_values(vertex); + if (input_values.empty()) { return 0.; } + return m_bdt->GetClassification(input_values); +} + +const xAOD::Vertex* BDTVertexWeightCalculator::getVertex( + const xAOD::VertexContainer& vertices) const { + float best_score = std::numeric_limits<float>::lowest(); + const xAOD::Vertex* vertex = nullptr; + for (const auto& v : vertices) { + if (v == nullptr) { + ATH_MSG_WARNING("Null vertex in container"); + continue; + } + if (v->vertexType() == xAOD::VxType::PriVtx) { + const float score = v->isAvailable<float>("score") + ? v->auxdata<float>("score") + : estimateSignalCompatibility(*v); + if (score > best_score) { + best_score = score; + vertex = v; + } + } + } + if (vertex == nullptr) { + ATH_MSG_WARNING("No vertex found"); + } + return vertex; +} diff --git a/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/src/DecorateVertexScoreAlg.cxx b/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/src/DecorateVertexScoreAlg.cxx new file mode 100644 index 0000000000000000000000000000000000000000..259a5da38ba9f5e50b0a50c23d1f7521a8cf7773 --- /dev/null +++ b/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/src/DecorateVertexScoreAlg.cxx @@ -0,0 +1,38 @@ +/* + Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration +*/ + +#include "TrkVertexWeightCalculators/DecorateVertexScoreAlg.h" + +#include "AsgDataHandles/WriteDecorHandle.h" +#include "xAODTracking/Vertex.h" +#include "xAODTracking/VertexContainer.h" + +DecorateVertexScoreAlg::DecorateVertexScoreAlg(const std::string& name, + ISvcLocator* svcLoc) + : EL::AnaAlgorithm(name, svcLoc) {} + +StatusCode DecorateVertexScoreAlg::initialize() { + ATH_CHECK(m_vertexContainer_key.initialize()); + ATH_CHECK(m_vertexScoreDecor_key.initialize()); + ATH_CHECK(m_vertexWeightCalculator.retrieve()); + return StatusCode::SUCCESS; +} + +StatusCode DecorateVertexScoreAlg::execute() { + SG::ReadHandle<xAOD::VertexContainer> vertices{m_vertexContainer_key}; + if (!vertices.isValid()) { + ATH_MSG_ERROR(m_vertexContainer_key.key() << " not found"); + return StatusCode::FAILURE; + } + + auto vertexScoreDecor = SG::makeHandle<float>(m_vertexScoreDecor_key); + + ATH_MSG_DEBUG("Decorating " << vertices->size() << " vertices with scores"); + for (const auto* vertex : *vertices) { + const float score = m_vertexWeightCalculator->estimateSignalCompatibility(*vertex); + vertexScoreDecor(*vertex) = score; + ATH_MSG_DEBUG("Decorated vertex with score: " << score); + } + return StatusCode::SUCCESS; +} diff --git a/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/src/components/TrkVertexWeightCalculators_entries.cxx b/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/src/components/TrkVertexWeightCalculators_entries.cxx index e6a77d6b2ab348bf6c33acb15e1b16c935a4b6a4..c507ee80f1a7abbd01b8535eacfa8342e2dccca6 100644 --- a/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/src/components/TrkVertexWeightCalculators_entries.cxx +++ b/Tracking/TrkVertexFitter/TrkVertexWeightCalculators/src/components/TrkVertexWeightCalculators_entries.cxx @@ -1,9 +1,11 @@ - #include "TrkVertexWeightCalculators/SumPtVertexWeightCalculator.h" - #include "TrkVertexWeightCalculators/TrueVertexDistanceWeightCalculator.h" +#include "TrkVertexWeightCalculators/SumPtVertexWeightCalculator.h" +#include "TrkVertexWeightCalculators/TrueVertexDistanceWeightCalculator.h" +#include "TrkVertexWeightCalculators/BDTVertexWeightCalculator.h" +#include "TrkVertexWeightCalculators/DecorateVertexScoreAlg.h" - - using namespace Trk; - - DECLARE_COMPONENT( SumPtVertexWeightCalculator ) - DECLARE_COMPONENT( TrueVertexDistanceWeightCalculator ) +using namespace Trk; +DECLARE_COMPONENT( SumPtVertexWeightCalculator ) +DECLARE_COMPONENT( TrueVertexDistanceWeightCalculator ) +DECLARE_COMPONENT( BDTVertexWeightCalculator ) +DECLARE_COMPONENT( DecorateVertexScoreAlg ) \ No newline at end of file