AsgElectronSelectorTool.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
*/

 
// Include this class's header
#include "ElectronPhotonSelectorTools/AsgElectronSelectorTool.h"
#include "ElectronPhotonSelectorTools/ElectronSelectorHelpers.h"
#include "EgammaAnalysisHelpers/AsgEGammaConfigHelper.h"
#include "EGSelectorConfigurationMapping.h"
#include "ElectronDNNCalculator.h"
 
 
//EDM includes
#include "xAODEgamma/Electron.h"
#include "xAODCaloEvent/CaloCluster.h"
#include "TEnv.h"
 
#include "AsgTools/CurrentContext.h"
#include "PathResolver/PathResolver.h"
// STL includes
#include <string>
#include <cstdint>
#include <cmath>
#include <unordered_map>
 
namespace AllowedVariables
{
  enum e {
    eta,
    et,
    f3,
    Rhad,
    Rhad1,
    Reta,
    weta2,
    f1,
    Eratio,
    deltaEta1,
    d0,
    qd0,
    d0significance,
    Rphi,
    dPOverP,
    deltaPhiRescaled2,
    trans_TRTPID,
    wtots1,
    EoverP,
    nPixHitsPlusDeadSensors,
    nSCTHitsPlusDeadSensors,
    SCTWeightedCharge
  };
 
  const std::unordered_map<std::string, int> variableMap = {
    {"eta", eta},
    {"et", et},
    {"f3", f3},
    {"Rhad", Rhad},
    {"Rhad1", Rhad1},
    {"Reta", Reta},
    {"weta2", weta2},
    {"f1", f1},
    {"Eratio", Eratio},
    {"deltaEta1", deltaEta1},
    {"d0", d0},
    {"qd0", qd0},
    {"d0significance", d0significance},
    {"Rphi", Rphi},
    {"dPOverP", dPOverP},
    {"deltaPhiRescaled2", deltaPhiRescaled2},
    {"trans_TRTPID", trans_TRTPID},
    {"wtots1", wtots1},
    {"EoverP", EoverP},
    {"nPixHitsPlusDeadSensors", nPixHitsPlusDeadSensors},
    {"nSCTHitsPlusDeadSensors", nSCTHitsPlusDeadSensors},
    {"SCTWeightedCharge", SCTWeightedCharge},
  };
};
 
  bool AsgElectronSelectorTool::SelectionLayer::initialize(
    TEnv& env,
    const std::string& prefixName,
    const std::string& layerName,
    const std::string& userModelFileName,
    const std::string& userQuantileFileName,
    asg::AsgTool* parentTool
) {
    prefix = prefixName;
    name = layerName;
    unsigned int numberOfExpectedBinCombinedMVA = AsgElectronSelectorTool::s_fnDiscEtBins * AsgElectronSelectorTool::s_fnDiscEtaBins;
    unsigned int numberOfExpectedEtaBins = AsgElectronSelectorTool::s_fnDiscEtBins;
    // 1. Getting quantiles and model files 
    std::string modelFile = !userModelFileName.empty()
        ? userModelFileName
        : env.GetValue((prefix + "inputModelFileName").c_str(), "");
 
    std::string quantileFile = !userQuantileFileName.empty()
        ? userQuantileFileName
        : env.GetValue((prefix + "inputQuantileFileName").c_str(), "");
 
    std::string modelPath = PathResolverFindCalibFile(modelFile);
    std::string quantilePath = PathResolverFindCalibFile(quantileFile);
 
    if (modelPath.empty() || quantilePath.empty()) {
        return false;
    }
 
    // 2. 
    m_multiClass = env.GetValue((prefix + "multiClass").c_str(), false);
    m_nOutputs = env.GetValue((prefix + "n_outputs").c_str(), m_multiClass ? 6 : 1);
    // Include cf node in numerator or denominator when combining different outputs
    m_cfSignal = env.GetValue((prefix + "cfSignal").c_str(), false);
 
    // 3. 
    std::string varStr = env.GetValue((prefix + "Variables").c_str(), "");
    std::stringstream varsStream(varStr);
    std::string substr;
 
    std::vector<std::string> tempVariables;
    m_enum_variables.clear();
 
    while (std::getline(varsStream, substr, ',')) {
        if (substr.empty()) continue;
        if (!prefix.empty() && substr.rfind(prefix, 0) == 0) {
            substr = substr.substr(prefix.length());
        }
        if (!AllowedVariables::variableMap.contains(substr)) {
            return false;
        }
        tempVariables.push_back(substr);
        m_enum_variables.push_back(AllowedVariables::variableMap.at(substr));
    }
 
    if (tempVariables.empty()) return false;
 
    m_cutAmbiguity = AsgConfigHelper::HelperInt((prefix + "CutAmbiguity").c_str(), env);
    m_cutBL        = AsgConfigHelper::HelperInt((prefix + "CutBL").c_str(), env);
    m_cutPi        = AsgConfigHelper::HelperInt((prefix + "CutPi").c_str(), env);
    m_cutSCT       = AsgConfigHelper::HelperInt((prefix + "CutSCT").c_str(), env);
    m_fractions    = AsgConfigHelper::HelperDouble((prefix + "Fractions").c_str(), env);
    m_cutSelector  = AsgConfigHelper::HelperDouble((prefix + "CutSelector").c_str(), env);
    m_cutSelectorCF = AsgConfigHelper::HelperDouble((prefix + "CutSelectorCF").c_str(), env);
    
 
 
    if (m_cutSelector.size() != numberOfExpectedBinCombinedMVA){
      parentTool->msg(MSG::ERROR) <<"Configuration issue :  cutSelector expected size " << numberOfExpectedBinCombinedMVA <<
                    " input size " << m_cutSelector.size();
      return false;
    }
    if (!m_cutSelectorCF.empty()){
      m_CFReject = true;
      if (m_cutSelectorCF.size() != numberOfExpectedBinCombinedMVA){
        parentTool->msg(MSG::ERROR) <<"Configuration issue :  cutSelectorCF expected size " << numberOfExpectedBinCombinedMVA <<
                      " input size " << m_cutSelectorCF.size();
        return false;
      }
      if(!m_multiClass){
        parentTool->msg(MSG::ERROR) <<"Configuration issue :  CF rejection is only defined "
                      "for multiClass: TRUE";
        return false;
      }
      } else {
      m_CFReject = false;
      }
 
    if (m_multiClass){
      unsigned int nFractionsPerBin = m_nOutputs - 1;
      unsigned int expectedFractionsSize = numberOfExpectedEtaBins * nFractionsPerBin;
 
      if (m_fractions.size() != expectedFractionsSize) {
          parentTool->msg(MSG::ERROR) <<"Configuration issue: multiclass is TRUE, but the number of fractions (" 
                        << m_fractions.size() << ") does not match " << m_nOutputs << " classes (" 
                        << expectedFractionsSize << " fractions expected).";
          return false;
      }
    }
      
    m_mvaTool = std::make_unique<ElectronDNNCalculator>(
        dynamic_cast<AsgElectronSelectorTool*>(parentTool),
        modelPath.c_str(),
        quantilePath.c_str(),
        tempVariables
    );
 
    m_defaultOutputs = std::vector<float>(m_nOutputs, -999.0f);
    return true;
}
 
