JetProducers/plugins/PileupJetIdProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    PileupJetIdProducer
0004 // Class:      PileupJetIdProducer
0005 //
0006 /**\class PileupJetIdProducer PileupJetIdProducer.cc CMGTools/PileupJetIdProducer/src/PileupJetIdProducer.cc
0007
0008 Description: [one line class summary]
0009
0010 Implementation:
0011 [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Pasquale Musella,40 2-A12,+41227671706,
0015 //         Created:  Wed Apr 18 15:48:47 CEST 2012
0016 //
0017 //
0018
0019 #include "RecoJets/JetProducers/plugins/PileupJetIdProducer.h"
0020
0021 #include <memory>
0022
0023 GBRForestsAndConstants::GBRForestsAndConstants(edm::ParameterSet const& iConfig)
0024     : runMvas_(iConfig.getParameter<bool>("runMvas")),
0025       produceJetIds_(iConfig.getParameter<bool>("produceJetIds")),
0026       inputIsCorrected_(iConfig.getParameter<bool>("inputIsCorrected")),
0027       applyJec_(iConfig.getParameter<bool>("applyJec")),
0028       jec_(iConfig.getParameter<std::string>("jec")),
0029       residualsFromTxt_(iConfig.getParameter<bool>("residualsFromTxt")),
0030       applyConstituentWeight_(false) {
0031   if (residualsFromTxt_) {
0032     residualsTxt_ = iConfig.getParameter<edm::FileInPath>("residualsTxt");
0033   }
0034
0035   std::vector<edm::ParameterSet> algos = iConfig.getParameter<std::vector<edm::ParameterSet>>("algos");
0036   for (auto const& algoPset : algos) {
0037     vAlgoGBRForestsAndConstants_.emplace_back(algoPset, runMvas_);
0038   }
0039
0040   if (!runMvas_) {
0041     assert(algos.size() == 1);
0042   }
0043
0044   edm::InputTag srcConstituentWeights = iConfig.getParameter<edm::InputTag>("srcConstituentWeights");
0045   if (!srcConstituentWeights.label().empty()) {
0046     applyConstituentWeight_ = true;
0047   }
0048 }
0049
0050 // ------------------------------------------------------------------------------------------
0051 PileupJetIdProducer::PileupJetIdProducer(const edm::ParameterSet& iConfig, GBRForestsAndConstants const* globalCache) {
0052   if (globalCache->produceJetIds()) {
0053     produces<edm::ValueMap<StoredPileupJetIdentifier>>("");
0054   }
0055   for (auto const& algoGBRForestsAndConstants : globalCache->vAlgoGBRForestsAndConstants()) {
0056     std::string const& label = algoGBRForestsAndConstants.label();
0057     algos_.emplace_back(label, std::make_unique<PileupJetIdAlgo>(&algoGBRForestsAndConstants));
0058     if (globalCache->runMvas()) {
0059       produces<edm::ValueMap<float>>(label + "Discriminant");
0060       produces<edm::ValueMap<int>>(label + "Id");
0061     }
0062   }
0063
0064   input_jet_token_ = consumes<edm::View<reco::Jet>>(iConfig.getParameter<edm::InputTag>("jets"));
0065   input_vertex_token_ = consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertexes"));
0066   input_vm_pujetid_token_ =
0067       consumes<edm::ValueMap<StoredPileupJetIdentifier>>(iConfig.getParameter<edm::InputTag>("jetids"));
0068   input_rho_token_ = consumes<double>(iConfig.getParameter<edm::InputTag>("rho"));
0069   parameters_token_ = esConsumes(edm::ESInputTag("", globalCache->jec()));
0070
0071   edm::InputTag srcConstituentWeights = iConfig.getParameter<edm::InputTag>("srcConstituentWeights");
0072   if (!srcConstituentWeights.label().empty()) {
0073     input_constituent_weights_token_ = consumes<edm::ValueMap<float>>(srcConstituentWeights);
0074   }
0075 }
0076
0077 // ------------------------------------------------------------------------------------------
0078 PileupJetIdProducer::~PileupJetIdProducer() {}
0079
0080 // ------------------------------------------------------------------------------------------
0081 void PileupJetIdProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0082   GBRForestsAndConstants const* gc = globalCache();
0083
0084   using namespace edm;
0085   using namespace std;
0086   using namespace reco;
0087
0088   // Input jets
0089   Handle<View<Jet>> jetHandle;
0090   iEvent.getByToken(input_jet_token_, jetHandle);
0091   const View<Jet>& jets = *jetHandle;
0092
0093   // Constituent weight (e.g PUPPI) Value Map
0094   edm::ValueMap<float> constituentWeights;
0095   if (!input_constituent_weights_token_.isUninitialized()) {
0096     constituentWeights = iEvent.get(input_constituent_weights_token_);
0097   }
0098
0099   // input variables
0100   Handle<ValueMap<StoredPileupJetIdentifier>> vmap;
0101   if (!gc->produceJetIds()) {
0102     iEvent.getByToken(input_vm_pujetid_token_, vmap);
0103   }
0104   // rho
0105   edm::Handle<double> rhoH;
0106   double rho = 0.;
0107
0108   // products
0109   vector<StoredPileupJetIdentifier> ids;
0110   map<string, vector<float>> mvas;
0111   map<string, vector<int>> idflags;
0112
0113   const VertexCollection* vertexes = nullptr;
0114   VertexCollection::const_iterator vtx;
0115   if (gc->produceJetIds()) {
0116     // vertexes
0117     Handle<VertexCollection> vertexHandle;
0118     iEvent.getByToken(input_vertex_token_, vertexHandle);
0119
0120     vertexes = vertexHandle.product();
0121
0122     // require basic quality cuts on the vertexes
0123     vtx = vertexes->begin();
0124     while (vtx != vertexes->end() && (vtx->isFake() || vtx->ndof() < 4)) {
0125       ++vtx;
0126     }
0127     if (vtx == vertexes->end()) {
0128       vtx = vertexes->begin();
0129     }
0130   }
0131
0132   // Loop over input jets
0133   bool ispat = true;
0134   for (unsigned int i = 0; i < jets.size(); ++i) {
0135     // Pick the first algo to compute the input variables
0136     auto algoi = algos_.begin();
0137     PileupJetIdAlgo* ialgo = algoi->second.get();
0138
0139     const Jet& jet = jets.at(i);
0140     const pat::Jet* patjet = nullptr;
0141     if (ispat) {
0142       patjet = dynamic_cast<const pat::Jet*>(&jet);
0143       ispat = patjet != nullptr;
0144     }
0145
0146     // Get jet energy correction
0147     float jec = 0.;
0148     if (gc->applyJec()) {
0149       // If haven't done it get rho from the event
0150       if (rho == 0.) {
0151         iEvent.getByToken(input_rho_token_, rhoH);
0152         rho = *rhoH;
0153       }
0154       // jet corrector
0155       if (not jecCor_) {
0156         initJetEnergyCorrector(iSetup, iEvent.isRealData());
0157       }
0158       if (ispat) {
0159         jecCor_->setJetPt(patjet->correctedJet(0).pt());
0160       } else {
0161         jecCor_->setJetPt(jet.pt());
0162       }
0163       jecCor_->setJetEta(jet.eta());
0164       jecCor_->setJetPhi(jet.phi());
0165       jecCor_->setJetA(jet.jetArea());
0166       jecCor_->setRho(rho);
0167       jec = jecCor_->getCorrection();
0168     }
0169     // If it was requested AND the input is an uncorrected jet apply the JEC
0170     bool applyJec = gc->applyJec() && (ispat || !gc->inputIsCorrected());
0171     std::unique_ptr<reco::Jet> corrJet;
0172
0173     if (applyJec) {
0174       float scale = jec;
0175       if (ispat) {
0176         corrJet = std::make_unique<pat::Jet>(patjet->correctedJet(0));
0177       } else {
0178         corrJet.reset(dynamic_cast<reco::Jet*>(jet.clone()));
0179       }
0180       corrJet->scaleEnergy(scale);
0181     }
0182     const reco::Jet* theJet = (applyJec ? corrJet.get() : &jet);
0183
0184     PileupJetIdentifier puIdentifier;
0185     if (gc->produceJetIds()) {
0186       // Compute the input variables
0187       ////////////////////////////// added PUPPI weight Value Map
0188       puIdentifier = ialgo->computeIdVariables(
0189           theJet, jec, &(*vtx), *vertexes, rho, constituentWeights, gc->applyConstituentWeight());
0190       ids.push_back(puIdentifier);
0191     } else {
0192       // Or read it from the value map
0193       puIdentifier = (*vmap)[jets.refAt(i)];
0194       puIdentifier.jetPt(theJet->pt());  // make sure JEC is applied when computing the MVA
0195       puIdentifier.jetEta(theJet->eta());
0196       puIdentifier.jetPhi(theJet->phi());
0197       ialgo->set(puIdentifier);
0198       puIdentifier = ialgo->computeMva();
0199     }
0200
0201     if (gc->runMvas()) {
0202       // Compute the MVA and WP
0203       mvas[algoi->first].push_back(puIdentifier.mva());
0204       idflags[algoi->first].push_back(puIdentifier.idFlag());
0205       for (++algoi; algoi != algos_.end(); ++algoi) {
0206         ialgo = algoi->second.get();
0207         ialgo->set(puIdentifier);
0208         PileupJetIdentifier id = ialgo->computeMva();
0209         mvas[algoi->first].push_back(id.mva());
0210         idflags[algoi->first].push_back(id.idFlag());
0211       }
0212     }
0213   }
0214
0215   // Produce the output value maps
0216   if (gc->runMvas()) {
0217     for (const auto& ialgo : algos_) {
0218       // MVA
0219       vector<float>& mva = mvas[ialgo.first];
0220       auto mvaout = std::make_unique<ValueMap<float>>();
0221       ValueMap<float>::Filler mvafiller(*mvaout);
0222       mvafiller.insert(jetHandle, mva.begin(), mva.end());
0223       mvafiller.fill();
0224       iEvent.put(std::move(mvaout), ialgo.first + "Discriminant");
0225
0226       // WP
0227       vector<int>& idflag = idflags[ialgo.first];
0228       auto idflagout = std::make_unique<ValueMap<int>>();
0229       ValueMap<int>::Filler idflagfiller(*idflagout);
0230       idflagfiller.insert(jetHandle, idflag.begin(), idflag.end());
0231       idflagfiller.fill();
0232       iEvent.put(std::move(idflagout), ialgo.first + "Id");
0233     }
0234   }
0235   // input variables
0236   if (gc->produceJetIds()) {
0237     assert(jetHandle->size() == ids.size());
0238     auto idsout = std::make_unique<ValueMap<StoredPileupJetIdentifier>>();
0239     ValueMap<StoredPileupJetIdentifier>::Filler idsfiller(*idsout);
0240     idsfiller.insert(jetHandle, ids.begin(), ids.end());
0241     idsfiller.fill();
0242     iEvent.put(std::move(idsout));
0243   }
0244 }
0245
0246 // ------------------------------------------------------------------------------------------
0247 void PileupJetIdProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0248   //The following says we do not know what parameters are allowed so do no validation
0249   // Please change this to state exactly what you do use, even if it is no parameters
0250   edm::ParameterSetDescription desc;
0251   desc.setUnknown();
0252   descriptions.addDefault(desc);
0253 }
0254
0255 // ------------------------------------------------------------------------------------------
0256 void PileupJetIdProducer::initJetEnergyCorrector(const edm::EventSetup& iSetup, bool isData) {
0257   GBRForestsAndConstants const* gc = globalCache();
0258
0259   //jet energy correction levels to apply on raw jet
0260   std::vector<std::string> jecLevels;
0261   jecLevels.push_back("L1FastJet");
0262   jecLevels.push_back("L2Relative");
0263   jecLevels.push_back("L3Absolute");
0264   if (isData && !gc->residualsFromTxt())
0265     jecLevels.push_back("L2L3Residual");
0266
0267   //check the corrector parameters needed according to the correction levels
0268   auto const& parameters = iSetup.getData(parameters_token_);
0269   for (std::vector<std::string>::const_iterator ll = jecLevels.begin(); ll != jecLevels.end(); ++ll) {
0270     const JetCorrectorParameters& ip = parameters[*ll];
0271     jetCorPars_.push_back(ip);
0272   }
0273   if (isData && gc->residualsFromTxt()) {
0274     jetCorPars_.push_back(JetCorrectorParameters(gc->residualsTxt().fullPath()));
0275   }
0276
0277   //instantiate the jet corrector
0278   jecCor_ = std::make_unique<FactorizedJetCorrector>(jetCorPars_);
0279 }
0280 //define this as a plug-in
0281 DEFINE_FWK_MODULE(PileupJetIdProducer);