Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:23:30

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/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::EDAnalyzer {
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   // set up output
0094   tree_ = fs_->make<TTree>("L1GenTree", "L1GenTree");
0095   tree_->Branch("Generator", "L1Analysis::L1AnalysisGeneratorDataFormat", l1GenData_.get(), 32000, 3);
0096 }
0097 
0098 L1GenTreeProducer::~L1GenTreeProducer() {
0099   // do anything here that needs to be done at desctruction time
0100   // (e.g. close files, deallocate resources etc.)
0101 }
0102 
0103 //
0104 // member functions
0105 //
0106 
0107 // ------------ method called to for each event  ------------
0108 void L1GenTreeProducer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0109   edm::Handle<GenEventInfoProduct> genInfo;
0110   iEvent.getByToken(genInfoToken_, genInfo);
0111 
0112   l1GenData_->Reset();
0113 
0114   if (genInfo.isValid()) {
0115     l1GenData_->weight = genInfo->weight();
0116     l1GenData_->pthat = genInfo->hasBinningValues() ? (genInfo->binningValues())[0] : 0.0;
0117   }
0118 
0119   edm::Handle<reco::GenJetCollection> genJets;
0120   iEvent.getByToken(genJetToken_, genJets);
0121 
0122   if (genJets.isValid()) {
0123     reco::GenJetCollection::const_iterator jetItr = genJets->begin();
0124     reco::GenJetCollection::const_iterator jetEnd = genJets->end();
0125     for (; jetItr != jetEnd; ++jetItr) {
0126       l1GenData_->jetPt.push_back(jetItr->pt());
0127       l1GenData_->jetEta.push_back(jetItr->eta());
0128       l1GenData_->jetPhi.push_back(jetItr->phi());
0129       l1GenData_->nJet++;
0130     }
0131 
0132   } else {
0133     edm::LogWarning("MissingProduct") << "Gen jets not found. Branch will not be filled" << std::endl;
0134   }
0135 
0136   edm::Handle<reco::GenParticleCollection> genParticles;
0137   iEvent.getByToken(genParticleToken_, genParticles);
0138 
0139   if (genParticles.isValid()) {
0140     int nPart{0};
0141 
0142     for (size_t i = 0; i < genParticles->size(); ++i) {
0143       const reco::GenParticle& p = (*genParticles)[i];
0144       int id = p.pdgId();
0145 
0146       // See if the parent was interesting
0147       int parentID = -10000;
0148       unsigned int nMo = p.numberOfMothers();
0149       for (unsigned int i = 0; i < nMo; ++i) {
0150         int thisParentID = dynamic_cast<const reco::GenParticle*>(p.mother(i))->pdgId();
0151         //
0152         // Is this a bottom hadron?
0153         int hundredsIndex = abs(thisParentID) / 100;
0154         int thousandsIndex = abs(thisParentID) / 1000;
0155         if (((abs(thisParentID) >= 23) && (abs(thisParentID) <= 25)) || (abs(thisParentID) == 6) ||
0156             (hundredsIndex == 5) || (hundredsIndex == 4) || (thousandsIndex == 5) || (thousandsIndex == 4))
0157           parentID = thisParentID;
0158       }
0159       if ((parentID == -10000) && (nMo > 0))
0160         parentID = dynamic_cast<const reco::GenParticle*>(p.mother(0))->pdgId();
0161       //
0162       // If the parent of this particle is interesting, store all of the info
0163       if ((parentID != p.pdgId()) &&
0164           ((parentID > -9999) || (abs(id) == 11) || (abs(id) == 13) || (abs(id) == 23) || (abs(id) == 24) ||
0165            (abs(id) == 25) || (abs(id) == 4) || (abs(id) == 5) || (abs(id) == 6))) {
0166         l1GenData_->partId.push_back(p.pdgId());
0167         l1GenData_->partStat.push_back(p.status());
0168         l1GenData_->partPt.push_back(p.pt());
0169         l1GenData_->partEta.push_back(p.eta());
0170         l1GenData_->partPhi.push_back(p.phi());
0171         l1GenData_->partE.push_back(p.energy());
0172         l1GenData_->partParent.push_back(parentID);
0173         l1GenData_->partCh.push_back(p.charge());
0174         ++nPart;
0175       }
0176     }
0177     l1GenData_->nPart = nPart;
0178   }
0179 
0180   edm::Handle<std::vector<PileupSummaryInfo>> puInfoCollection;
0181   iEvent.getByToken(pileupInfoToken_, puInfoCollection);
0182 
0183   if (!puInfoCollection.isValid()) {
0184     throw cms::Exception("ProductNotValid") << "pileupInfoSource not valid";
0185   }
0186 
0187   // Loop over vector, find in-time entry, then store the relevant info
0188   std::vector<PileupSummaryInfo>::const_iterator puItr = puInfoCollection->begin();
0189   std::vector<PileupSummaryInfo>::const_iterator puEnd = puInfoCollection->end();
0190   for (; puItr != puEnd; ++puItr) {
0191     int bx = puItr->getBunchCrossing();
0192     if (bx == 0) {
0193       l1GenData_->nMeanPU = puItr->getTrueNumInteractions();
0194       l1GenData_->nVtx = puItr->getPU_NumInteractions();
0195       break;
0196     }
0197   }
0198 
0199   tree_->Fill();
0200 }
0201 
0202 // ------------ method called once each job just before starting event loop  ------------
0203 void L1GenTreeProducer::beginJob(void) {}
0204 
0205 // ------------ method called once each job just after ending the event loop  ------------
0206 void L1GenTreeProducer::endJob() {}
0207 
0208 //define this as a plug-in
0209 DEFINE_FWK_MODULE(L1GenTreeProducer);