L1TNtuples/plugins/L1GenTreeProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    L1TriggerDPG/L1Ntuples
0004 // Class:      L1GenTreeProducer
0005 //
0006 /**\class L1GenTreeProducer L1GenTreeProducer.cc L1TriggerDPG/L1Ntuples/src/L1GenTreeProducer.cc
0007
0008 Description: Produce L1 Extra tree
0009
0010 Implementation:
0011
0012 */
0013 //
0014 // Original Author:
0015 //         Created:
0016 // $Id: L1GenTreeProducer.cc,v 1.8 2012/08/29 12:44:03 jbrooke Exp $
0017 //
0018 //
0019
0020 // system include files
0021 #include <memory>
0022
0023 // framework
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0030
0031 // input data formats
0032 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0033 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0034
0035 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0036 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0037 #include "HepMC/GenParticle.h"
0038 #include "HepMC/GenVertex.h"
0039 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0040
0041 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
0042
0043 // ROOT output stuff
0044 #include "FWCore/ServiceRegistry/interface/Service.h"
0045 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0046 #include "TTree.h"
0047
0048 #include "L1Trigger/L1TNtuples/interface/L1AnalysisGeneratorDataFormat.h"
0049
0050 //
0051 // class declaration
0052 //
0053
0054 class L1GenTreeProducer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0055 public:
0056   explicit L1GenTreeProducer(const edm::ParameterSet&);
0057   ~L1GenTreeProducer() override;
0058
0059 private:
0060   void beginJob(void) override;
0061   void analyze(const edm::Event&, const edm::EventSetup&) override;
0062   void endJob() override;
0063
0064 private:
0065   unsigned maxL1Upgrade_;
0066
0067   // output file
0068   edm::Service<TFileService> fs_;
0069
0070   // tree
0071   TTree* tree_;
0072
0073   // data format
0074   std::unique_ptr<L1Analysis::L1AnalysisGeneratorDataFormat> l1GenData_;
0075
0076   // EDM input tags
0077   edm::EDGetTokenT<reco::GenJetCollection> genJetToken_;
0078   edm::EDGetTokenT<reco::GenParticleCollection> genParticleToken_;
0079   edm::EDGetTokenT<std::vector<PileupSummaryInfo>> pileupInfoToken_;
0080   edm::EDGetTokenT<GenEventInfoProduct> genInfoToken_;
0081 };
0082
0083 L1GenTreeProducer::L1GenTreeProducer(const edm::ParameterSet& iConfig) {
0084   genJetToken_ = consumes<reco::GenJetCollection>(iConfig.getUntrackedParameter<edm::InputTag>("genJetToken"));
0085   genParticleToken_ =
0086       consumes<reco::GenParticleCollection>(iConfig.getUntrackedParameter<edm::InputTag>("genParticleToken"));
0087   pileupInfoToken_ =
0088       consumes<std::vector<PileupSummaryInfo>>(iConfig.getUntrackedParameter<edm::InputTag>("pileupInfoToken"));
0089   genInfoToken_ = consumes<GenEventInfoProduct>(iConfig.getParameter<edm::InputTag>("genInfoToken"));
0090
0091   l1GenData_ = std::make_unique<L1Analysis::L1AnalysisGeneratorDataFormat>();
0092
0093   usesResource(TFileService::kSharedResource);
0094   // set up output
0095   tree_ = fs_->make<TTree>("L1GenTree", "L1GenTree");
0096   tree_->Branch("Generator", "L1Analysis::L1AnalysisGeneratorDataFormat", l1GenData_.get(), 32000, 3);
0097 }
0098
0099 L1GenTreeProducer::~L1GenTreeProducer() {
0100   // do anything here that needs to be done at desctruction time
0101   // (e.g. close files, deallocate resources etc.)
0102 }
0103
0104 //
0105 // member functions
0106 //
0107
0108 // ------------ method called to for each event  ------------
0109 void L1GenTreeProducer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0110   edm::Handle<GenEventInfoProduct> genInfo;
0111   iEvent.getByToken(genInfoToken_, genInfo);
0112
0113   l1GenData_->Reset();
0114
0115   if (genInfo.isValid()) {
0116     l1GenData_->weight = genInfo->weight();
0117     l1GenData_->pthat = genInfo->hasBinningValues() ? (genInfo->binningValues())[0] : 0.0;
0118   }
0119
0120   edm::Handle<reco::GenJetCollection> genJets;
0121   iEvent.getByToken(genJetToken_, genJets);
0122
0123   if (genJets.isValid()) {
0124     reco::GenJetCollection::const_iterator jetItr = genJets->begin();
0125     reco::GenJetCollection::const_iterator jetEnd = genJets->end();
0126     for (; jetItr != jetEnd; ++jetItr) {
0127       l1GenData_->jetPt.push_back(jetItr->pt());
0128       l1GenData_->jetEta.push_back(jetItr->eta());
0129       l1GenData_->jetPhi.push_back(jetItr->phi());
0130       l1GenData_->nJet++;
0131     }
0132
0133   } else {
0134     edm::LogWarning("MissingProduct") << "Gen jets not found. Branch will not be filled" << std::endl;
0135   }
0136
0137   edm::Handle<reco::GenParticleCollection> genParticles;
0138   iEvent.getByToken(genParticleToken_, genParticles);
0139
0140   if (genParticles.isValid()) {
0141     int nPart{0};
0142
0143     for (size_t i = 0; i < genParticles->size(); ++i) {
0144       const reco::GenParticle& p = (*genParticles)[i];
0145       int id = p.pdgId();
0146
0147       // See if the parent was interesting
0148       int parentID = -10000;
0149       unsigned int nMo = p.numberOfMothers();
0150       for (unsigned int i = 0; i < nMo; ++i) {
0151         int thisParentID = dynamic_cast<const reco::GenParticle*>(p.mother(i))->pdgId();
0152         //
0153         // Is this a bottom hadron?
0154         int hundredsIndex = abs(thisParentID) / 100;
0155         int thousandsIndex = abs(thisParentID) / 1000;
0156         if (((abs(thisParentID) >= 23) && (abs(thisParentID) <= 25)) || (abs(thisParentID) == 6) ||
0157             (hundredsIndex == 5) || (hundredsIndex == 4) || (thousandsIndex == 5) || (thousandsIndex == 4))
0158           parentID = thisParentID;
0159       }
0160       if ((parentID == -10000) && (nMo > 0))
0161         parentID = dynamic_cast<const reco::GenParticle*>(p.mother(0))->pdgId();
0162       //
0163       // If the parent of this particle is interesting, store all of the info
0164       if ((parentID != p.pdgId()) &&
0165           ((parentID > -9999) || (abs(id) == 11) || (abs(id) == 13) || (abs(id) == 23) || (abs(id) == 24) ||
0166            (abs(id) == 25) || (abs(id) == 4) || (abs(id) == 5) || (abs(id) == 6))) {
0167         l1GenData_->partId.push_back(p.pdgId());
0168         l1GenData_->partStat.push_back(p.status());
0169         l1GenData_->partPt.push_back(p.pt());
0170         l1GenData_->partEta.push_back(p.eta());
0171         l1GenData_->partPhi.push_back(p.phi());
0172         l1GenData_->partE.push_back(p.energy());
0173         l1GenData_->partParent.push_back(parentID);
0174         l1GenData_->partCh.push_back(p.charge());
0175         ++nPart;
0176       }
0177     }
0178     l1GenData_->nPart = nPart;
0179   }
0180
0181   edm::Handle<std::vector<PileupSummaryInfo>> puInfoCollection;
0182   iEvent.getByToken(pileupInfoToken_, puInfoCollection);
0183
0184   if (!puInfoCollection.isValid()) {
0185     throw cms::Exception("ProductNotValid") << "pileupInfoSource not valid";
0186   }
0187
0188   // Loop over vector, find in-time entry, then store the relevant info
0189   std::vector<PileupSummaryInfo>::const_iterator puItr = puInfoCollection->begin();
0190   std::vector<PileupSummaryInfo>::const_iterator puEnd = puInfoCollection->end();
0191   for (; puItr != puEnd; ++puItr) {
0192     int bx = puItr->getBunchCrossing();
0193     if (bx == 0) {
0194       l1GenData_->nMeanPU = puItr->getTrueNumInteractions();
0195       l1GenData_->nVtx = puItr->getPU_NumInteractions();
0196       break;
0197     }
0198   }
0199
0200   tree_->Fill();
0201 }
0202
0203 // ------------ method called once each job just before starting event loop  ------------
0204 void L1GenTreeProducer::beginJob(void) {}
0205
0206 // ------------ method called once each job just after ending the event loop  ------------
0207 void L1GenTreeProducer::endJob() {}
0208
0209 //define this as a plug-in
0210 DEFINE_FWK_MODULE(L1GenTreeProducer);