Go to the documentation of this file.00001 #include <iostream>
00002 #include <sstream>
00003 #include <istream>
00004 #include <fstream>
00005 #include <iomanip>
00006 #include <string>
00007 #include <cmath>
00008 #include <functional>
00009
00010 #include "JetMETCorrections/MCJet/plugins/CaloMCTruthTreeProducer.h"
00011 #include "JetMETCorrections/MCJet/plugins/JetUtilMC.h"
00012 #include "FWCore/Framework/interface/Event.h"
00013 #include "DataFormats/Common/interface/Handle.h"
00014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00015 #include "FWCore/Framework/interface/EventSetup.h"
00016 #include "FWCore/Framework/interface/ESHandle.h"
00017 #include "Geometry/Records/interface/IdealGeometryRecord.h"
00018 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
00019 #include "DataFormats/JetReco/interface/CaloJet.h"
00020 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
00021 #include "DataFormats/JetReco/interface/GenJet.h"
00022 #include "DataFormats/JetReco/interface/GenJetCollection.h"
00023 using namespace edm;
00024 using namespace reco;
00025 using namespace std;
00026
00027
00028
00029 CaloMCTruthTreeProducer::CaloMCTruthTreeProducer(edm::ParameterSet const& cfg)
00030 {
00031 jets_ = cfg.getParameter<std::string> ("jets");
00032 genjets_ = cfg.getParameter<std::string> ("genjets");
00033 histogramFile_ = cfg.getParameter<std::string> ("histogramFile");
00034 }
00036 void CaloMCTruthTreeProducer::beginJob()
00037 {
00038 file_ = new TFile(histogramFile_.c_str(),"RECREATE");
00039 mcTruthTree_ = new TTree("mcTruthTree","mcTruthTree");
00040
00041 mcTruthTree_->Branch("ptJet", &ptJet_, "ptJet_/F");
00042 mcTruthTree_->Branch("ptGen", &ptGen_, "ptGen_/F");
00043 mcTruthTree_->Branch("ptHat", &ptHat_, "ptHat_/F");
00044 mcTruthTree_->Branch("emfJet", &emfJet_, "emfJet_/F");
00045 mcTruthTree_->Branch("etaJet", &etaJet_, "etaJet_/F");
00046 mcTruthTree_->Branch("etaGen", &etaGen_, "etaGen_/F");
00047 mcTruthTree_->Branch("phiJet", &phiJet_, "phiJet_/F");
00048 mcTruthTree_->Branch("phiGen", &phiGen_, "phiGen_/F");
00049 mcTruthTree_->Branch("dR", &dR_, "dR_/F");
00050 mcTruthTree_->Branch("rank", &rank_, "rank_/I");
00051 }
00053 void CaloMCTruthTreeProducer::endJob()
00054 {
00055 if (file_ !=0)
00056 {
00057 file_->cd();
00058 mcTruthTree_->Write();
00059 }
00060 file_ = 0;
00061 }
00063 void CaloMCTruthTreeProducer::analyze(edm::Event const& event, edm::EventSetup const& iSetup)
00064 {
00065 edm::Handle<GenJetCollection> genjets;
00066 edm::Handle<CaloJetCollection> jets;
00067 edm::Handle<GenEventInfoProduct> hEventInfo;
00068 CaloJetCollection::const_iterator i_jet,i_matched;
00069 GenJetCollection::const_iterator i_genjet;
00070 event.getByLabel (genjets_,genjets);
00071 event.getByLabel (jets_,jets);
00072 event.getByLabel("generator",hEventInfo);
00073 ptHat_ = hEventInfo->binningValues()[0];
00074 float rr;
00075 int njet(0);
00076 if (jets->size()>0 && genjets->size()>0)
00077 {
00078 for (i_genjet = genjets->begin(); i_genjet != genjets->end(); i_genjet++)
00079 {
00080 float rmin(99);
00081 for(i_jet = jets->begin();i_jet != jets->end(); i_jet++)
00082 {
00083 rr=radius(i_genjet,i_jet);
00084 if (rr<rmin)
00085 {
00086 rmin = rr;
00087 i_matched = i_jet;
00088 }
00089 }
00090 ptGen_ = i_genjet->pt();
00091 etaGen_ = i_genjet->eta();
00092 phiGen_ = i_genjet->phi();
00093 ptJet_ = i_matched->pt();
00094 etaJet_ = i_matched->eta();
00095 phiJet_ = i_matched->phi();
00096 emfJet_ = i_matched->emEnergyFraction();
00097 dR_ = rmin;
00098 rank_ = njet;
00099 mcTruthTree_->Fill();
00100 njet++;
00101 }
00102 }
00103 }
00105 CaloMCTruthTreeProducer::~CaloMCTruthTreeProducer()
00106 {
00107 delete file_;
00108 delete mcTruthTree_;
00109 }
00110