Tool_FeatureExtractor.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
*/
 
#include "PanTauAlgs/Tool_FeatureExtractor.h"
#include "PanTauAlgs/TauConstituent.h"
#include "PanTauAlgs/PanTauSeed.h"
#include "PanTauAlgs/TauFeature.h"
#include "PanTauAlgs/HelperFunctions.h"
 
 
bool sortTauConstituentMVA(const PanTau::TauConstituent* u, const PanTau::TauConstituent* v) {
  return u->getBDTValue() > v->getBDTValue();
}
 
 
bool sortTauConstituentEt(const PanTau::TauConstituent* u, const PanTau::TauConstituent* v) {
  return u->p4().Et() > v->p4().Et();
}
 
 
PanTau::Tool_FeatureExtractor::Tool_FeatureExtractor(const std::string& name) :
  asg::AsgTool(name)
{}
 
 
StatusCode PanTau::Tool_FeatureExtractor::initialize() {
 
  ATH_MSG_INFO(" initialize()");
  m_init=true;
    
  ATH_CHECK( HelperFunctions::bindToolHandle( m_Tool_InformationStore, m_Tool_InformationStoreName ) );
  ATH_CHECK( m_Tool_InformationStore.retrieve() );
    
  ATH_CHECK( m_Tool_InformationStore->getInfo_String("FeatureExtractor_VarTypeName_varTypeName_Ratio",        m_varTypeName_Ratio) );
  ATH_CHECK( m_Tool_InformationStore->getInfo_String("FeatureExtractor_VarTypeName_varTypeName_StdDev",       m_varTypeName_StdDev) );
  ATH_CHECK( m_Tool_InformationStore->getInfo_String("FeatureExtractor_VarTypeName_varTypeName_HLV",          m_varTypeName_HLV) );
  ATH_CHECK( m_Tool_InformationStore->getInfo_String("FeatureExtractor_VarTypeName_varTypeName_JetMoment",    m_varTypeName_JetMoment) );
  ATH_CHECK( m_Tool_InformationStore->getInfo_String("FeatureExtractor_VarTypeName_varTypeName_Combined",     m_varTypeName_Combined) );
  ATH_CHECK( m_Tool_InformationStore->getInfo_String("FeatureExtractor_VarTypeName_varTypeName_Basic",        m_varTypeName_Basic) );
  ATH_CHECK( m_Tool_InformationStore->getInfo_String("FeatureExtractor_VarTypeName_varTypeName_PID",          m_varTypeName_PID) );
    
  ATH_CHECK( m_Tool_InformationStore->getInfo_Int("FeatureExtractor_UseEmptySeeds",                           m_Config_UseEmptySeeds) );
    
  ATH_CHECK( m_Tool_InformationStore->getInfo_VecDouble("Common_BinEdges_Eta",               m_Config_BinEdges_Eta) );
  ATH_CHECK( m_Tool_InformationStore->getInfo_VecDouble("CellBased_EtaBinned_Pi0MVACut_1prong", m_Config_CellBased_EtaBinned_Pi0MVACut_1prong) );
  ATH_CHECK( m_Tool_InformationStore->getInfo_VecDouble("CellBased_EtaBinned_Pi0MVACut_3prong", m_Config_CellBased_EtaBinned_Pi0MVACut_3prong) );
    
  return StatusCode::SUCCESS;
}
 
 
void PanTau::Tool_FeatureExtractor::fillVariantsSeedEt(const std::vector<PanTau::TauConstituent*>& tauConstituents,
                               std::map<std::string, double>& variants_SeedEt) {
 
  //use different approaches to calculate total energy of seed:
  variants_SeedEt["EtAllConsts"] = 0.0;
    
  //loop over all constituents in seed
  for (unsigned int iConst = 0; iConst < tauConstituents.size(); iConst++) {
        
    //get current constituents
    PanTau::TauConstituent* curConstituent  = tauConstituents.at(iConst);
    double                  curEt           = curConstituent->p4().Et();
        
    //update the different Et definitions
    if (curConstituent->isOfType(PanTau::TauConstituent::t_Charged) ||
        curConstituent->isOfType(PanTau::TauConstituent::t_Neutral) ) {
      variants_SeedEt["EtAllConsts"]    += curEt;
    }
  }//end loop over constituents in seed    
}
 
 
void PanTau::Tool_FeatureExtractor::addFeatureWrtSeedEnergy(PanTau::TauFeature* targetMap,
                                const std::string& featName,
                                double numerator,
                                const std::map<std::string, double>& denominatorMap) {
  std::map<std::string, double>::const_iterator it = denominatorMap.begin();
  for(; it!=denominatorMap.end(); ++it) {
    std::string FullName = featName + it->first;
    float       value    = static_cast<float>(it->second);
    if(value <= 0. || std::isnan(value) || std::isinf(value) ) continue; 
    targetMap->addFeature(FullName, numerator / it->second);
  }
}
 
 
StatusCode PanTau::Tool_FeatureExtractor::execute(PanTau::PanTauSeed* inSeed) const {
    
  bool noAnyConstituents           = inSeed->isOfTechnicalQuality(PanTau::PanTauSeed::t_NoConstituentsAtAll);
  bool noSelConstituents           = inSeed->isOfTechnicalQuality(PanTau::PanTauSeed::t_NoSelectedConstituents);
  bool noValidInputTau             = inSeed->isOfTechnicalQuality(PanTau::PanTauSeed::t_NoValidInputTau);
  bool isBadSeed                   = (noAnyConstituents || noSelConstituents || noValidInputTau);
  if(static_cast<bool>(m_Config_UseEmptySeeds)) isBadSeed = noValidInputTau;
    
  if (isBadSeed) {
    ATH_MSG_DEBUG("Seed is not valid for feature extraction (no constituents or no valid input tau) - just fill isPanTauCandidate feature");
    inSeed->getFeatures()->addFeature("CellBased_" + m_varTypeName_Basic + "_isPanTauCandidate", 0);
    return StatusCode::SUCCESS;
  }
  inSeed->getFeatures()->addFeature("CellBased_" + m_varTypeName_Basic + "_isPanTauCandidate", 1);
    
    
  ATH_CHECK( calculateBasicFeatures(inSeed) );
    
  ATH_CHECK( addConstituentMomenta(inSeed) );
    
  // map containing different methods to calc seed et
  std::map<std::string, double> variants_SeedEt;
  // calculate the Et variants for the seed
  fillVariantsSeedEt(inSeed->getConstituentsAsList_All(), variants_SeedEt);
    
  //loop through all types of Constituents in tau and calculate type features for them
  //baseline
  ATH_CHECK( calculateFeatures(inSeed, PanTau::TauConstituent::t_NoType, variants_SeedEt) );  //=> all constituents
  ATH_CHECK( calculateFeatures(inSeed, PanTau::TauConstituent::t_Charged, variants_SeedEt) ); //=> charged ones in core
  ATH_CHECK( calculateFeatures(inSeed, PanTau::TauConstituent::t_Neutral, variants_SeedEt) ); //=> neutral ones in core
  ATH_CHECK( calculateFeatures(inSeed, PanTau::TauConstituent::t_Pi0Neut, variants_SeedEt) ); //=> pi0 tagged ones in core
    
  //fill the combined features
  ATH_CHECK( addCombinedFeatures(inSeed) );
    
  return StatusCode::SUCCESS;
}
 
 
StatusCode PanTau::Tool_FeatureExtractor::calculateBasicFeatures(PanTau::PanTauSeed* inSeed) const {
    
  PanTau::TauFeature* featureMap = inSeed->getFeatures();
    
  std::string featurePrefix = m_varTypeName_Basic;
  std::string name = "CellBased_" + featurePrefix;
    
  for(int iType=0; iType<PanTau::TauConstituent::t_nTypes; iType++) {
        
    bool foundIt = false;
    std::vector<TauConstituent*> curList = inSeed->getConstituentsOfType(iType, foundIt);
    if(!foundIt) continue;
        
    //store multiplicity of current type
    std::string     typeName        = PanTau::TauConstituent::getTypeName((PanTau::TauConstituent::Type)iType);
    unsigned int    nConstituents   = curList.size();
    featureMap->addFeature(name + "_N" + typeName + "Consts", nConstituents);        
  }
 
  return StatusCode::SUCCESS;
}
 
 
StatusCode PanTau::Tool_FeatureExtractor::addConstituentMomenta(PanTau::PanTauSeed* inSeed) const {
  TauFeature* tauFeatureMap = inSeed->getFeatures();
  std::string prefixVARType = m_varTypeName_HLV;
  for(int iType=0; iType<static_cast<int>(PanTau::TauConstituent::t_nTypes); iType++) {
    bool isOK;
    std::vector<PanTau::TauConstituent*>    list_TypeConstituents   = inSeed->getConstituentsOfType(iType, isOK); // list of constituents of current type
    TLorentzVector                          tlv_TypeConstituents    = inSeed->getSubsystemHLV(iType, isOK); // summed hlv of objects of current type
    std::string                             curTypeName             = PanTau::TauConstituent::getTypeName((PanTau::TauConstituent::Type)iType);
    std::string name = "CellBased_" + curTypeName + "_" + prefixVARType;
        
    if(!list_TypeConstituents.empty()) {
      tauFeatureMap->addFeature(name + "_SumM",   tlv_TypeConstituents.M());
    }
   
  } //end loop over constituent types
    
  return StatusCode::SUCCESS;
}
 
 
StatusCode PanTau::Tool_FeatureExtractor::calculateFeatures(PanTau::PanTauSeed* inSeed,
                                                            int tauConstituentType,
                                const std::map<std::string, double>& variants_SeedEt) const {
    
  std::string                             curTypeName                     = PanTau::TauConstituent::getTypeName( (PanTau::TauConstituent::Type)tauConstituentType );
  std::string                             curTypeName_All                 = PanTau::TauConstituent::AllConstituentsName();
  PanTau::TauFeature*                     tauFeatureMap                   = inSeed->getFeatures();
  TLorentzVector                          tlv_Reference                   = inSeed->getProtoMomentumCore();
    
  std::vector<PanTau::TauConstituent*>    list_AllConstituents            = inSeed->getConstituentsAsList_Core();
    
  bool                                    foundIt                         = false;
  std::vector<PanTau::TauConstituent*>    list_TypeConstituents;
  if (tauConstituentType != PanTau::TauConstituent::t_NoType) list_TypeConstituents = inSeed->getConstituentsOfType(tauConstituentType, foundIt);
  if (tauConstituentType == PanTau::TauConstituent::t_NoType) list_TypeConstituents = list_AllConstituents;
  if (!foundIt) return StatusCode::SUCCESS;
 
  unsigned int                            n_Constituents_Type             = list_TypeConstituents.size();
        
  //sort the lists by Et
  std::sort(list_AllConstituents.begin(),     list_AllConstituents.end(),     sortTauConstituentEt);
  std::sort(list_TypeConstituents.begin(),    list_TypeConstituents.end(),    sortTauConstituentEt);
  
  //make an additional list of constituents, but now ordered by BDT value
  std::vector<PanTau::TauConstituent*>    list_TypeConstituents_SortBDT = list_TypeConstituents;
  std::sort(list_TypeConstituents_SortBDT.begin(), list_TypeConstituents_SortBDT.end(), sortTauConstituentMVA);
    
  TLorentzVector  tlv_1st_BDT;
    
  if(!list_TypeConstituents_SortBDT.empty()) tlv_1st_BDT = list_TypeConstituents_SortBDT[0]->p4();
    
    
  // ===> Sum of Et, Et^2, E and E^2
  double                      sum_Et                      = 0;
  double                      sum_Et2                     = 0;
  // ===> Sum of Et (and E) times DeltaR, DeltaR', Angle 
  double                      sum_EtxDR                   = 0;
  // ===> Multiplicities
  unsigned int                num_EFOs                    = 0;
    
  for(unsigned int iTypeConst=0; iTypeConst<list_TypeConstituents.size(); iTypeConst++) {
        
    //get hep lorentz vector
    TLorentzVector tlv_curConst = list_TypeConstituents.at(iTypeConst)->p4();
        
    //final check (nan & inf)
    if (std::isnan(tlv_curConst.Pt()) || std::isinf(tlv_curConst.Pt())) continue;
        
    //ready to calc stuff
    //basically update each of the prepared sum_* and num_* variables above,
    // the sum HLV and the pointers to 1st, 2nd, 3rd leading constituents of current type
        
    //helpers to reduce function calls
    double hlp_Et               = tlv_curConst.Et();
    double hlp_Et2              = hlp_Et * hlp_Et;
    double hlp_DeltaR           = tlv_Reference.DeltaR(tlv_curConst);
 
    // update Sum of Et, Et^2, E and E^2
    sum_Et                      += hlp_Et;
    sum_Et2                     += hlp_Et2;
    // update Sum of Et (and E) times DeltaR, DeltaR', Angle 
    sum_EtxDR                   += hlp_Et * hlp_DeltaR;
    // update Multiplicities
    num_EFOs++;
  }//end loop over selected EFOs
    
    
    
    
  // naming string
  std::string prefixVARType = "";
    
  //if there are no EFOs of this type in the seed, all other variables cannot be calculated to have any reasonable value
  if( num_EFOs == 0 ) {
    return StatusCode::SUCCESS;
  }
    
  prefixVARType = m_varTypeName_PID;
    
  // Sorted by highest BDT score
  for(unsigned int iTypeConst=0; iTypeConst<n_Constituents_Type; iTypeConst++) {
        
    double value_BDT = list_TypeConstituents_SortBDT[iTypeConst]->getBDTValue();
    if( std::isnan(value_BDT) || std::isinf(value_BDT) ) continue;
        
    //correct BDT value based on BDT cut
    if(tauConstituentType != PanTau::TauConstituent::t_Charged) {
      double mvaCorrection = 0.0;
      double  etaCurConst = list_TypeConstituents_SortBDT[iTypeConst]->p4().Eta();
      int     etaBinIndex = m_HelperFunctions.getBinIndex(m_Config_BinEdges_Eta, std::abs(etaCurConst));
      int     numTrack    = inSeed->getTauJet()->nTracks();
      if(numTrack == 1) { mvaCorrection = m_Config_CellBased_EtaBinned_Pi0MVACut_1prong.at(etaBinIndex); }
      else              { mvaCorrection = m_Config_CellBased_EtaBinned_Pi0MVACut_3prong.at(etaBinIndex); }
 
      value_BDT = value_BDT - mvaCorrection;
    }
        
    std::string iConst = m_HelperFunctions.convertNumberToString(static_cast<double>(iTypeConst+1));
    tauFeatureMap->addFeature("CellBased_" + curTypeName + "_" + prefixVARType + "_BDTValues_BDTSort_" + iConst, value_BDT);
 
  }
    
  prefixVARType = PanTau::Tool_FeatureExtractor::varTypeName_Shots();
  //only execute if the constituent type is neutral
  if(PanTau::TauConstituent::isNeutralType(tauConstituentType)) {
        
    unsigned int            totalPhotonsInSeed      = 0;
        
    for(unsigned int iConst=0; iConst<list_TypeConstituents_SortBDT.size(); iConst++) {
            
      PanTau::TauConstituent*                 curConst            = list_TypeConstituents_SortBDT.at(iConst);
      std::vector<PanTau::TauConstituent*>    shotConstituents    = curConst->getShots();
      unsigned int                            nShots              = shotConstituents.size();
            
      unsigned int totalPhotonsInNeutral = 0;
            
      for(unsigned int iShot=0; iShot<nShots; iShot++) {
    PanTau::TauConstituent* curShot = shotConstituents.at(iShot);
    totalPhotonsInNeutral += curShot->getNPhotonsInShot();
      }//end loop over shots
      totalPhotonsInSeed  += totalPhotonsInNeutral;
            
    }//end loop over constituents in tau
    
    //number of shots in seed
    tauFeatureMap->addFeature("CellBased_" + curTypeName + "_" + prefixVARType + "_NPhotonsInSeed", totalPhotonsInSeed);
 
  }//end if check for shot info dumping
    
    
  prefixVARType = m_varTypeName_Ratio;
    
  if(curTypeName != curTypeName_All) addFeatureWrtSeedEnergy(tauFeatureMap, "CellBased_" + curTypeName + "_" + prefixVARType + "_EtOver", sum_Et, variants_SeedEt);
    
  if(tlv_1st_BDT.Pt() != 0) addFeatureWrtSeedEnergy(tauFeatureMap, "CellBased_" + curTypeName + "_" + prefixVARType + "_1stBDTEtOver",   tlv_1st_BDT.Et(), variants_SeedEt);
    
  prefixVARType =  m_varTypeName_StdDev;
 
  double stddev_Et            = m_HelperFunctions.stddev(sum_Et2, sum_Et, num_EFOs);
    
  if(stddev_Et > 0.)          addFeatureWrtSeedEnergy(tauFeatureMap, "CellBased_" + curTypeName + "_" + prefixVARType + "_Et_Wrt", stddev_Et, variants_SeedEt);
    
  //  Moment wrt X = Sum( Et * X ) / Sum(Et)
  //  ==> Named as Eflow_EFOType_JetMoment_EtX"   when using Et as weight
  //  ==> Named as Eflow_EFOType_JetMoment_EX"    when using E  as weight
  prefixVARType = m_varTypeName_JetMoment;
    
  if(sum_Et > 0.) {
    //Using transverse energy
    tauFeatureMap->addFeature("CellBased_" + curTypeName + "_" + prefixVARType + "_EtDRxTotalEt",    (sum_EtxDR / sum_Et) * variants_SeedEt.at("EtAllConsts"));
  }
 
  return StatusCode::SUCCESS;
}
 
 
StatusCode PanTau::Tool_FeatureExtractor::addCombinedFeatures(PanTau::PanTauSeed* inSeed) const {
    
  // Prepare some short notations for variables
  PanTau::TauFeature* tauFeatures     = inSeed->getFeatures();
    
  // Prepare the list of names for EFO Types & 
  // the 4 momenta of the different sub systems 
  // (ie. charged, neutral subsystem, etc...)
  TLorentzVector tlv_System[PanTau::TauConstituent::t_nTypes];
  TLorentzVector tlv_1stEFO[PanTau::TauConstituent::t_nTypes];
    
  // get input objects to calc combined features
  bool tlv_Sys_OK[PanTau::TauConstituent::t_nTypes];
  bool tlv_1st_OK[PanTau::TauConstituent::t_nTypes];
    
  //initialize arrays with default values
  for(unsigned int iType=0; iType<(unsigned int)PanTau::TauConstituent::t_nTypes; iType++) {
    tlv_System[iType]   = TLorentzVector();
    tlv_1stEFO[iType]   = TLorentzVector();
    tlv_Sys_OK[iType]   = false;
    tlv_1st_OK[iType]   = false;
  }
    
  for(int iType=0; iType<static_cast<int>(PanTau::TauConstituent::t_nTypes); iType++) {
        
    tlv_System[iType] = inSeed->getSubsystemHLV(iType, tlv_Sys_OK[iType]);
    if (!tlv_Sys_OK[iType]) continue;
        
    std::vector<TauConstituent*> typeConstituents = inSeed->getConstituentsOfType(iType, tlv_Sys_OK[iType]);
    if (typeConstituents.empty()) tlv_Sys_OK[iType] = false;
    if (!tlv_Sys_OK[iType]) continue;
        
    if (!typeConstituents.empty()) {
      tlv_1stEFO[iType] = typeConstituents.at(0)->p4();
      tlv_1st_OK[iType] = true;
    } else {
      tlv_1st_OK[iType] = false;
    }
  }    
    
  // From the extracted input, calc combined features
  std::string prefixVARType = m_varTypeName_Combined;
    
  // Combined Type-vs-Type Features
  // Loop over all EFO types...
  // for every type, loop over all other types (but neglect permutations, i.e.   A,B   is the same as   B,A)
  // calculate ratios, angles etc...
    
  for(int iType=0; iType<PanTau::TauConstituent::t_nTypes; iType++) {
        
    if(iType == static_cast<int>(PanTau::TauConstituent::t_NoType)) continue;
    int type_Denom = iType;
        
    for(int jType=0; jType<PanTau::TauConstituent::t_nTypes; jType++) {
            
      if(jType == static_cast<int>(PanTau::TauConstituent::t_NoType)) continue;
      int type_Nom = jType;
            
      if(jType == iType) continue;
            
      std::string typeName_Nom    = PanTau::TauConstituent::getTypeName((PanTau::TauConstituent::Type)type_Nom);
      std::string typeName_Denom  = PanTau::TauConstituent::getTypeName((PanTau::TauConstituent::Type)type_Denom);
            
      double sum_Et_Nom   = 0.0;
      double sum_Et_Denom = 0.0;
      if(tlv_Sys_OK[type_Nom] && tlv_Sys_OK[type_Denom]) {
    sum_Et_Nom   = tlv_System[type_Nom].Et();
    sum_Et_Denom = tlv_System[type_Denom].Et();
      }
            
      //following features symmetric in system A-B, hence skip multi calculation
      if(jType <= iType) continue;
            
      //Energy ratios between systems
      if(sum_Et_Denom > 0. && sum_Et_Nom > 0.) {
    tauFeatures->addFeature("CellBased_" + prefixVARType + "_Log" + typeName_Nom + "EtOver" + typeName_Denom + "Et",   std::log10(sum_Et_Nom / sum_Et_Denom) );
      }//end check for div by zero
      
      if(tlv_1st_OK[type_Nom] && tlv_1st_OK[type_Denom]) {
    //Delta R between 1st and 1st EFO
    tauFeatures->addFeature("CellBased_" + prefixVARType + "_DeltaR1st" + typeName_Nom + "To1st" + typeName_Denom, tlv_1stEFO[type_Nom].DeltaR( tlv_1stEFO[type_Denom] ) );
    
      } //end check for valid leading efo
            
    }//end loop over system B
  }//end loop over system A
    
  return StatusCode::SUCCESS;
}