PatExamples/src/AnalysisTasksAnalyzerJEC.cc

0001 #include "DataFormats/Common/interface/Handle.h"
0002 #include "DataFormats/JetReco/interface/Jet.h"
0003 #include "DataFormats/PatCandidates/interface/Jet.h"
0004 #include "PhysicsTools/PatExamples/interface/AnalysisTasksAnalyzerJEC.h"
0005 #include "TH2.h"
0006 #include "TProfile.h"
0007 /// default constructor
0008 AnalysisTasksAnalyzerJEC::AnalysisTasksAnalyzerJEC(const edm::ParameterSet& cfg, TFileDirectory& fs)
0009     : edm::BasicAnalyzer::BasicAnalyzer(cfg, fs),
0010       Jets_(cfg.getParameter<edm::InputTag>("Jets")),
0011       jecLevel_(cfg.getParameter<std::string>("jecLevel")),
0012       patJetCorrFactors_(cfg.getParameter<std::string>("patJetCorrFactors")),
0013       help_(cfg.getParameter<bool>("help")) {
0014   hists_["Response"] = fs.make<TH2F>("Response", "response; #eta;p_{T}(reco)/p_{T}(gen)", 5, -3.3, 3.3, 100, 0.4, 1.6);
0015   jetInEvents_ = 0;
0016 }
0017 AnalysisTasksAnalyzerJEC::AnalysisTasksAnalyzerJEC(const edm::ParameterSet& cfg,
0018                                                    TFileDirectory& fs,
0019                                                    edm::ConsumesCollector&& iC)
0020     : edm::BasicAnalyzer::BasicAnalyzer(cfg, fs),
0021       Jets_(cfg.getParameter<edm::InputTag>("Jets")),
0022       JetsToken_(iC.consumes<std::vector<pat::Jet> >(Jets_)),
0023       jecLevel_(cfg.getParameter<std::string>("jecLevel")),
0024       patJetCorrFactors_(cfg.getParameter<std::string>("patJetCorrFactors")),
0025       help_(cfg.getParameter<bool>("help")) {
0026   hists_["Response"] = fs.make<TH2F>("Response", "response; #eta;p_{T}(reco)/p_{T}(gen)", 5, -3.3, 3.3, 100, 0.4, 1.6);
0027   jetInEvents_ = 0;
0028 }
0029 /// deconstructor
0030 AnalysisTasksAnalyzerJEC::~AnalysisTasksAnalyzerJEC() {
0031   /// A Response function is nicer in a TProfile:
0032   TFile* file = new TFile("myResponse" + TString(jecLevel_) + ".root", "RECREATE");
0033   TProfile* prof = hists_["Response"]->ProfileX();
0034   prof->Write();
0035   file->Write();
0036   file->Close();
0037 }
0038 /// everything that needs to be done during the event loop
0039 void AnalysisTasksAnalyzerJEC::analyze(const edm::EventBase& event) {
0040   // define what Jet you are using; this is necessary as FWLite is not
0041   // capable of reading edm::Views
0042   using pat::Jet;
0043
0044   // Handle to the Jet collection
0045   edm::Handle<std::vector<Jet> > Jets;
0046   event.getByLabel(Jets_, Jets);
0047
0048   // loop Jet collection and fill histograms
0049   for (std::vector<Jet>::const_iterator jet_it = Jets->begin(); jet_it != Jets->end(); ++jet_it) {
0050     ///// You can get some help if you switch help to True in your config file
0051     if (help_ == true) {
0052       if (jetInEvents_ == 0) {
0053         std::cout << "\n \n number of available JEC sets: " << jet_it->availableJECSets().size() << std::endl;
0054         for (unsigned int k = 0; k < jet_it->availableJECSets().size(); ++k) {
0055           std::cout << "\n \n available JEC Set: " << jet_it->availableJECSets()[k] << std::endl;
0056         }
0057         std::cout << "\n \n*** You found out which JEC Sets exist! Now correct it in your config file." << std::endl;
0058         std::cout << "\n \n number of available JEC levels " << jet_it->availableJECLevels().size() << std::endl;
0059         for (unsigned int k = 0; k < jet_it->availableJECLevels().size(); ++k) {
0060           std::cout << "\n \n available JEC: " << jet_it->availableJECLevels(patJetCorrFactors_)[k] << std::endl;
0061         }
0062         std::cout << "\n \n**** You did it correctly congratulations!!!!  And you found out which JEC levels are saved "
0063                      "within the jets. Correct this in your configuration file."
0064                   << std::endl;
0065       }
0066     }
0067
0068     if (jet_it->genParticlesSize() > 0) {
0069       hists_["Response"]->Fill(jet_it->correctedJet(jecLevel_).eta(),
0070                                jet_it->correctedJet(jecLevel_).pt() / jet_it->genParticle(0)->pt());
0071       std::cout << "\n \n**** You did it correctly congratulations!!!! " << std::endl;
0072     }
0073     jetInEvents_ += 1;
0074   }
0075 }