ElectronMVAEstimatorRun2.cc

Go to the documentation of this file.

#include "RecoEgamma/ElectronIdentification/interface/ElectronMVAEstimatorRun2.h"

ElectronMVAEstimatorRun2::ElectronMVAEstimatorRun2(const edm::ParameterSet& conf):
  AnyMVAEstimatorRun2Base(conf),
  mvaVarMngr_(conf.getParameter<std::string>("variableDefinition"))
{

  const auto weightFileNames = conf.getParameter<std::vector<std::string> >("weightFileNames");
  const auto categoryCutStrings = conf.getParameter<std::vector<std::string> >("categoryCuts");

  if( (int)(categoryCutStrings.size()) != getNCategories() )
    throw cms::Exception("MVA config failure: ")
      << "wrong number of category cuts in ElectronMVAEstimatorRun2" << getTag() << std::endl;

  for (int i = 0; i < getNCategories(); ++i) {
      categoryFunctions_.emplace_back(categoryCutStrings[i]);
  }

  // Initialize GBRForests from weight files
  init(weightFileNames);
}

ElectronMVAEstimatorRun2::ElectronMVAEstimatorRun2( const std::string& mvaTag,
                                                    const std::string& mvaName,
                                                    int nCategories,
                                                    const std::string& variableDefinition,
                                                    const std::vector<std::string>& categoryCutStrings,
                                                    const std::vector<std::string> &weightFileNames,
                                                    bool debug)
  : AnyMVAEstimatorRun2Base( mvaName, mvaTag, nCategories, debug )
  , mvaVarMngr_            (variableDefinition)
{
  
  if( (int)(categoryCutStrings.size()) != getNCategories() )
    throw cms::Exception("MVA config failure: ")
      << "wrong number of category cuts in " << getName() << getTag() << std::endl;
  
  for (auto const& cut : categoryCutStrings) categoryFunctions_.emplace_back(cut);
  init(weightFileNames);
}

void ElectronMVAEstimatorRun2::init(const std::vector<std::string> &weightFileNames) {

  if(isDebug()) {
    std::cout << " *** Inside ElectronMVAEstimatorRun2" << getTag() << std::endl;
  }

  // Initialize GBRForests
  if( (int)(weightFileNames.size()) != getNCategories() )
    throw cms::Exception("MVA config failure: ")
      << "wrong number of weightfiles in ElectronMVAEstimatorRun2" << getTag() << std::endl;

  gbrForests_.clear();
  // Create a TMVA reader object for each category
  for(int i=0; i<getNCategories(); i++){

    std::vector<int> variablesInCategory;

    // Use unique_ptr so that all readers are properly cleaned up
    // when the vector clear() is called in the destructor

    std::vector<std::string> variableNamesInCategory;
    gbrForests_.push_back(createGBRForest(weightFileNames[i], variableNamesInCategory));

    nVariables_.push_back(variableNamesInCategory.size());

    variables_.push_back(variablesInCategory);

    if(isDebug()) {
      std::cout << " *** Inside ElectronMVAEstimatorRun2" << getTag() << std::endl;
      std::cout << " category " << i << " with nVariables " << nVariables_[i] << std::endl;
    }

    for (int j=0; j<nVariables_[i];++j) {
        int index = mvaVarMngr_.getVarIndex(variableNamesInCategory[j]);
        if(index == -1) {
            throw cms::Exception("MVA config failure: ")
               << "Concerning ElectronMVAEstimatorRun2" << getTag() << std::endl
               << "Variable " << variableNamesInCategory[j]
               << " not found in variable definition file!" << std::endl;
        }
        variables_[i].push_back(index);

    }
  }
}

float ElectronMVAEstimatorRun2::
mvaValue( const reco::Candidate* candidate, const std::vector<float>& auxVariables, int &iCategory) const {

  const reco::GsfElectron* electron = dynamic_cast<const reco::GsfElectron*>(candidate);

  if( electron == nullptr ) {
    throw cms::Exception("MVA failure: ")
      << " given particle is expected to be reco::GsfElectron or pat::Electron," << std::endl
      << " but appears to be neither" << std::endl;
  }

  iCategory = findCategory( electron );

  if (iCategory < 0) return -999;

  std::vector<float> vars;

  for (int i = 0; i < nVariables_[iCategory]; ++i) {
      vars.push_back(mvaVarMngr_.getValue(variables_[iCategory][i], *electron, auxVariables));
  }

  if(isDebug()) {
    std::cout << " *** Inside ElectronMVAEstimatorRun2" << getTag() << std::endl;
    std::cout << " category " << iCategory << std::endl;
    for (int i = 0; i < nVariables_[iCategory]; ++i) {
        std::cout << " " << mvaVarMngr_.getName(variables_[iCategory][i]) << " " << vars[i] << std::endl;
    }
  }
  const float response = gbrForests_.at(iCategory)->GetResponse(vars.data()); // The BDT score

  if(isDebug()) {
    std::cout << " ### MVA " << response << std::endl << std::endl;
  }

  return response;
}

int ElectronMVAEstimatorRun2::findCategory( const reco::Candidate* candidate) const {

  const reco::GsfElectron* electron = dynamic_cast<reco::GsfElectron const*>(candidate);

  if( electron == nullptr ) {
    throw cms::Exception("MVA failure: ")
      << " given particle is expected to be reco::GsfElectron or pat::Electron," << std::endl
      << " but appears to be neither" << std::endl;
  }

  return findCategory(*electron);

}

int ElectronMVAEstimatorRun2::findCategory(reco::GsfElectron const& electron) const {

  for (int i = 0; i < getNCategories(); ++i) {
      if (categoryFunctions_[i](electron)) return i;
  }

  edm::LogWarning  ("MVA warning") <<
      "category not defined for particle with pt " << electron.pt() << " GeV, eta " <<
          electron.superCluster()->eta() << " in ElectronMVAEstimatorRun2" << getTag();

  return -1;

}