std::vector<float> AsgElectronSelectorTool::SelectionLayer::evaluate(const std::vector<double>& inputs) const {
    if (!m_mvaTool) return m_defaultOutputs;
    auto scores = m_mvaTool->calculate(inputs);
    std::vector<float> res(scores.rows());
    for (int i = 0; i < scores.rows(); ++i) {
        res[i] = scores(i, 0);
    }
    return res;
}
 
 
//=============================================================================
// Standard constructor
//=============================================================================
AsgElectronSelectorTool::AsgElectronSelectorTool( const std::string& myname ) :
  AsgTool(myname),
  m_configFile{""}
{
 
  // Declare the needed properties
  declareProperty("WorkingPoint", m_workingPoint="", "The Working Point");
  declareProperty("ConfigFile", m_configFile="", "The config file to use");
 
  // model file name. Managed in the ElectronDNNCalculator.
  declareProperty("inputModelFileName", m_modelFileName="", "The input file name that holds the model" );
  // QuantileTransformer file name ( required for preprocessing ). Managed in the ElectronDNNCalculator.
  declareProperty("quantileFileName", m_quantileFileName="", "The input file name that holds the QuantileTransformer");
  // especially for  trigger electron
  declareProperty("skipDeltaPoverP",m_skipDeltaPoverP = false,"If true, it will skip the check of deltaPoverP");
 
  declareProperty("skipAmbiguityCut",m_skipAmbiguityCut = false,"If true, it will skip the ambiguity cut");
}
 
 
//=============================================================================
// Standard destructor
//=============================================================================
AsgElectronSelectorTool::~AsgElectronSelectorTool()
= default;
 
 
//=============================================================================
// Asgena initialize method
//=============================================================================
StatusCode AsgElectronSelectorTool::initialize()
{
  
if (m_configFile.empty()) {
  if (!m_workingPoint.empty()){
    m_configFile = AsgConfigHelper::findConfigFile(m_workingPoint, EgammaSelectors::ElectronDNNPointToConfFile);
    ATH_MSG_INFO("operating point : " << this->getOperatingPointName());
  }
}
else {
  ATH_MSG_INFO("Using USER provided config file: " << m_configFile);
}
  if (m_configFile.empty()){
    ATH_MSG_ERROR("Could not find configuration file " << m_configFile);
    return StatusCode::FAILURE;
  }
 
  std::string configFile = PathResolverFindCalibFile(m_configFile);
 ATH_MSG_INFO("CONFIG FILE (logical): " << m_configFile);
 ATH_MSG_INFO("CONFIG FILE (resolved full path): " << configFile);
 
 if (configFile.empty()){
    ATH_MSG_ERROR("Could not locate " << m_configFile);
    return StatusCode::FAILURE;
  }
 
  ATH_MSG_DEBUG("Configfile to use: " << m_configFile);
  TEnv env;
  env.ReadFile(configFile.c_str(), kEnvLocal);
 
  int nLayers = env.GetValue("nSelectionLayers", 1);
 
    const std::vector<std::string> defaultPrefixes = {
        "",                     // Layer 0: No prefix 
        "SecondSelection.",     // Layer 1
        "ThirdSelection.",      // Layer 2
        "FourthSelection."      // Layer 3
    };
 
    m_layers.clear();  
    for (int i = 0; i < nLayers; ++i) {
        std::string prefix = (i < static_cast<int>(defaultPrefixes.size()))
            ? defaultPrefixes[i]
            : ("Layer" + std::to_string(i + 1) + ".");
 
        std::string layerName = (i == 0) ? "FirstSelection"
            : (i == 1) ? "SecondSelection"
            : ("Layer" + std::to_string(i + 1));
 
        SelectionLayer layer;
        std::string userModel = (i == 0) ? m_modelFileName : "";
        std::string userQuantile = (i == 0) ? m_quantileFileName : "";
 
        if (!layer.initialize(env, prefix, layerName, userModel, userQuantile, this)) {
            ATH_MSG_ERROR("Error initializing layer " << i << " with prefix: " << prefix);
            return StatusCode::FAILURE;
        }
        m_layers.push_back(std::move(layer));
    }
 
  // do smooth interpolation between bins
  m_doSmoothBinInterpolation = env.GetValue("doSmoothBinInterpolation", false);
 
  // --------------------------------------------------------------------------
  // Register the cuts and check that the registration worked:
  // NOTE: THE ORDER IS IMPORTANT!!! Cut0 corresponds to bit 0, Cut1 to bit 1,...
  // use an int as a StatusCode
  int sc(1);
 
  // Cut position for the kineatic pre-selection
  m_cutPosition_kinematic = m_acceptMVA.addCut("kinematic", "pass kinematic");
  if (m_cutPosition_kinematic < 0) sc = 0;
 
  // NSilicon
  m_cutPosition_NSilicon = m_acceptMVA.addCut("NSCT", "pass NSCT");
  if (m_cutPosition_NSilicon < 0) sc = 0;
 
  // NPixel
  m_cutPosition_NPixel = m_acceptMVA.addCut("NPixel", "pass NPixel");
  if (m_cutPosition_NPixel < 0) sc = 0;
 
  // NBlayer
  m_cutPosition_NBlayer = m_acceptMVA.addCut("NBlayer", "pass NBlayer");
  if (m_cutPosition_NBlayer < 0) sc = 0;
 
  // Ambiguity
  m_cutPosition_ambiguity = m_acceptMVA.addCut("ambiguity", "pass ambiguity");
  if (m_cutPosition_ambiguity < 0) sc = 0;
 
 
  // Cut position for the likelihood selection - DO NOT CHANGE ORDER!
  m_cutPosition_MVA = m_acceptMVA.addCut("passMVA", "pass MVA");
  if (m_cutPosition_MVA < 0) sc = 0;
 
  // Check that we got everything OK
  if (sc == 0){
    ATH_MSG_ERROR("ERROR: Something went wrong with the setup of the decision objects...");
    return StatusCode::FAILURE;
  }
 
  ATH_MSG_INFO("Initialization SUCCESS");
  return StatusCode::SUCCESS;
}
 
 
//=============================================================================
// return the accept info object
//=============================================================================
 
