ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MPIAnalyzer/src/MPIntuple.cc
Revision: 1.10
Committed: Thu Jul 8 09:36:20 2010 UTC (14 years, 9 months ago) by naodell
Content type: text/plain
Branch: MAIN
Changes since 1.9: +33 -12 lines
Log Message:
Added genParticle P4.

File Contents

# User Rev Content
1 naodell 1.1 // Original Author: "Nathaniel Odell"
2     // Secondary Author: Steven Won
3 andrey 1.2 // With contributions from: Andrey Pozdnyakov
4 naodell 1.1 // Created: Thurs April 22 2010
5 naodell 1.10 // $Id: MPIntuple.cc,v 1.9 2010/06/14 00:08:41 naodell Exp $
6 naodell 1.1
7     // system include files
8     #include <memory>
9     #include <string>
10    
11     // user include files
12     #include "FWCore/Framework/interface/Frameworkfwd.h"
13     #include "FWCore/Framework/interface/EDAnalyzer.h"
14     #include "FWCore/Framework/interface/ESHandle.h"
15     #include "FWCore/Framework/interface/EventSetup.h"
16     #include "FWCore/Framework/interface/Event.h"
17     #include "FWCore/Framework/interface/MakerMacros.h"
18     #include "FWCore/ParameterSet/interface/ParameterSet.h"
19     #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
20     #include "Geometry/Records/interface/CaloGeometryRecord.h"
21 naodell 1.10 #include "DataFormats/Math/interface/deltaPhi.h"
22 naodell 1.1
23 naodell 1.4 // Jet and vertex functions
24 naodell 1.1 #include "DataFormats/JetReco/interface/CaloJet.h"
25     #include "DataFormats/JetReco/interface/CaloJetCollection.h"
26     #include "DataFormats/JetReco/interface/PFJet.h"
27     #include "DataFormats/JetReco/interface/PFJetCollection.h"
28     #include "DataFormats/JetReco/interface/GenJet.h"
29     #include "DataFormats/JetReco/interface/GenJetCollection.h"
30     #include "DataFormats/JetReco/interface/Jet.h"
31     #include "DataFormats/VertexReco/interface/Vertex.h"
32     #include "DataFormats/VertexReco/interface/VertexFwd.h"
33 naodell 1.7 #include "DataFormats/BTauReco/interface/JetTag.h"
34    
35 naodell 1.10 //GenParticles
36     #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
37 naodell 1.1
38 naodell 1.4 // Need to get correctors
39     #include "JetMETCorrections/Objects/interface/JetCorrector.h"
40    
41     // ntuple storage classes
42     #include "TCJet.h"
43     #include "TCPrimaryVtx.h"
44    
45 andrey 1.2 // Need for HLT trigger info:
46     #include "FWCore/Common/interface/TriggerNames.h"
47     #include "DataFormats/Common/interface/TriggerResults.h"
48     #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
49    
50 naodell 1.4 //Root stuff
51 naodell 1.1 #include "TROOT.h"
52     #include "TFile.h"
53     #include "TTree.h"
54     #include "TString.h"
55     #include "TObject.h"
56     #include "TObjArray.h"
57     #include "TClonesArray.h"
58     #include "TRefArray.h"
59     #include "TLorentzVector.h"
60     #include "TVector3.h"
61    
62     using namespace edm;
63     using namespace std;
64     using namespace reco;
65    
66     //
67     // class declaration
68     //
69    
70     class MPIntuple : public edm::EDAnalyzer {
71     public:
72     explicit MPIntuple(const edm::ParameterSet&);
73     ~MPIntuple();
74    
75    
76     private:
77     virtual void beginJob() ;
78     virtual void analyze(const edm::Event&, const edm::EventSetup&);
79     virtual void endJob() ;
80    
81 andrey 1.2 bool triggerDecision(edm::Handle<edm::TriggerResults> &hltR, int iTrigger);
82    
83 naodell 1.1 // ----------member data ---------------------------
84    
85     edm::InputTag PFJetHandle_;
86     edm::InputTag GenJetHandle_;
87     edm::InputTag PrimaryVtxHandle_;
88 naodell 1.4 edm::InputTag triggerResultsTag_;
89 naodell 1.10 edm::InputTag genParticles_;
90 naodell 1.1
91     // Counting variables //
92    
93 naodell 1.8 int eventNumber, runNumber, lumiSection, bunchCross;
94 naodell 1.1
95 naodell 1.6 TTree* sTree;
96 naodell 1.1 TFile* ntupleFile;
97    
98 naodell 1.4 TClonesArray* jet_ak5PF;
99 naodell 1.1 TClonesArray* jetP4_ak5Gen;
100 naodell 1.10 TClonesArray* genParticleP4;
101 naodell 1.5 TClonesArray* primaryVtx;
102 naodell 1.1
103     bool doGenJets_;
104     bool doPFJets_;
105 naodell 1.7 bool triggerHLT_;
106 naodell 1.8 bool isRealData;
107 naodell 1.1 string rootfilename;
108    
109 andrey 1.2 //Triggers
110 naodell 1.4 string hlTriggerResults_, hltName_, triggerName_;
111     TriggerNames triggerNames;
112 andrey 1.2 HLTConfigProvider hltConfig_;
113 naodell 1.4 vector<string> hlNames;
114 naodell 1.5 unsigned int triggerStatus;
115 andrey 1.2
116    
117 naodell 1.1 };
118    
119     MPIntuple::MPIntuple(const edm::ParameterSet& iConfig):
120    
121     PFJetHandle_(iConfig.getUntrackedParameter<edm::InputTag>("PFJetTag")),
122     GenJetHandle_(iConfig.getUntrackedParameter<edm::InputTag>("GenJetTag")),
123     PrimaryVtxHandle_(iConfig.getUntrackedParameter<edm::InputTag>("PrimaryVtxTag")),
124 naodell 1.10 genParticles_(iConfig.getUntrackedParameter<edm::InputTag>("genParticleTag")),
125 naodell 1.1 doGenJets_(iConfig.getUntrackedParameter<bool>("doGenJets")),
126     doPFJets_(iConfig.getUntrackedParameter<bool>("doPFJets")),
127 naodell 1.7 triggerHLT_(iConfig.getUntrackedParameter<bool>("triggerHLT")),
128 andrey 1.2 rootfilename(iConfig.getUntrackedParameter<string>("rootfilename")),
129     hlTriggerResults_(iConfig.getUntrackedParameter<string>("HLTriggerResults","TriggerResults")),
130     hltName_(iConfig.getUntrackedParameter<string>("hltName","HLT"))
131 naodell 1.1 {
132 naodell 1.10
133 naodell 1.1
134     }
135    
136    
137     MPIntuple::~MPIntuple()
138     {
139    
140     }
141    
142     //
143     // member functions
144     //
145    
146     // ------------ method called to for each event ------------
147     void MPIntuple::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
148     {
149    
150 naodell 1.6 eventNumber = iEvent.id().event();
151     runNumber = iEvent.id().run();
152 naodell 1.1 lumiSection = (unsigned int)iEvent.getLuminosityBlock().luminosityBlock();
153 naodell 1.8 bunchCross = (unsigned int)iEvent.bunchCrossing();
154     isRealData = iEvent.isRealData();
155    
156 naodell 1.10 int PFcount = 0;
157     int Gencount = 0;
158     int GenPcount = 0;
159     int Vtxcount = 0;
160 naodell 1.1
161     if(doPFJets_){
162 naodell 1.4
163 naodell 1.7 edm::Handle<reco::JetTagCollection> bTagHandle1;
164     iEvent.getByLabel("trackCountingHighEffBJetTags", bTagHandle1);
165     const reco::JetTagCollection & bTags1 = *(bTagHandle1.product());
166     reco::JetTagCollection::const_iterator jet_it_1;
167    
168     edm::Handle<reco::JetTagCollection> bTagHandle2;
169     iEvent.getByLabel("trackCountingHighPurBJetTags", bTagHandle2);
170     const reco::JetTagCollection & bTags2 = *(bTagHandle2.product());
171     reco::JetTagCollection::const_iterator jet_it_2;
172    
173 naodell 1.4 const JetCorrector* correctorL2 = JetCorrector::getJetCorrector ("ak5PFL2Relative",iSetup);
174     const JetCorrector* correctorL3 = JetCorrector::getJetCorrector ("ak5PFL3Absolute",iSetup);
175    
176 naodell 1.1 Handle<reco::PFJetCollection> PFJets;
177     iEvent.getByLabel(PFJetHandle_, PFJets);
178    
179     PFcount = 0;
180    
181     for(PFJetCollection::const_iterator jet_iter = PFJets->begin(); jet_iter!= PFJets->end(); ++jet_iter){
182    
183     reco::PFJet myJet = reco::PFJet(*jet_iter);
184 naodell 1.4
185 naodell 1.7 // float scale2 = correctorL2->correction(myJet);
186     // float scale3 = correctorL3->correction(myJet);
187 naodell 1.4
188 naodell 1.1 if(myJet.et() > 5){
189 naodell 1.7
190     for(jet_it_1 = bTags1.begin(); jet_it_1 != bTags1.end(); jet_it_1++) {
191     if( (fabs(jet_it_1->first->eta() - myJet.eta()) < .005) && (deltaPhi(jet_it_1->first->phi(),myJet.phi()) < .005)) break;
192     }
193    
194     for(jet_it_2 = bTags2.begin(); jet_it_2 != bTags2.end(); jet_it_2++) {
195     if( (fabs(jet_it_2->first->eta() - myJet.eta()) < .005) && (deltaPhi(jet_it_2->first->phi(),myJet.phi()) < .005)) break;
196     }
197 naodell 1.6
198     TCJet* jetCon = new ((*jet_ak5PF)[PFcount]) TCJet;
199    
200     jetCon->SetP4(myJet.px(), myJet.py(), myJet.pz(), myJet.energy());
201     jetCon->SetVtx(myJet.vx(), myJet.vy(), myJet.vz());
202 naodell 1.4
203 naodell 1.6 jetCon->SetChHadFrac(myJet.chargedHadronEnergyFraction());
204     jetCon->SetNeuHadFrac(myJet.neutralHadronEnergyFraction());
205     jetCon->SetChEmFrac(myJet.chargedEmEnergyFraction());
206     jetCon->SetNeuEmFrac(myJet.neutralEmEnergyFraction());
207 naodell 1.4
208 naodell 1.6 jetCon->SetNumConstit(myJet.chargedMultiplicity() + myJet.neutralMultiplicity());
209     jetCon->SetNumChPart(myJet.chargedMultiplicity());
210 naodell 1.4
211 naodell 1.6 jetCon->SetNumChPart(myJet.chargedMultiplicity());
212 naodell 1.4
213 naodell 1.7 if(jet_it_2 != bTags2.end()) jetCon->SetBDiscrTrkCountHiPure(jet_it_2->second);
214     if(jet_it_1 != bTags1.end()) jetCon->SetBDiscrTrkCountHiEff(jet_it_1->second);
215    
216     float scale2 = correctorL2->correction(myJet);
217     myJet.scaleEnergy(scale2);
218     float scale3 = correctorL3->correction(myJet);
219     myJet.scaleEnergy(scale3);
220    
221     //more corrections?
222    
223 naodell 1.6 jetCon->SetJetCorr(2, scale2);
224     jetCon->SetJetCorr(3, scale3);
225 naodell 1.7
226 naodell 1.1 ++PFcount;
227 naodell 1.4 }
228     }
229 naodell 1.1 }
230    
231 naodell 1.10 if(!isRealData){
232 naodell 1.1
233 naodell 1.10 Handle<reco::GenParticleCollection > genpcoll;
234     iEvent.getByLabel(genParticles_, genpcoll);
235    
236 naodell 1.1 Handle<reco::GenJetCollection> GenJets;
237     iEvent.getByLabel(GenJetHandle_, GenJets);
238    
239    
240     for(GenJetCollection::const_iterator jet_iter = GenJets->begin(); jet_iter!= GenJets->end(); ++jet_iter){
241    
242     reco::GenJet myJet = reco::GenJet(*jet_iter);
243    
244 naodell 1.4 if(myJet.pt() > 5){
245 naodell 1.1
246     new ((*jetP4_ak5Gen)[Gencount]) TLorentzVector(myJet.px(), myJet.py(), myJet.pz(), myJet.energy());
247    
248     ++Gencount;
249    
250     }
251     }
252 naodell 1.10
253     for(reco::GenParticleCollection::const_iterator iGen = genpcoll->begin(); iGen!=genpcoll->end(); ++iGen) {
254    
255     reco::GenParticle p = reco::GenParticle(*iGen);
256    
257     new ((*genParticleP4)[GenPcount]) TLorentzVector(p.px(), p.py(), p.pz(), p.energy());
258    
259     ++GenPcount;
260    
261     }
262    
263    
264 naodell 1.1 }
265 naodell 1.10
266    
267 naodell 1.5 Handle<reco::VertexCollection> primaryVtcs;
268     iEvent.getByLabel(PrimaryVtxHandle_, primaryVtcs);
269 naodell 1.10
270 naodell 1.5 for(VertexCollection::const_iterator vtx_iter = primaryVtcs->begin(); vtx_iter!= primaryVtcs->end(); ++vtx_iter){
271 naodell 1.1
272     reco::Vertex myVtx = reco::Vertex(*vtx_iter);
273    
274 naodell 1.6 TCPrimaryVtx* vtxCon = new ((*primaryVtx)[Vtxcount]) TCPrimaryVtx;
275 naodell 1.4
276 naodell 1.6 vtxCon->SetPosition(myVtx.x(), myVtx.y(), myVtx.z());
277     vtxCon->SetNDof(myVtx.ndof());
278     vtxCon->SetChi2(myVtx.chi2());
279    
280 naodell 1.1 ++Vtxcount;
281    
282     }
283    
284 andrey 1.2
285 naodell 1.7 if(triggerHLT_){
286 andrey 1.2
287 naodell 1.7 //---------- Filling HLT trigger bits! ------------
288 andrey 1.2
289 naodell 1.7 edm::Handle<TriggerResults> hltR;
290     triggerResultsTag_ = InputTag(hlTriggerResults_,"",hltName_);
291     iEvent.getByLabel(triggerResultsTag_,hltR);
292    
293     const TriggerNames & triggerNames = iEvent.triggerNames(*hltR);
294     hlNames=triggerNames.triggerNames();
295    
296 naodell 1.9 string MPI_TriggerNames[] = {"HLT_PixelTracks_Multiplicity70", "HLT_MinBiasBSC_NoBPTX", "HLT_PixelTracks_Multiplicity40","HLT_L1Tech_HCAL_HF", "HLT_IsoTrackHB_8E29", "HLT_IsoTrackHE_8E29", "HLT_L1Tech_RPC_TTU_RBst1_collisions", "HLT_L1_BscMinBiasOR_BptxPlusORMinus", "HLT_L1Tech_BSC_halo_forPhysicsBackground", "HLT_L1Tech_BSC_HighMultiplicity", "HLT_MinBiasPixel_DoubleIsoTrack5", "HLT_MinBiasPixel_DoubleTrack", "HLT_MinBiasPixel_SingleTrack", "HLT_ZeroBiasPixel_SingleTrack", "HLT_MinBiasBSC", "HLT_StoppedHSCP_8E29", "HLT_Jet15U_HcalNoiseFiltered", "HLT_QuadJet15U", "HLT_DiJetAve30U_8E29", "HLT_DiJetAve15U_8E29", "HLT_FwdJet20U", "HLT_Jet50U", "HLT_Jet30U", "HLT_Jet15U", "HLT_BTagMu_Jet10U", "HLT_DoubleJet15U_ForwardBackward", "HLT_BTagIP_Jet50U", "HLT_DoubleLooseIsoTau15", "HLT_SingleLooseIsoTau20", "HLT_HT100U", "HLT_MET100", "HLT_MET45"};
297 naodell 1.7
298     bool triggerPassed = false;
299     triggerStatus = 0x0;
300 naodell 1.4
301 naodell 1.7 for (uint i=0; i<hlNames.size(); ++i) {
302    
303     triggerPassed = triggerDecision(hltR, i);
304 naodell 1.4
305 naodell 1.7 if(triggerPassed){
306 andrey 1.3
307 naodell 1.7 for (uint j = 0; j != 32; ++j){
308    
309     if (hlNames[i] == MPI_TriggerNames[j])
310     {
311     // cout<<"trigger name: "<<hlNames[i]<<"list: "<<dec<<j+1<<endl;
312     triggerStatus |= 0x01 << j;
313    
314     }
315     }
316 naodell 1.4 }
317 naodell 1.7 }
318     }
319     // cout<< "total status: "<<hex<<triggerStatus<<endl;
320 naodell 1.4
321     if(PFcount > 0 || Gencount > 0 || Vtxcount > 0); sTree -> Fill();
322    
323     jet_ak5PF->Clear();
324 naodell 1.5 primaryVtx->Clear();
325 naodell 1.4 jetP4_ak5Gen->Clear();
326 andrey 1.2
327 naodell 1.1 }
328    
329    
330     // ------------ method called once each job just before starting event loop ------------
331     void MPIntuple::beginJob()
332     {
333    
334     ntupleFile = new TFile(rootfilename.c_str(), "RECREATE");
335 naodell 1.5 sTree = new TTree("mpiTree", "Tree for Jets");
336 naodell 1.1
337 naodell 1.4 jet_ak5PF = new TClonesArray("TCJet");
338 naodell 1.1 jetP4_ak5Gen = new TClonesArray("TLorentzVector");
339 naodell 1.10 genParticleP4 = new TClonesArray("TLorentzVector");
340 naodell 1.5 primaryVtx = new TClonesArray("TCPrimaryVtx");
341 naodell 1.1
342 naodell 1.4 sTree->Branch("jet_ak5PF",&jet_ak5PF, 6400, 0);
343 naodell 1.1 sTree->Branch("jetP4_ak5Gen",&jetP4_ak5Gen, 6400, 0);
344 naodell 1.10 sTree->Branch("genParticleP4",&genParticleP4, 6400, 0);
345 naodell 1.5 sTree->Branch("primaryVtx",&primaryVtx, 6400, 0);
346 naodell 1.1
347     sTree->Branch("eventNumber",&eventNumber, "eventNumber/I");
348     sTree->Branch("runNumber",&runNumber, "runNumber/I");
349     sTree->Branch("lumiSection",&lumiSection, "lumiSection/I");
350 naodell 1.6 sTree->Branch("triggerStatus",&triggerStatus, "triggerStatus/i");
351 naodell 1.8 sTree->Branch("isRealData",&isRealData, "isRealData/i");
352     sTree->Branch("bunchCross",&bunchCross, "bunchCross/i");
353 naodell 1.1
354     }
355    
356     // ------------ method called once each job just after ending the event loop ------------
357     void MPIntuple::endJob()
358     {
359     ntupleFile->Write();
360     ntupleFile->Close();
361    
362     }
363    
364    
365 andrey 1.2 bool MPIntuple::triggerDecision(edm::Handle<edm::TriggerResults> &hltR, int iTrigger)
366     {
367     bool triggerPassed = false;
368     if(hltR->wasrun(iTrigger) &&
369     hltR->accept(iTrigger) &&
370     !hltR->error(iTrigger) ){
371     triggerPassed = true;
372     }
373     return triggerPassed;
374     }
375    
376    
377 naodell 1.1 //define this as a plug-in
378     DEFINE_FWK_MODULE(MPIntuple);