PatAlgos/plugins/JetCorrFactorsProducer.cc

0001 /**
0002   \class    pat::JetCorrFactorsProducer JetCorrFactorsProducer.h "PhysicsTools/PatAlgos/interface/JetCorrFactorsProducer.h"
0003   \brief    Produces a ValueMap between JetCorrFactors and the to the originating reco jets
0004
0005    The JetCorrFactorsProducer produces a set of correction factors, defined in the class pat::JetCorrFactors. This vector
0006    is linked to the originating reco jets through an edm::ValueMap. The initializing parameters of the module can be found
0007    in the recoLayer1/jetCorrFactors_cfi.py of the PatAlgos package. In the standard PAT workflow the module has to be run
0008    before the creation of the pat::Jet. The edm::ValueMap will then be embedded into the pat::Jet.
0009
0010    Jets corrected up to a given correction level can then be accessed via the pat::Jet member function correctedJet. For
0011    more details have a look into the class description of the pat::Jet.
0012
0013    ATTENTION: available options for flavor corrections are
0014     * L5Flavor_gJ        L7Parton_gJ         gluon   from dijets
0015     * L5Flavor_qJ/_qT    L7Parton_qJ/_qT     quark   from dijets/top
0016     * L5Flavor_cJ/_cT    L7Parton_cJ/_cT     charm   from dijets/top
0017     * L5Flavor_bJ/_bT    L7Parton_bJ/_bT     beauty  from dijets/top
0018     *                    L7Parton_jJ/_tT     mixture from dijets/top
0019
0020    where mixture refers to the flavor mixture as determined from the MC sample the flavor dependent corrections have been
0021    derived from. 'J' and 'T' stand for a typical dijet (ttbar) sample.
0022
0023    L1Offset corrections require the collection of _offlinePrimaryVertices_, which are supposed to be added as an additional
0024    optional parameter _primaryVertices_ in the jetCorrFactors_cfi.py file.
0025
0026    L1FastJet corrections, which are an alternative to the standard L1Offset correction as recommended by the JetMET PAG the
0027    energy density parameter _rho_ is supposed to be added as an additional optional parameter _rho_ in the
0028    jetCorrFactors_cfi.py file.
0029
0030    NOTE:
0031     * the mixed mode (mc input mixture from dijets/ttbar) only exists for parton level corrections.
0032     * jJ and tT are not covered in this implementation of the JetCorrFactorsProducer
0033     * there are no gluon corrections available from the top sample on the L7Parton level.
0034 */
0035
0036 #include "CondFormats/JetMETObjects/interface/FactorizedJetCorrector.h"
0037 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
0038 #include "DataFormats/Common/interface/ValueMap.h"
0039 #include "DataFormats/Common/interface/View.h"
0040 #include "DataFormats/JetReco/interface/CaloJet.h"
0041 #include "DataFormats/JetReco/interface/JPTJet.h"
0042 #include "DataFormats/JetReco/interface/Jet.h"
0043 #include "DataFormats/PatCandidates/interface/Jet.h"
0044 #include "DataFormats/PatCandidates/interface/JetCorrFactors.h"
0045 #include "DataFormats/VertexReco/interface/Vertex.h"
0046 #include "FWCore/Framework/interface/ESHandle.h"
0047 #include "FWCore/Framework/interface/ESWatcher.h"
0048 #include "FWCore/Framework/interface/Event.h"
0049 #include "FWCore/Framework/interface/EventSetup.h"
0050 #include "FWCore/Framework/interface/stream/EDProducer.h"
0051 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0052 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0053 #include "FWCore/ParameterSet/interface/ParameterDescription.h"
0054 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0055 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0056 #include "FWCore/Utilities/interface/InputTag.h"
0057 #include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"
0058
0059 namespace pat {
0060
0061   class JetCorrFactorsProducer : public edm::stream::EDProducer<> {
0062   public:
0063     /// value map for JetCorrFactors (to be written into the event)
0064     typedef edm::ValueMap<pat::JetCorrFactors> JetCorrFactorsMap;
0065     /// map of correction levels to different flavors
0066     typedef std::map<JetCorrFactors::Flavor, std::vector<std::string> > FlavorCorrLevelMap;
0067
0068   public:
0069     /// default constructor
0070     explicit JetCorrFactorsProducer(const edm::ParameterSet& cfg);
0071     /// default destructor
0072     ~JetCorrFactorsProducer() override{};
0073     /// everything that needs to be done per event
0074     void produce(edm::Event& event, const edm::EventSetup& setup) override;
0075     /// description of configuration file parameters
0076     static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0077
0078   private:
0079     /// return true if the jec levels contain at least one flavor dependent correction level
0080     bool flavorDependent() const { return (levels_.size() > 1); };
0081     /// return the jec parameters as input to the FactorizedJetCorrector for different flavors
0082     std::vector<JetCorrectorParameters> params(const JetCorrectorParametersCollection& parameters,
0083                                                const std::vector<std::string>& levels) const;
0084     /// return an expanded version of correction levels for different flavors; the result should
0085     /// be of type ['L2Relative', 'L3Absolute', 'L5FLavor_gJ', 'L7Parton_gJ']; L7Parton_gT will
0086     /// result in an empty string as this correction level is not available
0087     std::vector<std::string> expand(const std::vector<std::string>& levels, const JetCorrFactors::Flavor& flavor);
0088     /// evaluate jet correction factor up to a given level
0089     float evaluate(edm::View<reco::Jet>::const_iterator& jet, const JetCorrFactors::Flavor& flavor, int level);
0090     /// determines the number of valid primary vertices for the standard L1Offset correction of JetMET
0091     int numberOf(const edm::Handle<std::vector<reco::Vertex> >& primaryVertices);
0092     /// map jet algorithm to payload in DB
0093
0094   private:
0095     /// use electromagnetic fraction for jet energy corrections or not (will only have an effect for jets CaloJets)
0096     bool emf_;
0097     /// input jet collection
0098     edm::EDGetTokenT<edm::View<reco::Jet> > srcToken_;
0099     /// type of flavor dependent JEC factors (only 'J' and 'T' are allowed)
0100     std::string type_;
0101     /// label of jec factors
0102     std::string label_;
0103     /// label of additional L1Offset corrector for JPT jets; for format reasons this string is
0104     /// kept in a vector of strings
0105     std::vector<std::string> extraJPTOffset_;
0106     /// label for L1Offset primaryVertex collection
0107     edm::InputTag primaryVertices_;
0108     edm::EDGetTokenT<std::vector<reco::Vertex> > primaryVerticesToken_;
0109     /// label for L1FastJet energy density parameter rho
0110     edm::InputTag rho_;
0111     edm::EDGetTokenT<double> rhoToken_;
0112     const edm::ESGetToken<JetCorrectorParametersCollection, JetCorrectionsRecord> parametersToken_;
0113     edm::ESWatcher<JetCorrectionsRecord> parametersWatcher_;
0114     /// use the NPV and rho with the JEC? (used for L1Offset/L1FastJet and L1FastJet, resp.)
0115     bool useNPV_;
0116     bool useRho_;
0117     /// jec levels for different flavors. In the default configuration
0118     /// this map would look like this:
0119     /// GLUON  : 'L2Relative', 'L3Absolute', 'L5FLavor_jg', L7Parton_jg'
0120     /// UDS    : 'L2Relative', 'L3Absolute', 'L5FLavor_jq', L7Parton_jq'
0121     /// CHARM  : 'L2Relative', 'L3Absolute', 'L5FLavor_jc', L7Parton_jc'
0122     /// BOTTOM : 'L2Relative', 'L3Absolute', 'L5FLavor_jb', L7Parton_jb'
0123     /// or just like this:
0124     /// NONE   : 'L2Relative', 'L3Absolute', 'L2L3Residual'
0125     /// per definition the vectors for all elements in this map should
0126     /// have the same size
0127     FlavorCorrLevelMap levels_;
0128     /// cache container for jet corrections
0129     std::map<JetCorrFactors::Flavor, std::unique_ptr<FactorizedJetCorrector> > correctors_;
0130     /// cache container for JPTOffset jet corrections
0131     std::unique_ptr<FactorizedJetCorrector> extraJPTOffsetCorrector_;
0132   };
0133
0134   inline int JetCorrFactorsProducer::numberOf(const edm::Handle<std::vector<reco::Vertex> >& primaryVertices) {
0135     int npv = 0;
0136     for (auto const& pv : *primaryVertices) {
0137       if (pv.ndof() >= 4)
0138         ++npv;
0139     }
0140     return npv;
0141   }
0142 }  // namespace pat
0143
0144 using namespace pat;
0145
0146 JetCorrFactorsProducer::JetCorrFactorsProducer(const edm::ParameterSet& cfg)
0147     : emf_(cfg.getParameter<bool>("emf")),
0148       srcToken_(consumes(cfg.getParameter<edm::InputTag>("src"))),
0149       type_(cfg.getParameter<std::string>("flavorType")),
0150       label_(cfg.getParameter<std::string>("@module_label")),
0151       parametersToken_{esConsumes(edm::ESInputTag("", cfg.getParameter<std::string>("payload")))},
0152       useNPV_(cfg.getParameter<bool>("useNPV")),
0153       useRho_(cfg.getParameter<bool>("useRho")) {
0154   std::vector<std::string> levels = cfg.getParameter<std::vector<std::string> >("levels");
0155   // fill the std::map for levels_, which might be flavor dependent or not;
0156   // flavor dependency is determined from the fact whether the std::string
0157   // L5Flavor or L7Parton can be found in levels; if flavor dependent four
0158   // vectors of strings will be filled into the map corresponding to GLUON,
0159   // UDS, CHARM and BOTTOM (according to JetCorrFactors::Flavor), 'L5Flavor'
0160   // and 'L7Parton' will be expanded accordingly; if not levels_ is filled
0161   // with only one vector of strings according to NONE. This vector will be
0162   // equivalent to the original vector of strings.
0163   if (std::find(levels.begin(), levels.end(), "L5Flavor") != levels.end() ||
0164       std::find(levels.begin(), levels.end(), "L7Parton") != levels.end()) {
0165     levels_[JetCorrFactors::GLUON] = expand(levels, JetCorrFactors::GLUON);
0166     levels_[JetCorrFactors::UDS] = expand(levels, JetCorrFactors::UDS);
0167     levels_[JetCorrFactors::CHARM] = expand(levels, JetCorrFactors::CHARM);
0168     levels_[JetCorrFactors::BOTTOM] = expand(levels, JetCorrFactors::BOTTOM);
0169   } else {
0170     levels_[JetCorrFactors::NONE] = levels;
0171   }
0172   // if the std::string L1JPTOffset can be found in levels an additional
0173   // parameter extraJPTOffset is needed, which should pass on the the usual
0174   // L1Offset correction, which is an additional input to the L1JPTOffset
0175   // corrector
0176   if (std::find(levels.begin(), levels.end(), "L1JPTOffset") != levels.end()) {
0177     extraJPTOffset_.push_back(cfg.getParameter<std::string>("extraJPTOffset"));
0178   }
0179   // if the std::string L1Offset can be found in levels an additional para-
0180   // meter primaryVertices is needed, which should pass on the offline pri-
0181   // mary vertex collection. The size of this collection is needed for the
0182   // L1Offset correction.
0183   if (useNPV_) {
0184     primaryVertices_ = cfg.getParameter<edm::InputTag>("primaryVertices");
0185     primaryVerticesToken_ = mayConsume<std::vector<reco::Vertex> >(primaryVertices_);
0186   }
0187   // if the std::string L1FastJet can be found in levels an additional
0188   // parameter rho is needed, which should pass on the energy density
0189   // parameter for the corresponding jet collection.
0190   if (useRho_) {
0191     if ((!extraJPTOffset_.empty() && extraJPTOffset_.front() == std::string("L1FastJet")) ||
0192         std::find(levels.begin(), levels.end(), "L1FastJet") != levels.end()) {
0193       rho_ = cfg.getParameter<edm::InputTag>("rho");
0194       rhoToken_ = mayConsume<double>(rho_);
0195     } else {
0196       edm::LogInfo message("Parameter rho not used");
0197       message << "Module is configured to use the parameter rho, but rho is only used     \n"
0198               << "for L1FastJet corrections. The configuration of levels does not contain \n"
0199               << "L1FastJet corrections though, so rho will not be used by this module.   \n";
0200     }
0201   }
0202   produces<JetCorrFactorsMap>();
0203 }
0204
0205 std::vector<std::string> JetCorrFactorsProducer::expand(const std::vector<std::string>& levels,
0206                                                         const JetCorrFactors::Flavor& flavor) {
0207   std::vector<std::string> expand;
0208   for (std::vector<std::string>::const_iterator level = levels.begin(); level != levels.end(); ++level) {
0209     if ((*level) == "L5Flavor" || (*level) == "L7Parton") {
0210       if (flavor == JetCorrFactors::GLUON) {
0211         if (*level == "L7Parton" && type_ == "T") {
0212           edm::LogWarning message("L7Parton::GLUON not available");
0213           message << "Jet energy corrections requested for level: L7Parton and type: 'T'. \n"
0214                   << "For this combination there is no GLUON correction available. The    \n"
0215                   << "correction for this flavor type will be taken from 'J'.";
0216         }
0217         expand.push_back(std::string(*level).append("_").append("g").append("J"));
0218       }
0219       if (flavor == JetCorrFactors::UDS)
0220         expand.push_back(std::string(*level).append("_").append("q").append(type_));
0221       if (flavor == JetCorrFactors::CHARM)
0222         expand.push_back(std::string(*level).append("_").append("c").append(type_));
0223       if (flavor == JetCorrFactors::BOTTOM)
0224         expand.push_back(std::string(*level).append("_").append("b").append(type_));
0225     } else {
0226       expand.push_back(*level);
0227     }
0228   }
0229   return expand;
0230 }
0231
0232 std::vector<JetCorrectorParameters> JetCorrFactorsProducer::params(const JetCorrectorParametersCollection& parameters,
0233                                                                    const std::vector<std::string>& levels) const {
0234   std::vector<JetCorrectorParameters> params;
0235   for (std::vector<std::string>::const_iterator level = levels.begin(); level != levels.end(); ++level) {
0236     const JetCorrectorParameters& ip = parameters[*level];  //ip.printScreen();
0237     params.push_back(ip);
0238   }
0239   return params;
0240 }
0241
0242 float JetCorrFactorsProducer::evaluate(edm::View<reco::Jet>::const_iterator& jet,
0243                                        const JetCorrFactors::Flavor& flavor,
0244                                        int level) {
0245   std::unique_ptr<FactorizedJetCorrector>& corrector = correctors_.find(flavor)->second;
0246   // add parameters for JPT corrections
0247   const reco::JPTJet* jpt = dynamic_cast<reco::JPTJet const*>(&*jet);
0248   if (jpt) {
0249     TLorentzVector p4;
0250     p4.SetPtEtaPhiE(jpt->getCaloJetRef()->pt(),
0251                     jpt->getCaloJetRef()->eta(),
0252                     jpt->getCaloJetRef()->phi(),
0253                     jpt->getCaloJetRef()->energy());
0254     if (extraJPTOffsetCorrector_) {
0255       extraJPTOffsetCorrector_->setJPTrawP4(p4);
0256       corrector->setJPTrawOff(extraJPTOffsetCorrector_->getSubCorrections()[0]);
0257     }
0258     corrector->setJPTrawP4(p4);
0259   }
0260   //For PAT jets undo previous jet energy corrections
0261   const Jet* patjet = dynamic_cast<Jet const*>(&*jet);
0262   if (patjet) {
0263     corrector->setJetEta(patjet->correctedP4(0).eta());
0264     corrector->setJetPt(patjet->correctedP4(0).pt());
0265     corrector->setJetPhi(patjet->correctedP4(0).phi());
0266     corrector->setJetE(patjet->correctedP4(0).energy());
0267   } else {
0268     corrector->setJetEta(jet->eta());
0269     corrector->setJetPt(jet->pt());
0270     corrector->setJetPhi(jet->phi());
0271     corrector->setJetE(jet->energy());
0272   }
0273   if (emf_ && dynamic_cast<const reco::CaloJet*>(&*jet)) {
0274     corrector->setJetEMF(dynamic_cast<const reco::CaloJet*>(&*jet)->emEnergyFraction());
0275   }
0276   return corrector->getSubCorrections()[level];
0277 }
0278
0279 void JetCorrFactorsProducer::produce(edm::Event& event, const edm::EventSetup& setup) {
0280   // get jet collection from the event
0281   edm::Handle<edm::View<reco::Jet> > jets;
0282   event.getByToken(srcToken_, jets);
0283
0284   // get primary vertices for L1Offset correction level if needed
0285   edm::Handle<std::vector<reco::Vertex> > primaryVertices;
0286   if (!primaryVertices_.label().empty())
0287     event.getByToken(primaryVerticesToken_, primaryVertices);
0288
0289   // get parameter rho for L1FastJet correction level if needed
0290   edm::Handle<double> rho;
0291   if (!rho_.label().empty())
0292     event.getByToken(rhoToken_, rho);
0293
0294   if (parametersWatcher_.check(setup)) {
0295     // retreive parameters from the DB
0296     auto const& parameters = setup.getData(parametersToken_);
0297     // initialize jet correctors
0298     for (auto const& flavor : levels_) {
0299       correctors_[flavor.first] = std::make_unique<FactorizedJetCorrector>(params(parameters, flavor.second));
0300     }
0301     // initialize extra jet corrector for jpt if needed
0302     if (!extraJPTOffset_.empty()) {
0303       extraJPTOffsetCorrector_ = std::make_unique<FactorizedJetCorrector>(params(parameters, extraJPTOffset_));
0304     }
0305   }
0306
0307   // fill the jetCorrFactors
0308   std::vector<JetCorrFactors> jcfs;
0309   for (edm::View<reco::Jet>::const_iterator jet = jets->begin(); jet != jets->end(); ++jet) {
0310     // the JetCorrFactors::CorrectionFactor is a std::pair<std::string, std::vector<float> >
0311     // the string corresponds to the label of the correction level, the vector contains four
0312     // floats if flavor dependent and one float else. Per construction jet energy corrections
0313     // will be flavor independent up to the first flavor dependent correction and flavor de-
0314     // pendent afterwards. The first correction level is predefined with label 'Uncorrected'.
0315     // Per definition it is flavor independent. The correction factor is 1.
0316     std::vector<JetCorrFactors::CorrectionFactor> jec;
0317     jec.emplace_back("Uncorrected", std::vector<float>(1, 1));
0318
0319     // pick the first element in the map (which could be the only one) and loop all jec
0320     // levels listed for that element. If this is not the only element all jec levels, which
0321     // are flavor independent will give the same correction factors until the first flavor
0322     // dependent correction level is reached. So the first element is still a good choice.
0323     FlavorCorrLevelMap::const_iterator corrLevel = levels_.begin();
0324     if (corrLevel == levels_.end()) {
0325       throw cms::Exception("No JECFactors")
0326           << "You request to create a jetCorrFactors object with no JEC Levels indicated. \n"
0327           << "This makes no sense, either you should correct this or drop the module from \n"
0328           << "the sequence.";
0329     }
0330     for (unsigned int idx = 0; idx < corrLevel->second.size(); ++idx) {
0331       std::vector<float> factors;
0332       if (corrLevel->second[idx].find("L5Flavor") != std::string::npos ||
0333           corrLevel->second[idx].find("L7Parton") != std::string::npos) {
0334         for (FlavorCorrLevelMap::const_iterator flavor = corrLevel; flavor != levels_.end(); ++flavor) {
0335           if (!primaryVertices_.label().empty()) {
0336             // if primaryVerticesToken_ has a value the number of primary vertices needs to be
0337             // specified
0338             correctors_.find(flavor->first)->second->setNPV(numberOf(primaryVertices));
0339           }
0340           if (!rho_.label().empty()) {
0341             // if rhoToken_ has a value the energy density parameter rho and the jet area need
0342             //  to be specified
0343             correctors_.find(flavor->first)->second->setRho(*rho);
0344             correctors_.find(flavor->first)->second->setJetA(jet->jetArea());
0345           }
0346           factors.push_back(evaluate(jet, flavor->first, idx));
0347         }
0348       } else {
0349         if (!primaryVertices_.label().empty()) {
0350           // if primaryVerticesToken_ has a value the number of primary vertices needs to be
0351           // specified
0352           correctors_.find(corrLevel->first)->second->setNPV(numberOf(primaryVertices));
0353         }
0354         if (!rho_.label().empty()) {
0355           // if rhoToken_ has a value the energy density parameter rho and the jet area need
0356           // to be specified
0357           correctors_.find(corrLevel->first)->second->setRho(*rho);
0358           correctors_.find(corrLevel->first)->second->setJetA(jet->jetArea());
0359         }
0360         factors.push_back(evaluate(jet, corrLevel->first, idx));
0361       }
0362       // push back the set of JetCorrFactors: the first entry corresponds to the label
0363       // of the correction level, which is taken from the first element in levels_. For
0364       // L5Flavor and L7Parton the part including the first '_' indicating the flavor
0365       // of the first element in levels_ is chopped of from the label to avoid confusion
0366       // of the correction levels. The second parameter corresponds to the set of jec
0367       // factors, which might be flavor dependent or not. In the default configuration
0368       // the CorrectionFactor will look like this: 'Uncorrected': 1 ; 'L2Relative': x ;
0369       // 'L3Absolute': x ; 'L5Flavor': v, x, y, z ; 'L7Parton': v, x, y, z
0370       jec.emplace_back((corrLevel->second[idx]).substr(0, (corrLevel->second[idx]).find('_')), factors);
0371     }
0372     // create the actual object with the scale factors we want the valuemap to refer to
0373     // label_ corresponds to the label of the module instance
0374     jcfs.emplace_back(label_, jec);
0375   }
0376   // build the value map
0377   auto jetCorrsMap = std::make_unique<JetCorrFactorsMap>();
0378   JetCorrFactorsMap::Filler filler(*jetCorrsMap);
0379   // jets and jetCorrs have their indices aligned by construction
0380   filler.insert(jets, jcfs.begin(), jcfs.end());
0381   filler.fill();  // do the actual filling
0382   // put our produced stuff in the event
0383   event.put(std::move(jetCorrsMap));
0384 }
0385
0386 void JetCorrFactorsProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0387   edm::ParameterSetDescription iDesc;
0388   iDesc.add<bool>("emf", false);
0389   iDesc.add<std::string>("flavorType", "J");
0390   iDesc.add<edm::InputTag>("src", edm::InputTag("ak5CaloJets"));
0391   iDesc.add<std::string>("payload", "AK5Calo");
0392   iDesc.add<bool>("useNPV", true);
0393   iDesc.add<edm::InputTag>("primaryVertices", edm::InputTag("offlinePrimaryVertices"));
0394   iDesc.add<bool>("useRho", true);
0395   iDesc.add<edm::InputTag>("rho", edm::InputTag("fixedGridRhoFastjetAllCalo"));
0396   iDesc.add<std::string>("extraJPTOffset", "L1Offset");
0397
0398   iDesc.add<std::vector<std::string> >("levels",
0399                                        {
0400                                            "L1Offset",
0401                                            "L2Relative",
0402                                            "L3Absolute",
0403                                            "L2L3Residual",
0404                                            "L5Flavor",
0405                                            "L7Parton",
0406                                        });
0407   descriptions.add("JetCorrFactorsProducer", iDesc);
0408 }
0409
0410 #include "FWCore/Framework/interface/MakerMacros.h"
0411 DEFINE_FWK_MODULE(JetCorrFactorsProducer);