Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-08-23 03:25:27

0001 #include "RecoMET/METPUSubtraction/interface/PFMETAlgorithmMVA.h"
0002 
0003 #include "FWCore/ParameterSet/interface/FileInPath.h"
0004 #include "FWCore/Utilities/interface/Exception.h"
0005 
0006 #include "DataFormats/METReco/interface/CommonMETData.h"
0007 
0008 #include <TFile.h>
0009 
0010 #include <iomanip>
0011 
0012 enum MVAType { kBaseline = 0 };
0013 
0014 const double Pi = std::cos(-1);
0015 
0016 const std::string PFMETAlgorithmMVA::updateVariableNames(std::string input) {
0017   if (input == "sumet")
0018     return "particleFlow_SumET";
0019   if (input == "npv")
0020     return "nPV";
0021   if (input == "pfu")
0022     return "particleFlow_U";
0023   if (input == "pfuphi")
0024     return "particleFlow_UPhi";
0025   if (input == "tksumet")
0026     return "track_SumET";
0027   if (input == "tku")
0028     return "track_U";
0029   if (input == "tkuphi")
0030     return "track_UPhi";
0031   if (input == "nopusumet")
0032     return "noPileUp_SumET";
0033   if (input == "nopuu")
0034     return "noPileUp_U";
0035   if (input == "nopuuphi")
0036     return "noPileUp_UPhi";
0037   if (input == "pusumet")
0038     return "pileUp_SumET";
0039   if (input == "pumet")
0040     return "pileUp_MET";
0041   if (input == "pumetphi")
0042     return "pileUp_METPhi";
0043   if (input == "pucsumet")
0044     return "pileUpCorrected_SumET";
0045   if (input == "pucu")
0046     return "pileUpCorrected_U";
0047   if (input == "pucuphi")
0048     return "pileUpCorrected_UPhi";
0049   if (input == "jetpt1")
0050     return "jet1_pT";
0051   if (input == "jeteta1")
0052     return "jet1_eta";
0053   if (input == "jetphi1")
0054     return "jet1_Phi";
0055   if (input == "jetpt2")
0056     return "jet2_pT";
0057   if (input == "jeteta2")
0058     return "jet2_eta";
0059   if (input == "jetphi2")
0060     return "jet2_Phi";
0061   if (input == "nalljet")
0062     return "nJets";
0063   if (input == "njet")
0064     return "numJetsPtGt30";
0065   if (input == "uphi_mva")
0066     return "PhiCor_UPhi";
0067   if (input == "uphix_mva")
0068     return "PhiCor_UPhi";
0069   if (input == "ux_mva")
0070     return "RecoilCor_U";
0071   return input;
0072 }
0073 
0074 const GBRForest* PFMETAlgorithmMVA::loadMVAfromFile(const edm::FileInPath& inputFileName, const std::string& mvaName) {
0075   if (inputFileName.location() == edm::FileInPath::Unknown)
0076     throw cms::Exception("PFMETAlgorithmMVA::loadMVA") << " Failed to find File = " << inputFileName << " !!\n";
0077   std::unique_ptr<TFile> inputFile(new TFile(inputFileName.fullPath().data()));
0078 
0079   std::vector<std::string>* lVec = (std::vector<std::string>*)inputFile->Get("varlist");
0080 
0081   if (lVec == nullptr) {
0082     throw cms::Exception("PFMETAlgorithmMVA::loadMVA")
0083         << " Failed to load mva file : " << inputFileName.fullPath().data() << " is not a proper file !!\n";
0084   }
0085 
0086   std::vector<std::string> variableNames;
0087   for (unsigned int i = 0; i < lVec->size(); ++i) {
0088     variableNames.push_back(updateVariableNames(lVec->at(i)));
0089   }
0090 
0091   if (mvaName.find(mvaNameU_) != std::string::npos)
0092     varForU_ = variableNames;
0093   else if (mvaName.find(mvaNameDPhi_) != std::string::npos)
0094     varForDPhi_ = variableNames;
0095   else if (mvaName.find(mvaNameCovU1_) != std::string::npos)
0096     varForCovU1_ = variableNames;
0097   else if (mvaName.find(mvaNameCovU2_) != std::string::npos)
0098     varForCovU2_ = variableNames;
0099   else
0100     throw cms::Exception("PFMETAlgorithmMVA::loadMVA")
0101         << "MVA MET weight file tree names do not match specified inputs" << std::endl;
0102 
0103   const GBRForest* mva = (GBRForest*)inputFile->Get(mvaName.data());
0104   if (!mva)
0105     throw cms::Exception("PFMETAlgorithmMVA::loadMVA")
0106         << " Failed to load MVA = " << mvaName.data() << " from file = " << inputFileName.fullPath().data() << " !!\n";
0107 
0108   return mva;
0109 }
0110 
0111 PFMETAlgorithmMVA::PFMETAlgorithmMVA(const edm::ParameterSet& cfg, edm::ConsumesCollector iC)
0112     : utils_(cfg),
0113       mvaReaderU_(nullptr),
0114       mvaReaderDPhi_(nullptr),
0115       mvaReaderCovU1_(nullptr),
0116       mvaReaderCovU2_(nullptr),
0117       cfg_(cfg) {
0118   mvaType_ = kBaseline;
0119 
0120   edm::ParameterSet cfgInputRecords = cfg_.getParameter<edm::ParameterSet>("inputRecords");
0121   mvaNameU_ = cfgInputRecords.getParameter<std::string>("U");
0122   mvaNameDPhi_ = cfgInputRecords.getParameter<std::string>("DPhi");
0123   mvaNameCovU1_ = cfgInputRecords.getParameter<std::string>("CovU1");
0124   mvaNameCovU2_ = cfgInputRecords.getParameter<std::string>("CovU2");
0125 
0126   loadMVAfromDB_ = cfg.getParameter<bool>("loadMVAfromDB");
0127   if (loadMVAfromDB_) {
0128     mvaTokenU_ = iC.esConsumes(edm::ESInputTag("", mvaNameU_));
0129     mvaTokenDPhi_ = iC.esConsumes(edm::ESInputTag("", mvaNameDPhi_));
0130     mvaTokenCovU1_ = iC.esConsumes(edm::ESInputTag("", mvaNameCovU1_));
0131     mvaTokenCovU2_ = iC.esConsumes(edm::ESInputTag("", mvaNameCovU2_));
0132   }
0133 }
0134 
0135 PFMETAlgorithmMVA::~PFMETAlgorithmMVA() {
0136   if (!loadMVAfromDB_) {
0137     delete mvaReaderU_;
0138     delete mvaReaderDPhi_;
0139     delete mvaReaderCovU1_;
0140     delete mvaReaderCovU2_;
0141   }
0142 }
0143 
0144 //-------------------------------------------------------------------------------
0145 void PFMETAlgorithmMVA::initialize(const edm::EventSetup& es) {
0146   if (loadMVAfromDB_) {
0147     mvaReaderU_ = &es.getData(mvaTokenU_);
0148     mvaReaderDPhi_ = &es.getData(mvaTokenDPhi_);
0149     mvaReaderCovU1_ = &es.getData(mvaTokenCovU1_);
0150     mvaReaderCovU2_ = &es.getData(mvaTokenCovU2_);
0151   } else {
0152     edm::ParameterSet cfgInputFileNames = cfg_.getParameter<edm::ParameterSet>("inputFileNames");
0153 
0154     edm::FileInPath inputFileNameU = cfgInputFileNames.getParameter<edm::FileInPath>("U");
0155     mvaReaderU_ = loadMVAfromFile(inputFileNameU, mvaNameU_);
0156     edm::FileInPath inputFileNameDPhi = cfgInputFileNames.getParameter<edm::FileInPath>("DPhi");
0157     mvaReaderDPhi_ = loadMVAfromFile(inputFileNameDPhi, mvaNameDPhi_);
0158     edm::FileInPath inputFileNameCovU1 = cfgInputFileNames.getParameter<edm::FileInPath>("CovU1");
0159     mvaReaderCovU1_ = loadMVAfromFile(inputFileNameCovU1, mvaNameCovU1_);
0160     edm::FileInPath inputFileNameCovU2 = cfgInputFileNames.getParameter<edm::FileInPath>("CovU2");
0161     mvaReaderCovU2_ = loadMVAfromFile(inputFileNameCovU2, mvaNameCovU2_);
0162   }
0163 }
0164 
0165 //-------------------------------------------------------------------------------
0166 void PFMETAlgorithmMVA::setInput(const std::vector<reco::PUSubMETCandInfo>& leptons,
0167                                  const std::vector<reco::PUSubMETCandInfo>& jets,
0168                                  const std::vector<reco::PUSubMETCandInfo>& pfCandidates,
0169                                  const std::vector<reco::Vertex::Point>& vertices) {
0170   utils_.computeAllSums(jets, leptons, pfCandidates);
0171 
0172   sumLeptonPx_ = utils_.getLeptonsSumMEX();
0173   sumLeptonPy_ = utils_.getLeptonsSumMEY();
0174 
0175   chargedSumLeptonPx_ = utils_.getLeptonsChSumMEX();
0176   chargedSumLeptonPy_ = utils_.getLeptonsChSumMEY();
0177 
0178   const std::vector<reco::PUSubMETCandInfo> jets_cleaned = utils_.getCleanedJets();
0179 
0180   CommonMETData pfRecoil_data = utils_.computeRecoil(MvaMEtUtilities::kPF);
0181   CommonMETData chHSRecoil_data = utils_.computeRecoil(MvaMEtUtilities::kChHS);
0182   CommonMETData hsRecoil_data = utils_.computeRecoil(MvaMEtUtilities::kHS);
0183   CommonMETData puRecoil_data = utils_.computeRecoil(MvaMEtUtilities::kPU);
0184   CommonMETData hsMinusNeutralPUMEt_data = utils_.computeRecoil(MvaMEtUtilities::kHSMinusNeutralPU);
0185 
0186   reco::Candidate::LorentzVector jet1P4 = utils_.leadJetP4(jets_cleaned);
0187   reco::Candidate::LorentzVector jet2P4 = utils_.subleadJetP4(jets_cleaned);
0188 
0189   var_["particleFlow_U"] = pfRecoil_data.met;
0190   var_["particleFlow_SumET"] = pfRecoil_data.sumet;
0191   var_["particleFlow_UPhi"] = pfRecoil_data.phi;
0192 
0193   var_["track_SumET"] = chHSRecoil_data.sumet / var_["particleFlow_SumET"];
0194   var_["track_U"] = chHSRecoil_data.met;
0195   var_["track_UPhi"] = chHSRecoil_data.phi;
0196 
0197   var_["noPileUp_SumET"] = hsRecoil_data.sumet / var_["particleFlow_SumET"];
0198   var_["noPileUp_U"] = hsRecoil_data.met;
0199   var_["noPileUp_UPhi"] = hsRecoil_data.phi;
0200 
0201   var_["pileUp_SumET"] = puRecoil_data.sumet / var_["particleFlow_SumET"];
0202   var_["pileUp_MET"] = puRecoil_data.met;
0203   var_["pileUp_METPhi"] = puRecoil_data.phi;
0204 
0205   var_["pileUpCorrected_SumET"] = hsMinusNeutralPUMEt_data.sumet / var_["particleFlow_SumET"];
0206   var_["pileUpCorrected_U"] = hsMinusNeutralPUMEt_data.met;
0207   var_["pileUpCorrected_UPhi"] = hsMinusNeutralPUMEt_data.phi;
0208 
0209   var_["jet1_pT"] = jet1P4.pt();
0210   var_["jet1_eta"] = jet1P4.eta();
0211   var_["jet1_Phi"] = jet1P4.phi();
0212   var_["jet2_pT"] = jet2P4.pt();
0213   var_["jet2_eta"] = jet2P4.eta();
0214   var_["jet2_Phi"] = jet2P4.phi();
0215 
0216   var_["numJetsPtGt30"] = utils_.numJetsAboveThreshold(jets_cleaned, 30.);
0217   var_["nJets"] = jets_cleaned.size();
0218   var_["nPV"] = vertices.size();
0219 }
0220 
0221 //-------------------------------------------------------------------------------
0222 std::unique_ptr<float[]> PFMETAlgorithmMVA::createFloatVector(std::vector<std::string> variableNames) {
0223   std::unique_ptr<float[]> floatVector(new float[variableNames.size()]);
0224   int i = 0;
0225   for (const auto& variableName : variableNames) {
0226     floatVector[i++] = var_[variableName];
0227   }
0228   return floatVector;
0229 }
0230 
0231 //-------------------------------------------------------------------------------
0232 void PFMETAlgorithmMVA::evaluateMVA() {
0233   // CV: MVAs needs to be evaluated in order { DPhi, U1, CovU1, CovU2 }
0234   //     as MVA for U1 (CovU1, CovU2) uses output of DPhi (DPhi and U1) MVA
0235   mvaOutputDPhi_ = GetResponse(mvaReaderDPhi_, varForDPhi_);
0236   var_["PhiCor_UPhi"] = var_["particleFlow_UPhi"] + mvaOutputDPhi_;
0237   mvaOutputU_ = GetResponse(mvaReaderU_, varForU_);
0238   var_["RecoilCor_U"] = var_["particleFlow_U"] * mvaOutputU_;
0239   var_["RecoilCor_UPhi"] = var_["PhiCor_UPhi"];
0240   mvaOutputCovU1_ = std::pow(GetResponse(mvaReaderCovU1_, varForCovU1_) * mvaOutputU_ * var_["particleFlow_U"], 2);
0241   mvaOutputCovU2_ = std::pow(GetResponse(mvaReaderCovU2_, varForCovU2_) * mvaOutputU_ * var_["particleFlow_U"], 2);
0242 
0243   // compute MET(Photon check)
0244   if (hasPhotons_) {
0245     //Fix events with unphysical properties
0246     double sumLeptonPt = std::max(sqrt(sumLeptonPx_ * sumLeptonPx_ + sumLeptonPy_ * sumLeptonPy_), 1.);
0247     if (var_["track_U"] / sumLeptonPt < 0.1 || var_["noPileUp_U"] / sumLeptonPt < 0.1) {
0248       mvaOutputU_ = 1.;
0249       mvaOutputDPhi_ = 0.;
0250     }
0251   }
0252   computeMET();
0253 }
0254 //-------------------------------------------------------------------------------
0255 
0256 void PFMETAlgorithmMVA::computeMET() {
0257   double U = var_["RecoilCor_U"];
0258   double Phi = var_["PhiCor_UPhi"];
0259   if (U < 0.)
0260     Phi += Pi;  //RF: No sign flip for U necessary in that case?
0261   double cosPhi = std::cos(Phi);
0262   double sinPhi = std::sin(Phi);
0263   double metPx = U * cosPhi - sumLeptonPx_;
0264   double metPy = U * sinPhi - sumLeptonPy_;
0265   double metPt = sqrt(metPx * metPx + metPy * metPy);
0266   mvaMEt_.SetCoordinates(metPx, metPy, 0., metPt);
0267   // compute MET uncertainties in dirrections parallel and perpendicular to hadronic recoil
0268   // (neglecting uncertainties on lepton momenta)
0269   mvaMEtCov_(0, 0) = mvaOutputCovU1_ * cosPhi * cosPhi + mvaOutputCovU2_ * sinPhi * sinPhi;
0270   mvaMEtCov_(0, 1) = -mvaOutputCovU1_ * sinPhi * cosPhi + mvaOutputCovU2_ * sinPhi * cosPhi;
0271   mvaMEtCov_(1, 0) = mvaMEtCov_(0, 1);
0272   mvaMEtCov_(1, 1) = mvaOutputCovU1_ * sinPhi * sinPhi + mvaOutputCovU2_ * cosPhi * cosPhi;
0273 }
0274 
0275 //-------------------------------------------------------------------------------
0276 const float PFMETAlgorithmMVA::GetResponse(const GBRForest* Reader, std::vector<std::string>& variableNames) {
0277   std::unique_ptr<float[]> mvaInputVector = createFloatVector(variableNames);
0278   double result = Reader->GetResponse(mvaInputVector.get());
0279   return result;
0280 }
0281 
0282 //-------------------------------------------------------------------------------
0283 void PFMETAlgorithmMVA::print(std::ostream& stream) const {
0284   stream << "<PFMETAlgorithmMVA::print>:" << std::endl;
0285   for (const auto& entry : var_)
0286     stream << entry.first << " = " << entry.second << std::endl;
0287   stream << " covU1 = " << mvaOutputCovU1_ << ", covU2 = " << mvaOutputCovU2_ << std::endl;
0288   stream << " sum(leptons): Pt = " << sqrt(sumLeptonPx_ * sumLeptonPx_ + sumLeptonPy_ * sumLeptonPy_) << ","
0289          << " phi = " << atan2(sumLeptonPy_, sumLeptonPx_) << " "
0290          << "(Px = " << sumLeptonPx_ << ", Py = " << sumLeptonPy_ << ")";
0291   stream << " MVA output: U = " << mvaOutputU_ << ", dPhi = " << mvaOutputDPhi_ << ","
0292          << " covU1 = " << mvaOutputCovU1_ << ", covU2 = " << mvaOutputCovU2_ << std::endl;
0293   stream << std::endl;
0294 }