DeltaRTool.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
*/
 
// Author: James Catmore (james.catmore@cern.ch)
//
 
#include "DerivationFrameworkTools/DeltaRTool.h"
#include <utility> //for std::pair
 
namespace DerivationFramework {
 
  StatusCode DeltaRTool::initialize()
  {
    if (m_sgName.key().empty()) {
      ATH_MSG_ERROR("No SG name provided for the output of invariant mass tool!");
      return StatusCode::FAILURE;
    }
    ATH_CHECK(m_sgName.initialize());
    ATH_CHECK(m_containerName.initialize());
 
    if (!m_containerName2.key().empty()) {
      ATH_CHECK(m_containerName2.initialize());
    }
 
    ATH_CHECK(initializeParser( {m_expression, m_2ndExpression} ));
    return StatusCode::SUCCESS;
  }
 
  StatusCode DeltaRTool::finalize()
  {
    ATH_CHECK(finalizeParser());
    return StatusCode::SUCCESS;
  }
 
  StatusCode DeltaRTool::addBranches(const EventContext& ctx) const
  {
    // Write deltaRs to SG for access by downstream algs
    if (evtStore()->contains<std::vector<float> >(m_sgName.key())) { // FIXME Use Handles
      ATH_MSG_ERROR("Tool is attempting to write a StoreGate key " << m_sgName << " which already exists. Please use a different key");
      return StatusCode::FAILURE;
    }
    std::unique_ptr<std::vector<float> > deltaRs(new std::vector<float>());
    ATH_CHECK(getDeltaRs(deltaRs.get(), ctx ));
 
    SG::WriteHandle<std::vector<float> > writeHandle(m_sgName, ctx);
    ATH_CHECK(writeHandle.record(std::move(deltaRs)));
 
    return StatusCode::SUCCESS;
  }
 
  StatusCode DeltaRTool::getDeltaRs(std::vector<float>* deltaRs, const EventContext& ctx) const
  {
 
    // check the relevant information is available
    if (m_containerName.key().empty()) {
      ATH_MSG_WARNING("Input container missing - returning zero");
      deltaRs->push_back(0.0);
      return StatusCode::FAILURE;
    }
    bool secondContainer(false);
    if (!m_containerName2.key().empty()) secondContainer=true;
 
    // get the relevant branches
    SG::ReadHandle<xAOD::IParticleContainer> particles{m_containerName, ctx};
 
    const xAOD::IParticleContainer* secondParticles(nullptr);
    if (secondContainer) {
      SG::ReadHandle<xAOD::IParticleContainer> particleHdl2{m_containerName2, ctx};
      secondParticles=particleHdl2.cptr();
    }
 
    // get the positions of the elements which pass the requirement
    std::vector<int> entries, entries2;
    if (!m_expression.empty()) {entries = m_parser[kDeltaRToolParser1]->evaluateAsVector();}
    else {entries.assign(particles->size(),1);} // default: include all elements
    unsigned int nEntries = entries.size();
    // check the sizes are compatible
    if (particles->size() != nEntries ) {
      ATH_MSG_FATAL("Branch sizes incompatible");
      return StatusCode::FAILURE;
    }
    unsigned int nEntries2(0);
    if (secondContainer) {
      if (!m_2ndExpression.empty()) {entries2 =  m_parser[kDeltaRToolParser2]->evaluateAsVector();}
      else {entries2.assign(secondParticles->size(),1);} // default: include all elements
      nEntries2 = entries2.size();
      // check the sizes are compatible
      if (secondParticles->size() != nEntries2 ) {
        ATH_MSG_FATAL("Branch sizes incompatible - returning zero");
        return StatusCode::FAILURE;
      }
    }
 
    // Double loop to get the pairs for which the mass should be calculated
    std::vector<std::pair<unsigned, unsigned>> pairs;
    if (!secondContainer) {
      for (unsigned outerIt=0; outerIt<nEntries; ++outerIt) {
        for (unsigned innerIt=outerIt+1; innerIt<nEntries; ++innerIt) {
          if (entries[outerIt]==1 && entries[innerIt]==1) {
            pairs.emplace_back(outerIt, innerIt);
          }
        }
      }
    }
 
    if (secondContainer) {
      for (unsigned coll1It=0; coll1It<nEntries; ++coll1It) {
        for (unsigned coll2It=0; coll2It<nEntries2; ++coll2It) {
          //coverity[copy_paste_error]
          if (entries[coll1It]==1 && entries2[coll2It]==1) {
            pairs.emplace_back(coll1It, coll2It);
          }
        }
      }
    }
 
    // Loop over the pairs; calculate the mass; put into vector and return
    for (const auto & [first, second] : pairs) {
      if (!secondContainer) {
        float phi1f = ((*particles)[first])->p4().Phi(); float phi2f = ((*particles)[second])->p4().Phi();
        float eta1f = ((*particles)[first])->p4().Eta(); float eta2f = ((*particles)[second])->p4().Eta();
        float deltaR = calculateDeltaR(phi1f,phi2f,eta1f,eta2f);
        deltaRs->push_back(deltaR);
      }
      if (secondContainer) {
        float phi1f = ((*particles)[first])->p4().Phi(); float phi2f = ((*secondParticles)[second])->p4().Phi();
        float eta1f = ((*particles)[first])->p4().Eta(); float eta2f = ((*secondParticles)[second])->p4().Eta();
        float deltaR = calculateDeltaR(phi1f,phi2f,eta1f,eta2f);
        deltaRs->push_back(deltaR);
      }
    }
 
    return StatusCode::SUCCESS;
 
  }
 
  float DeltaRTool::calculateDeltaR(float phi1, float phi2, float eta1, float eta2)
  {
    float deltaPhi = fabs(phi1-phi2);
    if (deltaPhi>TMath::Pi()) deltaPhi = 2.0*TMath::Pi() - deltaPhi;
    float deltaPhiSq = deltaPhi*deltaPhi;
    float deltaEtaSq = (eta1-eta2)*(eta1-eta2);
    float deltaR = sqrt(deltaPhiSq+deltaEtaSq);
    return deltaR;
  }
}