AnalysisTasksAnalyzerJEC.cc

Go to the documentation of this file.

#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/JetReco/interface/Jet.h"
#include "DataFormats/PatCandidates/interface/Jet.h"
#include "PhysicsTools/PatExamples/interface/AnalysisTasksAnalyzerJEC.h"
#include "TH2.h"
#include "TProfile.h"
AnalysisTasksAnalyzerJEC::AnalysisTasksAnalyzerJEC(const edm::ParameterSet& cfg, TFileDirectory& fs)
    : edm::BasicAnalyzer::BasicAnalyzer(cfg, fs),
      Jets_(cfg.getParameter<edm::InputTag>("Jets")),
      jecLevel_(cfg.getParameter<std::string>("jecLevel")),
      patJetCorrFactors_(cfg.getParameter<std::string>("patJetCorrFactors")),
      help_(cfg.getParameter<bool>("help")) {
  hists_["Response"] = fs.make<TH2F>("Response", "response; #eta;p_{T}(reco)/p_{T}(gen)", 5, -3.3, 3.3, 100, 0.4, 1.6);
  jetInEvents_ = 0;
}
AnalysisTasksAnalyzerJEC::AnalysisTasksAnalyzerJEC(const edm::ParameterSet& cfg,
                                                   TFileDirectory& fs,
                                                   edm::ConsumesCollector&& iC)
    : edm::BasicAnalyzer::BasicAnalyzer(cfg, fs),
      Jets_(cfg.getParameter<edm::InputTag>("Jets")),
      JetsToken_(iC.consumes<std::vector<pat::Jet> >(Jets_)),
      jecLevel_(cfg.getParameter<std::string>("jecLevel")),
      patJetCorrFactors_(cfg.getParameter<std::string>("patJetCorrFactors")),
      help_(cfg.getParameter<bool>("help")) {
  hists_["Response"] = fs.make<TH2F>("Response", "response; #eta;p_{T}(reco)/p_{T}(gen)", 5, -3.3, 3.3, 100, 0.4, 1.6);
  jetInEvents_ = 0;
}
AnalysisTasksAnalyzerJEC::~AnalysisTasksAnalyzerJEC() {
  TFile* file = new TFile("myResponse" + TString(jecLevel_) + ".root", "RECREATE");
  TProfile* prof = hists_["Response"]->ProfileX();
  prof->Write();
  file->Write();
  file->Close();
}
void AnalysisTasksAnalyzerJEC::analyze(const edm::EventBase& event) {
  // define what Jet you are using; this is necessary as FWLite is not
  // capable of reading edm::Views
  using pat::Jet;
 
  // Handle to the Jet collection
  edm::Handle<std::vector<Jet> > Jets;
  event.getByLabel(Jets_, Jets);
 
  // loop Jet collection and fill histograms
  for (std::vector<Jet>::const_iterator jet_it = Jets->begin(); jet_it != Jets->end(); ++jet_it) {
    if (help_ == true) {
      if (jetInEvents_ == 0) {
        std::cout << "\n \n number of available JEC sets: " << jet_it->availableJECSets().size() << std::endl;
        for (unsigned int k = 0; k < jet_it->availableJECSets().size(); ++k) {
          std::cout << "\n \n available JEC Set: " << jet_it->availableJECSets()[k] << std::endl;
        }
        std::cout << "\n \n*** You found out which JEC Sets exist! Now correct it in your config file." << std::endl;
        std::cout << "\n \n number of available JEC levels " << jet_it->availableJECLevels().size() << std::endl;
        for (unsigned int k = 0; k < jet_it->availableJECLevels().size(); ++k) {
          std::cout << "\n \n available JEC: " << jet_it->availableJECLevels(patJetCorrFactors_)[k] << std::endl;
        }
        std::cout << "\n \n**** You did it correctly congratulations!!!!  And you found out which JEC levels are saved "
                     "within the jets. Correct this in your configuration file."
                  << std::endl;
      }
    }
 
    if (jet_it->genParticlesSize() > 0) {
      hists_["Response"]->Fill(jet_it->correctedJet(jecLevel_).eta(),
                               jet_it->correctedJet(jecLevel_).pt() / jet_it->genParticle(0)->pt());
      std::cout << "\n \n**** You did it correctly congratulations!!!! " << std::endl;
    }
    jetInEvents_ += 1;
  }
}