Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-08-23 03:25:23

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       usePuppi_(iConfig.getParameter<bool>("usePuppi")) {
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 
0045 // ------------------------------------------------------------------------------------------
0046 PileupJetIdProducer::PileupJetIdProducer(const edm::ParameterSet& iConfig, GBRForestsAndConstants const* globalCache) {
0047   if (globalCache->produceJetIds()) {
0048     produces<edm::ValueMap<StoredPileupJetIdentifier>>("");
0049   }
0050   for (auto const& algoGBRForestsAndConstants : globalCache->vAlgoGBRForestsAndConstants()) {
0051     std::string const& label = algoGBRForestsAndConstants.label();
0052     algos_.emplace_back(label, std::make_unique<PileupJetIdAlgo>(&algoGBRForestsAndConstants));
0053     if (globalCache->runMvas()) {
0054       produces<edm::ValueMap<float>>(label + "Discriminant");
0055       produces<edm::ValueMap<int>>(label + "Id");
0056     }
0057   }
0058 
0059   input_jet_token_ = consumes<edm::View<reco::Jet>>(iConfig.getParameter<edm::InputTag>("jets"));
0060   input_vertex_token_ = consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertexes"));
0061   input_vm_pujetid_token_ =
0062       consumes<edm::ValueMap<StoredPileupJetIdentifier>>(iConfig.getParameter<edm::InputTag>("jetids"));
0063   input_rho_token_ = consumes<double>(iConfig.getParameter<edm::InputTag>("rho"));
0064   parameters_token_ = esConsumes(edm::ESInputTag("", globalCache->jec()));
0065 }
0066 
0067 // ------------------------------------------------------------------------------------------
0068 PileupJetIdProducer::~PileupJetIdProducer() {}
0069 
0070 // ------------------------------------------------------------------------------------------
0071 void PileupJetIdProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0072   GBRForestsAndConstants const* gc = globalCache();
0073 
0074   using namespace edm;
0075   using namespace std;
0076   using namespace reco;
0077 
0078   // Input jets
0079   Handle<View<Jet>> jetHandle;
0080   iEvent.getByToken(input_jet_token_, jetHandle);
0081   const View<Jet>& jets = *jetHandle;
0082 
0083   // input variables
0084   Handle<ValueMap<StoredPileupJetIdentifier>> vmap;
0085   if (!gc->produceJetIds()) {
0086     iEvent.getByToken(input_vm_pujetid_token_, vmap);
0087   }
0088   // rho
0089   edm::Handle<double> rhoH;
0090   double rho = 0.;
0091 
0092   // products
0093   vector<StoredPileupJetIdentifier> ids;
0094   map<string, vector<float>> mvas;
0095   map<string, vector<int>> idflags;
0096 
0097   const VertexCollection* vertexes = nullptr;
0098   VertexCollection::const_iterator vtx;
0099   if (gc->produceJetIds()) {
0100     // vertexes
0101     Handle<VertexCollection> vertexHandle;
0102     iEvent.getByToken(input_vertex_token_, vertexHandle);
0103 
0104     vertexes = vertexHandle.product();
0105 
0106     // require basic quality cuts on the vertexes
0107     vtx = vertexes->begin();
0108     while (vtx != vertexes->end() && (vtx->isFake() || vtx->ndof() < 4)) {
0109       ++vtx;
0110     }
0111     if (vtx == vertexes->end()) {
0112       vtx = vertexes->begin();
0113     }
0114   }
0115 
0116   // Loop over input jets
0117   bool ispat = true;
0118   for (unsigned int i = 0; i < jets.size(); ++i) {
0119     // Pick the first algo to compute the input variables
0120     auto algoi = algos_.begin();
0121     PileupJetIdAlgo* ialgo = algoi->second.get();
0122 
0123     const Jet& jet = jets.at(i);
0124     const pat::Jet* patjet = nullptr;
0125     if (ispat) {
0126       patjet = dynamic_cast<const pat::Jet*>(&jet);
0127       ispat = patjet != nullptr;
0128     }
0129 
0130     // Get jet energy correction
0131     float jec = 0.;
0132     if (gc->applyJec()) {
0133       // If haven't done it get rho from the event
0134       if (rho == 0.) {
0135         iEvent.getByToken(input_rho_token_, rhoH);
0136         rho = *rhoH;
0137       }
0138       // jet corrector
0139       if (not jecCor_) {
0140         initJetEnergyCorrector(iSetup, iEvent.isRealData());
0141       }
0142       if (ispat) {
0143         jecCor_->setJetPt(patjet->correctedJet(0).pt());
0144       } else {
0145         jecCor_->setJetPt(jet.pt());
0146       }
0147       jecCor_->setJetEta(jet.eta());
0148       jecCor_->setJetA(jet.jetArea());
0149       jecCor_->setRho(rho);
0150       jec = jecCor_->getCorrection();
0151     }
0152     // If it was requested AND the input is an uncorrected jet apply the JEC
0153     bool applyJec = gc->applyJec() && (ispat || !gc->inputIsCorrected());
0154     std::unique_ptr<reco::Jet> corrJet;
0155 
0156     if (applyJec) {
0157       float scale = jec;
0158       if (ispat) {
0159         corrJet = std::make_unique<pat::Jet>(patjet->correctedJet(0));
0160       } else {
0161         corrJet.reset(dynamic_cast<reco::Jet*>(jet.clone()));
0162       }
0163       corrJet->scaleEnergy(scale);
0164     }
0165     const reco::Jet* theJet = (applyJec ? corrJet.get() : &jet);
0166 
0167     PileupJetIdentifier puIdentifier;
0168     if (gc->produceJetIds()) {
0169       // Compute the input variables
0170       puIdentifier = ialgo->computeIdVariables(theJet, jec, &(*vtx), *vertexes, rho, gc->usePuppi());
0171       ids.push_back(puIdentifier);
0172     } else {
0173       // Or read it from the value map
0174       puIdentifier = (*vmap)[jets.refAt(i)];
0175       puIdentifier.jetPt(theJet->pt());  // make sure JEC is applied when computing the MVA
0176       puIdentifier.jetEta(theJet->eta());
0177       puIdentifier.jetPhi(theJet->phi());
0178       ialgo->set(puIdentifier);
0179       puIdentifier = ialgo->computeMva();
0180     }
0181 
0182     if (gc->runMvas()) {
0183       // Compute the MVA and WP
0184       mvas[algoi->first].push_back(puIdentifier.mva());
0185       idflags[algoi->first].push_back(puIdentifier.idFlag());
0186       for (++algoi; algoi != algos_.end(); ++algoi) {
0187         ialgo = algoi->second.get();
0188         ialgo->set(puIdentifier);
0189         PileupJetIdentifier id = ialgo->computeMva();
0190         mvas[algoi->first].push_back(id.mva());
0191         idflags[algoi->first].push_back(id.idFlag());
0192       }
0193     }
0194   }
0195 
0196   // Produce the output value maps
0197   if (gc->runMvas()) {
0198     for (const auto& ialgo : algos_) {
0199       // MVA
0200       vector<float>& mva = mvas[ialgo.first];
0201       auto mvaout = std::make_unique<ValueMap<float>>();
0202       ValueMap<float>::Filler mvafiller(*mvaout);
0203       mvafiller.insert(jetHandle, mva.begin(), mva.end());
0204       mvafiller.fill();
0205       iEvent.put(std::move(mvaout), ialgo.first + "Discriminant");
0206 
0207       // WP
0208       vector<int>& idflag = idflags[ialgo.first];
0209       auto idflagout = std::make_unique<ValueMap<int>>();
0210       ValueMap<int>::Filler idflagfiller(*idflagout);
0211       idflagfiller.insert(jetHandle, idflag.begin(), idflag.end());
0212       idflagfiller.fill();
0213       iEvent.put(std::move(idflagout), ialgo.first + "Id");
0214     }
0215   }
0216   // input variables
0217   if (gc->produceJetIds()) {
0218     assert(jetHandle->size() == ids.size());
0219     auto idsout = std::make_unique<ValueMap<StoredPileupJetIdentifier>>();
0220     ValueMap<StoredPileupJetIdentifier>::Filler idsfiller(*idsout);
0221     idsfiller.insert(jetHandle, ids.begin(), ids.end());
0222     idsfiller.fill();
0223     iEvent.put(std::move(idsout));
0224   }
0225 }
0226 
0227 // ------------------------------------------------------------------------------------------
0228 void PileupJetIdProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0229   //The following says we do not know what parameters are allowed so do no validation
0230   // Please change this to state exactly what you do use, even if it is no parameters
0231   edm::ParameterSetDescription desc;
0232   desc.setUnknown();
0233   descriptions.addDefault(desc);
0234 }
0235 
0236 // ------------------------------------------------------------------------------------------
0237 void PileupJetIdProducer::initJetEnergyCorrector(const edm::EventSetup& iSetup, bool isData) {
0238   GBRForestsAndConstants const* gc = globalCache();
0239 
0240   //jet energy correction levels to apply on raw jet
0241   std::vector<std::string> jecLevels;
0242   jecLevels.push_back("L1FastJet");
0243   jecLevels.push_back("L2Relative");
0244   jecLevels.push_back("L3Absolute");
0245   if (isData && !gc->residualsFromTxt())
0246     jecLevels.push_back("L2L3Residual");
0247 
0248   //check the corrector parameters needed according to the correction levels
0249   auto const& parameters = iSetup.getData(parameters_token_);
0250   for (std::vector<std::string>::const_iterator ll = jecLevels.begin(); ll != jecLevels.end(); ++ll) {
0251     const JetCorrectorParameters& ip = parameters[*ll];
0252     jetCorPars_.push_back(ip);
0253   }
0254   if (isData && gc->residualsFromTxt()) {
0255     jetCorPars_.push_back(JetCorrectorParameters(gc->residualsTxt().fullPath()));
0256   }
0257 
0258   //instantiate the jet corrector
0259   jecCor_ = std::make_unique<FactorizedJetCorrector>(jetCorPars_);
0260 }
0261 //define this as a plug-in
0262 DEFINE_FWK_MODULE(PileupJetIdProducer);