File indexing completed on 2023-03-17 11:18:32
0001 #include <memory>
0002 #include "FWCore/Framework/interface/Frameworkfwd.h"
0003 #include "FWCore/Framework/interface/stream/EDProducer.h"
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/MakerMacros.h"
0006 #include "FWCore/Framework/interface/ESHandle.h"
0007 #include "FWCore/Framework/interface/EventSetup.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009
0010 #include "RecoJets/JetProducers/interface/MVAJetPuId.h"
0011
0012 #include "DataFormats/JetReco/interface/Jet.h"
0013 #include "DataFormats/Common/interface/ValueMap.h"
0014 #include "DataFormats/JetReco/interface/PileupJetIdentifier.h"
0015 #include "DataFormats/VertexReco/interface/Vertex.h"
0016
0017 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
0018 #include "CondFormats/JetMETObjects/interface/FactorizedJetCorrector.h"
0019 #include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"
0020
0021 #include "FWCore/ParameterSet/interface/FileInPath.h"
0022
0023 class MVAJetPuIdProducer : public edm::stream::EDProducer<> {
0024 public:
0025 explicit MVAJetPuIdProducer(const edm::ParameterSet &);
0026
0027 static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0028
0029 private:
0030 void produce(edm::Event &, const edm::EventSetup &) override;
0031
0032 void initJetEnergyCorrector(const edm::EventSetup &iSetup, bool isData);
0033
0034 edm::InputTag jets_, vertexes_, jetids_, rho_;
0035 std::string jec_;
0036 bool runMvas_, produceJetIds_, inputIsCorrected_, applyJec_;
0037 std::vector<std::pair<std::string, MVAJetPuId *>> algos_;
0038
0039 bool residualsFromTxt_;
0040 edm::FileInPath residualsTxt_;
0041 FactorizedJetCorrector *jecCor_;
0042 std::vector<JetCorrectorParameters> jetCorPars_;
0043
0044 edm::EDGetTokenT<edm::View<reco::Jet>> input_jet_token_;
0045 edm::EDGetTokenT<reco::VertexCollection> input_vertex_token_;
0046 edm::EDGetTokenT<edm::ValueMap<StoredPileupJetIdentifier>> input_vm_pujetid_token_;
0047 edm::EDGetTokenT<double> input_rho_token_;
0048 };
0049
0050 MVAJetPuIdProducer::MVAJetPuIdProducer(const edm::ParameterSet &iConfig) {
0051 runMvas_ = iConfig.getParameter<bool>("runMvas");
0052 produceJetIds_ = iConfig.getParameter<bool>("produceJetIds");
0053 jets_ = iConfig.getParameter<edm::InputTag>("jets");
0054 vertexes_ = iConfig.getParameter<edm::InputTag>("vertexes");
0055 jetids_ = iConfig.getParameter<edm::InputTag>("jetids");
0056 inputIsCorrected_ = iConfig.getParameter<bool>("inputIsCorrected");
0057 applyJec_ = iConfig.getParameter<bool>("applyJec");
0058 jec_ = iConfig.getParameter<std::string>("jec");
0059 rho_ = iConfig.getParameter<edm::InputTag>("rho");
0060 residualsFromTxt_ = iConfig.getParameter<bool>("residualsFromTxt");
0061 if (residualsFromTxt_)
0062 residualsTxt_ = iConfig.getParameter<edm::FileInPath>("residualsTxt");
0063 std::vector<edm::ParameterSet> algos = iConfig.getParameter<std::vector<edm::ParameterSet>>("algos");
0064
0065 jecCor_ = nullptr;
0066
0067 if (!runMvas_)
0068 assert(algos.size() == 1);
0069
0070 if (produceJetIds_) {
0071 produces<edm::ValueMap<StoredPileupJetIdentifier>>("");
0072 }
0073 for (std::vector<edm::ParameterSet>::iterator it = algos.begin(); it != algos.end(); ++it) {
0074 std::string label = it->getParameter<std::string>("label");
0075 algos_.push_back(std::make_pair(label, new MVAJetPuId(*it)));
0076 if (runMvas_) {
0077 produces<edm::ValueMap<float>>(label + "Discriminant");
0078 produces<edm::ValueMap<int>>(label + "Id");
0079 }
0080 }
0081
0082 input_jet_token_ = consumes<edm::View<reco::Jet>>(jets_);
0083 input_vertex_token_ = consumes<reco::VertexCollection>(vertexes_);
0084 input_vm_pujetid_token_ = consumes<edm::ValueMap<StoredPileupJetIdentifier>>(jetids_);
0085 input_rho_token_ = consumes<double>(rho_);
0086 }
0087
0088 void MVAJetPuIdProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0089 using namespace edm;
0090 using namespace std;
0091 using namespace reco;
0092
0093 Handle<View<Jet>> jetHandle;
0094 iEvent.getByToken(input_jet_token_, jetHandle);
0095 const View<Jet> &jets = *jetHandle;
0096 Handle<VertexCollection> vertexHandle;
0097 if (produceJetIds_) {
0098 iEvent.getByToken(input_vertex_token_, vertexHandle);
0099 }
0100 const VertexCollection &vertexes = *(vertexHandle.product());
0101 Handle<ValueMap<StoredPileupJetIdentifier>> vmap;
0102 if (!produceJetIds_) {
0103 iEvent.getByToken(input_vm_pujetid_token_, vmap);
0104 }
0105 edm::Handle<double> rhoH;
0106 double rho = 0.;
0107
0108 vector<StoredPileupJetIdentifier> ids;
0109 map<string, vector<float>> mvas;
0110 map<string, vector<int>> idflags;
0111
0112 VertexCollection::const_iterator vtx;
0113 if (produceJetIds_) {
0114 vtx = vertexes.begin();
0115 while (vtx != vertexes.end() && (vtx->isFake() || vtx->ndof() < 4)) {
0116 ++vtx;
0117 }
0118 if (vtx == vertexes.end()) {
0119 vtx = vertexes.begin();
0120 }
0121 }
0122
0123 for (unsigned int i = 0; i < jets.size(); ++i) {
0124 vector<pair<string, MVAJetPuId *>>::iterator algoi = algos_.begin();
0125 MVAJetPuId *ialgo = algoi->second;
0126
0127 const Jet &jet = jets[i];
0128
0129 float jec = 0.;
0130 if (applyJec_) {
0131 if (rho == 0.) {
0132 iEvent.getByToken(input_rho_token_, rhoH);
0133 rho = *rhoH;
0134 }
0135 jecCor_->setJetPt(jet.pt());
0136 jecCor_->setJetEta(jet.eta());
0137 jecCor_->setJetA(jet.jetArea());
0138 jecCor_->setRho(rho);
0139 jec = jecCor_->getCorrection();
0140 }
0141
0142 bool applyJec = applyJec_ || !inputIsCorrected_;
0143 reco::Jet *corrJet = nullptr;
0144 if (applyJec) {
0145 float scale = jec;
0146 corrJet = dynamic_cast<reco::Jet *>(jet.clone());
0147 corrJet->scaleEnergy(scale);
0148 }
0149 const reco::Jet *theJet = (applyJec ? corrJet : &jet);
0150 PileupJetIdentifier puIdentifier;
0151 if (produceJetIds_) {
0152 puIdentifier = ialgo->computeIdVariables(theJet, jec, &(*vtx), vertexes, runMvas_);
0153 ids.push_back(puIdentifier);
0154 } else {
0155 puIdentifier = (*vmap)[jets.refAt(i)];
0156 puIdentifier.jetPt(theJet->pt());
0157 puIdentifier.jetEta(theJet->eta());
0158 puIdentifier.jetPhi(theJet->phi());
0159 ialgo->set(puIdentifier);
0160 puIdentifier = ialgo->computeMva();
0161 }
0162
0163 if (runMvas_) {
0164 for (; algoi != algos_.end(); ++algoi) {
0165 ialgo = algoi->second;
0166 ialgo->set(puIdentifier);
0167 PileupJetIdentifier id = ialgo->computeMva();
0168 mvas[algoi->first].push_back(id.mva());
0169 idflags[algoi->first].push_back(id.idFlag());
0170 }
0171 }
0172
0173 if (corrJet) {
0174 delete corrJet;
0175 }
0176 }
0177
0178 if (runMvas_) {
0179 for (vector<pair<string, MVAJetPuId *>>::iterator ialgo = algos_.begin(); ialgo != algos_.end(); ++ialgo) {
0180 vector<float> &mva = mvas[ialgo->first];
0181 auto mvaout = std::make_unique<ValueMap<float>>();
0182 ValueMap<float>::Filler mvafiller(*mvaout);
0183 mvafiller.insert(jetHandle, mva.begin(), mva.end());
0184 mvafiller.fill();
0185 iEvent.put(std::move(mvaout), ialgo->first + "Discriminant");
0186
0187 vector<int> &idflag = idflags[ialgo->first];
0188 auto idflagout = std::make_unique<ValueMap<int>>();
0189 ValueMap<int>::Filler idflagfiller(*idflagout);
0190 idflagfiller.insert(jetHandle, idflag.begin(), idflag.end());
0191 idflagfiller.fill();
0192 iEvent.put(std::move(idflagout), ialgo->first + "Id");
0193 }
0194 }
0195 if (produceJetIds_) {
0196 assert(jetHandle->size() == ids.size());
0197 auto idsout = std::make_unique<ValueMap<StoredPileupJetIdentifier>>();
0198 ValueMap<StoredPileupJetIdentifier>::Filler idsfiller(*idsout);
0199 idsfiller.insert(jetHandle, ids.begin(), ids.end());
0200 idsfiller.fill();
0201 iEvent.put(std::move(idsout));
0202 }
0203 }
0204
0205 void MVAJetPuIdProducer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0206 edm::ParameterSetDescription desc;
0207 desc.add<bool>("runMvas", true);
0208 desc.add<bool>("inputIsCorrected", true);
0209 desc.add<edm::InputTag>("vertexes", edm::InputTag("hltPixelVertices"));
0210 desc.add<bool>("produceJetIds", true);
0211 desc.add<std::string>("jec", "AK4PF");
0212 desc.add<bool>("residualsFromTxt", false);
0213 desc.add<bool>("applyJec", false);
0214 desc.add<edm::InputTag>("jetids", edm::InputTag(""));
0215 desc.add<edm::InputTag>("rho", edm::InputTag("hltFixedGridRhoFastjetAll"));
0216 desc.add<edm::InputTag>("jets", edm::InputTag("hltAK4PFJetsCorrected"));
0217 edm::ParameterSetDescription vpsd1;
0218 vpsd1.add<std::vector<std::string>>("tmvaVariables",
0219 {
0220 "rho",
0221 "nParticles",
0222 "nCharged",
0223 "majW",
0224 "minW",
0225 "frac01",
0226 "frac02",
0227 "frac03",
0228 "frac04",
0229 "ptD",
0230 "beta",
0231 "betaStar",
0232 "dR2Mean",
0233 "pull",
0234 "jetR",
0235 "jetRchg",
0236 });
0237 vpsd1.add<std::string>("tmvaMethod", "JetID");
0238 vpsd1.add<bool>("cutBased", false);
0239 vpsd1.add<std::string>("tmvaWeights", "RecoJets/JetProducers/data/MVAJetPuID.weights.xml.gz");
0240 vpsd1.add<std::vector<std::string>>("tmvaSpectators",
0241 {
0242 "jetEta",
0243 "jetPt",
0244 });
0245 vpsd1.add<std::string>("label", "CATEv0");
0246 vpsd1.add<int>("version", -1);
0247 {
0248 edm::ParameterSetDescription psd0;
0249 psd0.add<std::vector<double>>("Pt2030_Tight",
0250 {
0251 0.73,
0252 0.05,
0253 -0.26,
0254 -0.42,
0255 });
0256 psd0.add<std::vector<double>>("Pt2030_Loose",
0257 {
0258 -0.63,
0259 -0.6,
0260 -0.55,
0261 -0.45,
0262 });
0263 psd0.add<std::vector<double>>("Pt3050_Medium",
0264 {
0265 0.1,
0266 -0.36,
0267 -0.54,
0268 -0.54,
0269 });
0270 psd0.add<std::vector<double>>("Pt1020_Tight",
0271 {
0272 -0.83,
0273 -0.81,
0274 -0.74,
0275 -0.81,
0276 });
0277 psd0.add<std::vector<double>>("Pt2030_Medium",
0278 {
0279 0.1,
0280 -0.36,
0281 -0.54,
0282 -0.54,
0283 });
0284 psd0.add<std::vector<double>>("Pt010_Tight",
0285 {
0286 -0.83,
0287 -0.81,
0288 -0.74,
0289 -0.81,
0290 });
0291 psd0.add<std::vector<double>>("Pt1020_Loose",
0292 {
0293 -0.95,
0294 -0.96,
0295 -0.94,
0296 -0.95,
0297 });
0298 psd0.add<std::vector<double>>("Pt010_Medium",
0299 {
0300 -0.83,
0301 -0.92,
0302 -0.9,
0303 -0.92,
0304 });
0305 psd0.add<std::vector<double>>("Pt1020_Medium",
0306 {
0307 -0.83,
0308 -0.92,
0309 -0.9,
0310 -0.92,
0311 });
0312 psd0.add<std::vector<double>>("Pt010_Loose",
0313 {
0314 -0.95,
0315 -0.96,
0316 -0.94,
0317 -0.95,
0318 });
0319 psd0.add<std::vector<double>>("Pt3050_Loose",
0320 {
0321 -0.63,
0322 -0.6,
0323 -0.55,
0324 -0.45,
0325 });
0326 psd0.add<std::vector<double>>("Pt3050_Tight",
0327 {
0328 0.73,
0329 0.05,
0330 -0.26,
0331 -0.42,
0332 });
0333 vpsd1.add<edm::ParameterSetDescription>("JetIdParams", psd0);
0334 }
0335 vpsd1.add<double>("impactParTkThreshold", 1.0);
0336 std::vector<edm::ParameterSet> temp1;
0337 temp1.reserve(1);
0338
0339 desc.addVPSet("algos", vpsd1, temp1);
0340
0341 descriptions.add("MVAJetPuIdProducer", desc);
0342 }
0343
0344 DEFINE_FWK_MODULE(MVAJetPuIdProducer);
