ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MPIAnalyzer/src/MPIntuple.cc
Revision: 1.7
Committed: Sat May 29 05:20:11 2010 UTC (14 years, 11 months ago) by naodell
Content type: text/plain
Branch: MAIN
Changes since 1.6: +84 -44 lines
Log Message:
Corrected jet correction retrieval method.

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.7 // $Id: MPIntuple.cc,v 1.6 2010/05/23 21:56:11 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     int eventNumber, runNumber, lumiSection;
92    
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.1
104     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.6 // isMC = iEvent.???();
151 naodell 1.1
152 naodell 1.6 int PFcount = 0;
153 naodell 1.1 int Gencount = 0;
154     int Vtxcount = 0;
155    
156     if(doPFJets_){
157 naodell 1.4
158 naodell 1.7 edm::Handle<reco::JetTagCollection> bTagHandle1;
159     iEvent.getByLabel("trackCountingHighEffBJetTags", bTagHandle1);
160     const reco::JetTagCollection & bTags1 = *(bTagHandle1.product());
161     reco::JetTagCollection::const_iterator jet_it_1;
162    
163     edm::Handle<reco::JetTagCollection> bTagHandle2;
164     iEvent.getByLabel("trackCountingHighPurBJetTags", bTagHandle2);
165     const reco::JetTagCollection & bTags2 = *(bTagHandle2.product());
166     reco::JetTagCollection::const_iterator jet_it_2;
167    
168 naodell 1.4 const JetCorrector* correctorL2 = JetCorrector::getJetCorrector ("ak5PFL2Relative",iSetup);
169     const JetCorrector* correctorL3 = JetCorrector::getJetCorrector ("ak5PFL3Absolute",iSetup);
170    
171 naodell 1.1 Handle<reco::PFJetCollection> PFJets;
172     iEvent.getByLabel(PFJetHandle_, PFJets);
173    
174     PFcount = 0;
175    
176     for(PFJetCollection::const_iterator jet_iter = PFJets->begin(); jet_iter!= PFJets->end(); ++jet_iter){
177    
178     reco::PFJet myJet = reco::PFJet(*jet_iter);
179 naodell 1.4
180 naodell 1.7 // float scale2 = correctorL2->correction(myJet);
181     // float scale3 = correctorL3->correction(myJet);
182 naodell 1.4
183 naodell 1.1 if(myJet.et() > 5){
184 naodell 1.7
185     for(jet_it_1 = bTags1.begin(); jet_it_1 != bTags1.end(); jet_it_1++) {
186     if( (fabs(jet_it_1->first->eta() - myJet.eta()) < .005) && (deltaPhi(jet_it_1->first->phi(),myJet.phi()) < .005)) break;
187     }
188    
189     for(jet_it_2 = bTags2.begin(); jet_it_2 != bTags2.end(); jet_it_2++) {
190     if( (fabs(jet_it_2->first->eta() - myJet.eta()) < .005) && (deltaPhi(jet_it_2->first->phi(),myJet.phi()) < .005)) break;
191     }
192 naodell 1.6
193     TCJet* jetCon = new ((*jet_ak5PF)[PFcount]) TCJet;
194    
195     jetCon->SetP4(myJet.px(), myJet.py(), myJet.pz(), myJet.energy());
196     jetCon->SetVtx(myJet.vx(), myJet.vy(), myJet.vz());
197 naodell 1.4
198 naodell 1.6 jetCon->SetChHadFrac(myJet.chargedHadronEnergyFraction());
199     jetCon->SetNeuHadFrac(myJet.neutralHadronEnergyFraction());
200     jetCon->SetChEmFrac(myJet.chargedEmEnergyFraction());
201     jetCon->SetNeuEmFrac(myJet.neutralEmEnergyFraction());
202 naodell 1.4
203 naodell 1.6 jetCon->SetNumConstit(myJet.chargedMultiplicity() + myJet.neutralMultiplicity());
204     jetCon->SetNumChPart(myJet.chargedMultiplicity());
205 naodell 1.4
206 naodell 1.6 jetCon->SetNumChPart(myJet.chargedMultiplicity());
207 naodell 1.4
208 naodell 1.7 if(jet_it_2 != bTags2.end()) jetCon->SetBDiscrTrkCountHiPure(jet_it_2->second);
209     if(jet_it_1 != bTags1.end()) jetCon->SetBDiscrTrkCountHiEff(jet_it_1->second);
210    
211     float scale2 = correctorL2->correction(myJet);
212     myJet.scaleEnergy(scale2);
213     float scale3 = correctorL3->correction(myJet);
214     myJet.scaleEnergy(scale3);
215    
216     //more corrections?
217    
218 naodell 1.6 jetCon->SetJetCorr(2, scale2);
219     jetCon->SetJetCorr(3, scale3);
220 naodell 1.7
221     //if(( (lumiSection==666690) && (eventNumber==2439)) ||
222     // ( (lumiSection==666702) && (eventNumber==3629)) ||
223     // ( (lumiSection==666702) && (eventNumber==3681)) ||
224     // ( (lumiSection==666709) && (eventNumber==4320)) ||
225     // ( (lumiSection==666709) && (eventNumber==4328)) ||
226     // ( (lumiSection==666724) && (eventNumber==5855)) ||
227     // ( (lumiSection==666742) && (eventNumber==7642)) ||
228     // ( (lumiSection==666754) && (eventNumber==8869)) ||
229     // ( (lumiSection==666670) && (eventNumber==466)) ||
230     // ( (lumiSection==666675) && (eventNumber==919)) ||
231     // ( (lumiSection==666718) && (eventNumber==5235)))
232     //{
233     // std::cout << runNumber << " " << lumiSection << " " << eventNumber << "\n";
234     // std::cout << "(" << myJet.pt() << "," << myJet.eta() << "," << myJet.phi() << ", [" << scale2 << ", " << scale3 << "]) " << std::endl;
235     // std::cout << "(" << jetCon->Pt(3) << "," << jetCon->P4(3).Eta() << "," << jetCon->P4(3).Phi() << "," << ", [" << jetCon->JetCorr(2) << ", " << jetCon->JetCorr(3) << "]) " << std::endl;
236     //}
237    
238 naodell 1.6
239 naodell 1.1 ++PFcount;
240 naodell 1.4 }
241     }
242 naodell 1.1 }
243    
244     if(doGenJets_){
245    
246     Handle<reco::GenJetCollection> GenJets;
247     iEvent.getByLabel(GenJetHandle_, GenJets);
248    
249    
250     for(GenJetCollection::const_iterator jet_iter = GenJets->begin(); jet_iter!= GenJets->end(); ++jet_iter){
251    
252     reco::GenJet myJet = reco::GenJet(*jet_iter);
253    
254 naodell 1.4 if(myJet.pt() > 5){
255 naodell 1.1
256     new ((*jetP4_ak5Gen)[Gencount]) TLorentzVector(myJet.px(), myJet.py(), myJet.pz(), myJet.energy());
257    
258     ++Gencount;
259    
260     }
261     }
262     }
263 naodell 1.4
264    
265 naodell 1.5 Handle<reco::VertexCollection> primaryVtcs;
266     iEvent.getByLabel(PrimaryVtxHandle_, primaryVtcs);
267 naodell 1.1
268 naodell 1.5 for(VertexCollection::const_iterator vtx_iter = primaryVtcs->begin(); vtx_iter!= primaryVtcs->end(); ++vtx_iter){
269 naodell 1.1
270     reco::Vertex myVtx = reco::Vertex(*vtx_iter);
271    
272 naodell 1.6 TCPrimaryVtx* vtxCon = new ((*primaryVtx)[Vtxcount]) TCPrimaryVtx;
273 naodell 1.4
274 naodell 1.6 vtxCon->SetPosition(myVtx.x(), myVtx.y(), myVtx.z());
275     vtxCon->SetNDof(myVtx.ndof());
276     vtxCon->SetChi2(myVtx.chi2());
277    
278 naodell 1.1 ++Vtxcount;
279    
280     }
281    
282 andrey 1.2
283 naodell 1.7 if(triggerHLT_){
284 andrey 1.2
285 naodell 1.7 //---------- Filling HLT trigger bits! ------------
286 andrey 1.2
287 naodell 1.7 edm::Handle<TriggerResults> hltR;
288     triggerResultsTag_ = InputTag(hlTriggerResults_,"",hltName_);
289     iEvent.getByLabel(triggerResultsTag_,hltR);
290    
291     const TriggerNames & triggerNames = iEvent.triggerNames(*hltR);
292     hlNames=triggerNames.triggerNames();
293    
294     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"};
295    
296     bool triggerPassed = false;
297     triggerStatus = 0x0;
298 naodell 1.4
299 naodell 1.7 for (uint i=0; i<hlNames.size(); ++i) {
300    
301     triggerPassed = triggerDecision(hltR, i);
302 naodell 1.4
303 naodell 1.7 if(triggerPassed){
304 andrey 1.3
305 naodell 1.7 for (uint j = 0; j != 32; ++j){
306    
307     if (hlNames[i] == MPI_TriggerNames[j])
308     {
309     // cout<<"trigger name: "<<hlNames[i]<<"list: "<<dec<<j+1<<endl;
310     triggerStatus |= 0x01 << j;
311    
312     }
313     }
314 naodell 1.4 }
315 naodell 1.7 }
316     }
317     // cout<< "total status: "<<hex<<triggerStatus<<endl;
318 naodell 1.4
319     if(PFcount > 0 || Gencount > 0 || Vtxcount > 0); sTree -> Fill();
320    
321     jet_ak5PF->Clear();
322 naodell 1.5 primaryVtx->Clear();
323 naodell 1.4 jetP4_ak5Gen->Clear();
324 andrey 1.2
325 naodell 1.1 }
326    
327    
328     // ------------ method called once each job just before starting event loop ------------
329     void MPIntuple::beginJob()
330     {
331    
332     ntupleFile = new TFile(rootfilename.c_str(), "RECREATE");
333 naodell 1.5 sTree = new TTree("mpiTree", "Tree for Jets");
334 naodell 1.1
335 naodell 1.4 jet_ak5PF = new TClonesArray("TCJet");
336 naodell 1.1 jetP4_ak5Gen = new TClonesArray("TLorentzVector");
337 naodell 1.5 primaryVtx = new TClonesArray("TCPrimaryVtx");
338 naodell 1.1
339 naodell 1.4 sTree->Branch("jet_ak5PF",&jet_ak5PF, 6400, 0);
340 naodell 1.1 sTree->Branch("jetP4_ak5Gen",&jetP4_ak5Gen, 6400, 0);
341 naodell 1.5 sTree->Branch("primaryVtx",&primaryVtx, 6400, 0);
342 naodell 1.1
343     sTree->Branch("eventNumber",&eventNumber, "eventNumber/I");
344     sTree->Branch("runNumber",&runNumber, "runNumber/I");
345     sTree->Branch("lumiSection",&lumiSection, "lumiSection/I");
346 naodell 1.6 sTree->Branch("triggerStatus",&triggerStatus, "triggerStatus/i");
347 naodell 1.7 // sTree->Branch("isMC",&isMC,"isMC/O");
348 naodell 1.5
349 naodell 1.1
350     }
351    
352     // ------------ method called once each job just after ending the event loop ------------
353     void MPIntuple::endJob()
354     {
355     ntupleFile->Write();
356     ntupleFile->Close();
357    
358     }
359    
360    
361 andrey 1.2 bool MPIntuple::triggerDecision(edm::Handle<edm::TriggerResults> &hltR, int iTrigger)
362     {
363     bool triggerPassed = false;
364     if(hltR->wasrun(iTrigger) &&
365     hltR->accept(iTrigger) &&
366     !hltR->error(iTrigger) ){
367     triggerPassed = true;
368     }
369     return triggerPassed;
370     }
371    
372    
373 naodell 1.1 //define this as a plug-in
374     DEFINE_FWK_MODULE(MPIntuple);