ParticleSelectionAlg.icc

Go to the documentation of this file.

// Dear emacs, this is -*- c++ -*-
 
/*
  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
 
 
#include "xAODBase/IParticle.h"
#include "xAODCore/AuxContainerBase.h"
#include "AthContainers/ConstDataVector.h"
#include "xAODEventInfo/EventInfo.h"
#include "xAODCutFlow/CutBookkeeperContainer.h"
#include "PATCore/AcceptInfo.h"
#include "PATCore/AcceptData.h"
 
 
 
template<class CONT, class AUXCONT>
StatusCode ParticleSelectionAlg::selectParticles(const xAOD::IParticleContainer* inContainer,
                                                 const std::vector<int>& resultVec) const
{
  ATH_MSG_DEBUG ("selectParticlesStepTwo<CONT,AUXCONT>(...) " << name() << "...");
 
  // Get the type of particle that the current container CONT has
  typedef typename CONT::base_value_type PART;
 
  // Create a new output container and its associated auxiliary container
  CONT* outCont = nullptr;
  ConstDataVector<CONT>* outContConst = nullptr;
 
  if ( m_outOwnPolicy == SG::OWN_ELEMENTS ) {
    outCont = new CONT( m_outOwnPolicy );
    ATH_CHECK( evtStore()->record ( outCont, m_outCollKey.value() ) );
    if ( m_writeSplitAux.value() ) {
      xAOD::AuxContainerBase* outAuxContainer = new xAOD::AuxContainerBase();
      ATH_CHECK( evtStore()->record( outAuxContainer, m_outCollKey.value() + "Aux." ) );
      outCont->setStore( outAuxContainer );
      ATH_MSG_VERBOSE("Recorded xAOD::AuxContainerBase for container: " << m_outCollKey.value() );
    }
    else {
      AUXCONT* outAuxContainer = new AUXCONT();
      ATH_CHECK( evtStore()->record( outAuxContainer, m_outCollKey.value() + "Aux." ) );
      outCont->setStore( outAuxContainer );
      ATH_MSG_VERBOSE("Recorded AuxContainer for container: " << m_outCollKey.value() );
    }
  }
  else {
    outContConst = new ConstDataVector<CONT>( m_outOwnPolicy );
    ATH_CHECK( evtStore()->record ( outContConst, m_outCollKey.value() ) );
    ATH_MSG_VERBOSE("Recorded ConstDataVector for container: " << m_outCollKey.value() );
  }
  // ATH_MSG_VERBOSE("StoreGate dump: " << evtStore()->dump() );
 
  // Get the CutBookkeeperContainer, if cut-flow bookkeeping was requested
  double eventWeight = 1.0;
  xAOD::CutBookkeeperContainer* cutBKCont = nullptr;
  const xAOD::EventInfo* evtInfo = nullptr;
  if (m_doCutFlow){
    ATH_CHECK( outputMetaStore()->retrieve(cutBKCont, m_cutBKCName.value() ) );
    // Also, get the event info as we need to get the MC-weight for this event
    ATH_CHECK( evtStore()->retrieve( evtInfo, m_evtInfoName.value() ));
    const bool isSim = evtInfo->eventType(xAOD::EventInfo::EventType::IS_SIMULATION);
    if (isSim){
      eventWeight = static_cast<double>( evtInfo->mcEventWeight() );
      ATH_MSG_VERBOSE("Got MC event weight of " << eventWeight );
    }
  }
 
 
  // Now, loop over the input container and check which particles to write out
  ATH_MSG_DEBUG("Input container has size " << inContainer->size() );
  for ( std::size_t i=0; i<inContainer->size(); ++i ) {
    const xAOD::IParticle* partBase = inContainer->at(i);
    const PART* part = static_cast<const PART*>(partBase);
 
    // ================================
    // Apply the selection tools
    // ================================
    // The default object pass value. As soon as one cut is not passed, then
    // the remainder of the cuts will not even be tried
    bool passEverything = true;
 
    //------------- for the AsgSelectionTools --------------
    // Loop over all selection tools
    ATH_MSG_VERBOSE("Loop over all selection tools");
    for ( std::size_t toolIdx=0; toolIdx < m_selTools.size(); ++toolIdx ){
      if (passEverything){
        ATH_MSG_VERBOSE("Now going to try AsgSelectionTools number " << toolIdx );
        const asg::AcceptData acceptData = m_selTools[toolIdx]->accept(part);
        if (!m_doCutFlow){ passEverything &= static_cast<bool>(acceptData); }
        else {
          const std::size_t cbkStartIdx = m_selToolIdxOffset[toolIdx];
          const unsigned int nCuts = acceptData.getNCuts();
          for ( unsigned int iCut=0; iCut<nCuts; ++iCut ){
            passEverything &= acceptData.getCutResult(iCut);
            if (passEverything){
              const std::size_t currentCBKIdx = cbkStartIdx + iCut;
              xAOD::CutBookkeeper* cutBK = cutBKCont->at(currentCBKIdx);
              cutBK->addNAcceptedEvents(1);
              cutBK->addSumOfEventWeights(eventWeight);
              cutBK->addSumOfEventWeightsSquared(eventWeight*eventWeight);
            }
          }
        } // Done doing detailed particle cut-flow for this tool
        ATH_MSG_VERBOSE("AsgSelectionTools number " << toolIdx << " passed/failed: " << passEverything );
      }
    }
 
    //------------- for the ExpressionParsing in this algorithm --------------
    ATH_MSG_VERBOSE("Looking at expression parser result: " << passEverything);
    if ( passEverything && !(resultVec.empty()) ){
      // If this particle is not accepted by the expression parser, go to the next one
      ATH_MSG_VERBOSE("Now going to try expression '" << m_selection.value() << "'" );
      passEverything &= static_cast<bool>(resultVec.at(i));
      ATH_MSG_VERBOSE("Expression '" << m_selection.value() << "' passed/failed: " << passEverything
                      << ", particle index=" << part->index() <<", pt=" << 0.001*(part->pt()) << " GeV, eta="
                      << part->eta() << ", phi=" << part->phi() );
      if (passEverything && m_doCutFlow){
        xAOD::CutBookkeeper* cutBK = cutBKCont->at(m_idxSelParster);
        cutBK->addNAcceptedEvents(1);
        cutBK->addSumOfEventWeights(eventWeight);
        cutBK->addSumOfEventWeightsSquared(eventWeight*eventWeight);
      }
    }
 
    // Now, if all cuts are passed, write out the current particle to the
    // output container
    if (passEverything){
      ATH_MSG_VERBOSE("Going to fill output container with particle with pt=" << 0.001*(part->pt())
                      << " GeV, eta=" << part->eta() );
      if ( m_outOwnPolicy == SG::OWN_ELEMENTS ) {
        PART* newPart = new PART();
        outCont->push_back(newPart);
        *newPart = *part;
      }
      else {
        outContConst->push_back(part);
      }
    }
 
  } // End: Loop over input particles
 
  return StatusCode::SUCCESS;
}