diff --git a/RecoMET/METPUSubtraction/interface/PFMETAlgorithmMVA.h b/RecoMET/METPUSubtraction/interface/PFMETAlgorithmMVA.h index 2f5a081bb3d99..140546bf052ae 100644 --- a/RecoMET/METPUSubtraction/interface/PFMETAlgorithmMVA.h +++ b/RecoMET/METPUSubtraction/interface/PFMETAlgorithmMVA.h @@ -58,21 +58,14 @@ class PFMETAlgorithmMVA void print(std::ostream&) const; private: + const std::string updateVariableNames(std::string input); + const GBRForest* loadMVAfromFile(const edm::FileInPath& inputFileName, const std::string& mvaName); + const GBRForest* loadMVAfromDB(const edm::EventSetup& es, const std::string& mvaName); - void setInput(double, double, double, - double, double, double, - double, double, double, - double, double, double, - double, double, double, - double, double, double, - double, double, double, - double, double, - double); - - void evaluateU(); - void evaluateDPhi(); - void evaluateCovU1(); - void evaluateCovU2(); + const float evaluateU(); + const float evaluateDPhi(); + const float evaluateCovU1(); + const float evaluateCovU2(); MvaMEtUtilities utils_; @@ -83,41 +76,29 @@ class PFMETAlgorithmMVA int mvaType_; bool hasPhotons_; - - Float_t pfSumEt_; - Float_t pfU_; - Float_t pfPhi_; - Float_t tkSumEt_; - Float_t tkU_; - Float_t tkPhi_; - Float_t npuSumEt_; - Float_t npuU_; - Float_t npuPhi_; - Float_t puSumEt_; - Float_t puMEt_; - Float_t puPhi_; - Float_t pucSumEt_; - Float_t pucU_; - Float_t pucPhi_; - Float_t jet1Pt_; - Float_t jet1Eta_; - Float_t jet1Phi_; - Float_t jet2Pt_; - Float_t jet2Eta_; - Float_t jet2Phi_; - Float_t numJetsPtGt30_; - Float_t numJets_; - Float_t numVertices_; - - Float_t* mvaInputU_; - Float_t* mvaInputDPhi_; - Float_t* mvaInputCovU1_; - Float_t* mvaInputCovU2_; + + double dZcut_; + std::unique_ptr createFloatVector(std::vector variableNames); + const float GetResponse(const GBRForest *Reader, std::vector &variableNames); + void computeMET(); + std::map var_; + + + float* mvaInputU_; + float* mvaInputDPhi_; + float* mvaInputCovU1_; + float* mvaInputCovU2_; - Float_t mvaOutputU_; - Float_t mvaOutputDPhi_; - Float_t mvaOutputCovU1_; - Float_t mvaOutputCovU2_; + float mvaOutputU_; + float mvaOutputDPhi_; + float mvaOutputCovU1_; + float mvaOutputCovU2_; + + std::vector varForU_; + std::vector varForDPhi_; + std::vector varForCovU1_; + std::vector varForCovU2_; + double sumLeptonPx_; double sumLeptonPy_; diff --git a/RecoMET/METPUSubtraction/plugins/PFMETProducerMVA.h b/RecoMET/METPUSubtraction/plugins/PFMETProducerMVA.h index b7d179de568aa..9a7bb2226d4dc 100644 --- a/RecoMET/METPUSubtraction/plugins/PFMETProducerMVA.h +++ b/RecoMET/METPUSubtraction/plugins/PFMETProducerMVA.h @@ -106,7 +106,6 @@ namespace reco edm::EDGetTokenT > srcRho_; std::string correctorLabel_; - bool isOld42_ ; bool useType1_; double globalThreshold_; diff --git a/RecoMET/METPUSubtraction/python/mvaPFMET_cff.py b/RecoMET/METPUSubtraction/python/mvaPFMET_cff.py index 686457533f264..2d94f158565ad 100644 --- a/RecoMET/METPUSubtraction/python/mvaPFMET_cff.py +++ b/RecoMET/METPUSubtraction/python/mvaPFMET_cff.py @@ -72,11 +72,17 @@ globalThreshold = cms.double(-1.),#pfMet.globalThreshold, minCorrJetPt = cms.double(-1.), inputFileNames = cms.PSet( - U = cms.FileInPath('RecoMET/METPUSubtraction/data/gbrmet_53_Sep2013_type1.root'), - DPhi = cms.FileInPath('RecoMET/METPUSubtraction/data/gbrmetphi_53_June2013_type1.root'), - CovU1 = cms.FileInPath('RecoMET/METPUSubtraction/data/gbru1cov_53_Dec2012.root'), - CovU2 = cms.FileInPath('RecoMET/METPUSubtraction/data/gbru2cov_53_Dec2012.root') - ), + U = cms.FileInPath('RecoMET/METPUSubtraction/data/gbrmet_7_2_X_MINIAOD_BX25PU20_Mar2015.root'), + DPhi = cms.FileInPath('RecoMET/METPUSubtraction/data/gbrphi_7_2_X_MINIAOD_BX25PU20_Mar2015.root'), + CovU1 = cms.FileInPath('RecoMET/METPUSubtraction/data/gbru1cov_7_2_X_MINIAOD_BX25PU20_Mar2015.root'), + CovU2 = cms.FileInPath('RecoMET/METPUSubtraction/data/gbru2cov_7_2_X_MINIAOD_BX25PU20_Mar2015.root') + ), + inputRecords = cms.PSet( + U = cms.string("RecoilCor"), + DPhi = cms.string("PhiCor"), + CovU1 = cms.string("CovU1"), + CovU2 = cms.string("CovU2") + ), loadMVAfromDB = cms.bool(False), corrector = cms.string("ak4PFL1Fastjet"), diff --git a/RecoMET/METPUSubtraction/src/PFMETAlgorithmMVA.cc b/RecoMET/METPUSubtraction/src/PFMETAlgorithmMVA.cc index 2b983812b355f..6ece3fb222d87 100644 --- a/RecoMET/METPUSubtraction/src/PFMETAlgorithmMVA.cc +++ b/RecoMET/METPUSubtraction/src/PFMETAlgorithmMVA.cc @@ -15,41 +15,83 @@ enum MVAType { kBaseline = 0 }; -const double Pi=cos(-1); +const double Pi=std::cos(-1); -namespace +const std::string PFMETAlgorithmMVA::updateVariableNames(std::string input) { - const GBRForest* 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"; - TFile* inputFile = new TFile(inputFileName.fullPath().data()); - - //const GBRForest* mva = dynamic_cast(inputFile->Get(mvaName.data())); // CV: dynamic_cast fails for some reason ?! - 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"; - - delete inputFile; + 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; +} - return mva; +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 inputFile(new TFile(inputFileName.fullPath().data()) ); + + std::vector *lVec = (std::vector*)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"; } - const GBRForest* loadMVAfromDB(const edm::EventSetup& es, const std::string& mvaName) + std::vector variableNames; + for(unsigned int i=0; i< lVec->size();++i) { - edm::ESHandle mva; - es.get().get(mvaName, mva); - return mva.product(); + 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 mva; + es.get().get(mvaName, mva); + return mva.product(); } PFMETAlgorithmMVA::PFMETAlgorithmMVA(const edm::ParameterSet& cfg) : utils_(cfg), - mvaInputU_(nullptr), - mvaInputDPhi_(nullptr), - mvaInputCovU1_(nullptr), - mvaInputCovU2_(nullptr), mvaReaderU_(nullptr), mvaReaderDPhi_(nullptr), mvaReaderCovU1_(nullptr), @@ -60,19 +102,10 @@ PFMETAlgorithmMVA::PFMETAlgorithmMVA(const edm::ParameterSet& cfg) loadMVAfromDB_ = cfg.getParameter("loadMVAfromDB"); - mvaInputU_ = new Float_t[25]; - mvaInputDPhi_ = new Float_t[23]; - mvaInputCovU1_ = new Float_t[26]; - mvaInputCovU2_ = new Float_t[26]; } PFMETAlgorithmMVA::~PFMETAlgorithmMVA() { - delete mvaInputU_; - delete mvaInputDPhi_; - delete mvaInputCovU1_; - delete mvaInputCovU2_; - if ( !loadMVAfromDB_ ) { delete mvaReaderU_; delete mvaReaderDPhi_; @@ -81,26 +114,23 @@ PFMETAlgorithmMVA::~PFMETAlgorithmMVA() } } +//------------------------------------------------------------------------------- void PFMETAlgorithmMVA::initialize(const edm::EventSetup& es) { + edm::ParameterSet cfgInputRecords = cfg_.getParameter("inputRecords"); + mvaNameU_ = cfgInputRecords.getParameter("U"); + mvaNameDPhi_ = cfgInputRecords.getParameter("DPhi"); + mvaNameCovU1_ = cfgInputRecords.getParameter("CovU1"); + mvaNameCovU2_ = cfgInputRecords.getParameter("CovU2"); + if ( loadMVAfromDB_ ) { - edm::ParameterSet cfgInputRecords = cfg_.getParameter("inputRecords"); - mvaNameU_ = cfgInputRecords.getParameter("U"); mvaReaderU_ = loadMVAfromDB(es, mvaNameU_); - mvaNameDPhi_ = cfgInputRecords.getParameter("DPhi"); mvaReaderDPhi_ = loadMVAfromDB(es, mvaNameDPhi_); - mvaNameCovU1_ = cfgInputRecords.getParameter("CovU1"); mvaReaderCovU1_ = loadMVAfromDB(es, mvaNameCovU1_); - mvaNameCovU2_ = cfgInputRecords.getParameter("CovU2"); mvaReaderCovU2_ = loadMVAfromDB(es, mvaNameCovU2_); } else { edm::ParameterSet cfgInputFileNames = cfg_.getParameter("inputFileNames"); - mvaNameU_ = "U1Correction"; - mvaNameDPhi_ = "PhiCorrection"; - mvaNameCovU1_ = "CovU1"; - mvaNameCovU2_ = "CovU2"; - edm::FileInPath inputFileNameU = cfgInputFileNames.getParameter("U"); mvaReaderU_ = loadMVAfromFile(inputFileNameU, mvaNameU_); edm::FileInPath inputFileNameDPhi = cfgInputFileNames.getParameter("DPhi"); @@ -139,255 +169,116 @@ void PFMETAlgorithmMVA::setInput(const std::vector& lept reco::Candidate::LorentzVector jet1P4 = utils_.leadJetP4(jets_cleaned); reco::Candidate::LorentzVector jet2P4 = utils_.subleadJetP4(jets_cleaned); - double pfSumEt = pfRecoil_data.sumet; - double pfU = pfRecoil_data.met; - double pfPhi = pfRecoil_data.phi; - double tkSumEt = chHSRecoil_data.sumet; - double tkU = chHSRecoil_data.met; - double tkPhi = chHSRecoil_data.phi; - double npuSumEt = hsRecoil_data.sumet; - double npuU = hsRecoil_data.met; - double npuPhi = hsRecoil_data.phi; - double puSumEt = puRecoil_data.sumet; - double puMEt = puRecoil_data.met; - double puPhi = puRecoil_data.phi; - double pucSumEt = hsMinusNeutralPUMEt_data.sumet; - double pucU = hsMinusNeutralPUMEt_data.met; - double pucPhi = hsMinusNeutralPUMEt_data.phi; - double jet1Pt = jet1P4.pt(); - double jet1Eta = jet1P4.eta(); - double jet1Phi = jet1P4.phi(); - double jet2Pt = jet2P4.pt(); - double jet2Eta = jet2P4.eta(); - double jet2Phi = jet2P4.phi(); - double numJetsPtGt30 = utils_.numJetsAboveThreshold(jets_cleaned, 30.); - double numJets = jets_cleaned.size(); - double numVertices = vertices.size(); - - setInput(pfSumEt, pfU, pfPhi, - tkSumEt, tkU, tkPhi, - npuSumEt, npuU, npuPhi, - puSumEt, puMEt, puPhi, - pucSumEt, pucU, pucPhi, - jet1Pt, jet1Eta, jet1Phi, - jet2Pt, jet2Eta, jet2Phi, - numJetsPtGt30, numJets, - numVertices); + 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(); } -void PFMETAlgorithmMVA::setInput(double pfSumEt, double pfU, double pfPhi, - double tkSumEt, double tkU, double tkPhi, - double npuSumEt, double npuU, double npuPhi, - double puSumEt, double puMEt, double puPhi, - double pucSumEt, double pucU, double pucPhi, - double jet1Pt, double jet1Eta, double jet1Phi, - double jet2Pt, double jet2Eta, double jet2Phi, - double numJetsPtGt30, double numJets, - double numVertices) +//------------------------------------------------------------------------------- +std::unique_ptr PFMETAlgorithmMVA::createFloatVector(std::vector variableNames) { - // protection against "empty events" - if ( pfSumEt < 1. ) pfSumEt = 1.; - - pfSumEt_ = pfSumEt; - pfU_ = pfU; - pfPhi_ = pfPhi; - tkSumEt_ = tkSumEt/pfSumEt_; - tkU_ = tkU; - tkPhi_ = tkPhi; - npuSumEt_ = npuSumEt/pfSumEt_; - npuU_ = npuU; - npuPhi_ = npuPhi; - puSumEt_ = puSumEt/pfSumEt_; - puMEt_ = puMEt; - puPhi_ = puPhi; - pucSumEt_ = pucSumEt/pfSumEt_; - pucU_ = pucU; - pucPhi_ = pucPhi; - jet1Pt_ = jet1Pt; - jet1Eta_ = jet1Eta; - jet1Phi_ = jet1Phi; - jet2Pt_ = jet2Pt; - jet2Eta_ = jet2Eta; - jet2Phi_ = jet2Phi; - numJetsPtGt30_ = numJetsPtGt30; - numJets_ = numJets; - numVertices_ = numVertices; + std::unique_ptr 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 - evaluateDPhi(); - evaluateU(); - evaluateCovU1(); - evaluateCovU2(); + 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_ = GetResponse(mvaReaderCovU1_, varForCovU1_)* mvaOutputU_ * var_["particleFlow_U"]; + mvaOutputCovU2_ = GetResponse(mvaReaderCovU2_, varForCovU2_)* mvaOutputU_ * var_["particleFlow_U"]; + // compute MET(Photon check) if(hasPhotons_) { //Fix events with unphysical properties double sumLeptonPt = std::max(sqrt(sumLeptonPx_*sumLeptonPx_+sumLeptonPy_*sumLeptonPy_),1.); - if(tkU_/sumLeptonPt < 0.1 || npuU_/sumLeptonPt < 0.1 ) mvaOutputU_ = 1.; - if(tkU_/sumLeptonPt < 0.1 || npuU_/sumLeptonPt < 0.1 ) mvaOutputDPhi_ = 0.; + if(var_["track_U"]/sumLeptonPt < 0.1 || var_["noPileUp_U"]/sumLeptonPt < 0.1 ) { + mvaOutputU_ = 1.; + mvaOutputDPhi_ = 0.; + } } - double U = pfU_*mvaOutputU_; - double Phi = pfPhi_ + mvaOutputDPhi_; - if ( U < 0. ) Phi += Pi; - double cosPhi = cos(Phi); - double sinPhi = 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; + computeMET(); } //------------------------------------------------------------------------------- -//------------------------------------------------------------------------------- -void PFMETAlgorithmMVA::evaluateU() -{ - mvaInputU_[0] = pfSumEt_; // PH: helps flattens response vs. Nvtx - mvaInputU_[1] = numVertices_; - mvaInputU_[2] = pfU_; - mvaInputU_[3] = pfPhi_; - mvaInputU_[4] = tkSumEt_; - mvaInputU_[5] = tkU_; - mvaInputU_[6] = tkPhi_; - mvaInputU_[7] = npuSumEt_; - mvaInputU_[8] = npuU_; - mvaInputU_[9] = npuPhi_; - mvaInputU_[10] = puSumEt_; - mvaInputU_[11] = puMEt_; - mvaInputU_[12] = puPhi_; - mvaInputU_[13] = pucSumEt_; - mvaInputU_[14] = pucU_; - mvaInputU_[15] = pucPhi_; - mvaInputU_[16] = jet1Pt_; - mvaInputU_[17] = jet1Eta_; - mvaInputU_[18] = jet1Phi_; - mvaInputU_[19] = jet2Pt_; - mvaInputU_[20] = jet2Eta_; - mvaInputU_[21] = jet2Phi_; - mvaInputU_[22] = numJets_; - mvaInputU_[23] = numJetsPtGt30_; - mvaInputU_[24] = pfPhi_ + mvaOutputDPhi_; - mvaOutputU_ = mvaReaderU_->GetResponse(mvaInputU_); -} - -void PFMETAlgorithmMVA::evaluateDPhi() -{ - mvaInputDPhi_[0] = numVertices_; - mvaInputDPhi_[1] = pfU_; - mvaInputDPhi_[2] = pfPhi_; - mvaInputDPhi_[3] = tkSumEt_; - mvaInputDPhi_[4] = tkU_; - mvaInputDPhi_[5] = tkPhi_; - mvaInputDPhi_[6] = npuSumEt_; - mvaInputDPhi_[7] = npuU_; - mvaInputDPhi_[8] = npuPhi_; - mvaInputDPhi_[9] = puSumEt_; - mvaInputDPhi_[10] = puMEt_; - mvaInputDPhi_[11] = puPhi_; - mvaInputDPhi_[12] = pucSumEt_; - mvaInputDPhi_[13] = pucU_; - mvaInputDPhi_[14] = pucPhi_; - mvaInputDPhi_[15] = jet1Pt_; - mvaInputDPhi_[16] = jet1Eta_; - mvaInputDPhi_[17] = jet1Phi_; - mvaInputDPhi_[18] = jet2Pt_; - mvaInputDPhi_[19] = jet2Eta_; - mvaInputDPhi_[20] = jet2Phi_; - mvaInputDPhi_[21] = numJets_; - mvaInputDPhi_[22] = numJetsPtGt30_; - mvaOutputDPhi_ = mvaReaderDPhi_->GetResponse(mvaInputDPhi_); +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; } -void PFMETAlgorithmMVA::evaluateCovU1() -{ - mvaInputCovU1_[0] = pfSumEt_; // PH: helps flattens response vs. Nvtx - mvaInputCovU1_[1] = numVertices_; - mvaInputCovU1_[2] = pfU_; - mvaInputCovU1_[3] = pfPhi_; - mvaInputCovU1_[4] = tkSumEt_; - mvaInputCovU1_[5] = tkU_; - mvaInputCovU1_[6] = tkPhi_; - mvaInputCovU1_[7] = npuSumEt_; - mvaInputCovU1_[8] = npuU_; - mvaInputCovU1_[9] = npuPhi_; - mvaInputCovU1_[10] = puSumEt_; - mvaInputCovU1_[11] = puMEt_; - mvaInputCovU1_[12] = puPhi_; - mvaInputCovU1_[13] = pucSumEt_; - mvaInputCovU1_[14] = pucU_; - mvaInputCovU1_[15] = pucPhi_; - mvaInputCovU1_[16] = jet1Pt_; - mvaInputCovU1_[17] = jet1Eta_; - mvaInputCovU1_[18] = jet1Phi_; - mvaInputCovU1_[19] = jet2Pt_; - mvaInputCovU1_[20] = jet2Eta_; - mvaInputCovU1_[21] = jet2Phi_; - mvaInputCovU1_[22] = numJets_; - mvaInputCovU1_[23] = numJetsPtGt30_; - mvaInputCovU1_[24] = pfPhi_ + mvaOutputDPhi_; - mvaInputCovU1_[25] = mvaOutputU_*pfU_; - mvaOutputCovU1_ = mvaReaderCovU1_->GetResponse(mvaInputCovU1_)*mvaOutputU_*pfU_; +//------------------------------------------------------------------------------- +const float PFMETAlgorithmMVA::GetResponse(const GBRForest * Reader, std::vector &variableNames ) +{ + std::unique_ptr mvaInputVector = createFloatVector(variableNames); + double result = Reader->GetResponse( mvaInputVector.get() ); + return result; } -void PFMETAlgorithmMVA::evaluateCovU2() -{ - mvaInputCovU2_[0] = pfSumEt_; // PH: helps flattens response vs. Nvtx - mvaInputCovU2_[1] = numVertices_; - mvaInputCovU2_[2] = pfU_; - mvaInputCovU2_[3] = pfPhi_; - mvaInputCovU2_[4] = tkSumEt_; - mvaInputCovU2_[5] = tkU_; - mvaInputCovU2_[6] = tkPhi_; - mvaInputCovU2_[7] = npuSumEt_; - mvaInputCovU2_[8] = npuU_; - mvaInputCovU2_[9] = npuPhi_; - mvaInputCovU2_[10] = puSumEt_; - mvaInputCovU2_[11] = puMEt_; - mvaInputCovU2_[12] = puPhi_; - mvaInputCovU2_[13] = pucSumEt_; - mvaInputCovU2_[14] = pucU_; - mvaInputCovU2_[15] = pucPhi_; - mvaInputCovU2_[16] = jet1Pt_; - mvaInputCovU2_[17] = jet1Eta_; - mvaInputCovU2_[18] = jet1Phi_; - mvaInputCovU2_[19] = jet2Pt_; - mvaInputCovU2_[20] = jet2Eta_; - mvaInputCovU2_[21] = jet2Phi_; - mvaInputCovU2_[22] = numJets_; - mvaInputCovU2_[23] = numJetsPtGt30_; - mvaInputCovU2_[24] = pfPhi_ + mvaOutputDPhi_; - mvaInputCovU2_[25] = mvaOutputU_*pfU_; - mvaOutputCovU2_ = mvaReaderCovU2_->GetResponse(mvaInputCovU2_)*mvaOutputU_*pfU_; -} +//------------------------------------------------------------------------------- void PFMETAlgorithmMVA::print(std::ostream& stream) const { stream << ":" << std::endl; - stream << " PF: sumEt = " << pfSumEt_ << ", U = " << pfU_ << ", phi = " << pfPhi_ << std::endl; - stream << " TK: sumEt = " << tkSumEt_ << ", U = " << tkU_ << ", phi = " << tkPhi_ << std::endl; - stream << " NPU: sumEt = " << npuSumEt_ << ", U = " << npuU_ << ", phi = " << npuPhi_ << std::endl; - stream << " PU: sumEt = " << puSumEt_ << ", MEt = " << puMEt_ << ", phi = " << puPhi_ << std::endl; - stream << " PUC: sumEt = " << pucSumEt_ << ", U = " << pucU_ << ", phi = " << pucPhi_ << std::endl; - stream << " jet1: Pt = " << jet1Pt_ << ", eta = " << jet1Eta_ << ", phi = " << jet1Phi_ << std::endl; - stream << " jet2: Pt = " << jet2Pt_ << ", eta = " << jet2Eta_ << ", phi = " << jet2Phi_ << std::endl; - stream << " num. jets = " << numJets_ << " (" << numJetsPtGt30_ << " with Pt > 30 GeV)" << std::endl; - stream << " num. vertices = " << numVertices_ << std::endl; - stream << " MVA output: U = " << mvaOutputU_ << ", dPhi = " << mvaOutputDPhi_ << "," - << " covU1 = " << mvaOutputCovU1_ << ", covU2 = " << mvaOutputCovU2_ << 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_ << ")" << std::endl; + << " phi = " << atan2(sumLeptonPy_, sumLeptonPx_) << " " + << "(Px = " << sumLeptonPx_ << ", Py = " << sumLeptonPy_ << ")" ; + stream << " MVA output: U = " << mvaOutputU_ << ", dPhi = " << mvaOutputDPhi_ << "," << " covU1 = " << mvaOutputCovU1_ << ", covU2 = " << mvaOutputCovU2_ << std::endl; + stream << std::endl; } diff --git a/RecoMET/METPUSubtraction/test/mvaMETOnMiniAOD_cfg.py b/RecoMET/METPUSubtraction/test/mvaMETOnMiniAOD_cfg.py index 40b7f0453ccb3..8a138d427f4b6 100644 --- a/RecoMET/METPUSubtraction/test/mvaMETOnMiniAOD_cfg.py +++ b/RecoMET/METPUSubtraction/test/mvaMETOnMiniAOD_cfg.py @@ -23,7 +23,7 @@ # Input source process.source = cms.Source("PoolSource", secondaryFileNames = cms.untracked.vstring(), - fileNames = cms.untracked.vstring('root://eoscms//eos/cms/store/relval/CMSSW_7_2_1/RelValZEE_13/MINIAODSIM/PU25ns_PHYS14_25_V1_Phys14-v1/00000/C6B792AB-D15E-E411-A787-02163E00F4FB.root'), + fileNames = cms.untracked.vstring('root://eoscms//eos/cms/store/relval/CMSSW_7_5_0_pre4/RelValZMM_13/MINIAODSIM/PU50ns_MCRUN2_75_V0-v1/00000/66CB72C0-70F8-E411-81DC-002590596468.root'), skipEvents = cms.untracked.uint32(0) ) diff --git a/RecoMET/METPUSubtraction/test/testMVAMetProducer.py b/RecoMET/METPUSubtraction/test/testMVAMetProducer.py index 482e19108104a..4f4b7f4fff56b 100644 --- a/RecoMET/METPUSubtraction/test/testMVAMetProducer.py +++ b/RecoMET/METPUSubtraction/test/testMVAMetProducer.py @@ -17,6 +17,7 @@ process.source = cms.Source("PoolSource", fileNames = cms.untracked.vstring( +"root://eoscms//eos/cms/store/relval/CMSSW_7_5_0_pre4/RelValZMM_13/GEN-SIM-RECO/PU50ns_MCRUN2_75_V0-v1/00000/24B4BCD7-22F8-E411-A740-0025905A48F2.root" ), skipEvents = cms.untracked.uint32(0) )