PFMETAlgorithmMVA.cc

Go to the documentation of this file.

#include "RecoMET/METPUSubtraction/interface/PFMETAlgorithmMVA.h"

#include "FWCore/ParameterSet/interface/FileInPath.h"
#include "FWCore/Utilities/interface/Exception.h"

#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "CondFormats/DataRecord/interface/GBRWrapperRcd.h"

#include "DataFormats/METReco/interface/CommonMETData.h"

#include <TFile.h>

#include <iomanip>

enum MVAType { kBaseline = 0 };

const double Pi=std::cos(-1);

const std::string PFMETAlgorithmMVA::updateVariableNames(std::string input)
{
  if(input=="sumet")     return "particleFlow_SumET";
  if(input=="npv")       return "nPV";
  if(input=="pfu")       return "particleFlow_U";
  if(input=="pfuphi")    return "particleFlow_UPhi";
  if(input=="tksumet")   return "track_SumET";
  if(input=="tku")       return "track_U";
  if(input=="tkuphi")    return "track_UPhi";
  if(input=="nopusumet") return "noPileUp_SumET";
  if(input=="nopuu")     return "noPileUp_U";
  if(input=="nopuuphi")  return "noPileUp_UPhi";
  if(input=="pusumet")   return "pileUp_SumET";
  if(input=="pumet")     return "pileUp_MET";
  if(input=="pumetphi")  return "pileUp_METPhi";
  if(input=="pucsumet")  return "pileUpCorrected_SumET";
  if(input=="pucu")      return "pileUpCorrected_U";
  if(input=="pucuphi")   return "pileUpCorrected_UPhi";
  if(input=="jetpt1")    return "jet1_pT";
  if(input=="jeteta1")   return "jet1_eta";
  if(input=="jetphi1")   return "jet1_Phi";
  if(input=="jetpt2")    return "jet2_pT";
  if(input=="jeteta2")   return "jet2_eta";
  if(input=="jetphi2")   return "jet2_Phi";
  if(input=="nalljet")   return "nJets";
  if(input=="njet")      return "numJetsPtGt30";
  if(input=="uphi_mva")  return "PhiCor_UPhi";
  if(input=="uphix_mva") return "PhiCor_UPhi";
  if(input=="ux_mva")    return "RecoilCor_U";
  return input;
}

const GBRForest* PFMETAlgorithmMVA::loadMVAfromFile(const edm::FileInPath& inputFileName, const std::string& mvaName)
{
  if ( inputFileName.location()==edm::FileInPath::Unknown ) throw cms::Exception("PFMETAlgorithmMVA::loadMVA") 
    << " Failed to find File = " << inputFileName << " !!\n";
  std::unique_ptr<TFile> inputFile(new TFile(inputFileName.fullPath().data()) );

  std::vector<std::string> *lVec = (std::vector<std::string>*)inputFile->Get("varlist");

  if(lVec==nullptr) {
    throw cms::Exception("PFMETAlgorithmMVA::loadMVA")
      << " Failed to load mva file : " << inputFileName.fullPath().data() << " is not a proper file !!\n";
  }

  std::vector<std::string> variableNames;
  for(unsigned int i=0; i< lVec->size();++i)
  {
      variableNames.push_back(updateVariableNames(lVec->at(i)));
  }

  if(mvaName.find(mvaNameU_)      != std::string::npos) varForU_    = variableNames;
  else if(mvaName.find(mvaNameDPhi_)   != std::string::npos) varForDPhi_ = variableNames;
  else if(mvaName.find(mvaNameCovU1_)  != std::string::npos) varForCovU1_ = variableNames;
  else if(mvaName.find(mvaNameCovU2_)  != std::string::npos) varForCovU2_ = variableNames;
  else throw cms::Exception("PFMETAlgorithmMVA::loadMVA") << "MVA MET weight file tree names do not match specified inputs" << std::endl;


  const GBRForest* mva = (GBRForest*)inputFile->Get(mvaName.data());
  if ( !mva )
    throw cms::Exception("PFMETAlgorithmMVA::loadMVA")
      << " Failed to load MVA = " << mvaName.data() << " from file = " << inputFileName.fullPath().data() << " !!\n";

  return mva;
}

