JetCorrFactorsProducer.cc

Go to the documentation of this file.

#include <vector>
#include <string>
#include <iostream>

#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "DataFormats/JetReco/interface/JPTJet.h"
#include "DataFormats/JetReco/interface/CaloJet.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/PatCandidates/interface/Jet.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "PhysicsTools/PatAlgos/plugins/JetCorrFactorsProducer.h"
#include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"

#include "FWCore/ParameterSet/interface/ParameterDescription.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"

using namespace pat;

JetCorrFactorsProducer::JetCorrFactorsProducer(const edm::ParameterSet& cfg)
    : emf_(cfg.getParameter<bool>("emf")),
      srcToken_(consumes<edm::View<reco::Jet> >(cfg.getParameter<edm::InputTag>("src"))),
      type_(cfg.getParameter<std::string>("flavorType")),
      label_(cfg.getParameter<std::string>("@module_label")),
      payload_(cfg.getParameter<std::string>("payload")),
      useNPV_(cfg.getParameter<bool>("useNPV")),
      useRho_(cfg.getParameter<bool>("useRho")),
      cacheId_(0) {
  std::vector<std::string> levels = cfg.getParameter<std::vector<std::string> >("levels");
  // fill the std::map for levels_, which might be flavor dependent or not;
  // flavor dependency is determined from the fact whether the std::string
  // L5Flavor or L7Parton can be found in levels; if flavor dependent four
  // vectors of strings will be filled into the map corresponding to GLUON,
  // UDS, CHARM and BOTTOM (according to JetCorrFactors::Flavor), 'L5Flavor'
  // and 'L7Parton' will be expanded accordingly; if not levels_ is filled
  // with only one vector of strings according to NONE. This vector will be
  // equivalent to the original vector of strings.
  if (std::find(levels.begin(), levels.end(), "L5Flavor") != levels.end() ||
      std::find(levels.begin(), levels.end(), "L7Parton") != levels.end()) {
    levels_[JetCorrFactors::GLUON] = expand(levels, JetCorrFactors::GLUON);
    levels_[JetCorrFactors::UDS] = expand(levels, JetCorrFactors::UDS);
    levels_[JetCorrFactors::CHARM] = expand(levels, JetCorrFactors::CHARM);
    levels_[JetCorrFactors::BOTTOM] = expand(levels, JetCorrFactors::BOTTOM);
  } else {
    levels_[JetCorrFactors::NONE] = levels;
  }
  // if the std::string L1JPTOffset can be found in levels an additional
  // parameter extraJPTOffset is needed, which should pass on the the usual
  // L1Offset correction, which is an additional input to the L1JPTOffset
  // corrector
  if (std::find(levels.begin(), levels.end(), "L1JPTOffset") != levels.end()) {
    if (cfg.existsAs<std::string>("extraJPTOffset")) {
      extraJPTOffset_.push_back(cfg.getParameter<std::string>("extraJPTOffset"));
    } else {
      throw cms::Exception("No parameter extraJPTOffset specified")
          << "The configured correction levels contain a L1JPTOffset correction, which re- \n"
          << "quires the additional parameter extraJPTOffset or type std::string. This     \n"
          << "string should correspond to the L1Offset corrections that should be applied  \n"
          << "together with the JPTL1Offset corrections. These corrections can be of type  \n"
          << "L1Offset or L1FastJet.                                                       \n";
    }
  }
  // if the std::string L1Offset can be found in levels an additional para-
  // meter primaryVertices is needed, which should pass on the offline pri-
  // mary vertex collection. The size of this collection is needed for the
  // L1Offset correction.
  if (useNPV_) {
    if (cfg.existsAs<edm::InputTag>("primaryVertices")) {
      primaryVertices_ = cfg.getParameter<edm::InputTag>("primaryVertices");
      primaryVerticesToken_ = mayConsume<std::vector<reco::Vertex> >(primaryVertices_);
    } else {
      throw cms::Exception("No primaryVertices specified")
          << "The configured correction levels contain an L1Offset or L1FastJet correction, \n"
          << "which requires the number of offlinePrimaryVertices. Please specify this col- \n"
          << "lection as additional optional parameter primaryVertices of type edm::InputTag\n"
          << "in the jetCorrFactors module.                                                 \n";
    }
  }
  // if the std::string L1FastJet can be found in levels an additional
  // parameter rho is needed, which should pass on the energy density
  // parameter for the corresponding jet collection.
  if (useRho_) {
    if ((!extraJPTOffset_.empty() && extraJPTOffset_.front() == std::string("L1FastJet")) ||
        std::find(levels.begin(), levels.end(), "L1FastJet") != levels.end()) {
      if (cfg.existsAs<edm::InputTag>("rho")) {
        rho_ = cfg.getParameter<edm::InputTag>("rho");
        rhoToken_ = mayConsume<double>(rho_);
      } else {
        throw cms::Exception("No parameter rho specified")
            << "The configured correction levels contain a L1FastJet correction, which re- \n"
            << "quires the energy density parameter rho. Please specify this collection as \n"
            << "additional optional parameter rho of type edm::InputTag in the jetCorrFac- \n"
            << "tors module.                                                               \n";
      }
    } else {
      edm::LogInfo message("Parameter rho not used");
      message << "Module is configured to use the parameter rho, but rho is only used     \n"
              << "for L1FastJet corrections. The configuration of levels does not contain \n"
              << "L1FastJet corrections though, so rho will not be used by this module.   \n";
    }
  }
  produces<JetCorrFactorsMap>();
}

