L1GenTreeProducer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    L1TriggerDPG/L1Ntuples
// Class:      L1GenTreeProducer
// 
//
// Original Author:  
//         Created:  
// $Id: L1GenTreeProducer.cc,v 1.8 2012/08/29 12:44:03 jbrooke Exp $
//
//


// system include files
#include <memory>

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

// input data formats
#include "DataFormats/JetReco/interface/GenJetCollection.h"

#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "HepMC/GenParticle.h"
#include "HepMC/GenVertex.h"
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"

#include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"


// ROOT output stuff
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "TTree.h"

#include "L1Trigger/L1TNtuples/interface/L1AnalysisGeneratorDataFormat.h"

//
// class declaration
//

class L1GenTreeProducer : public edm::EDAnalyzer {
public:
  explicit L1GenTreeProducer(const edm::ParameterSet&);
  ~L1GenTreeProducer();
  
  
private:
  virtual void beginJob(void) ;
  virtual void analyze(const edm::Event&, const edm::EventSetup&);
  virtual void endJob();  

private:

  unsigned maxL1Upgrade_;

  // output file
  edm::Service<TFileService> fs_;
  
  // tree
  TTree * tree_;
 
  // data format
  L1Analysis::L1AnalysisGeneratorDataFormat * l1GenData_;

  // EDM input tags
  edm::EDGetTokenT<reco::GenJetCollection> genJetToken_;
  edm::EDGetTokenT<reco::GenParticleCollection> genParticleToken_;
  edm::EDGetTokenT<std::vector<PileupSummaryInfo> > pileupInfoToken_;

};



L1GenTreeProducer::L1GenTreeProducer(const edm::ParameterSet& iConfig)
{

  genJetToken_ = consumes<reco::GenJetCollection>(iConfig.getUntrackedParameter<edm::InputTag>("genJetToken"));
  genParticleToken_ = consumes<reco::GenParticleCollection>(iConfig.getUntrackedParameter<edm::InputTag>("genParticleToken"));
  pileupInfoToken_ = consumes<std::vector<PileupSummaryInfo> >(iConfig.getUntrackedParameter<edm::InputTag>("pileupInfoToken"));
  
  // set up output
  tree_=fs_->make<TTree>("L1GenTree", "L1GenTree");
  tree_->Branch("Generator", "L1Analysis::L1AnalysisGeneratorDataFormat", &l1GenData_, 32000, 3);

}


L1GenTreeProducer::~L1GenTreeProducer()
{
 
  // do anything here that needs to be done at desctruction time
  // (e.g. close files, deallocate resources etc.)

}


//
// member functions
//

// ------------ method called to for each event  ------------
void
L1GenTreeProducer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
  
  l1GenData_->Reset();

  edm::Handle<reco::GenJetCollection> genJets;
  iEvent.getByToken(genJetToken_, genJets);


  if (genJets.isValid()){ 

    reco::GenJetCollection::const_iterator jetItr = genJets->begin();
    reco::GenJetCollection::const_iterator jetEnd = genJets->end();
    for( ; jetItr != jetEnd ; ++jetItr) {
      l1GenData_->jetPt.push_back( jetItr->pt() );
      l1GenData_->jetEta.push_back( jetItr->eta() );
      l1GenData_->jetPhi.push_back( jetItr->phi() );
      l1GenData_->nJet++;
    }

  } else {
    edm::LogWarning("MissingProduct") << "Gen jets not found. Branch will not be filled" << std::endl;
  }

  edm::Handle<reco::GenParticleCollection> genParticles;
  iEvent.getByToken(genParticleToken_, genParticles);

  for(size_t i = 0; i < genParticles->size(); ++ i) {
    const reco::GenParticle & p = (*genParticles)[i];
    int id = p.pdgId();
    //int st = p.status();  
    if (abs(id) == 13) {
      unsigned int nMo=p.numberOfMothers();
      //std::cout << "id " << id << "; st " << st 
      //<< "; nMo " << nMo << std::endl;
      for(unsigned int i=0;i<nMo;++i){
    //int thisParentID = dynamic_cast<const reco::GenParticle*>(p.mother(i))->pdgId();
    //std::cout << "   mother ID " << thisParentID << std::endl;
      }
    }
    
    //
    // See if the parent was interesting
    int parentID = -10000;
    unsigned int nMo=p.numberOfMothers();
    for(unsigned int i=0;i<nMo;++i){
      int thisParentID = dynamic_cast
    <const reco::GenParticle*>(p.mother(i))->pdgId();
      //
      // Is this a bottom hadron?
      int hundredsIndex = abs(thisParentID)/100;
      int thousandsIndex = abs(thisParentID)/1000;
      if ( ((abs(thisParentID) >= 23) && 
        (abs(thisParentID) <= 25)) ||
       (abs(thisParentID) == 6) ||
       (hundredsIndex == 5) ||
       (hundredsIndex == 4) ||
       (thousandsIndex == 5) ||
       (thousandsIndex == 4) 
       )
    parentID = thisParentID;
    }
    if ((parentID == -10000) && (nMo > 0)) 
      parentID = dynamic_cast
    <const reco::GenParticle*>(p.mother(0))->pdgId();
    //
    // If the parent of this particle is interesting, store all of the info
    if ((parentID != p.pdgId()) &&
    ((parentID > -9999) 
     || (abs(id) == 11)
     || (abs(id) == 13)
     || (abs(id) == 23)
     || (abs(id) == 24)
     || (abs(id) == 25)
     || (abs(id) == 4)
     || (abs(id) == 5)
     || (abs(id) == 6))
    )
      {
    l1GenData_->partId.push_back(p.pdgId());
    l1GenData_->partStat.push_back(p.status());
    l1GenData_->partPt.push_back(p.pt());
    l1GenData_->partEta.push_back(p.eta());
    l1GenData_->partPhi.push_back(p.phi());
    l1GenData_->partE.push_back(p.energy());
    l1GenData_->partParent.push_back(parentID);
      }
  }


  edm::Handle<std::vector<PileupSummaryInfo>> puInfoCollection;
  iEvent.getByToken(pileupInfoToken_, puInfoCollection);

  if (!puInfoCollection.isValid()) {
    throw cms::Exception("ProductNotValid") << "pileupInfoSource not valid";
  }

  // Loop over vector, find in-time entry, then store the relevant info
  std::vector<PileupSummaryInfo>::const_iterator puItr = puInfoCollection->begin();
  std::vector<PileupSummaryInfo>::const_iterator puEnd = puInfoCollection->end();
  for( ; puItr != puEnd; ++puItr) {
    int bx = puItr->getBunchCrossing();
    if (bx == 0) {
      l1GenData_->nPUPoissonMean = puItr->getTrueNumInteractions();
      l1GenData_->nVtx           = puItr->getPU_NumInteractions();
      break;
    }
  }



  tree_->Fill();

}

// ------------ method called once each job just before starting event loop  ------------
void 
L1GenTreeProducer::beginJob(void)
{
}

// ------------ method called once each job just after ending the event loop  ------------
void 
L1GenTreeProducer::endJob() {
}

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