const GBRForest* PFMETAlgorithmMVA::loadMVAfromDB(const edm::EventSetup& es, const std::string& mvaName)
{
  edm::ESHandle<GBRForest> mva;
  es.get<GBRWrapperRcd>().get(mvaName, mva);
  return mva.product();
}

PFMETAlgorithmMVA::PFMETAlgorithmMVA(const edm::ParameterSet& cfg) 
  : utils_(cfg),
    mvaReaderU_(nullptr),
    mvaReaderDPhi_(nullptr),
    mvaReaderCovU1_(nullptr),
    mvaReaderCovU2_(nullptr),
    cfg_(cfg)
{
  mvaType_ = kBaseline;

  loadMVAfromDB_ = cfg.getParameter<bool>("loadMVAfromDB");
  
}

PFMETAlgorithmMVA::~PFMETAlgorithmMVA()
{
  if ( !loadMVAfromDB_ ) {
    delete mvaReaderU_;
    delete mvaReaderDPhi_;
    delete mvaReaderCovU1_;
    delete mvaReaderCovU2_;
  }
}

//-------------------------------------------------------------------------------
void PFMETAlgorithmMVA::initialize(const edm::EventSetup& es)
{
  edm::ParameterSet cfgInputRecords = cfg_.getParameter<edm::ParameterSet>("inputRecords");
  mvaNameU_       = cfgInputRecords.getParameter<std::string>("U");
  mvaNameDPhi_    = cfgInputRecords.getParameter<std::string>("DPhi");
  mvaNameCovU1_   = cfgInputRecords.getParameter<std::string>("CovU1");
  mvaNameCovU2_   = cfgInputRecords.getParameter<std::string>("CovU2");

  if ( loadMVAfromDB_ ) {
    mvaReaderU_     = loadMVAfromDB(es, mvaNameU_);
    mvaReaderDPhi_  = loadMVAfromDB(es, mvaNameDPhi_);
    mvaReaderCovU1_ = loadMVAfromDB(es, mvaNameCovU1_);
    mvaReaderCovU2_ = loadMVAfromDB(es, mvaNameCovU2_);
  } else {
    edm::ParameterSet cfgInputFileNames = cfg_.getParameter<edm::ParameterSet>("inputFileNames");
    
    edm::FileInPath inputFileNameU = cfgInputFileNames.getParameter<edm::FileInPath>("U");
    mvaReaderU_     = loadMVAfromFile(inputFileNameU, mvaNameU_);
    edm::FileInPath inputFileNameDPhi = cfgInputFileNames.getParameter<edm::FileInPath>("DPhi");
    mvaReaderDPhi_  = loadMVAfromFile(inputFileNameDPhi, mvaNameDPhi_);
    edm::FileInPath inputFileNameCovU1 = cfgInputFileNames.getParameter<edm::FileInPath>("CovU1");
    mvaReaderCovU1_ = loadMVAfromFile(inputFileNameCovU1, mvaNameCovU1_);
    edm::FileInPath inputFileNameCovU2 = cfgInputFileNames.getParameter<edm::FileInPath>("CovU2");
    mvaReaderCovU2_ = loadMVAfromFile(inputFileNameCovU2, mvaNameCovU2_);
  }
}