const asg::AcceptInfo& AsgElectronSelectorTool::getAcceptInfo() const
{
  return m_acceptMVA;
}
 
//=============================================================================
// The main accept method: the actual cuts are applied here
//=============================================================================
asg::AcceptData AsgElectronSelectorTool::accept( const EventContext& ctx, const xAOD::Electron* eg, double mu ) const
{
  ATH_MSG_VERBOSE("\t AsgElectronSelectorTool::accept( &ctx, *eg, mu= "<<(&ctx)<<", "<<eg<<", "<<mu<<" )");
 
  // Setup return accept with AcceptInfo
  asg::AcceptData acceptData(&m_acceptMVA);
 
  if (!eg){
    throw std::runtime_error("AsgElectronSelectorTool: Failed, no electron object was passed");
  }
 
  const xAOD::CaloCluster* cluster = eg->caloCluster();
  if (!cluster){
    ATH_MSG_DEBUG("exiting because cluster is NULL " << cluster);
    return acceptData;
  }
 
  if(!cluster->hasSampling(CaloSampling::CaloSample::EMB2) && !cluster->hasSampling(CaloSampling::CaloSample::EME2)){
    ATH_MSG_DEBUG("Failed, cluster is missing samplings EMB2 and EME2");
    return acceptData;
  }
 
  const double energy = cluster->e();
  const float eta = cluster->etaBE(2);
 
  if(isForwardElectron(eg, eta)){
    ATH_MSG_DEBUG("Failed, this is a forward electron! The AsgElectronSelectorTool is only suitable for central electrons!");
    return acceptData;
  }
 
  const xAOD::TrackParticle* track = eg->trackParticle();
  if (!track){
    ATH_MSG_DEBUG("exiting because track is NULL " << track);
    return acceptData;
  }
 
  // transverse energy of the electron (using the track eta)
  double et = (std::cosh(track->eta()) != 0.) ? energy / std::cosh(track->eta()) : 0.;
 
  // number of track hits
  uint8_t nSiHitsPlusDeadSensors(0);
  uint8_t nPixHitsPlusDeadSensors(0);
  bool passBLayerRequirement(false);
  uint8_t ambiguityBit(0);
 
  bool allFound = true;
  std::string notFoundList = "";
 
  // get the ambiguity type from the decoration
  if (!m_skipAmbiguityCut){
    static const SG::Accessor<uint8_t> ambiguityTypeAcc("ambiguityType");
    if (ambiguityTypeAcc.isAvailable(*eg)) {
      ambiguityBit = ambiguityTypeAcc(*eg);
    }
    else {
      allFound = false;
      notFoundList += "ambiguityType ";
    }
  }
 
  nSiHitsPlusDeadSensors = ElectronSelectorHelpers::numberOfSiliconHitsAndDeadSensors(*track);
  nPixHitsPlusDeadSensors = ElectronSelectorHelpers::numberOfPixelHitsAndDeadSensors(*track);
  passBLayerRequirement = ElectronSelectorHelpers::passBLayerRequirement(*track);
 
  // calculate the output of the selector tool
  SelectorOutputs mvaOutputs = runSelections(ctx, eg, mu);
 
  if (!allFound){
    throw std::runtime_error("AsgElectronSelectorTool: Not all variables needed for the decision are found. The following variables are missing: " + notFoundList );
  }
 
  // Set up the individual cuts
  bool passKine(true);
  bool passNSilicon(true);
  bool passNPixel(true);
  bool passNBlayer(true);
  bool passAmbiguity(true);
  bool passMVA(true);
 
 
  if (std::abs(eta) > 2.47){
    ATH_MSG_DEBUG("This electron is fabs(eta)>2.47 Returning False.");
    passKine = false;
  }
 
  unsigned int etBin = getDiscEtBin(et);
  unsigned int etaBin = getDiscEtaBin(eta);
 
  // sanity
  if (etBin >= s_fnDiscEtBins){
    ATH_MSG_DEBUG("Cannot evaluate model for Et " << et << ". Returning false..");
    passKine = false;
  }
 
  // Return if the kinematic requirements are not fulfilled
  acceptData.setCutResult(m_cutPosition_kinematic, passKine);
  if (!passKine){return acceptData;}
 
for (size_t idx = 0; idx < m_layers.size(); ++idx) {
    const auto& layer = m_layers[idx];
    const auto& scores = mvaOutputs.layerOutputs[idx];
 
  // ambiguity bit
  if (!layer.m_cutAmbiguity.empty()){
    if (!ElectronSelectorHelpers::passAmbiguity((xAOD::AmbiguityTool::AmbiguityType)ambiguityBit, layer.m_cutAmbiguity[etaBin])){
      ATH_MSG_DEBUG("MVA macro: ambiguity Bit cut failed.");
      passAmbiguity = false;
    }
  }
 
  // blayer cut
  if (!layer.m_cutBL.empty()) {
    if(layer.m_cutBL[etaBin] == 1 && !passBLayerRequirement){
      ATH_MSG_DEBUG("MVA macro: Blayer cut failed.");
      passNBlayer = false;
    }
  }
  // pixel cut
  if (!layer.m_cutPi.empty()){
    if (nPixHitsPlusDeadSensors < layer.m_cutPi[etaBin]){
      ATH_MSG_DEBUG("MVA macro: Pixels Failed.");
      passNPixel = false;
    }
  }
  // SCT cut
  if (!layer.m_cutSCT.empty()){
    if (nSiHitsPlusDeadSensors < layer.m_cutSCT[etaBin]){
      ATH_MSG_DEBUG( "MVA macro: Silicon Failed.");
      passNSilicon = false;
    }
  }
  unsigned int ibin_combinedMVA = etBin*s_fnDiscEtaBins+etaBin; // Must change if number of eta bins changes!.
  
  // First cut on the CF discriminant
  // If empty, continue only with prompt ID
  if (!layer.m_cutSelectorCF.empty()){
    double cutDiscriminantCF;
    // To protect against a binning mismatch, which should never happen
    if (ibin_combinedMVA >= layer.m_cutSelectorCF.size()){
      throw std::runtime_error("AsgElectronSelectorTool: The desired eta/pt bin is outside of the range specified by the input. This should never happen! This indicates a mismatch between the binning in the configuration file and the tool implementation." );
    }
    if (m_doSmoothBinInterpolation){
      cutDiscriminantCF = interpolateCuts(layer.m_cutSelectorCF, et, eta);
    }
    else{
      cutDiscriminantCF = layer.m_cutSelectorCF.at(ibin_combinedMVA);
    }
    // Determine if the calculated mva score value passes the combined cut
    ATH_MSG_DEBUG("MVA macro: CF Discriminant: ");
    double mvaScoreCF = combineOutputsCF(scores);
    if (mvaScoreCF < cutDiscriminantCF){
      ATH_MSG_DEBUG("MVA macro: CF cut failed.");
      passMVA = false;
    }
  }
  // (Second) cut on prompt discriminant
  // Cut on Score MVA First Selection
  if (!layer.m_cutSelector.empty()) {
    double cutDiscriminant;
    // To protect against a binning mismatch, which should never happen
    if (ibin_combinedMVA >= layer.m_cutSelector.size()){
      throw std::runtime_error("AsgElectronSelectorTool: The desired eta/pt bin is outside of the range specified by the input. This should never happen! This indicates a mismatch between the binning in the configuration file and the tool implementation." );
    }
    if (m_doSmoothBinInterpolation){
      cutDiscriminant = interpolateCuts(layer.m_cutSelector, et, eta);
    }
    else{
      cutDiscriminant = layer.m_cutSelector.at(ibin_combinedMVA);
    }
    
    // Determine if the calculated mva score value passes the combined cut
    ATH_MSG_DEBUG("MVA macro: Prompt Discriminant: ");
    double disc = layer.m_multiClass
            ? combineOutputs(scores, eta, layer.m_fractions, layer.m_cfSignal)
            : transformMLOutput(scores.at(0));
 
        if (disc < cutDiscriminant) {
            ATH_MSG_DEBUG("MVA macro: Prompt cut failed in layer " << layer.name);
            passMVA = false;
        }
  }
    if (!passMVA || !passNSilicon || !passNPixel || !passNBlayer || !passAmbiguity) {
        break;
    }
} 
  // Set the individual cut bits in the return object
  acceptData.setCutResult(m_cutPosition_NSilicon, passNSilicon);
  acceptData.setCutResult(m_cutPosition_NPixel, passNPixel);
  acceptData.setCutResult(m_cutPosition_NBlayer, passNBlayer);
  acceptData.setCutResult(m_cutPosition_ambiguity, passAmbiguity);
  acceptData.setCutResult(m_cutPosition_MVA, passMVA);
 
  return acceptData;
 
}
 
