IsolatedParticles/plugins/IsolatedParticlesGeneratedJets.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    IsolatedParticles
0004 // Class:      IsolatedParticlesGeneratedJets
0005 //
0006 /**\class IsolatedParticlesGeneratedJets IsolatedParticlesGeneratedJets.cc Calibration/IsolatedParticles/plugins/IsolatedParticlesGeneratedJets.cc
0007
0008  Description: Studies properties of jets at generator level in context of
0009               isolated particles
0010
0011  Implementation:
0012      <Notes on implementation>
0013 */
0014 //
0015 // Original Author:  Seema Sharma
0016 //         Created:  Thu Mar  4 18:52:02 CST 2010
0017 //
0018 //
0019
0020 // user include files
0021 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0022 #include "DataFormats/JetReco/interface/GenJet.h"
0023 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0024 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
0025
0026 #include "FWCore/Framework/interface/Frameworkfwd.h"
0027 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0028 #include "FWCore/Framework/interface/Event.h"
0029 #include "FWCore/Framework/interface/MakerMacros.h"
0030 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032
0033 //TFile Service
0034 #include "FWCore/ServiceRegistry/interface/Service.h"
0035 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0036
0037 #include "RecoJets/JetProducers/interface/JetMatchingTools.h"
0038
0039 // root objects
0040 #include "TROOT.h"
0041 #include "TSystem.h"
0042 #include "TFile.h"
0043 #include "TH1F.h"
0044 #include "TH2F.h"
0045 #include "TProfile.h"
0046 #include "TDirectory.h"
0047 #include "TTree.h"
0048
0049 class IsolatedParticlesGeneratedJets : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0050 public:
0051   explicit IsolatedParticlesGeneratedJets(const edm::ParameterSet &);
0052   ~IsolatedParticlesGeneratedJets() override {}
0053
0054   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0055
0056 private:
0057   void beginJob() override;
0058   void analyze(const edm::Event &, const edm::EventSetup &) override;
0059   void endJob() override {}
0060
0061   void bookHistograms();
0062   void clearTreeVectors();
0063
0064   const bool debug_;
0065   TTree *tree_;
0066
0067   const edm::EDGetTokenT<reco::GenJetCollection> tok_jets_;
0068   const edm::EDGetTokenT<reco::GenParticleCollection> tok_parts_;
0069
0070   std::vector<int> *t_gjetN;
0071   std::vector<double> *t_gjetE, *t_gjetPt, *t_gjetEta, *t_gjetPhi;
0072   std::vector<std::vector<double> > *t_jetTrkP;
0073   std::vector<std::vector<double> > *t_jetTrkPt;
0074   std::vector<std::vector<double> > *t_jetTrkEta;
0075   std::vector<std::vector<double> > *t_jetTrkPhi;
0076   std::vector<std::vector<double> > *t_jetTrkPdg;
0077   std::vector<std::vector<double> > *t_jetTrkCharge;
0078 };
0079
0080 IsolatedParticlesGeneratedJets::IsolatedParticlesGeneratedJets(const edm::ParameterSet &iConfig)
0081     : debug_(iConfig.getUntrackedParameter<bool>("Debug", false)),
0082       tok_jets_(consumes<reco::GenJetCollection>(iConfig.getParameter<edm::InputTag>("JetSource"))),
0083       tok_parts_(consumes<reco::GenParticleCollection>(iConfig.getParameter<edm::InputTag>("ParticleSource"))) {
0084   usesResource(TFileService::kSharedResource);
0085 }
0086
0087 void IsolatedParticlesGeneratedJets::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0088   edm::ParameterSetDescription desc;
0089   desc.addUntracked<bool>("Debug", true);
0090   desc.add<edm::InputTag>("JetSource", edm::InputTag("ak5GenJets"));
0091   desc.add<edm::InputTag>("ParticleSource", edm::InputTag("genParticles"));
0092   descriptions.add("isolatedParticlesGeneratedJets", desc);
0093 }
0094
0095 void IsolatedParticlesGeneratedJets::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0096   //using namespace edm;
0097   clearTreeVectors();
0098
0099   //=== genJet information
0100   edm::Handle<reco::GenJetCollection> genJets;
0101   iEvent.getByToken(tok_jets_, genJets);
0102
0103   //=== genJet information
0104   edm::Handle<reco::GenParticleCollection> genParticles;
0105   iEvent.getByToken(tok_parts_, genParticles);
0106
0107   JetMatchingTools jetMatching(iEvent, consumesCollector());
0108   std::vector<std::vector<const reco::GenParticle *> > genJetConstituents(genJets->size());
0109
0110   int njets = 0;
0111   for (unsigned iGenJet = 0; iGenJet < genJets->size(); ++iGenJet) {
0112     const reco::GenJet &genJet = (*genJets)[iGenJet];
0113
0114     double genJetE = genJet.energy();
0115     double genJetPt = genJet.pt();
0116     double genJetEta = genJet.eta();
0117     double genJetPhi = genJet.phi();
0118
0119     if (genJetPt > 30.0 && std::abs(genJetEta) < 3.0) {
0120       njets++;
0121
0122       std::vector<const reco::GenParticle *> genJetConstituents = jetMatching.getGenParticles((*genJets)[iGenJet]);
0123       std::vector<double> v_trkP, v_trkPt, v_trkEta, v_trkPhi, v_trkPdg, v_trkCharge;
0124
0125       if (debug_)
0126         edm::LogVerbatim("IsoTrack") << "Jet(pt,Eta,Phi) " << genJetPt << " " << genJetEta << " " << genJetPhi;
0127       for (unsigned int ic = 0; ic < genJetConstituents.size(); ic++) {
0128         if (debug_)
0129           edm::LogVerbatim("IsoTrack") << "p,pt,eta,phi " << genJetConstituents[ic]->p() << " "
0130                                        << genJetConstituents[ic]->pt() << " " << genJetConstituents[ic]->eta() << " "
0131                                        << genJetConstituents[ic]->phi();
0132
0133         v_trkP.push_back(genJetConstituents[ic]->p());
0134         v_trkPt.push_back(genJetConstituents[ic]->pt());
0135         v_trkEta.push_back(genJetConstituents[ic]->eta());
0136         v_trkPhi.push_back(genJetConstituents[ic]->phi());
0137         v_trkPdg.push_back(genJetConstituents[ic]->pdgId());
0138         v_trkCharge.push_back(genJetConstituents[ic]->charge());
0139
0140       }  //loop over genjet constituents
0141
0142       t_gjetE->push_back(genJetE);
0143       t_gjetPt->push_back(genJetPt);
0144       t_gjetEta->push_back(genJetEta);
0145       t_gjetPhi->push_back(genJetPhi);
0146
0147       t_jetTrkP->push_back(v_trkP);
0148       t_jetTrkPt->push_back(v_trkPt);
0149       t_jetTrkEta->push_back(v_trkEta);
0150       t_jetTrkPhi->push_back(v_trkPhi);
0151       t_jetTrkPdg->push_back(v_trkPdg);
0152       t_jetTrkCharge->push_back(v_trkCharge);
0153
0154     }  // if jetPt>30
0155
0156   }  //loop over genjets
0157
0158   t_gjetN->push_back(njets);
0159
0160   if (debug_) {
0161     unsigned int indx = 0;
0162     reco::GenParticleCollection::const_iterator ig = genParticles->begin();
0163     for (; ig != genParticles->end(); ++ig, ++indx) {
0164       edm::LogVerbatim("IsoTrack") << "Track " << indx << " Status " << ig->status() << " charge " << ig->charge()
0165                                    << " pdgId " << ig->pdgId() << " mass " << ig->mass() << " P " << ig->momentum()
0166                                    << " E " << ig->energy() << " Origin " << ig->vertex();
0167     }
0168   }
0169
0170   tree_->Fill();
0171 }
0172
0173 void IsolatedParticlesGeneratedJets::beginJob() { bookHistograms(); }
0174
0175 void IsolatedParticlesGeneratedJets::clearTreeVectors() {
0176   t_gjetN->clear();
0177   t_gjetE->clear();
0178   t_gjetPt->clear();
0179   t_gjetEta->clear();
0180   t_gjetPhi->clear();
0181
0182   t_jetTrkP->clear();
0183   t_jetTrkPt->clear();
0184   t_jetTrkEta->clear();
0185   t_jetTrkPhi->clear();
0186   t_jetTrkPdg->clear();
0187   t_jetTrkCharge->clear();
0188 }
0189
0190 void IsolatedParticlesGeneratedJets::bookHistograms() {
0191   edm::Service<TFileService> fs;
0192   tree_ = fs->make<TTree>("tree", "tree");
0193
0194   t_gjetN = new std::vector<int>();
0195   t_gjetE = new std::vector<double>();
0196   t_gjetPt = new std::vector<double>();
0197   t_gjetEta = new std::vector<double>();
0198   t_gjetPhi = new std::vector<double>();
0199
0200   t_jetTrkP = new std::vector<std::vector<double> >();
0201   t_jetTrkPt = new std::vector<std::vector<double> >();
0202   t_jetTrkEta = new std::vector<std::vector<double> >();
0203   t_jetTrkPhi = new std::vector<std::vector<double> >();
0204   t_jetTrkPdg = new std::vector<std::vector<double> >();
0205   t_jetTrkCharge = new std::vector<std::vector<double> >();
0206
0207   tree_->Branch("t_gjetN", "std::vector<int>", &t_gjetN);
0208   tree_->Branch("t_gjetE", "std::vector<double>", &t_gjetE);
0209   tree_->Branch("t_gjetPt", "std::vector<double>", &t_gjetPt);
0210   tree_->Branch("t_gjetEta", "std::vector<double>", &t_gjetEta);
0211   tree_->Branch("t_gjetPhi", "std::vector<double>", &t_gjetPhi);
0212
0213   tree_->Branch("t_jetTrkP", "std::vector<vector<double> >", &t_jetTrkP);
0214   tree_->Branch("t_jetTrkPt", "std::vector<vector<double> >", &t_jetTrkPt);
0215   tree_->Branch("t_jetTrkEta", "std::vector<vector<double> >", &t_jetTrkEta);
0216   tree_->Branch("t_jetTrkPhi", "std::vector<vector<double> >", &t_jetTrkPhi);
0217   tree_->Branch("t_jetTrkPdg", "std::vector<vector<double> >", &t_jetTrkPdg);
0218   tree_->Branch("t_jetTrkCharge", "std::vector<vector<double> >", &t_jetTrkCharge);
0219 }
0220
0221 //define this as a plug-in
0222 DEFINE_FWK_MODULE(IsolatedParticlesGeneratedJets);