MCJet/plugins/PFMCTruthTreeProducer.cc

0001 #include <iostream>
0002 #include <sstream>
0003 #include <istream>
0004 #include <fstream>
0005 #include <iomanip>
0006 #include <string>
0007 #include <cmath>
0008 #include <functional>
0009
0010 #include "JetMETCorrections/MCJet/plugins/PFMCTruthTreeProducer.h"
0011 #include "JetMETCorrections/MCJet/plugins/JetUtilMC.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "DataFormats/Common/interface/Handle.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/Framework/interface/ESHandle.h"
0017 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0018 #include "DataFormats/JetReco/interface/PFJet.h"
0019 #include "DataFormats/JetReco/interface/GenJet.h"
0020 using namespace edm;
0021 using namespace reco;
0022 using namespace std;
0023 //namespace cms
0024 //{
0025
0026 PFMCTruthTreeProducer::PFMCTruthTreeProducer(edm::ParameterSet const& cfg) {
0027   jets_ = consumes<PFJetCollection>(edm::InputTag(cfg.getParameter<std::string>("jets")));
0028   genjets_ = consumes<GenJetCollection>(edm::InputTag(cfg.getParameter<std::string>("genjets")));
0029   gen_ = consumes<GenEventInfoProduct>(edm::InputTag("generator"));
0030   histogramFile_ = cfg.getParameter<std::string>("histogramFile");
0031 }
0032 //////////////////////////////////////////////////////////////////////////////////////////
0033 void PFMCTruthTreeProducer::beginJob() {
0034   file_ = new TFile(histogramFile_.c_str(), "RECREATE");
0035   mcTruthTree_ = new TTree("mcTruthTree", "mcTruthTree");
0036
0037   mcTruthTree_->Branch("ptJet", &ptJet_, "ptJet_/F");
0038   mcTruthTree_->Branch("ptGen", &ptGen_, "ptGen_/F");
0039   mcTruthTree_->Branch("ptHat", &ptHat_, "ptHat_/F");
0040   mcTruthTree_->Branch("chfJet", &chfJet_, "chfJet_/F");
0041   mcTruthTree_->Branch("nhfJet", &nhfJet_, "nhfJet_/F");
0042   mcTruthTree_->Branch("cemfJet", &cemfJet_, "cemfJet_/F");
0043   mcTruthTree_->Branch("nemfJet", &nemfJet_, "nemfJet_/F");
0044   mcTruthTree_->Branch("cmultiJet", &cmultiJet_, "cmultiJet_/I");
0045   mcTruthTree_->Branch("nmultiJet", &nmultiJet_, "nmultiJet_/I");
0046   mcTruthTree_->Branch("etaJet", &etaJet_, "etaJet_/F");
0047   mcTruthTree_->Branch("etaGen", &etaGen_, "etaGen_/F");
0048   mcTruthTree_->Branch("phiJet", &phiJet_, "phiJet_/F");
0049   mcTruthTree_->Branch("phiGen", &phiGen_, "phiGen_/F");
0050   mcTruthTree_->Branch("dR", &dR_, "dR_/F");
0051   mcTruthTree_->Branch("rank", &rank_, "rank_/I");
0052 }
0053 //////////////////////////////////////////////////////////////////////////////////////////
0054 void PFMCTruthTreeProducer::endJob() {
0055   if (file_ != nullptr) {
0056     file_->cd();
0057     mcTruthTree_->Write();
0058   }
0059   file_ = nullptr;
0060 }
0061 //////////////////////////////////////////////////////////////////////////////////////////
0062 void PFMCTruthTreeProducer::analyze(edm::Event const& event, edm::EventSetup const& iSetup) {
0063   edm::Handle<GenJetCollection> genjets;
0064   edm::Handle<PFJetCollection> jets;
0065   edm::Handle<GenEventInfoProduct> hEventInfo;
0066   PFJetCollection::const_iterator i_jet, i_matched;
0067   GenJetCollection::const_iterator i_genjet;
0068   event.getByToken(genjets_, genjets);
0069   event.getByToken(jets_, jets);
0070   event.getByToken(gen_, hEventInfo);
0071   ptHat_ = hEventInfo->binningValues()[0];
0072   float rr;
0073   int njet(0);
0074   if (!jets->empty() && !genjets->empty()) {
0075     for (i_genjet = genjets->begin(); i_genjet != genjets->end(); i_genjet++) {
0076       float rmin(99);
0077       for (i_jet = jets->begin(); i_jet != jets->end(); i_jet++) {
0078         rr = radius(i_genjet, i_jet);
0079         if (rr < rmin) {
0080           rmin = rr;
0081           i_matched = i_jet;
0082         }
0083       }
0084       ptGen_ = i_genjet->pt();
0085       etaGen_ = i_genjet->eta();
0086       phiGen_ = i_genjet->phi();
0087       ptJet_ = i_matched->pt();
0088       etaJet_ = i_matched->eta();
0089       phiJet_ = i_matched->phi();
0090       chfJet_ = i_matched->chargedHadronEnergyFraction();
0091       nhfJet_ = i_matched->neutralHadronEnergyFraction();
0092       cemfJet_ = i_matched->chargedEmEnergyFraction();
0093       nemfJet_ = i_matched->neutralEmEnergyFraction();
0094       cmultiJet_ = i_matched->chargedMultiplicity();
0095       nmultiJet_ = i_matched->neutralMultiplicity();
0096       dR_ = rmin;
0097       rank_ = njet;
0098       mcTruthTree_->Fill();
0099       njet++;
0100     }
0101   }
0102 }
0103 //////////////////////////////////////////////////////////////////////////////////////////
0104 PFMCTruthTreeProducer::~PFMCTruthTreeProducer() {
0105   delete file_;
0106   delete mcTruthTree_;
0107 }
0108 //}