//=============================================================================
// The main result method: the actual mvaScore is calculated here
//=============================================================================
double AsgElectronSelectorTool::calculate( const EventContext& ctx, const xAOD::Electron* eg, double mu ) const
{
  // Get all outputs of the mva tool
  SelectorOutputs mvaOutputs = runSelections(ctx, eg, mu);
  if (mvaOutputs.layerOutputs.empty()) return -999.0;
 
    auto discScores = getDiscriminant(mvaOutputs, eg);
    return discScores.empty() ? -999.0 : discScores.back(); //Returns final layer score
}
 
std::vector<double> AsgElectronSelectorTool::getDiscriminant(const SelectorOutputs& mvaOutputs, const xAOD::Electron* eg ) const
{
  std::vector<double> discriminants;
  const xAOD::CaloCluster* cluster = eg->caloCluster();
  const float eta = cluster->etaBE(2);
  // If a binary model is used, vector will have one entry
 
  for (size_t i = 0; i < mvaOutputs.layerOutputs.size(); ++i) {
        const auto& scores = mvaOutputs.layerOutputs[i];
        if (scores.empty()) continue;
 
        const auto& layer = m_layers[i];
        double disc = layer.m_multiClass
            ? combineOutputs(scores, eta, layer.m_fractions, layer.m_cfSignal)
            : transformMLOutput(scores.at(0));
 
        discriminants.push_back(disc);
    }
 
  return discriminants;
}
 
