ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MPIAnalyzer/src/MPIntuple.cc
Revision: 1.9
Committed: Mon Jun 14 00:08:41 2010 UTC (14 years, 10 months ago) by naodell
Content type: text/plain
Branch: MAIN
Changes since 1.8: +2 -20 lines
Log Message:
""

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