ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MPIAnalyzer/src/MPIntuple.cc
Revision: 1.5
Committed: Sun May 23 02:33:41 2010 UTC (14 years, 11 months ago) by naodell
Content type: text/plain
Branch: MAIN
Changes since 1.4: +16 -11 lines
Log Message:
Add MC/Data discriminator. Changed some object names.

File Contents

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