ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MPIAnalyzer/src/MPIntuple.cc
Revision: 1.14
Committed: Tue Oct 5 21:27:59 2010 UTC (14 years, 6 months ago) by naodell
Content type: text/plain
Branch: MAIN
Changes since 1.13: +6 -6 lines
Log Message:
DPS ntuple production

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