std::vector<float> AsgElectronSelectorTool::calculateMultipleOutputs(const EventContext &ctx, const xAOD::Electron *eg, double mu) const
{
  SelectorOutputs mvaOutputs = runSelections(ctx, eg, mu);
  
  if (!mvaOutputs.layerOutputs.empty()) {
      return mvaOutputs.layerOutputs.back();
  }
  return {};
}
AsgElectronSelectorTool::SelectorOutputs AsgElectronSelectorTool::runSelections(const EventContext &ctx, const xAOD::Electron *eg, double mu) const
{
  
  SelectorOutputs mvaOutputs;
  ATH_MSG_VERBOSE("\t AsgElectronSelectorTool::calculateMultipleOutputs( &ctx, *eg, mu= "<<(&ctx)<<", "<<eg<<", "<<mu<<" )");
  if (!eg){
    throw std::runtime_error("AsgElectronSelectorTool: Failed, no electron object was passed" );
  }
 
  const xAOD::CaloCluster* cluster = eg->caloCluster();
  if (!cluster){
    ATH_MSG_DEBUG("Failed, no cluster.");
    // Return a default value
    for (const auto& l : m_layers) mvaOutputs.layerOutputs.push_back(l.m_defaultOutputs);
    return mvaOutputs;    
  }
 
  if (!cluster->hasSampling(CaloSampling::CaloSample::EMB2) && !cluster->hasSampling(CaloSampling::CaloSample::EME2)){
    ATH_MSG_DEBUG("Failed, cluster is missing samplings EMB2 and EME2.");
    // Return a default value
    for (const auto& l : m_layers) mvaOutputs.layerOutputs.push_back(l.m_defaultOutputs);
    return mvaOutputs;
  }
 
  const double energy = cluster->e();
  const float eta = cluster->etaBE(2);
 
  if (isForwardElectron(eg, eta)){
    ATH_MSG_DEBUG("Failed, this is a forward electron! The AsgElectronSelectorTool is only suitable for central electrons!");
    // Return a default value
    for (const auto& l : m_layers) mvaOutputs.layerOutputs.push_back(l.m_defaultOutputs);
    return mvaOutputs;    
  }
 
  const xAOD::TrackParticle* track = eg->trackParticle();
  if (!track){
    ATH_MSG_DEBUG("Failed, no track.");
    // Return a default value
    for (const auto& l : m_layers) mvaOutputs.layerOutputs.push_back(l.m_defaultOutputs);
    return mvaOutputs;    
  }
 
  // transverse energy of the electron (using the track eta)
  const double et = energy / std::cosh(track->eta());
 
  // Variables used in the ML model
  // track quantities
  double  SCTWeightedCharge(0.0);
  uint8_t nSCTHitsPlusDeadSensors(0);
  uint8_t nPixHitsPlusDeadSensors(0);
  float d0(0.0), d0sigma(0.0), d0significance(0.0), qd0(0.0);
  float trackqoverp(0.0);
  double dPOverP(0.0);
  float TRT_PID(0.0);
  double trans_TRTPID(0.0);
 
  // Track Cluster matching
  float deltaEta1(0), deltaPhiRescaled2(0), EoverP(0);
 
  // Calorimeter
  float Reta(0), Rphi(0), Rhad1(0), Rhad(0), w2(0), f1(0), Eratio(0), f3(0), wtots1(0);
 
  bool allFound = true;
  std::string notFoundList = "";
 
  // retrieve track variables
  trackqoverp = track->qOverP();
  d0 = track->d0();
  qd0 = (eg->charge())*track->d0();
  float vard0 = track->definingParametersCovMatrix()(0, 0);
  if (vard0 > 0){
    d0sigma = std::sqrt(vard0);
  }
  d0significance = (d0sigma == 0.) ? -99999. : std::abs(d0 / d0sigma);
 
  const static SG::Accessor<float> trans_TRT_PID_acc("transformed_e_probability_ht");
  if (!trans_TRT_PID_acc.isAvailable(*eg)) {
    // most probable case, need to compute the variable
 
    if (!track->summaryValue(TRT_PID, xAOD::eProbabilityHT)) {
      allFound = false;
      notFoundList += "eProbabilityHT ";
    }
 
    // Transform the TRT PID output for use in the LH tool.
    const double tau = 15.0;
    const double fEpsilon = 1.0e-30; // to avoid zero division
    double pid_tmp = TRT_PID;
    if (pid_tmp >= 1.0)
      pid_tmp = 1.0 - 1.0e-15; // this number comes from TMVA
    else if (pid_tmp <= fEpsilon)
      pid_tmp = fEpsilon;
    trans_TRTPID = -std::log(1.0 / pid_tmp - 1.0) * (1. / tau);
  }
  else
  {
    // it means the variable have been already computed by another tool
    // usually this is the EGammaVariableCorrection, which means that
    // it is also fudged (only MC)
    trans_TRTPID = trans_TRT_PID_acc(*eg);
  }
 
  //Change default value of TRT PID to 0.15 instead of 0 when there is no information from the TRT
  if ((std::abs(trans_TRTPID) < 1.0e-6) && (std::abs(eta) > 2.01)){
      trans_TRTPID = 0.15;
  }
 
  unsigned int index;
  if (track->indexOfParameterAtPosition(index, xAOD::LastMeasurement)){
    double refittedTrack_LMqoverp = track->charge() / std::sqrt(std::pow(track->parameterPX(index), 2) +
                                                                std::pow(track->parameterPY(index), 2) +
                                                                std::pow(track->parameterPZ(index), 2));
 
    dPOverP = 1 - trackqoverp / (refittedTrack_LMqoverp);
  }
  else if (!m_skipDeltaPoverP) {
    allFound = false;
    notFoundList += "deltaPoverP ";
  }
 
  EoverP =  energy * std::abs(trackqoverp);
 
  nPixHitsPlusDeadSensors = ElectronSelectorHelpers::numberOfPixelHitsAndDeadSensors(*track);
  nSCTHitsPlusDeadSensors = ElectronSelectorHelpers::numberOfSCTHitsAndDeadSensors(*track);
 
  float charge = 0;
  uint8_t SCT = 0;
  for (unsigned TPit = 0; TPit < eg->nTrackParticles(); TPit++) {
    uint8_t temp_NSCTHits = 0;
    if (eg->trackParticle(TPit)) {
      eg->trackParticle(TPit)->summaryValue(temp_NSCTHits, xAOD::numberOfSCTHits);
      SCT += temp_NSCTHits;
      charge += temp_NSCTHits*(eg->trackParticle(TPit)->charge());
    }
  }
  if (SCT)
    SCTWeightedCharge = (eg->charge()*charge/SCT);
  else {
    ATH_MSG_WARNING("No SCT hit for any track associated to electron ! nTP = " << eg->nTrackParticles());
  }
 
  // retrieve Calorimeter variables
  // reta = e237/e277
  if (!eg->showerShapeValue(Reta, xAOD::EgammaParameters::Reta)){
    allFound = false;
    notFoundList += "Reta ";
  }
  // rphi e233/e237
  if (!eg->showerShapeValue(Rphi, xAOD::EgammaParameters::Rphi)){
    allFound = false;
    notFoundList += "Rphi ";
  }
  // rhad1 = ethad1/et
  if (!eg->showerShapeValue(Rhad1, xAOD::EgammaParameters::Rhad1)){
    allFound = false;
    notFoundList += "Rhad1 ";
  }
  // rhad = ethad/et
  if (!eg->showerShapeValue(Rhad, xAOD::EgammaParameters::Rhad)){
    allFound = false;
    notFoundList += "Rhad ";
  }
  // shower width in 2nd sampling
  if (!eg->showerShapeValue(w2, xAOD::EgammaParameters::weta2)){
    allFound = false;
    notFoundList += "weta2 ";
  }
  // fraction of energy reconstructed in the 1st sampling
  if (!eg->showerShapeValue(f1, xAOD::EgammaParameters::f1)){
    allFound = false;
    notFoundList += "f1 ";
  }
  // E of 2nd max between max and min in strips
  if (!eg->showerShapeValue(Eratio, xAOD::EgammaParameters::Eratio)){
    allFound = false;
    notFoundList += "Eratio ";
  }
  // fraction of energy reconstructed in the 3rd sampling
  if (!eg->showerShapeValue(f3, xAOD::EgammaParameters::f3)){
    allFound = false;
    notFoundList += "f3 ";
  }
 
  // Set f3 to default value in eta region where it is poorly modelled
  if (std::abs(eta) > 2.01) {
    f3 = 0.05;
  }
 
  // Shower width in first sampling of the calorimeter
  if (!eg->showerShapeValue(wtots1, xAOD::EgammaParameters::wtots1)){
    allFound = false;
    notFoundList += "wtots1 ";
  }
 
  // retrieve Track Cluster matching variables
  // difference between cluster eta (sampling 1) and the eta of the track
  if (!eg->trackCaloMatchValue(deltaEta1, xAOD::EgammaParameters::deltaEta1)){
    allFound = false;
    notFoundList += "deltaEta1 ";
  }
  // difference between the cluster phi (sampling 2) and the phi of the track extrapolated from the last measurement point.
  if (!eg->trackCaloMatchValue(deltaPhiRescaled2, xAOD::EgammaParameters::deltaPhiRescaled2)){
    allFound = false;
    notFoundList += "deltaPhiRescaled2 ";
  }
 
 
  ATH_MSG_VERBOSE(Form("Vars: eta=%8.5f, et=%8.5f, f3=%8.5f, rHad==%8.5f, rHad1=%8.5f, Reta=%8.5f, w2=%8.5f, f1=%8.5f, Emaxs1=%8.5f, deltaEta1=%8.5f, d0=%8.5f, qd0=%8.5f, d0significance=%8.5f, Rphi=%8.5f, dPOverP=%8.5f, deltaPhiRescaled2=%8.5f, TRT_PID=%8.5f, trans_TRTPID=%8.5f, mu=%8.5f, wtots1=%8.5f, EoverP=%8.5f, nPixHitsPlusDeadSensors=%2df, nSCTHitsPlusDeadSensors=%2df, SCTWeightedCharge=%8.5f",
                       eta, et, f3, Rhad, Rhad1, Reta,
                       w2, f1, Eratio,
                       deltaEta1, d0, qd0,
                       d0significance,
                       Rphi, dPOverP, deltaPhiRescaled2,
                       TRT_PID, trans_TRTPID,
                       mu,
                       wtots1, EoverP, int(nPixHitsPlusDeadSensors), int(nSCTHitsPlusDeadSensors), SCTWeightedCharge));
 
  if (!allFound){
    throw std::runtime_error("AsgElectronSelectorTool: Not all variables needed for MVA calculation are found. The following variables are missing: " + notFoundList );
  }
 
  
  mvaOutputs.layerOutputs.reserve(m_layers.size());
  for (const auto& layer : m_layers) {
      std::vector<double> variableValues;
      variableValues.reserve(layer.m_enum_variables.size());
      for (const auto varID : layer.m_enum_variables) {
          
        switch(varID)
      {
        case AllowedVariables::eta:
          variableValues.push_back(std::abs(eta));  break; // TODO - rename to abseta?
        case AllowedVariables::et:
          variableValues.push_back(et); break;
        case AllowedVariables::f3:
          variableValues.push_back(f3); break;
        case AllowedVariables::Rhad:
          variableValues.push_back(Rhad); break;
        case AllowedVariables::Rhad1:
          variableValues.push_back(Rhad1); break;
        case AllowedVariables::Reta:
          variableValues.push_back(Reta); break;
        case AllowedVariables::weta2:
          variableValues.push_back(w2); break;
        case AllowedVariables::f1:
          variableValues.push_back(f1); break;
        case AllowedVariables::Eratio:
          variableValues.push_back(Eratio); break;
        case AllowedVariables::deltaEta1:
          variableValues.push_back(deltaEta1); break;
        case AllowedVariables::d0:
          variableValues.push_back(d0); break;
        case AllowedVariables::qd0:
          variableValues.push_back(qd0); break;
        case AllowedVariables::d0significance:
          variableValues.push_back(d0significance); break;
        case AllowedVariables::Rphi:
          variableValues.push_back(Rphi); break;
        case AllowedVariables::dPOverP:
          variableValues.push_back(dPOverP); break;
        case AllowedVariables::deltaPhiRescaled2:
          variableValues.push_back(deltaPhiRescaled2); break;
        case AllowedVariables::trans_TRTPID:
          variableValues.push_back(trans_TRTPID); break;
        case AllowedVariables::wtots1:
          variableValues.push_back(wtots1); break;
        case AllowedVariables::EoverP:
          variableValues.push_back(EoverP); break;
        case AllowedVariables::nPixHitsPlusDeadSensors:
          variableValues.push_back(nPixHitsPlusDeadSensors); break;
        case AllowedVariables::nSCTHitsPlusDeadSensors:
          variableValues.push_back(nSCTHitsPlusDeadSensors); break;
        case AllowedVariables::SCTWeightedCharge:
          variableValues.push_back(SCTWeightedCharge); break;
        default:
          // Handle unknown varID or error case
          throw std::runtime_error("AsgElectronSelectorTool: unknown variable "
            "index, something went wrong in initialization!" );
          break;
 
      }
    }
 
      mvaOutputs.layerOutputs.push_back(layer.evaluate(variableValues));
  }  
    
return mvaOutputs;
}
 