//-------------------------------------------------------------------------------
void PFMETAlgorithmMVA::setInput(const std::vector<reco::PUSubMETCandInfo>& leptons,
                 const std::vector<reco::PUSubMETCandInfo>& jets,
                 const std::vector<reco::PUSubMETCandInfo>& pfCandidates,
                 const std::vector<reco::Vertex::Point>& vertices)
{

  
  utils_.computeAllSums( jets, leptons, pfCandidates);
  
  sumLeptonPx_        = utils_.getLeptonsSumMEX();
  sumLeptonPy_        = utils_.getLeptonsSumMEY();

  chargedSumLeptonPx_ = utils_.getLeptonsChSumMEX();
  chargedSumLeptonPy_ = utils_.getLeptonsChSumMEY();

  const std::vector<reco::PUSubMETCandInfo> jets_cleaned = utils_.getCleanedJets();

  CommonMETData pfRecoil_data  = utils_.computeRecoil( MvaMEtUtilities::kPF );
  CommonMETData chHSRecoil_data  = utils_.computeRecoil( MvaMEtUtilities::kChHS );
  CommonMETData hsRecoil_data = utils_.computeRecoil( MvaMEtUtilities::kHS );
  CommonMETData puRecoil_data = utils_.computeRecoil( MvaMEtUtilities::kPU );
  CommonMETData hsMinusNeutralPUMEt_data = utils_.computeRecoil( MvaMEtUtilities::kHSMinusNeutralPU );

  reco::Candidate::LorentzVector jet1P4 = utils_.leadJetP4(jets_cleaned);
  reco::Candidate::LorentzVector jet2P4 = utils_.subleadJetP4(jets_cleaned);


  var_["particleFlow_U"]        = pfRecoil_data.met;
  var_["particleFlow_SumET"]    = pfRecoil_data.sumet;
  var_["particleFlow_UPhi"]     = pfRecoil_data.phi;

  var_["track_SumET"]           = chHSRecoil_data.sumet/var_["particleFlow_SumET"];
  var_["track_U"]               = chHSRecoil_data.met;
  var_["track_UPhi"]            = chHSRecoil_data.phi;

  var_["noPileUp_SumET"]        = hsRecoil_data.sumet/var_["particleFlow_SumET"];
  var_["noPileUp_U"]            = hsRecoil_data.met;
  var_["noPileUp_UPhi"]         = hsRecoil_data.phi;

  var_["pileUp_SumET"]          = puRecoil_data.sumet/var_["particleFlow_SumET"];
  var_["pileUp_MET"]            = puRecoil_data.met;
  var_["pileUp_METPhi"]         = puRecoil_data.phi;

  var_["pileUpCorrected_SumET"] = hsMinusNeutralPUMEt_data.sumet/var_["particleFlow_SumET"];
  var_["pileUpCorrected_U"]     = hsMinusNeutralPUMEt_data.met;
  var_["pileUpCorrected_UPhi"]  = hsMinusNeutralPUMEt_data.phi;

  var_["jet1_pT"]               = jet1P4.pt();
  var_["jet1_eta"]              = jet1P4.eta();
  var_["jet1_Phi"]              = jet1P4.phi();
  var_["jet2_pT"]               = jet2P4.pt();
  var_["jet2_eta"]              = jet2P4.eta();
  var_["jet2_Phi"]              = jet2P4.phi();

  var_["numJetsPtGt30"]         = utils_.numJetsAboveThreshold(jets_cleaned, 30.);
  var_["nJets"]                 = jets_cleaned.size();
  var_["nPV"]                   = vertices.size();
}

//-------------------------------------------------------------------------------
std::unique_ptr<float[]> PFMETAlgorithmMVA::createFloatVector(std::vector<std::string> variableNames)
{
  std::unique_ptr<float[]> floatVector(new float[variableNames.size()]);
    int i = 0;
    for(auto variableName: variableNames)
    {
        floatVector[i++] = var_[variableName];
    }
    return floatVector;
}

