CMS 3D CMS Logo

L1GenTreeProducer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: L1TriggerDPG/L1Ntuples
4 // Class: L1GenTreeProducer
5 //
13 //
14 // Original Author:
15 // Created:
16 // $Id: L1GenTreeProducer.cc,v 1.8 2012/08/29 12:44:03 jbrooke Exp $
17 //
18 //
19 
20 
21 // system include files
22 #include <memory>
23 
24 // framework
31 
32 // input data formats
34 
37 #include "HepMC/GenParticle.h"
38 #include "HepMC/GenVertex.h"
40 
42 
43 
44 // ROOT output stuff
47 #include "TTree.h"
48 
50 
51 //
52 // class declaration
53 //
54 
56 public:
57  explicit L1GenTreeProducer(const edm::ParameterSet&);
59 
60 
61 private:
62  virtual void beginJob(void) ;
63  virtual void analyze(const edm::Event&, const edm::EventSetup&);
64  virtual void endJob();
65 
66 private:
67 
68  unsigned maxL1Upgrade_;
69 
70  // output file
72 
73  // tree
74  TTree * tree_;
75 
76  // data format
77  std::unique_ptr<L1Analysis::L1AnalysisGeneratorDataFormat> l1GenData_;
78 
79  // EDM input tags
83 
84 };
85 
86 
87 
89 {
90 
91  genJetToken_ = consumes<reco::GenJetCollection>(iConfig.getUntrackedParameter<edm::InputTag>("genJetToken"));
92  genParticleToken_ = consumes<reco::GenParticleCollection>(iConfig.getUntrackedParameter<edm::InputTag>("genParticleToken"));
93  pileupInfoToken_ = consumes<std::vector<PileupSummaryInfo> >(iConfig.getUntrackedParameter<edm::InputTag>("pileupInfoToken"));
94 
95  l1GenData_ = std::make_unique<L1Analysis::L1AnalysisGeneratorDataFormat>();
96 
97  // set up output
98  tree_=fs_->make<TTree>("L1GenTree", "L1GenTree");
99  tree_->Branch("Generator", "L1Analysis::L1AnalysisGeneratorDataFormat", l1GenData_.get(), 32000, 3);
100 }
101 
102 
104 {
105 
106  // do anything here that needs to be done at desctruction time
107  // (e.g. close files, deallocate resources etc.)
108 
109 }
110 
111 
112 //
113 // member functions
114 //
115 
116 // ------------ method called to for each event ------------
117 void
119 {
120 
121  l1GenData_->Reset();
122 
124  iEvent.getByToken(genJetToken_, genJets);
125 
126 
127  if (genJets.isValid()){
128 
129  reco::GenJetCollection::const_iterator jetItr = genJets->begin();
130  reco::GenJetCollection::const_iterator jetEnd = genJets->end();
131  for( ; jetItr != jetEnd ; ++jetItr) {
132  l1GenData_->jetPt.push_back( jetItr->pt() );
133  l1GenData_->jetEta.push_back( jetItr->eta() );
134  l1GenData_->jetPhi.push_back( jetItr->phi() );
135  l1GenData_->nJet++;
136  }
137 
138  } else {
139  edm::LogWarning("MissingProduct") << "Gen jets not found. Branch will not be filled" << std::endl;
140  }
141 
143  iEvent.getByToken(genParticleToken_, genParticles);
144 
145  if (genParticles.isValid()) {
146 
147  int nPart {0};
148 
149  for(size_t i = 0; i < genParticles->size(); ++ i) {
150  const reco::GenParticle & p = (*genParticles)[i];
151  int id = p.pdgId();
152 
153  // See if the parent was interesting
154  int parentID = -10000;
155  unsigned int nMo=p.numberOfMothers();
156  for(unsigned int i=0;i<nMo;++i){
157  int thisParentID = dynamic_cast
158  <const reco::GenParticle*>(p.mother(i))->pdgId();
159  //
160  // Is this a bottom hadron?
161  int hundredsIndex = abs(thisParentID)/100;
162  int thousandsIndex = abs(thisParentID)/1000;
163  if ( ((abs(thisParentID) >= 23) &&
164  (abs(thisParentID) <= 25)) ||
165  (abs(thisParentID) == 6) ||
166  (hundredsIndex == 5) ||
167  (hundredsIndex == 4) ||
168  (thousandsIndex == 5) ||
169  (thousandsIndex == 4)
170  )
171  parentID = thisParentID;
172  }
173  if ((parentID == -10000) && (nMo > 0))
174  parentID = dynamic_cast
175  <const reco::GenParticle*>(p.mother(0))->pdgId();
176  //
177  // If the parent of this particle is interesting, store all of the info
178  if ((parentID != p.pdgId()) &&
179  ((parentID > -9999)
180  || (abs(id) == 11)
181  || (abs(id) == 13)
182  || (abs(id) == 23)
183  || (abs(id) == 24)
184  || (abs(id) == 25)
185  || (abs(id) == 4)
186  || (abs(id) == 5)
187  || (abs(id) == 6))
188  )
189  {
190  l1GenData_->partId.push_back(p.pdgId());
191  l1GenData_->partStat.push_back(p.status());
192  l1GenData_->partPt.push_back(p.pt());
193  l1GenData_->partEta.push_back(p.eta());
194  l1GenData_->partPhi.push_back(p.phi());
195  l1GenData_->partE.push_back(p.energy());
196  l1GenData_->partParent.push_back(parentID);
197  l1GenData_->partCh.push_back(p.charge());
198  ++nPart;
199  }
200  }
201  l1GenData_->nPart = nPart;
202  }
203 
204 
206  iEvent.getByToken(pileupInfoToken_, puInfoCollection);
207 
208  if (!puInfoCollection.isValid()) {
209  throw cms::Exception("ProductNotValid") << "pileupInfoSource not valid";
210  }
211 
212  // Loop over vector, find in-time entry, then store the relevant info
213  std::vector<PileupSummaryInfo>::const_iterator puItr = puInfoCollection->begin();
214  std::vector<PileupSummaryInfo>::const_iterator puEnd = puInfoCollection->end();
215  for( ; puItr != puEnd; ++puItr) {
216  int bx = puItr->getBunchCrossing();
217  if (bx == 0) {
218  l1GenData_->nMeanPU = puItr->getTrueNumInteractions();
219  l1GenData_->nVtx = puItr->getPU_NumInteractions();
220  break;
221  }
222  }
223 
224 
225 
226  tree_->Fill();
227 
228 }
229 
230 // ------------ method called once each job just before starting event loop ------------
231 void
233 {
234 }
235 
236 // ------------ method called once each job just after ending the event loop ------------
237 void
239 }
240 
241 //define this as a plug-in
virtual double pt() const final
transverse momentum
T getUntrackedParameter(std::string const &, T const &) const
virtual double eta() const final
momentum pseudorapidity
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:460
edm::EDGetTokenT< std::vector< PileupSummaryInfo > > pileupInfoToken_
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
virtual int status() const final
status word
virtual double phi() const final
momentum azimuthal angle
virtual int charge() const final
electric charge
Definition: LeafCandidate.h:91
edm::Service< TFileService > fs_
int iEvent
Definition: GenABIO.cc:230
virtual size_t numberOfMothers() const
number of mothers
edm::EDGetTokenT< reco::GenJetCollection > genJetToken_
edm::EDGetTokenT< reco::GenParticleCollection > genParticleToken_
virtual int pdgId() const final
PDG identifier.
virtual double energy() const final
energy
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
virtual void endJob()
bool isValid() const
Definition: HandleBase.h:74
TString nPart(Int_t part, TString string, TString delimit=";", Bool_t removerest=true)
std::unique_ptr< L1Analysis::L1AnalysisGeneratorDataFormat > l1GenData_
L1GenTreeProducer(const edm::ParameterSet &)
virtual void analyze(const edm::Event &, const edm::EventSetup &)
virtual void beginJob(void)
virtual const Candidate * mother(size_type=0) const
return mother at a given position, i = 0, ... numberOfMothers() - 1 (read only mode) ...