//=============================================================================
//=============================================================================
std::string AsgElectronSelectorTool::getOperatingPointName() const
{
  return m_workingPoint;
}
 
//=============================================================================
asg::AcceptData AsgElectronSelectorTool::accept( const xAOD::IParticle* part ) const
{
  return accept(Gaudi::Hive::currentContext(), part);
}
asg::AcceptData AsgElectronSelectorTool::accept( const EventContext& ctx, const xAOD::IParticle* part ) const
{
  ATH_MSG_VERBOSE("\t AsgElectronSelectorTool::accept( &ctx, *part= "<<(&ctx)<<", "<<part<<" )");
  const xAOD::Electron* eg = dynamic_cast<const xAOD::Electron*>(part);
  if (eg){
    return accept(ctx, eg);
  }
  else {
    ATH_MSG_DEBUG("AsgElectronSelectorTool::could not cast to const Electron");
    // Setup return accept with AcceptInfo
    asg::AcceptData acceptData(&m_acceptMVA);
    return acceptData;
  }
}
 
double AsgElectronSelectorTool::calculate( const EventContext& ctx, const xAOD::IParticle* part ) const
{
  ATH_MSG_VERBOSE("\t AsgElectronSelectorTool::calculate( &ctx, *part"<<(&ctx)<<", "<<part<<" )");
  const xAOD::Electron* eg = dynamic_cast<const xAOD::Electron*>(part);
  if (eg){
    return calculate(ctx, eg);
  }
  else {
    ATH_MSG_DEBUG("AsgElectronSelectorTool::could not cast to const Electron");
    // Return a default value
    return -999.;
  }
}
 
