IsolatedParticlesGeneratedJets.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    IsolatedParticles
// Class:      IsolatedParticlesGeneratedJets
//
//
// Original Author:  Seema Sharma
//         Created:  Thu Mar  4 18:52:02 CST 2010
//
//

// user include files
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
#include "DataFormats/JetReco/interface/GenJet.h"
#include "DataFormats/JetReco/interface/GenJetCollection.h"
#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/one/EDAnalyzer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

//TFile Service
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"

#include "RecoJets/JetProducers/interface/JetMatchingTools.h"

// root objects
#include "TROOT.h"
#include "TSystem.h"
#include "TFile.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TProfile.h"
#include "TDirectory.h"
#include "TTree.h"

class IsolatedParticlesGeneratedJets : public edm::one::EDAnalyzer<edm::one::SharedResources> {
public:
  explicit IsolatedParticlesGeneratedJets(const edm::ParameterSet &);
  ~IsolatedParticlesGeneratedJets() override {}

  static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);

private:
  void beginJob() override;
  void analyze(const edm::Event &, const edm::EventSetup &) override;
  void endJob() override {}

  void bookHistograms();
  void clearTreeVectors();

  const bool debug_;
  TTree *tree_;

  const edm::EDGetTokenT<reco::GenJetCollection> tok_jets_;
  const edm::EDGetTokenT<reco::GenParticleCollection> tok_parts_;

  std::vector<int> *t_gjetN;
  std::vector<double> *t_gjetE, *t_gjetPt, *t_gjetEta, *t_gjetPhi;
  std::vector<std::vector<double> > *t_jetTrkP;
  std::vector<std::vector<double> > *t_jetTrkPt;
  std::vector<std::vector<double> > *t_jetTrkEta;
  std::vector<std::vector<double> > *t_jetTrkPhi;
  std::vector<std::vector<double> > *t_jetTrkPdg;
  std::vector<std::vector<double> > *t_jetTrkCharge;
};

IsolatedParticlesGeneratedJets::IsolatedParticlesGeneratedJets(const edm::ParameterSet &iConfig)
    : debug_(iConfig.getUntrackedParameter<bool>("Debug", false)),
      tok_jets_(consumes<reco::GenJetCollection>(iConfig.getParameter<edm::InputTag>("JetSource"))),
      tok_parts_(consumes<reco::GenParticleCollection>(iConfig.getParameter<edm::InputTag>("ParticleSource"))) {
  usesResource(TFileService::kSharedResource);
}

void IsolatedParticlesGeneratedJets::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
  edm::ParameterSetDescription desc;
  desc.addUntracked<bool>("Debug", true);
  desc.add<edm::InputTag>("JetSource", edm::InputTag("ak5GenJets"));
  desc.add<edm::InputTag>("ParticleSource", edm::InputTag("genParticles"));
  descriptions.add("isolatedParticlesGeneratedJets", desc);
}

