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, variants_SeedEt) );
    
  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;
  
  double SumCharge = 0;
  double AbsCharge = 0;
    
  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);        
        
    //count charge, i.e. skip if not charged
    if(iType != static_cast<int>(PanTau::TauConstituent::t_Charged)) continue;
    for(unsigned int iConst=0; iConst<nConstituents; iConst++) {
      PanTau::TauConstituent* curConstituent = curList[iConst];
      SumCharge += curConstituent->getCharge();
      AbsCharge += std::abs(static_cast<double>(curConstituent->getCharge()));
    }
  }
    
  //store charge
  featureMap->addFeature(name + "_SumCharge", SumCharge);
  featureMap->addFeature(name + "_AbsCharge", AbsCharge);
    
  //proto 4-vector (just the sum of all constituents)
  // will have better four momentum after mode ID
  TLorentzVector tlv_ProtoMomentumCore = inSeed->getProtoMomentumCore();
  featureMap->addFeature(name + "_ProtoMomentumCore_pt", tlv_ProtoMomentumCore.Perp());
  featureMap->addFeature(name + "_ProtoMomentumCore_eta", tlv_ProtoMomentumCore.Eta());
  featureMap->addFeature(name + "_ProtoMomentumCore_phi", tlv_ProtoMomentumCore.Phi());
  featureMap->addFeature(name + "_ProtoMomentumCore_m", tlv_ProtoMomentumCore.M());
    
  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
    unsigned int                            n_Constituents_Type     = list_TypeConstituents.size();          // number of objects 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;
        
    std::vector<PanTau::TauConstituent*>    list_TypeConstituents_SortBDT = inSeed->getConstituentsOfType(iType, isOK);
    std::sort(list_TypeConstituents_SortBDT.begin(), list_TypeConstituents_SortBDT.end(), sortTauConstituentMVA);
        
    if(!list_TypeConstituents.empty()) {
      tauFeatureMap->addFeature(name + "_SumPt",  tlv_TypeConstituents.Perp());
      tauFeatureMap->addFeature(name + "_SumEta", tlv_TypeConstituents.Eta());
      tauFeatureMap->addFeature(name + "_SumPhi", tlv_TypeConstituents.Phi());
      tauFeatureMap->addFeature(name + "_SumM",   tlv_TypeConstituents.M());
    }
    
    //store 4-vectors of current type (bdt sort)
    std::vector<double> curConstsBDT_pt    = std::vector<double>(0);
    std::vector<double> curConstsBDT_eta   = std::vector<double>(0);
    std::vector<double> curConstsBDT_phi   = std::vector<double>(0);
    std::vector<double> curConstsBDT_m     = std::vector<double>(0);
    for(unsigned int iConst=0; iConst<n_Constituents_Type; iConst++) {
      TLorentzVector tlv_curConstBDT = list_TypeConstituents_SortBDT[iConst]->p4();
      curConstsBDT_pt.push_back(tlv_curConstBDT.Perp());
      curConstsBDT_eta.push_back(tlv_curConstBDT.Eta());
      curConstsBDT_phi.push_back(tlv_curConstBDT.Phi());
      curConstsBDT_m.push_back(tlv_curConstBDT.M());
    }
    tauFeatureMap->addVecFeature(name + "_BDTSort_Constituents_pt", curConstsBDT_pt);
    tauFeatureMap->addVecFeature(name + "_BDTSort_Constituents_eta", curConstsBDT_eta);
    tauFeatureMap->addVecFeature(name + "_BDTSort_Constituents_phi", curConstsBDT_phi);
    tauFeatureMap->addVecFeature(name + "_BDTSort_Constituents_m", curConstsBDT_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();
    
    
  // ===> hlv for the leading EFOs and the summed HLV
  TLorentzVector              tlv_TypeConstituents;
  // ===> 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
    tlv_TypeConstituents += tlv_curConst;
        
    //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)) {
        
    TLorentzVector          totalTLV_SumShots       = TLorentzVector(0., 0., 0., 0.);
    unsigned int            totalPhotonsInSeed      = 0;
        
    std::vector<TLorentzVector> allShotTLVs = std::vector<TLorentzVector>(0);
        
    for(unsigned int iConst=0; iConst<list_TypeConstituents_SortBDT.size(); iConst++) {
            
      PanTau::TauConstituent*                 curConst            = list_TypeConstituents_SortBDT.at(iConst);
      TLorentzVector                          tlv_CurConst        = curConst->p4();
      std::vector<PanTau::TauConstituent*>    shotConstituents    = curConst->getShots();
      unsigned int                            nShots              = shotConstituents.size();
            
      unsigned int totalPhotonsInNeutral = 0;
      TLorentzVector tlv_SumShots = TLorentzVector(0., 0., 0., 0.);
            
      for(unsigned int iShot=0; iShot<nShots; iShot++) {
    PanTau::TauConstituent* curShot = shotConstituents.at(iShot);
    totalPhotonsInNeutral += curShot->getNPhotonsInShot();
    tlv_SumShots += curShot->p4();
    allShotTLVs.push_back(curShot->p4());
      }//end loop over shots
      totalTLV_SumShots   += tlv_SumShots;
      totalPhotonsInSeed  += totalPhotonsInNeutral;
            
      std::string iConstStr = m_HelperFunctions.convertNumberToString(static_cast<double>(iConst+1));
                       
      tauFeatureMap->addFeature("CellBased_" + curTypeName + "_" + prefixVARType + "_nPhotons_BDTSort_" + iConstStr, totalPhotonsInNeutral);
 
      if(tlv_CurConst.Et() > 0.) tauFeatureMap->addFeature("CellBased_" + curTypeName + "_" + prefixVARType + "_EtSumShotsOverConstEt_BDTSort_" + iConstStr, tlv_SumShots.Et() / tlv_CurConst.Et());
            
    }//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);
    
  if(tlv_1st_BDT.Pt() != 0 && sum_Et > 0.) tauFeatureMap->addFeature("CellBased_" + curTypeName + "_" + prefixVARType + "_1stBDTEtOverTypeEt",    tlv_1st_BDT.Et() / sum_Et);    
    
  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 + "_EtDR",            sum_EtxDR / sum_Et);
    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 std::map<std::string, double>& variants_SeedEt) const {
    
  // Prepare some short notations for variables
  PanTau::TauFeature* tauFeatures     = inSeed->getFeatures();
        
  //et: EFO Type
  int et_Charged      = PanTau::TauConstituent::t_Charged;
  int et_Pi0Neut      = PanTau::TauConstituent::t_Pi0Neut;
  int et_Neutral      = PanTau::TauConstituent::t_Neutral;
    
  bool foundIt;
  std::vector<PanTau::TauConstituent*>    list_NeutralConstituents = inSeed->getConstituentsOfType(et_Neutral, foundIt);
    
  // Prepare the list of names for EFO Types & 
  // the 4 momenta of the different sub systems 
  // (ie. charged, neutral subsystem, etc...)
  std::string    name_EFOType[PanTau::TauConstituent::t_nTypes];
  double         num_EFOs[PanTau::TauConstituent::t_nTypes];
  TLorentzVector tlv_System[PanTau::TauConstituent::t_nTypes];
  TLorentzVector tlv_1stEFO[PanTau::TauConstituent::t_nTypes];
  TLorentzVector tlv_2ndEFO[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++) {
    name_EFOType[iType] = "";
    num_EFOs[iType]     = 0.;
    tlv_System[iType]   = TLorentzVector();
    tlv_1stEFO[iType]   = TLorentzVector();
    tlv_2ndEFO[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++) {
    name_EFOType[iType] = PanTau::TauConstituent::getTypeName((PanTau::TauConstituent::Type)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;
        
    num_EFOs[iType] = typeConstituents.size();
        
    if (!typeConstituents.empty()) {
      tlv_1stEFO[iType] = typeConstituents.at(0)->p4();
      tlv_1st_OK[iType] = true;
    } else {
      tlv_1st_OK[iType] = false;
    }
        
    if (typeConstituents.size() > 1) {
      tlv_2ndEFO[iType] = typeConstituents.at(1)->p4();
    } 
  }    
    
  // 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();
      }
            
      // Fraction of leading EFO of system A wrt complete system B
      // this is not symmetric wrt system A and B, hence do this for all combinations
      if(tlv_1st_OK[type_Nom]) {
    if(tlv_1stEFO[type_Nom].Et() > 0. && sum_Et_Denom > 0.) {
      tauFeatures->addFeature("CellBased_" + prefixVARType + "_Log1st" + typeName_Nom + "EtOver" + typeName_Denom + "Et", std::log10(tlv_1stEFO[type_Nom].Et() / sum_Et_Denom) );
    }
      }
            
      //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
            
      //Angles between systems
      if(tlv_Sys_OK[type_Nom] && tlv_Sys_OK[type_Denom]) {
    tauFeatures->addFeature("CellBased_" + prefixVARType + "_Angle" + typeName_Nom + "To" + typeName_Denom, tlv_System[type_Nom].Angle( tlv_System[type_Denom].Vect()) );
      }//end check for valid system pointer            
      
      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] ) );
    
    //Angles between 1st and 1st EFO
    tauFeatures->addFeature("CellBased_" + prefixVARType + "_Angle1st" + typeName_Nom + "To1st" + typeName_Denom, tlv_1stEFO[type_Nom].Angle( tlv_1stEFO[type_Denom].Vect() ) );
      } //end check for valid leading efo
            
    }//end loop over system B
  }//end loop over system A
    
    
  // Combined Selected-Type Features
 
  //setup arrays for combination of selected charged and neutral combinations
  const int cTypes = 1;
  const int nTypes = 2;
  int index_charged[cTypes] = {et_Charged};
  int index_neutral[nTypes] = {et_Pi0Neut, et_Neutral};
    
  for(int cType=0; cType<cTypes; cType++) {
    for(int nType=0; nType<nTypes; nType++) {
            
      int et_c = index_charged[cType];
      int et_n = index_neutral[nType];
            
      std::string name_cType = PanTau::TauConstituent::getTypeName((PanTau::TauConstituent::Type)et_c);
      std::string name_nType = PanTau::TauConstituent::getTypeName((PanTau::TauConstituent::Type)et_n);
            
      if(!tlv_Sys_OK[et_c] || !tlv_Sys_OK[et_n]) continue;
            
      //mean Et fraction of charged+neutral system wrt total ET
      if(num_EFOs[et_c] + num_EFOs[et_n] > 0.) {
    double mean_cTypenTypeEt = ( tlv_System[et_c].Et() + tlv_System[et_n].Et() ) / (num_EFOs[et_c] + num_EFOs[et_n]);
    addFeatureWrtSeedEnergy(tauFeatures, "CellBased_" + prefixVARType + "_Mean" + name_cType + name_nType + "Et_Wrt", mean_cTypenTypeEt, variants_SeedEt);
      }
            
      //invariant masses
      tauFeatures->addFeature("CellBased_" + prefixVARType + "_InvMass"  + name_cType + name_nType,   ( tlv_System[et_c]  +  tlv_System[et_n] ).M() );
            
    }//end loop neutral types
  }//end loop charged types
    
  return StatusCode::SUCCESS;
}