File indexing completed on 2023-03-17 10:43:32
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
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
0034 #include "FWCore/ServiceRegistry/interface/Service.h"
0035 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0036
0037 #include "RecoJets/JetProducers/interface/JetMatchingTools.h"
0038
0039
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
0097 clearTreeVectors();
0098
0099
0100 edm::Handle<reco::GenJetCollection> genJets;
0101 iEvent.getByToken(tok_jets_, genJets);
0102
0103
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 }
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 }
0155
0156 }
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
0222 DEFINE_FWK_MODULE(IsolatedParticlesGeneratedJets);