/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/PhysicsTools/PatAlgos/plugins/JetCorrFactorsProducer.cc

Go to the documentation of this file.

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

#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "DataFormats/JetReco/interface/CaloJet.h"
#include "DataFormats/VertexReco/interface/Vertex.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" )), 
  src_(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"))
{

  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 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");
    }
    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 in the jetCorrFactors\n"
        << "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(std::find(levels.begin(), levels.end(), "L1FastJet")!=levels.end()){
      if(cfg.existsAs<edm::InputTag>("rho")){
        rho_=cfg.getParameter<edm::InputTag>("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 in the jetCorrFactors module. \n";
      }
    }
    else{
      edm::LogWarning message( "Parameter rho not used" );
      message << "Module is configured to use the parameter rho, but but rho is only used \n"
              << "for L1FastJet corrections at the moment. The configuration of levels    \n"
              << "does not contain L1FastJet corrections though, so rho will not be used  \n"
              << "throughout 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, boost::shared_ptr<FactorizedJetCorrector>& corrector, int level)
{
  corrector->setJetEta(jet->eta()); corrector->setJetPt(jet->pt()); 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.getByLabel(src_, jets);

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

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

  // retreive parameters from the DB 
  edm::ESHandle<JetCorrectorParametersCollection> parameters;
  setup.get<JetCorrectionsRecord>().get(payload(), parameters); 

  // initialize jet correctors
  std::map<JetCorrFactors::Flavor, boost::shared_ptr<FactorizedJetCorrector> > corrector;
  for(FlavorCorrLevelMap::const_iterator flavor=levels_.begin(); flavor!=levels_.end(); ++flavor){
    corrector[flavor->first] = boost::shared_ptr<FactorizedJetCorrector>( new FactorizedJetCorrector(params(*parameters, flavor->second)) );
  }

  // 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){
      bool flavorDependent=false;
      std::vector<float> factors;
      if(flavorDependent || 
         corrLevel->second[idx].find("L5Flavor")!=std::string::npos || 
         corrLevel->second[idx].find("L7Parton")!=std::string::npos){
        flavorDependent=true;
        // after the first encounter all subsequent correction levels are flavor dependent
        for(FlavorCorrLevelMap::const_iterator flavor=corrLevel; flavor!=levels_.end(); ++flavor){
          if(!primaryVertices_.label().empty()){
            // if primaryVertices_ has a value the number of primary vertices needs to be 
            // specified
            corrector.find(flavor->first)->second->setNPV(numberOf(primaryVertices));
          }
          if(!rho_.label().empty()){
            // if rho_ has a value the energy density parameter rho and the jet area need
            //  to be specified
            corrector.find(flavor->first)->second->setRho(*rho);
            corrector.find(flavor->first)->second->setJetA(jet->jetArea());
          }
          factors.push_back(evaluate(jet, corrector.find(flavor->first)->second, idx)); 
        }
      }
      else{
        if(!primaryVertices_.label().empty()){
          // if primaryVertices_ has a value the number of primary vertices needs to be 
          // specified
          corrector.find(corrLevel->first)->second->setNPV(numberOf(primaryVertices));
        }
        if(!rho_.label().empty()){
          // if rho_ has a value the energy density parameter rho and the jet area need
          //  to be specified
          corrector.find(corrLevel->first)->second->setRho(*rho);
          corrector.find(corrLevel->first)->second->setJetA(jet->jetArea());
        }
        factors.push_back(evaluate(jet, corrector.find(corrLevel->first)->second, 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
  std::auto_ptr<JetCorrFactorsMap> jetCorrsMap(new 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(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", false);
  iDesc.add<edm::InputTag>("rho", edm::InputTag("kt6PFJets", "rho"));

  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("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);