Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:35:04

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()); /*make sure JEC is applied when computing the MVA*/
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);