diff --git a/Reconstruction/Jet/JetRec/JetRec/JetClusterer.h b/Reconstruction/Jet/JetRec/JetRec/JetClusterer.h index 32fbe46d8b8bed5843b1d39c548b81740c8b46c2..64b71f7bd2e950f6d0c316b40ef139c8b640341b 100644 --- a/Reconstruction/Jet/JetRec/JetRec/JetClusterer.h +++ b/Reconstruction/Jet/JetRec/JetRec/JetClusterer.h @@ -19,13 +19,14 @@ #include "xAODEventInfo/EventInfo.h" #include "StoreGate/ReadHandleKey.h" +#include "StoreGate/WriteHandleKey.h" #include "JetInterface/IJetProvider.h" #include "AsgTools/AsgTool.h" #include "JetRec/PseudoJetContainer.h" #include "JetRec/JetFromPseudojet.h" - +#include "JetEDM/PseudoJetVector.h" #include "fastjet/PseudoJet.hh" #include "fastjet/AreaDefinition.hh" @@ -61,7 +62,7 @@ protected: SG::ReadHandleKey<PseudoJetContainer> m_inputPseudoJets {this, "InputPseudoJets", "inputpseudojet", "input constituents"}; /// used to build the key under which the final PJ will be stored in evtStore() - std::string m_finalPseudoJets; + SG::WriteHandleKey<PseudoJetVector> m_finalPseudoJets {this, "FinalPseudoJets_DONOTSET", "", "output pseudojets -- autoconfigured name"}; // Job options. Gaudi::Property<std::string> m_jetalg {this, "JetAlgorithm", "AntiKt", "alg type : AntiKt, Kt, CA..."}; diff --git a/Reconstruction/Jet/JetRec/JetRec/JetGroomer.h b/Reconstruction/Jet/JetRec/JetRec/JetGroomer.h index 3cbff98faf653e31b81a69fc7cfc2fb7ec9e6ce0..60bca3c63ccdc57afbdca8c32a31ad5505fbdab6 100644 --- a/Reconstruction/Jet/JetRec/JetRec/JetGroomer.h +++ b/Reconstruction/Jet/JetRec/JetRec/JetGroomer.h @@ -9,16 +9,23 @@ /// /// \class JetGroomer /// -/// Creates a new JetContainer by grooming an input jet collection +/// Base class for tools that create a new JetContainer by grooming an input jet collection /// -/// This tool implements the IJetProvider interface. The JetContainer it returns is built by -/// running an IJetGroomer tool (e.g. JetTrimmer) on the input jets. +/// This tool implements the IJetProvider interface. Children should override (only) +/// the initialize and insertGroomedJet methods, to implement a specific instance of +/// jet grooming. /// #include "AsgTools/AsgTool.h" + #include "StoreGate/ReadHandleKey.h" +#include "StoreGate/WriteHandleKey.h" + #include "JetInterface/IJetProvider.h" -#include "JetInterface/IJetGroomer.h" + +#include "JetEDM/PseudoJetVector.h" +#include "JetRec/PseudoJetContainer.h" + #include "xAODJet/JetContainer.h" #include "xAODJet/JetAuxContainer.h" @@ -32,17 +39,23 @@ class JetGroomer using asg::AsgTool::AsgTool; virtual StatusCode initialize() override; - virtual std::pair<std::unique_ptr<xAOD::JetContainer>, std::unique_ptr<SG::IAuxStore> > getJets() const override; - private: - // Handle Input JetContainer - SG::ReadHandleKey<xAOD::JetContainer> m_inputJetsKey {this, "UngroomedJets", "", "Jet collection to be groomed"}; + // From IJetProvider + virtual std::pair<std::unique_ptr<xAOD::JetContainer>, std::unique_ptr<SG::IAuxStore> > getJets() const override final; + + // Implementation of grooming goes here + // The jet is inserted into the output container, which is necessary for speed + // in the xAOD container paradigm + virtual void insertGroomedJet(const xAOD::Jet&, const PseudoJetContainer&, xAOD::JetContainer&, PseudoJetVector&) const = 0; + + protected: + /// Handle Input JetContainer (this contains the parent ungroomed jets to be trimmed) + SG::ReadHandleKey<xAOD::JetContainer> m_inputJetContainer {this, "UngroomedJets", "ungroomedinput", "Input ungroomed jet container"}; - // Handle Input PseudoJetContainer - // Needed to extract the constituents - SG::ReadHandleKey<PseudoJetContainer> m_inputPseudoJetsKey {this, "ParentPseudoJets", "", "Input constituents"}; + /// This is the input to the parent JetContainer. It is needed in order to re-assign the ghost constituents from the final groomed PJ to the xAOD::Jet + SG::ReadHandleKey<PseudoJetContainer> m_inputPseudoJets {this, "ParentPseudoJets", "inputpseudojet", "input constituents of parent JetContainer"}; - ToolHandle<IJetGroomer> m_groomer ={this , "Groomer" , {} , "Tool grooming the jets (trim, prune, softdrop etc)"}; + SG::WriteHandleKey<PseudoJetVector> m_finalPseudoJets {this, "FinalPseudoJets_DONOTSET", "", "output pseudojets -- autoconfigured name"}; }; diff --git a/Reconstruction/Jet/JetRec/JetRec/PseudoJetTranslator.h b/Reconstruction/Jet/JetRec/JetRec/PseudoJetTranslator.h index 6053fd625d034a8709f3dd97407ee0b0517a6596..84cd8386558978265c2208b4f2ff98c3aa8a1d66 100644 --- a/Reconstruction/Jet/JetRec/JetRec/PseudoJetTranslator.h +++ b/Reconstruction/Jet/JetRec/JetRec/PseudoJetTranslator.h @@ -15,14 +15,14 @@ public: PseudoJetTranslator(bool saveArea, bool saveArea4Vec) : m_saveArea(saveArea), m_saveArea4Vec(saveArea4Vec) {} - xAOD::Jet * translate(const fastjet::PseudoJet& pj, - const PseudoJetContainer& pjCont, - xAOD::JetContainer& jetCont) const ; + xAOD::Jet& translate(const fastjet::PseudoJet& pj, + const PseudoJetContainer& pjCont, + xAOD::JetContainer& jetCont) const ; - xAOD::Jet * translate(const fastjet::PseudoJet& pj, - const PseudoJetContainer& pjCont, - xAOD::JetContainer& jetCont, - const xAOD::Jet &parent) const ; + xAOD::Jet& translate(const fastjet::PseudoJet& pj, + const PseudoJetContainer& pjCont, + xAOD::JetContainer& jetCont, + const xAOD::Jet &parent) const ; diff --git a/Reconstruction/Jet/JetRec/Root/JetClusterer.cxx b/Reconstruction/Jet/JetRec/Root/JetClusterer.cxx index e7ac8847ef1a8fc1dd0eccac20c9ff6cdd604428..fba227864d178fa7510fbb259ad40682dfbd8ca8 100644 --- a/Reconstruction/Jet/JetRec/Root/JetClusterer.cxx +++ b/Reconstruction/Jet/JetRec/Root/JetClusterer.cxx @@ -14,7 +14,6 @@ #include "xAODJet/JetAuxContainer.h" #include "JetEDM/FastJetUtils.h" -#include "JetEDM/PseudoJetVector.h" #include "JetRec/PseudoJetTranslator.h" @@ -25,15 +24,15 @@ namespace JetClustererHelper { //const xAOD::EventInfo* pevinfo = handle.cptr(); auto ievt = ei->eventNumber(); auto irun = ei->runNumber(); - + if ( ei->eventType(xAOD::EventInfo::IS_SIMULATION)) { // For MC, use the channel and MC event number ievt = ei->mcEventNumber(); irun = ei->mcChannelNumber(); - } + } seeds.push_back(ievt); seeds.push_back(irun); - } + } } StatusCode JetClusterer::initialize() { @@ -50,17 +49,22 @@ StatusCode JetClusterer::initialize() { ATH_MSG_ERROR("Invalid jet size parameter: " << m_jetrad); return StatusCode::FAILURE; } - + // buld an empty ClusterSequence, just for the fastjet splash screen to appear during initialization (?) fastjet::JetDefinition jetdef(m_fjalg, m_jetrad); PseudoJetVector empty; fastjet::ClusterSequence cs(empty, jetdef); cs.inclusive_jets(m_ptmin); - + // Input DataHandles + if( !m_finalPseudoJets.empty() ) { + ATH_MSG_WARNING("A non-empty value was found for the FinalPseudoJets WriteHandleKey -- this will be ignored!"); + } + ATH_CHECK( m_eventinfokey.initialize() ); ATH_CHECK( m_inputPseudoJets.initialize() ); - m_finalPseudoJets = m_inputPseudoJets.key() +"FinalPJ"; + m_finalPseudoJets = name() + "FinalPJ"; + ATH_CHECK( m_finalPseudoJets.initialize() ); return StatusCode::SUCCESS; } @@ -78,7 +82,7 @@ std::pair<std::unique_ptr<xAOD::JetContainer>, std::unique_ptr<SG::IAuxStore> > return nullreturn; } - // Build the container to be returned + // Build the container to be returned // Avoid memory leaks with unique_ptr auto jets = std::make_unique<xAOD::JetContainer>(); auto auxCont = std::make_unique<xAOD::JetAuxContainer>(); @@ -94,7 +98,7 @@ std::pair<std::unique_ptr<xAOD::JetContainer>, std::unique_ptr<SG::IAuxStore> > bool useArea = m_ghostarea <= 0 ; if ( useArea ) { ATH_MSG_DEBUG("Creating input cluster sequence"); - clSequence = new fastjet::ClusterSequence(*pseudoJetVector, jetdef); + clSequence = new fastjet::ClusterSequence(*pseudoJetVector, jetdef); } else { // Prepare ghost area specifications ------------- ATH_MSG_DEBUG("Creating input area cluster sequence"); @@ -104,7 +108,7 @@ std::pair<std::unique_ptr<xAOD::JetContainer>, std::unique_ptr<SG::IAuxStore> > else {return nullreturn;} } - + // ----------------------- // Build a new pointer to a PseudoJetVector containing the final PseudoJet // This allows us to own the vector of PseudoJet which we will put in the evt store. @@ -116,7 +120,7 @@ std::pair<std::unique_ptr<xAOD::JetContainer>, std::unique_ptr<SG::IAuxStore> > msg() << " Pseudojet with pt " << std::setprecision(4) << pj.Et()*1e-3 << " has " << pj.constituents().size() << " constituents" << endmsg; } } - + // Let fastjet deal with deletion of ClusterSequence, so we don't need to also put it in the EventStore. clSequence->delete_self_when_unused(); @@ -128,19 +132,19 @@ std::pair<std::unique_ptr<xAOD::JetContainer>, std::unique_ptr<SG::IAuxStore> > PseudoJetTranslator pjTranslator(useArea, useArea); for (const fastjet::PseudoJet & pj: *pjVector ) { // create the xAOD::Jet from the PseudoJet, doing the signal & ghost constituents extraction - xAOD::Jet* jet = pjTranslator.translate(pj, *pjContHandle, *jets); - + xAOD::Jet& jet = pjTranslator.translate(pj, *pjContHandle, *jets); + // Add the PseudoJet onto the xAOD jet. Maybe we should do it in the above JetFromPseudojet call ?? - pjAccessor(*jet) = &pj; + pjAccessor(jet) = &pj; - jet->setInputType( xAOD::JetInput::Type( (int) m_inputType) ) ; + jet.setInputType( xAOD::JetInput::Type( (int) m_inputType) ) ; xAOD::JetAlgorithmType::ID ialg = xAOD::JetAlgorithmType::algId(m_fjalg); - jet->setAlgorithmType(ialg); - jet->setSizeParameter((float)m_jetrad); - if(useArea) jet->setAttribute(xAOD::JetAttribute::JetGhostArea, (float)m_ghostarea); - - ATH_MSG_VERBOSE( " xAOD::Jet with pt " << std::setprecision(4) << jet->pt()*1e-3 << " has " << jet->getConstituents().size() << " constituents" ); - ATH_MSG_VERBOSE( " Leading constituent is of type " << jet->getConstituents()[0].rawConstituent()->type()); + jet.setAlgorithmType(ialg); + jet.setSizeParameter((float)m_jetrad); + if(useArea) jet.setAttribute(xAOD::JetAttribute::JetGhostArea, (float)m_ghostarea); + + ATH_MSG_VERBOSE( " xAOD::Jet with pt " << std::setprecision(4) << jet.pt()*1e-3 << " has " << jet.getConstituents().size() << " constituents" ); + ATH_MSG_VERBOSE( " Leading constituent is of type " << jet.getConstituents()[0].rawConstituent()->type()); } // ------------------------------------- @@ -163,9 +167,9 @@ fastjet::AreaDefinition JetClusterer::buildAreaDefinition(bool & seedsok) const fastjet::GhostedAreaSpec gspec(5.0, 1, m_ghostarea); seedsok = true; - + if ( m_ranopt == 1 ) { - // Use run/event number as random number seeds. + // Use run/event number as random number seeds. auto evtInfoHandle = SG::makeHandle(m_eventinfokey); if (!evtInfoHandle.isValid()){ ATH_MSG_ERROR("Unable to retrieve event info"); @@ -178,7 +182,7 @@ fastjet::AreaDefinition JetClusterer::buildAreaDefinition(bool & seedsok) const gspec.set_random_status(inseeds); } - + ATH_MSG_DEBUG("Active area specs:"); ATH_MSG_DEBUG(" Requested ghost area: " << m_ghostarea); ATH_MSG_DEBUG(" Actual ghost area: " << gspec.actual_ghost_area()); diff --git a/Reconstruction/Jet/JetRec/Root/JetGroomer.cxx b/Reconstruction/Jet/JetRec/Root/JetGroomer.cxx index 44aefcb547a5bd9e113ea5d13be5ddac2c67b8ab..65138a0b6222fda95d80f2f654a301bcaf3b9164 100644 --- a/Reconstruction/Jet/JetRec/Root/JetGroomer.cxx +++ b/Reconstruction/Jet/JetRec/Root/JetGroomer.cxx @@ -4,6 +4,7 @@ #include "JetRec/JetGroomer.h" +#include "fastjet/PseudoJet.hh" #include "JetRec/PseudoJetContainer.h" using xAOD::JetContainer; @@ -12,47 +13,73 @@ StatusCode JetGroomer::initialize() { ATH_MSG_DEBUG("Initializing..."); - if(m_inputJetsKey.empty()){ + if(m_inputJetContainer.empty()){ ATH_MSG_ERROR("Jet grooming requested with no input ungroomed jets"); return StatusCode::FAILURE; - } else if(m_inputPseudoJetsKey.empty()){ + } else if(m_inputPseudoJets.empty()){ ATH_MSG_ERROR("Jet grooming requested with no input pseudojets"); return StatusCode::FAILURE; } else{ - ATH_CHECK(m_inputJetsKey.initialize()); - ATH_CHECK(m_inputPseudoJetsKey.initialize()); - ATH_CHECK(m_groomer.retrieve()); + if(!m_finalPseudoJets.empty()){ + ATH_MSG_WARNING("A non-empty value was found for the FinalPseudoJets WriteHandleKey -- this will be ignored!"); + } + + ATH_CHECK(m_inputJetContainer.initialize()); + ATH_CHECK(m_inputPseudoJets.initialize()); + m_finalPseudoJets = name() + "FinalPJ"; + ATH_CHECK(m_finalPseudoJets.initialize()); } return StatusCode::SUCCESS; } +// Common operations for any jet groomer std::pair<std::unique_ptr<xAOD::JetContainer>,std::unique_ptr<SG::IAuxStore> > JetGroomer::getJets() const { + // Return this in case of any problems + auto nullreturn = std::make_pair(std::unique_ptr<xAOD::JetContainer>(nullptr), std::unique_ptr<SG::IAuxStore>(nullptr)); + // ----------------------- // retrieve input - SG::ReadHandle<JetContainer> inputJets(m_inputJetsKey); + SG::ReadHandle<xAOD::JetContainer> jetContHandle(m_inputJetContainer); + if(!jetContHandle.isValid()) { + ATH_MSG_ERROR("No valid JetContainer with key "<< m_inputJetContainer.key()); + return nullreturn; + } - if(inputJets.isValid()) { - ATH_MSG_DEBUG("Retrieval of ungroomed jets succeeded"); - } else { - ATH_MSG_ERROR("Retrieval of ungroomed jets failed"); - return std::make_pair(std::unique_ptr<xAOD::JetContainer>(nullptr),std::unique_ptr<SG::IAuxStore>(nullptr)); + SG::ReadHandle<PseudoJetContainer> pjContHandle(m_inputPseudoJets); + if(!pjContHandle.isValid()) { + ATH_MSG_ERROR("No valid PseudoJetContainer with key "<< m_inputPseudoJets.key()); + return nullreturn; } - ATH_MSG_DEBUG("Grooming jets"); - // Create the output containers - // Make sure that memory is managed safely - auto outjets = std::make_unique<xAOD::JetContainer>(); - auto outaux = std::make_unique<xAOD::JetAuxContainer>(); - outjets->setStore(outaux.get()); + // Build the container to be returned + // Avoid memory leaks with unique_ptr + auto groomedJets = std::make_unique<xAOD::JetContainer>(); + auto auxCont = std::make_unique<xAOD::JetAuxContainer>(); + groomedJets->setStore(auxCont.get()); - SG::ReadHandle<PseudoJetContainer> inputPseudoJets(m_inputPseudoJetsKey); + // ----------------------- + // Build a new pointer to a PseudoJetVector containing the final groomed PseudoJets + // This allows us to own the vector of PseudoJet which we will put in the evt store. + // Thus the contained PseudoJet will be kept frozen there and we can safely use pointer to them from the xAOD::Jet objects + auto groomPJVector = std::make_unique<PseudoJetVector>( ); + groomPJVector->resize( jetContHandle->size() ); + + // loop over input jets + for (const xAOD::Jet* parentJet: *jetContHandle){ + // Child will create a groomed jet and insert it into the output container + this->insertGroomedJet(*parentJet, *pjContHandle, *groomedJets, *groomPJVector); + } - for (const xAOD::Jet* jet : *inputJets){ - m_groomer->groom(*jet, *inputPseudoJets, *outjets); + // ------------------------------------- + // record final PseudoJetVector + SG::WriteHandle<PseudoJetVector> pjVectorHandle(m_finalPseudoJets); + if(!pjVectorHandle.record(std::move(groomPJVector))){ + ATH_MSG_ERROR("Can't record PseudoJetVector under key "<< m_finalPseudoJets); + return nullreturn; } - return std::make_pair(std::move(outjets),std::move(outaux)); + return std::make_pair(std::move(groomedJets),std::move(auxCont)); } diff --git a/Reconstruction/Jet/JetRec/Root/JetTrimmer.cxx b/Reconstruction/Jet/JetRec/Root/JetTrimmer.cxx index 507e04ecf793816c78453aade05f05810189d17c..5c52857d1e88f09c83153957663d9084b3d21e02 100644 --- a/Reconstruction/Jet/JetRec/Root/JetTrimmer.cxx +++ b/Reconstruction/Jet/JetRec/Root/JetTrimmer.cxx @@ -73,8 +73,8 @@ int JetTrimmer::groom(const xAOD::Jet& jin, // Add jet to collection. xAOD::Jet* pjet = m_bld->add(pjtrim, pjContainer, jets, &jin); pjet->setAttribute<int>("TransformType", xAOD::JetTransform::Trim); - pjet->setAttribute("RClus", m_rclus); - pjet->setAttribute("PtFrac", m_ptfrac); + pjet->setAttribute("RClus", float(m_rclus)); + pjet->setAttribute("PtFrac", float(m_ptfrac)); pjet->setAttribute<int>("NTrimSubjets", nptrim); ATH_MSG_DEBUG("Properties after trimming:"); ATH_MSG_DEBUG(" ncon: " << pjtrim.constituents().size() << "/" diff --git a/Reconstruction/Jet/JetRec/Root/PseudoJetTranslator.cxx b/Reconstruction/Jet/JetRec/Root/PseudoJetTranslator.cxx index c194e9044adf4a0b8898c392c6514862d25952ed..d5ae49d2096e4c0845a69add3e71b138fc58114f 100644 --- a/Reconstruction/Jet/JetRec/Root/PseudoJetTranslator.cxx +++ b/Reconstruction/Jet/JetRec/Root/PseudoJetTranslator.cxx @@ -4,53 +4,54 @@ #include "JetRec/PseudoJetTranslator.h" -xAOD::Jet* PseudoJetTranslator::translate(const fastjet::PseudoJet& pj, +xAOD::Jet& PseudoJetTranslator::translate(const fastjet::PseudoJet& pj, const PseudoJetContainer& pjCont, xAOD::JetContainer& jetCont) const { - xAOD::Jet* jet = new xAOD::Jet(); - jetCont.push_back(jet); - jet->setJetP4( xAOD::JetFourMom_t( pj.pt(), pj.eta(), pj.phi(), pj.m() ) ); + // Create a new jet in place at the end of the container + jetCont.emplace_back(new xAOD::Jet()); + xAOD::Jet& jet = *jetCont.back(); + jet.setJetP4( xAOD::JetFourMom_t( pj.pt(), pj.eta(), pj.phi(), pj.m() ) ); - static SG::AuxElement::Accessor<const fastjet::PseudoJet*> pjAccessor("PseudoJet"); - pjAccessor(*jet) = &pj; + const static SG::AuxElement::Accessor<const fastjet::PseudoJet*> pjAccessor("PseudoJet"); + pjAccessor(jet) = &pj; // Record the jet-finding momentum, i.e. the one used to find/groom the jet. - jet->setJetP4(xAOD::JetConstitScaleMomentum, jet->jetP4()); + jet.setJetP4(xAOD::JetConstitScaleMomentum, jet.jetP4()); // save area if needed --------- if( pj.has_area() ){ - if(m_saveArea) jet->setAttribute(xAOD::JetAttribute::ActiveArea,pj.area()); + if(m_saveArea) jet.setAttribute(xAOD::JetAttribute::ActiveArea,pj.area()); if(m_saveArea4Vec){ fastjet::PseudoJet pja = pj.area_4vector(); xAOD::JetFourMom_t fvarea(pja.pt(), pja.eta(), pja.phi(), pja.m()); - jet->setAttribute(xAOD::JetAttribute::ActiveArea4vec, fvarea); + jet.setAttribute(xAOD::JetAttribute::ActiveArea4vec, fvarea); } }// area ------------- - pjCont.extractConstituents(*jet, pj); + pjCont.extractConstituents(jet, pj); return jet; } -xAOD::Jet* PseudoJetTranslator::translate(const fastjet::PseudoJet& pj, +xAOD::Jet& PseudoJetTranslator::translate(const fastjet::PseudoJet& pj, const PseudoJetContainer& pjCont, xAOD::JetContainer& jetCont, const xAOD::Jet &parent ) const { - xAOD::Jet * jet = translate(pj, pjCont, jetCont); + xAOD::Jet& jet = translate(pj, pjCont, jetCont); const xAOD::JetContainer* parentCont = dynamic_cast<const xAOD::JetContainer*>(parent.container()); if ( parentCont == 0 ) { return jet ;} // can this happen? if so THIS IS an ERROR ! should do something ElementLink<xAOD::JetContainer> el(*parentCont, parent.index()); - jet->setAttribute("Parent", el); + jet.setAttribute("Parent", el); - jet->setInputType(parent.getInputType()); - jet->setAlgorithmType(parent.getAlgorithmType()); - jet->setSizeParameter(parent.getSizeParameter()); - jet->setConstituentsSignalState(parent.getConstituentsSignalState()); - jet->setAttribute(xAOD::JetAttribute::JetGhostArea, parent.getAttribute<float>(xAOD::JetAttribute::JetGhostArea)); + jet.setInputType(parent.getInputType()); + jet.setAlgorithmType(parent.getAlgorithmType()); + jet.setSizeParameter(parent.getSizeParameter()); + jet.setConstituentsSignalState(parent.getConstituentsSignalState()); + jet.setAttribute(xAOD::JetAttribute::JetGhostArea, parent.getAttribute<float>(xAOD::JetAttribute::JetGhostArea)); return jet; } diff --git a/Reconstruction/Jet/JetRec/share/JetRecAlgTestCfg.py b/Reconstruction/Jet/JetRec/share/JetRecAlgTestCfg.py index ad9f2e7cc7553c5f160a5cc769766a10c38b25a0..9232884db3d359d1aeb11019c4c070e5eaff954c 100755 --- a/Reconstruction/Jet/JetRec/share/JetRecAlgTestCfg.py +++ b/Reconstruction/Jet/JetRec/share/JetRecAlgTestCfg.py @@ -32,8 +32,8 @@ ConfigFlags.lock() from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator # Get a ComponentAccumulator setting up the fundamental Athena job -from AthenaConfiguration.MainServicesConfig import MainServicesCfg -cfg=MainServicesCfg(ConfigFlags) +from AthenaConfiguration.MainServicesConfig import MainServicesCfg +cfg=MainServicesCfg(ConfigFlags) # Add the components for reading in pool files from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg @@ -120,7 +120,7 @@ def JetBuildAlgCfg(ConfigFlags,buildjetsname): "pjmergealg_"+buildjetsname, InputPJContainers = pjcs, OutputContainer = "PseudoJetMerged_"+buildjetsname) - + buildcfg.addEventAlgo(mergepjalg) # Create the JetClusterer, set some standard options @@ -162,15 +162,14 @@ def JetGroomAlgCfg(ConfigFlags,buildjetsname,groomjetsname): # Create the JetGroomer, provide it with a JetTrimmer jtrim = CompFactory.JetTrimming("trimSmallR2Frac5",RClus=0.2,PtFrac=0.05) - jtrim.InputJetContainer = buildjetsname - jtrim.InputPseudoJets = "PseudoJetMerged_"+buildjetsname - jtrim.TrimmedOutputPseudoJets = "PseudoJet"+groomjetsname + jtrim.UngroomedJets = buildjetsname + jtrim.ParentPseudoJets = "PseudoJetMerged_"+buildjetsname # Create the JetRecAlg, configure it to use the builder # using constructor syntax instead # (equivalent to setting properties with "=") jra = CompFactory.JetRecAlg( - "JRA_groom", + "JRA_trim", Provider = jtrim, # Single ToolHandle Modifiers = [], # ToolHandleArray OutputContainer = groomjetsname) diff --git a/Reconstruction/Jet/JetRec/src/JetTrimming.cxx b/Reconstruction/Jet/JetRec/src/JetTrimming.cxx index 445674e9ffa5bb41b0eaa2e33b0cf4ad78c771b5..6379a4d8d0f6e41b49d2b966542f669eda932e13 100644 --- a/Reconstruction/Jet/JetRec/src/JetTrimming.cxx +++ b/Reconstruction/Jet/JetRec/src/JetTrimming.cxx @@ -7,124 +7,77 @@ #include "fastjet/PseudoJet.hh" #include "fastjet/JetDefinition.hh" #include "fastjet/Selector.hh" -#include "fastjet/tools/Filter.hh" #include "JetEDM/PseudoJetVector.h" #include "JetRec/PseudoJetTranslator.h" -using fastjet::PseudoJet; -using xAOD::JetContainer; - - + // tool implementing the operations to convert fastjet::PseudoJet -> xAOD::Jet +const static PseudoJetTranslator s_pjTranslator(false, false); // false => do not save jet areas //********************************************************************** StatusCode JetTrimming::initialize() { + ATH_CHECK( JetGroomer::initialize() ); + + // Unfortunately not possible to do this because there is no + // declareProperty overload for CheckedProperty in AsgTool + // + // Enforce upper/lower limits of Gaudi::Properties + // m_rclus.verifier().setBounds(0.,10.); + // m_ptfrac.verifier().setBounds(0.,1.); if ( m_rclus < 0.0 || m_rclus > 10.0 ) { - ATH_MSG_ERROR("Invalid value for RClus " << m_rclus); + ATH_MSG_ERROR("Invalid value " << m_rclus << "for RClus. Allowable range is 0-10"); return StatusCode::FAILURE; } if ( m_ptfrac < 0.0 || m_ptfrac > 1.0 ) { - ATH_MSG_ERROR("Invalid value for PtFrac " << m_ptfrac); + ATH_MSG_ERROR("Invalid value " << m_ptfrac << "for PtFrac. Allowable range is 0-1"); return StatusCode::FAILURE; } - ATH_CHECK( m_inputPseudoJets.initialize() ); - ATH_CHECK( m_inputJetContainer.initialize() ); + // Define the trimmer + m_trimmer = std::make_unique<fastjet::Filter>(fastjet::JetDefinition(fastjet::kt_algorithm, m_rclus), + fastjet::SelectorPtFractionMin(m_ptfrac)); ATH_MSG_INFO(" Recluster R: " << m_rclus); ATH_MSG_INFO(" pT faction min: " << m_ptfrac); + return StatusCode::SUCCESS; } //********************************************************************** -std::pair<std::unique_ptr<xAOD::JetContainer>, std::unique_ptr<SG::IAuxStore> > JetTrimming::getJets() const { - // Return this in case of any problems - auto nullreturn = std::make_pair(std::unique_ptr<xAOD::JetContainer>(nullptr), std::unique_ptr<SG::IAuxStore>(nullptr)); - - // ----------------------- - // retrieve input - SG::ReadHandle<xAOD::JetContainer> jetContHandle(m_inputJetContainer); - if(!jetContHandle.isValid()) { - ATH_MSG_ERROR("No valid JetContainer with key "<< m_inputJetContainer.key()); - return nullreturn; - } - - SG::ReadHandle<PseudoJetContainer> pjContHandle(m_inputPseudoJets); - if(!pjContHandle.isValid()) { - ATH_MSG_ERROR("No valid PseudoJetContainer with key "<< m_inputPseudoJets.key()); - return nullreturn; - } - - // Build the container to be returned - // Avoid memory leaks with unique_ptr - auto trimmedJets = std::make_unique<xAOD::JetContainer>(); - auto auxCont = std::make_unique<xAOD::JetAuxContainer>(); - trimmedJets->setStore(auxCont.get()); - - - // The trimmer - fastjet::Filter trimmer(fastjet::JetDefinition(fastjet::kt_algorithm, m_rclus), - fastjet::SelectorPtFractionMin(m_ptfrac)); - - // tool implementing the operations to convert fastjet::PseudoJet -> xAOD::Jet - PseudoJetTranslator pjTranslator(false, false); // false => do not save jet areas - - // ----------------------- - // Build a new pointer to a PseudoJetVector containing the final trimmed PseudoJets - // This allows us to own the vector of PseudoJet which we will put in the evt store. - // Thus the contained PseudoJet will be kept frozen there and we can safely use pointer to them from the xAOD::Jet objects - auto trimPJVector = std::make_unique<PseudoJetVector>( ); - trimPJVector->resize( jetContHandle->size() ); - // Accessors to retrieve and set the pointers to PseudoJet - static SG::AuxElement::ConstAccessor<const fastjet::PseudoJet*> pjConstAcc("PseudoJet"); - static SG::AuxElement::Accessor<const fastjet::PseudoJet*> pjAcc("PseudoJet"); - - - // loop over input jets - int count = 0; - for (const xAOD::Jet* parentJet: *jetContHandle){ - // retrieve the PseudoJet from the parent : - const PseudoJet * parentPJ = pjConstAcc(*parentJet); - - // Trim : - PseudoJet trimmedPJ = trimmer(*parentPJ) ; - (*trimPJVector)[count] = trimmedPJ; // save a *copy* of this trimmed PJ - - ATH_MSG_VERBOSE(" Input cluster sequence: " << parentPJ->associated_cluster_sequence()); - ATH_MSG_VERBOSE(" Trimmed cluster sequence: " << trimmedPJ.associated_cluster_sequence()); - - // build the xAOD::Jet from the PseudoJet, and put it in the container - xAOD::Jet* jet = pjTranslator.translate(trimmedPJ, *pjContHandle, *trimmedJets, *parentJet); - - // decorate with the pointer to the PJ we keep in the evt store. - pjAcc(*jet) = & (*trimPJVector)[count] ; - - int nptrim = trimmedPJ.pieces().size(); - jet->setAttribute<int>(xAOD::JetAttribute::TransformType, xAOD::JetTransform::Trim); - jet->setAttribute<int>(xAOD::JetAttribute::NTrimSubjets, nptrim); - // Need to convert from GaudiProperty - jet->setAttribute(xAOD::JetAttribute::RClus, float(m_rclus)); - jet->setAttribute(xAOD::JetAttribute::PtFrac, float(m_ptfrac)); - ATH_MSG_DEBUG("Properties after trimming:"); - ATH_MSG_DEBUG(" ncon: " << trimmedPJ.constituents().size() << "/" - << parentPJ->constituents().size()); - ATH_MSG_DEBUG(" nsub: " << nptrim); - count++; - } - - // ------------------------------------- - // record final PseudoJetVector - SG::WriteHandle<PseudoJetVector> pjVectorHandle(m_finalPseudoJets); - if(!pjVectorHandle.record(std::move(trimPJVector))){ - ATH_MSG_ERROR("Can't record PseudoJetVector under key "<< m_finalPseudoJets); - return nullreturn; - } - - // Return the jet container and aux, use move to transfer - // ownership of pointers to caller - return std::make_pair(std::move(trimmedJets), std::move(auxCont)); +void JetTrimming::insertGroomedJet(const xAOD::Jet& parentjet, const PseudoJetContainer& inpjcont, xAOD::JetContainer& outcont, PseudoJetVector& trimpjvec) const { + + const static SG::AuxElement::Accessor<const fastjet::PseudoJet*> s_pjAcc("PseudoJet"); + const static SG::AuxElement::ConstAccessor<const fastjet::PseudoJet*> s_pjConstAcc("PseudoJet"); + + // retrieve the PseudoJet from the parent : + const fastjet::PseudoJet& parentPJ = *s_pjConstAcc(parentjet); + + // Trim : + fastjet::PseudoJet trimmedPJ = m_trimmer->result(parentPJ) ; + ATH_MSG_VERBOSE(" Input cluster sequence: " << parentPJ.associated_cluster_sequence()); + ATH_MSG_VERBOSE(" Trimmed cluster sequence: " << trimmedPJ.associated_cluster_sequence()); + + // build the xAOD::Jet from the PseudoJet, and put it in the container + xAOD::Jet& jet = s_pjTranslator.translate(trimmedPJ, inpjcont, outcont, parentjet); + // The vector is resized externally to match the jet container size, + // so just fill the corresponding entry + trimpjvec[jet.index()] = trimmedPJ; // save a *copy* of this trimmed PJ + + // decorate with the pointer to the PJ we keep in the evt store. + s_pjAcc(jet) = & trimpjvec[jet.index()]; + + int nptrim = trimmedPJ.pieces().size(); + jet.setAttribute<int>(xAOD::JetAttribute::TransformType, xAOD::JetTransform::Trim); + jet.setAttribute<int>(xAOD::JetAttribute::NTrimSubjets, nptrim); + // Need to convert from GaudiProperty + jet.setAttribute(xAOD::JetAttribute::RClus, float(m_rclus)); + jet.setAttribute(xAOD::JetAttribute::PtFrac, float(m_ptfrac)); + ATH_MSG_VERBOSE("Properties after trimming:"); + ATH_MSG_VERBOSE(" ncon: " << trimmedPJ.constituents().size() << "/" + << parentPJ.constituents().size()); + ATH_MSG_VERBOSE(" nsub: " << nptrim); } diff --git a/Reconstruction/Jet/JetRec/src/JetTrimming.h b/Reconstruction/Jet/JetRec/src/JetTrimming.h index faa7d192ac173db6b800ade3607c9b25ed5701fe..4302969abe17fc70af7f247b3efafc6493f26855 100644 --- a/Reconstruction/Jet/JetRec/src/JetTrimming.h +++ b/Reconstruction/Jet/JetRec/src/JetTrimming.h @@ -15,47 +15,35 @@ /// running the trimming procedure on each jet of the parent JetContainer the tool is acting on. /// Obviously the parent JetContainer must be present in the evt store, but also the input PseudoJetContainer from which it has been built. -#include "xAODEventInfo/EventInfo.h" -#include "StoreGate/ReadHandleKey.h" -#include "xAODJet/JetAuxContainer.h" +#include "fastjet/PseudoJet.hh" +#include "fastjet/tools/Filter.hh" -#include "JetInterface/IJetProvider.h" -#include "AsgTools/AsgTool.h" +#include "xAODJet/JetContainer.h" +#include "JetInterface/IJetProvider.h" +#include "JetRec/JetGroomer.h" #include "JetRec/PseudoJetContainer.h" - class JetTrimming -: public asg::AsgTool, - virtual public JetProvider<xAOD::JetAuxContainer> { +: virtual public JetGroomer { ASG_TOOL_CLASS(JetTrimming, IJetProvider) public: - using asg::AsgTool::AsgTool; + using JetGroomer::JetGroomer; - StatusCode initialize() override; + StatusCode initialize() override final; - /// Return the final jets with their aux store. - /// Can return a pair of null pointers if an error occurs. - std::pair<std::unique_ptr<xAOD::JetContainer>, std::unique_ptr<SG::IAuxStore> > getJets() const override; + virtual void insertGroomedJet(const xAOD::Jet&, const PseudoJetContainer&, xAOD::JetContainer&, PseudoJetVector&) const override final; +private: -protected: + // The filter object that will apply the grooming + std::unique_ptr<fastjet::Filter> m_trimmer; - - /// Handle Input JetContainer (this contains the parent ungroomed jets to be trimmed) - SG::ReadHandleKey<xAOD::JetContainer> m_inputJetContainer {this, "InputJetContainer", "ungroomedinput", "Input ungroomed jet container"}; - - /// This is the input to the parent JetContainer. It is needed in order to re-assign the ghost constituents from the final trimmed PJ to the xAOD::Jet - SG::ReadHandleKey<PseudoJetContainer> m_inputPseudoJets {this, "InputPseudoJets", "inputpseudojet", "input constituents of parent JetContainer"}; - - // Job options. Gaudi::Property<float> m_rclus {this, "RClus", 0.3 , "R for reclustering (0 for none)"}; Gaudi::Property<float> m_ptfrac {this, "PtFrac", 0.03, "pT fraction for retaining subjets"}; - Gaudi::Property<std::string> m_finalPseudoJets {this, "TrimmedOutputPseudoJets", "undef", "Key to save the final trimmed pj. Necessary in order to ensure a valid PJ can be linked from xAOD::Jet"}; - }; diff --git a/Reconstruction/Jet/JetRec/src/components/JetRec_entries.cxx b/Reconstruction/Jet/JetRec/src/components/JetRec_entries.cxx index 87ea18b902ceac02071a4eaa10fbbc5250129603..bc1988eef2711cb477c5cedbee0cf4f921bb7184 100644 --- a/Reconstruction/Jet/JetRec/src/components/JetRec_entries.cxx +++ b/Reconstruction/Jet/JetRec/src/components/JetRec_entries.cxx @@ -27,7 +27,6 @@ #include "JetRec/JetConstitRemover.h" #include "JetRec/JetClusterer.h" #include "JetRec/JetCopier.h" -#include "JetRec/JetGroomer.h" DECLARE_COMPONENT( JetToolRunner ) DECLARE_COMPONENT( JetRecTool ) @@ -48,7 +47,6 @@ DECLARE_COMPONENT( JetPseudojetCopier ) DECLARE_COMPONENT( JetConstitRemover ) DECLARE_COMPONENT( JetClusterer ) DECLARE_COMPONENT( JetCopier ) -DECLARE_COMPONENT( JetGroomer ) DECLARE_COMPONENT( PseudoJetMerger ) DECLARE_COMPONENT( JetAlgorithm )