std::vector<std::string> JetCorrFactorsProducer::expand(const std::vector<std::string>& levels,
                                                        const JetCorrFactors::Flavor& flavor) {
  std::vector<std::string> expand;
  for (std::vector<std::string>::const_iterator level = levels.begin(); level != levels.end(); ++level) {
    if ((*level) == "L5Flavor" || (*level) == "L7Parton") {
      if (flavor == JetCorrFactors::GLUON) {
        if (*level == "L7Parton" && type_ == "T") {
          edm::LogWarning message("L7Parton::GLUON not available");
          message << "Jet energy corrections requested for level: L7Parton and type: 'T'. \n"
                  << "For this combination there is no GLUON correction available. The    \n"
                  << "correction for this flavor type will be taken from 'J'.";
        }
        expand.push_back(std::string(*level).append("_").append("g").append("J"));
      }
      if (flavor == JetCorrFactors::UDS)
        expand.push_back(std::string(*level).append("_").append("q").append(type_));
      if (flavor == JetCorrFactors::CHARM)
        expand.push_back(std::string(*level).append("_").append("c").append(type_));
      if (flavor == JetCorrFactors::BOTTOM)
        expand.push_back(std::string(*level).append("_").append("b").append(type_));
    } else {
      expand.push_back(*level);
    }
  }
  return expand;
}

std::vector<JetCorrectorParameters> JetCorrFactorsProducer::params(const JetCorrectorParametersCollection& parameters,
                                                                   const std::vector<std::string>& levels) const {
  std::vector<JetCorrectorParameters> params;
  for (std::vector<std::string>::const_iterator level = levels.begin(); level != levels.end(); ++level) {
    const JetCorrectorParameters& ip = parameters[*level];  //ip.printScreen();
    params.push_back(ip);
  }
  return params;
}

float JetCorrFactorsProducer::evaluate(edm::View<reco::Jet>::const_iterator& jet,
                                       const JetCorrFactors::Flavor& flavor,
                                       int level) {
  std::unique_ptr<FactorizedJetCorrector>& corrector = correctors_.find(flavor)->second;
  // add parameters for JPT corrections
  const reco::JPTJet* jpt = dynamic_cast<reco::JPTJet const*>(&*jet);
  if (jpt) {
    TLorentzVector p4;
    p4.SetPtEtaPhiE(jpt->getCaloJetRef()->pt(),
                    jpt->getCaloJetRef()->eta(),
                    jpt->getCaloJetRef()->phi(),
                    jpt->getCaloJetRef()->energy());
    if (extraJPTOffsetCorrector_) {
      extraJPTOffsetCorrector_->setJPTrawP4(p4);
      corrector->setJPTrawOff(extraJPTOffsetCorrector_->getSubCorrections()[0]);
    }
    corrector->setJPTrawP4(p4);
  }
  //For PAT jets undo previous jet energy corrections
  const Jet* patjet = dynamic_cast<Jet const*>(&*jet);
  if (patjet) {
    corrector->setJetEta(patjet->correctedP4(0).eta());
    corrector->setJetPt(patjet->correctedP4(0).pt());
    corrector->setJetPhi(patjet->correctedP4(0).phi());
    corrector->setJetE(patjet->correctedP4(0).energy());
  } else {
    corrector->setJetEta(jet->eta());
    corrector->setJetPt(jet->pt());
    corrector->setJetPhi(jet->phi());
    corrector->setJetE(jet->energy());
  }
  if (emf_ && dynamic_cast<const reco::CaloJet*>(&*jet)) {
    corrector->setJetEMF(dynamic_cast<const reco::CaloJet*>(&*jet)->emEnergyFraction());
  }
  return corrector->getSubCorrections()[level];
}

