File indexing completed on 2023-03-17 11:20:15
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
0234
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
0244 if (hasPhotons_) {
0245
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;
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
0268
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 }