asg::AcceptData AsgElectronSelectorTool::accept( const EventContext& ctx, const xAOD::Egamma* eg, double mu ) const
{
  ATH_MSG_VERBOSE("\t AsgElectronSelectorTool::accept( &ctx, *eg, mu= "<<(&ctx)<<", "<<eg<<", "<<mu<<" )");
  const xAOD::Electron* ele = dynamic_cast<const xAOD::Electron*>(eg);
  if (ele){
    return accept(ctx, ele, mu);
  }
  else {
    ATH_MSG_DEBUG("AsgElectronSelectorTool::could not cast to const Electron");
    // Setup return accept with AcceptInfo
    asg::AcceptData acceptData(&m_acceptMVA);
    return acceptData;
  }
}
 
double AsgElectronSelectorTool::calculate( const EventContext& ctx, const xAOD::Egamma* eg, double mu ) const
{
  ATH_MSG_VERBOSE("\t AsgElectronSelectorTool::calculate( &ctx, *eg, mu= "<<(&ctx)<<", "<<eg<<", "<<mu<<" )");
  const xAOD::Electron* ele = dynamic_cast<const xAOD::Electron*>(eg);
  if (ele){
    return calculate(ctx, ele, mu);
  }
  else {
    ATH_MSG_DEBUG("AsgElectronSelectorTool::could not cast to const Electron");
    return -999.;
  }
}
 