//-------------------------------------------------------------------------------
void PFMETAlgorithmMVA::evaluateMVA()
{
  // CV: MVAs needs to be evaluated in order { DPhi, U1, CovU1, CovU2 }
  //     as MVA for U1 (CovU1, CovU2) uses output of DPhi (DPhi and U1) MVA
  mvaOutputDPhi_  = GetResponse(mvaReaderDPhi_, varForDPhi_);
  var_["PhiCor_UPhi"] = var_["particleFlow_UPhi"] + mvaOutputDPhi_;
  mvaOutputU_     = GetResponse(mvaReaderU_, varForU_);
  var_["RecoilCor_U"] = var_["particleFlow_U"] * mvaOutputU_;
  var_["RecoilCor_UPhi"] = var_["PhiCor_UPhi"];
  mvaOutputCovU1_ = std::pow(GetResponse(mvaReaderCovU1_, varForCovU1_)* mvaOutputU_ * var_["particleFlow_U"], 2);
  mvaOutputCovU2_ = std::pow(GetResponse(mvaReaderCovU2_, varForCovU2_)* mvaOutputU_ * var_["particleFlow_U"], 2);


  // compute MET(Photon check)
  if(hasPhotons_) { 
    //Fix events with unphysical properties
    double sumLeptonPt = std::max(sqrt(sumLeptonPx_*sumLeptonPx_+sumLeptonPy_*sumLeptonPy_),1.);
    if(var_["track_U"]/sumLeptonPt < 0.1 || var_["noPileUp_U"]/sumLeptonPt <  0.1 ) {
      mvaOutputU_      = 1.;
      mvaOutputDPhi_   = 0.;
    }
  }
  computeMET();
}
//-------------------------------------------------------------------------------

void PFMETAlgorithmMVA::computeMET()
{
    double U      = var_["RecoilCor_U"];
    double Phi    = var_["PhiCor_UPhi"];
    if ( U < 0. ) Phi += Pi; //RF: No sign flip for U necessary in that case?
    double cosPhi = std::cos(Phi);
    double sinPhi = std::sin(Phi);
    double metPx  = U*cosPhi - sumLeptonPx_;
    double metPy  = U*sinPhi - sumLeptonPy_;
    double metPt  = sqrt(metPx*metPx + metPy*metPy);
    mvaMEt_.SetCoordinates(metPx, metPy, 0., metPt);
    // compute MET uncertainties in dirrections parallel and perpendicular to hadronic recoil
    // (neglecting uncertainties on lepton momenta)
    mvaMEtCov_(0, 0) =  mvaOutputCovU1_*cosPhi*cosPhi + mvaOutputCovU2_*sinPhi*sinPhi;
    mvaMEtCov_(0, 1) = -mvaOutputCovU1_*sinPhi*cosPhi + mvaOutputCovU2_*sinPhi*cosPhi;
    mvaMEtCov_(1, 0) =  mvaMEtCov_(0, 1);
    mvaMEtCov_(1, 1) =  mvaOutputCovU1_*sinPhi*sinPhi + mvaOutputCovU2_*cosPhi*cosPhi;
}

//-------------------------------------------------------------------------------
const float PFMETAlgorithmMVA::GetResponse(const GBRForest * Reader, std::vector<std::string> &variableNames )
{
    std::unique_ptr<float[]> mvaInputVector = createFloatVector(variableNames);
    double result = Reader->GetResponse( mvaInputVector.get() );
    return result;
}

//-------------------------------------------------------------------------------
void PFMETAlgorithmMVA::print(std::ostream& stream) const
{
  stream << "<PFMETAlgorithmMVA::print>:" << std::endl;
    for(auto entry: var_)
        stream << entry.first << " = " << entry.second << std::endl;
  stream << " covU1 = " << mvaOutputCovU1_ << ", covU2 = " << mvaOutputCovU2_ << std::endl;
  stream << " sum(leptons): Pt = " << sqrt(sumLeptonPx_*sumLeptonPx_ + sumLeptonPy_*sumLeptonPy_) << ","
  << " phi = " << atan2(sumLeptonPy_, sumLeptonPx_) << " "
  << "(Px = " << sumLeptonPx_ << ", Py = " << sumLeptonPy_ << ")" ;
  stream << " MVA output: U = " << mvaOutputU_ << ", dPhi = " << mvaOutputDPhi_ << "," << " covU1 = " << mvaOutputCovU1_ << ", covU2 = " << mvaOutputCovU2_ << std::endl;
  stream << std::endl;
}