void IsolatedParticlesGeneratedJets::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
  //using namespace edm;
  clearTreeVectors();

  //=== genJet information
  edm::Handle<reco::GenJetCollection> genJets;
  iEvent.getByToken(tok_jets_, genJets);

  //=== genJet information
  edm::Handle<reco::GenParticleCollection> genParticles;
  iEvent.getByToken(tok_parts_, genParticles);

  JetMatchingTools jetMatching(iEvent, consumesCollector());
  std::vector<std::vector<const reco::GenParticle *> > genJetConstituents(genJets->size());

  int njets = 0;
  for (unsigned iGenJet = 0; iGenJet < genJets->size(); ++iGenJet) {
    const reco::GenJet &genJet = (*genJets)[iGenJet];

    double genJetE = genJet.energy();
    double genJetPt = genJet.pt();
    double genJetEta = genJet.eta();
    double genJetPhi = genJet.phi();

    if (genJetPt > 30.0 && std::abs(genJetEta) < 3.0) {
      njets++;

      std::vector<const reco::GenParticle *> genJetConstituents = jetMatching.getGenParticles((*genJets)[iGenJet]);
      std::vector<double> v_trkP, v_trkPt, v_trkEta, v_trkPhi, v_trkPdg, v_trkCharge;

      if (debug_)
        edm::LogVerbatim("IsoTrack") << "Jet(pt,Eta,Phi) " << genJetPt << " " << genJetEta << " " << genJetPhi;
      for (unsigned int ic = 0; ic < genJetConstituents.size(); ic++) {
        if (debug_)
          edm::LogVerbatim("IsoTrack") << "p,pt,eta,phi " << genJetConstituents[ic]->p() << " "
                                       << genJetConstituents[ic]->pt() << " " << genJetConstituents[ic]->eta() << " "
                                       << genJetConstituents[ic]->phi();

        v_trkP.push_back(genJetConstituents[ic]->p());
        v_trkPt.push_back(genJetConstituents[ic]->pt());
        v_trkEta.push_back(genJetConstituents[ic]->eta());
        v_trkPhi.push_back(genJetConstituents[ic]->phi());
        v_trkPdg.push_back(genJetConstituents[ic]->pdgId());
        v_trkCharge.push_back(genJetConstituents[ic]->charge());

      }  //loop over genjet constituents

      t_gjetE->push_back(genJetE);
      t_gjetPt->push_back(genJetPt);
      t_gjetEta->push_back(genJetEta);
      t_gjetPhi->push_back(genJetPhi);

      t_jetTrkP->push_back(v_trkP);
      t_jetTrkPt->push_back(v_trkPt);
      t_jetTrkEta->push_back(v_trkEta);
      t_jetTrkPhi->push_back(v_trkPhi);
      t_jetTrkPdg->push_back(v_trkPdg);
      t_jetTrkCharge->push_back(v_trkCharge);

    }  // if jetPt>30

  }  //loop over genjets

  t_gjetN->push_back(njets);

  if (debug_) {
    unsigned int indx = 0;
    reco::GenParticleCollection::const_iterator ig = genParticles->begin();
    for (; ig != genParticles->end(); ++ig, ++indx) {
      edm::LogVerbatim("IsoTrack") << "Track " << indx << " Status " << ig->status() << " charge " << ig->charge()
                                   << " pdgId " << ig->pdgId() << " mass " << ig->mass() << " P " << ig->momentum()
                                   << " E " << ig->energy() << " Origin " << ig->vertex();
    }
  }

  tree_->Fill();
}

void IsolatedParticlesGeneratedJets::beginJob() { bookHistograms(); }

void IsolatedParticlesGeneratedJets::clearTreeVectors() {
  t_gjetN->clear();
  t_gjetE->clear();
  t_gjetPt->clear();
  t_gjetEta->clear();
  t_gjetPhi->clear();

  t_jetTrkP->clear();
  t_jetTrkPt->clear();
  t_jetTrkEta->clear();
  t_jetTrkPhi->clear();
  t_jetTrkPdg->clear();
  t_jetTrkCharge->clear();
}

void IsolatedParticlesGeneratedJets::bookHistograms() {
  edm::Service<TFileService> fs;
  tree_ = fs->make<TTree>("tree", "tree");

  t_gjetN = new std::vector<int>();
  t_gjetE = new std::vector<double>();
  t_gjetPt = new std::vector<double>();
  t_gjetEta = new std::vector<double>();
  t_gjetPhi = new std::vector<double>();

  t_jetTrkP = new std::vector<std::vector<double> >();
  t_jetTrkPt = new std::vector<std::vector<double> >();
  t_jetTrkEta = new std::vector<std::vector<double> >();
  t_jetTrkPhi = new std::vector<std::vector<double> >();
  t_jetTrkPdg = new std::vector<std::vector<double> >();
  t_jetTrkCharge = new std::vector<std::vector<double> >();

  tree_->Branch("t_gjetN", "std::vector<int>", &t_gjetN);
  tree_->Branch("t_gjetE", "std::vector<double>", &t_gjetE);
  tree_->Branch("t_gjetPt", "std::vector<double>", &t_gjetPt);
  tree_->Branch("t_gjetEta", "std::vector<double>", &t_gjetEta);
  tree_->Branch("t_gjetPhi", "std::vector<double>", &t_gjetPhi);

  tree_->Branch("t_jetTrkP", "std::vector<vector<double> >", &t_jetTrkP);
  tree_->Branch("t_jetTrkPt", "std::vector<vector<double> >", &t_jetTrkPt);
  tree_->Branch("t_jetTrkEta", "std::vector<vector<double> >", &t_jetTrkEta);
  tree_->Branch("t_jetTrkPhi", "std::vector<vector<double> >", &t_jetTrkPhi);
  tree_->Branch("t_jetTrkPdg", "std::vector<vector<double> >", &t_jetTrkPdg);
  tree_->Branch("t_jetTrkCharge", "std::vector<vector<double> >", &t_jetTrkCharge);
}

//define this as a plug-in
DEFINE_FWK_MODULE(IsolatedParticlesGeneratedJets);