bool AsgElectronSelectorTool::isForwardElectron( const xAOD::Egamma* eg, const float eta ) const
{
  static const SG::ConstAccessor< uint16_t > accAuthor( "author" );
 
  if (accAuthor.isAvailable(*eg)){
    // cannot just do eg->author() because it isn't always filled
    // at trigger level
    if (accAuthor(*eg) == xAOD::EgammaParameters::AuthorFwdElectron){
      ATH_MSG_DEBUG("Failed, this is a forward electron! The AsgElectronSelectorTool is only suitable for central electrons!");
      return true;
    }
  }
  else{
    //Check for fwd via eta range the old logic
    if (std::abs(eta) > 2.5){
      ATH_MSG_DEBUG("Failed, cluster->etaBE(2) range due to " << eta << " seems like a fwd electron" );
      return true;
    }
  }
 
  return false;
}
 
 
double AsgElectronSelectorTool::transformMLOutput( float score ) const
{
  // returns transformed or non-transformed output
  constexpr double oneOverTau = 1 ;
  constexpr double fEpsilon = 1.0e-30; // to avoid zero division
  if (score >= 1.0) score = 1.0 - 1.0e-15; // this number comes from TMVA
  else if (score <= fEpsilon) score = fEpsilon;
  //cppcheck-suppress invalidFunctionArg
  score = -std::log(1.0 / score - 1.0) * oneOverTau;
  ATH_MSG_DEBUG("score is " << score);
  return score;
}
 
 
double AsgElectronSelectorTool::combineOutputs( const std::vector<float>& mvaScores, double eta, const std::vector<double>& fractions, bool cfSignal ) const
{
  unsigned int etaBin = getDiscEtaBin(eta);
  double disc = 0;
  size_t nClasses = mvaScores.size();
  size_t nFractionsPerBin = nClasses - 1;
 
  if (cfSignal && nClasses >= 3) {
 
    double numerator = mvaScores.at(0) * (1.0 - fractions.at(nFractionsPerBin * etaBin + 0)) +
                       mvaScores.at(1) * fractions.at(nFractionsPerBin * etaBin + 0);
                       
    double denominator = 0.0;
    for (size_t i = 2; i < nClasses; ++i) {
      denominator += mvaScores.at(i) * fractions.at(nFractionsPerBin * etaBin + i - 1);
    }
    if (denominator != 0.0) {
      disc = numerator / denominator;
    }
  }
  else {
 
    double numerator = mvaScores.at(0);
    double denominator = 0.0;
    
    for (size_t i = 1; i < nClasses; ++i) {
      denominator += mvaScores.at(i) * fractions.at(nFractionsPerBin * etaBin + i - 1);
    }
    if (denominator != 0.0) {
      disc = numerator / denominator;
    }
  }
  // Log transform to have values in reasonable range
  return std::log(disc);
}
 
double AsgElectronSelectorTool::combineOutputsCF( const std::vector<float>& mvaScores)
{
  double disc = 0;
  disc = mvaScores.at(0) / mvaScores.at(1);
 
  return std::log(disc);
}
 
 
// Gets the Discriminant Eta bin [0,s_fnDiscEtaBins-1] given the eta
unsigned int AsgElectronSelectorTool::getDiscEtaBin( double eta )
{
  const unsigned int nEtaBins = s_fnDiscEtaBins;
  const double etaBins[nEtaBins] = {0.1, 0.6, 0.8, 1.15, 1.37, 1.52, 1.81, 2.01, 2.37, 2.47};
  for (unsigned int etaBin = 0; etaBin < nEtaBins; ++etaBin){
    if (std::abs(eta) < etaBins[etaBin]) return etaBin;
  }
  return (nEtaBins-1);
}
 
// Gets the Discriminant Et bin (MeV) [0,s_fnDiscEtBins-1]
unsigned int AsgElectronSelectorTool::getDiscEtBin( double et )
{
  static const double GeV = 1000;
  const unsigned int nEtBins = s_fnDiscEtBins;
  const double etBins[nEtBins] = {7*GeV,10*GeV,15*GeV,20*GeV,25*GeV,30*GeV,35*GeV,40*GeV,45*GeV,6000*GeV};
  for (unsigned int etBin = 0; etBin < nEtBins; ++etBin){
    if (et < etBins[etBin]) return etBin;
  }
  return (nEtBins-1);
}
 
 
// Note that this will only perform the cut interpolation up to ~45 GeV, so
// no smoothing is done above this for the high ET LH binning yet
double AsgElectronSelectorTool::interpolateCuts( const std::vector<double>& cuts,double et,double eta )
{
  const int etbin = getDiscEtBin(et);
  const int etabin = getDiscEtaBin(eta);
  unsigned int ibin_combinedML = etbin*s_fnDiscEtaBins+etabin;
  double cut = cuts.at(ibin_combinedML);
  const double GeV = 1000;
  const double eTBins[10] = {5.5*GeV,8.5*GeV,12.5*GeV,17.5*GeV,22.5*GeV,27.5*GeV,32.5*GeV,37.5*GeV,42.5*GeV,47.5*GeV};
 
  if (et >= eTBins[9]) return cut; // no interpolation for electrons above 47.5 GeV
  if (et <= eTBins[0]) return cut; // no interpolation for electrons below 5.5 GeV
 
  // find the bin where the value is smaller than the next bin
  // Start with bin = 1, since it always has to be at least in
  // bin 1 because of previous cut
  int bin = 1;
  while ( et > eTBins[bin] ) bin++;
 
  double etLow = eTBins[bin-1];
  double etUp = eTBins[bin];
  double discLow = cuts.at((bin-1) * s_fnDiscEtaBins+etabin);
  double discUp = cuts.at((bin) * s_fnDiscEtaBins+etabin);
 
  double gradient = ( discUp - discLow ) / ( etUp - etLow );
 
  return discLow + (et - etLow) * gradient;
}