std::string JetCorrFactorsProducer::payload() { return payload_; }

void JetCorrFactorsProducer::produce(edm::Event& event, const edm::EventSetup& setup) {
  // get jet collection from the event
  edm::Handle<edm::View<reco::Jet> > jets;
  event.getByToken(srcToken_, jets);

  // get primary vertices for L1Offset correction level if needed
  edm::Handle<std::vector<reco::Vertex> > primaryVertices;
  if (!primaryVertices_.label().empty())
    event.getByToken(primaryVerticesToken_, primaryVertices);

  // get parameter rho for L1FastJet correction level if needed
  edm::Handle<double> rho;
  if (!rho_.label().empty())
    event.getByToken(rhoToken_, rho);

  auto const& rec = setup.get<JetCorrectionsRecord>();
  if (cacheId_ != rec.cacheIdentifier()) {
    // retreive parameters from the DB
    edm::ESHandle<JetCorrectorParametersCollection> parameters;
    setup.get<JetCorrectionsRecord>().get(payload(), parameters);
    // initialize jet correctors
    for (FlavorCorrLevelMap::const_iterator flavor = levels_.begin(); flavor != levels_.end(); ++flavor) {
      correctors_[flavor->first].reset(new FactorizedJetCorrector(params(*parameters, flavor->second)));
    }
    // initialize extra jet corrector for jpt if needed
    if (!extraJPTOffset_.empty()) {
      extraJPTOffsetCorrector_.reset(new FactorizedJetCorrector(params(*parameters, extraJPTOffset_)));
    }
    cacheId_ = rec.cacheIdentifier();
  }

  // fill the jetCorrFactors
  std::vector<JetCorrFactors> jcfs;
  for (edm::View<reco::Jet>::const_iterator jet = jets->begin(); jet != jets->end(); ++jet) {
    // the JetCorrFactors::CorrectionFactor is a std::pair<std::string, std::vector<float> >
    // the string corresponds to the label of the correction level, the vector contains four
    // floats if flavor dependent and one float else. Per construction jet energy corrections
    // will be flavor independent up to the first flavor dependent correction and flavor de-
    // pendent afterwards. The first correction level is predefined with label 'Uncorrected'.
    // Per definition it is flavor independent. The correction factor is 1.
    std::vector<JetCorrFactors::CorrectionFactor> jec;
    jec.push_back(
        std::make_pair<std::string, std::vector<float> >(std::string("Uncorrected"), std::vector<float>(1, 1)));

    // pick the first element in the map (which could be the only one) and loop all jec
    // levels listed for that element. If this is not the only element all jec levels, which
    // are flavor independent will give the same correction factors until the first flavor
    // dependent correction level is reached. So the first element is still a good choice.
    FlavorCorrLevelMap::const_iterator corrLevel = levels_.begin();
    if (corrLevel == levels_.end()) {
      throw cms::Exception("No JECFactors")
          << "You request to create a jetCorrFactors object with no JEC Levels indicated. \n"
          << "This makes no sense, either you should correct this or drop the module from \n"
          << "the sequence.";
    }
    for (unsigned int idx = 0; idx < corrLevel->second.size(); ++idx) {
      std::vector<float> factors;
      if (corrLevel->second[idx].find("L5Flavor") != std::string::npos ||
          corrLevel->second[idx].find("L7Parton") != std::string::npos) {
        for (FlavorCorrLevelMap::const_iterator flavor = corrLevel; flavor != levels_.end(); ++flavor) {
          if (!primaryVertices_.label().empty()) {
            // if primaryVerticesToken_ has a value the number of primary vertices needs to be
            // specified
            correctors_.find(flavor->first)->second->setNPV(numberOf(primaryVertices));
          }
          if (!rho_.label().empty()) {
            // if rhoToken_ has a value the energy density parameter rho and the jet area need
            //  to be specified
            correctors_.find(flavor->first)->second->setRho(*rho);
            correctors_.find(flavor->first)->second->setJetA(jet->jetArea());
          }
          factors.push_back(evaluate(jet, flavor->first, idx));
        }
      } else {
        if (!primaryVertices_.label().empty()) {
          // if primaryVerticesToken_ has a value the number of primary vertices needs to be
          // specified
          correctors_.find(corrLevel->first)->second->setNPV(numberOf(primaryVertices));
        }
        if (!rho_.label().empty()) {
          // if rhoToken_ has a value the energy density parameter rho and the jet area need
          // to be specified
          correctors_.find(corrLevel->first)->second->setRho(*rho);
          correctors_.find(corrLevel->first)->second->setJetA(jet->jetArea());
        }
        factors.push_back(evaluate(jet, corrLevel->first, idx));
      }
      // push back the set of JetCorrFactors: the first entry corresponds to the label
      // of the correction level, which is taken from the first element in levels_. For
      // L5Flavor and L7Parton the part including the first '_' indicating the flavor
      // of the first element in levels_ is chopped of from the label to avoid confusion
      // of the correction levels. The second parameter corresponds to the set of jec
      // factors, which might be flavor dependent or not. In the default configuration
      // the CorrectionFactor will look like this: 'Uncorrected': 1 ; 'L2Relative': x ;
      // 'L3Absolute': x ; 'L5Flavor': v, x, y, z ; 'L7Parton': v, x, y, z
      jec.push_back(std::make_pair((corrLevel->second[idx]).substr(0, (corrLevel->second[idx]).find("_")), factors));
    }
    // create the actual object with the scale factors we want the valuemap to refer to
    // label_ corresponds to the label of the module instance
    JetCorrFactors corrFactors(label_, jec);
    jcfs.push_back(corrFactors);
  }
  // build the value map
  auto jetCorrsMap = std::make_unique<JetCorrFactorsMap>();
  JetCorrFactorsMap::Filler filler(*jetCorrsMap);
  // jets and jetCorrs have their indices aligned by construction
  filler.insert(jets, jcfs.begin(), jcfs.end());
  filler.fill();  // do the actual filling
  // put our produced stuff in the event
  event.put(std::move(jetCorrsMap));
}

void JetCorrFactorsProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  edm::ParameterSetDescription iDesc;
  iDesc.add<bool>("emf", false);
  iDesc.add<std::string>("flavorType", "J");
  iDesc.add<edm::InputTag>("src", edm::InputTag("ak5CaloJets"));
  iDesc.add<std::string>("payload", "AK5Calo");
  iDesc.add<bool>("useNPV", true);
  iDesc.add<edm::InputTag>("primaryVertices", edm::InputTag("offlinePrimaryVertices"));
  iDesc.add<bool>("useRho", true);
  iDesc.add<edm::InputTag>("rho", edm::InputTag("fixedGridRhoFastjetAllCalo"));
  iDesc.add<std::string>("extraJPTOffset", "L1Offset");

  std::vector<std::string> levels;
  levels.push_back(std::string("L1Offset"));
  levels.push_back(std::string("L2Relative"));
  levels.push_back(std::string("L3Absolute"));
  levels.push_back(std::string("L2L3Residual"));
  levels.push_back(std::string("L5Flavor"));
  levels.push_back(std::string("L7Parton"));
  iDesc.add<std::vector<std::string> >("levels", levels);
  descriptions.add("JetCorrFactorsProducer", iDesc);
}

#include "FWCore/Framework/interface/MakerMacros.h"

DEFINE_FWK_MODULE(JetCorrFactorsProducer);