ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MPIAnalyzer/src/MPIntuple.cc
Revision: 1.13
Committed: Fri Jul 23 12:45:32 2010 UTC (14 years, 9 months ago) by naodell
Content type: text/plain
Branch: MAIN
Changes since 1.12: +4 -25 lines
Log Message:
At least 4 jets in event.

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.13 // $Id: MPIntuple.cc,v 1.12 2010/07/09 17:30:39 devildog 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.6 eventNumber = iEvent.id().event();
162     runNumber = iEvent.id().run();
163 naodell 1.1 lumiSection = (unsigned int)iEvent.getLuminosityBlock().luminosityBlock();
164 naodell 1.8 bunchCross = (unsigned int)iEvent.bunchCrossing();
165     isRealData = iEvent.isRealData();
166    
167 naodell 1.10 int PFcount = 0;
168     int Gencount = 0;
169     int Vtxcount = 0;
170 naodell 1.1
171     if(doPFJets_){
172 naodell 1.4
173 naodell 1.7 edm::Handle<reco::JetTagCollection> bTagHandle1;
174     iEvent.getByLabel("trackCountingHighEffBJetTags", bTagHandle1);
175     const reco::JetTagCollection & bTags1 = *(bTagHandle1.product());
176     reco::JetTagCollection::const_iterator jet_it_1;
177    
178     edm::Handle<reco::JetTagCollection> bTagHandle2;
179     iEvent.getByLabel("trackCountingHighPurBJetTags", bTagHandle2);
180     const reco::JetTagCollection & bTags2 = *(bTagHandle2.product());
181     reco::JetTagCollection::const_iterator jet_it_2;
182    
183 naodell 1.4 const JetCorrector* correctorL2 = JetCorrector::getJetCorrector ("ak5PFL2Relative",iSetup);
184     const JetCorrector* correctorL3 = JetCorrector::getJetCorrector ("ak5PFL3Absolute",iSetup);
185    
186 naodell 1.1 Handle<reco::PFJetCollection> PFJets;
187     iEvent.getByLabel(PFJetHandle_, PFJets);
188    
189     PFcount = 0;
190    
191     for(PFJetCollection::const_iterator jet_iter = PFJets->begin(); jet_iter!= PFJets->end(); ++jet_iter){
192    
193     reco::PFJet myJet = reco::PFJet(*jet_iter);
194 naodell 1.4
195 naodell 1.7 // float scale2 = correctorL2->correction(myJet);
196     // float scale3 = correctorL3->correction(myJet);
197 naodell 1.4
198 naodell 1.1 if(myJet.et() > 5){
199 naodell 1.7
200     for(jet_it_1 = bTags1.begin(); jet_it_1 != bTags1.end(); jet_it_1++) {
201     if( (fabs(jet_it_1->first->eta() - myJet.eta()) < .005) && (deltaPhi(jet_it_1->first->phi(),myJet.phi()) < .005)) break;
202     }
203    
204     for(jet_it_2 = bTags2.begin(); jet_it_2 != bTags2.end(); jet_it_2++) {
205     if( (fabs(jet_it_2->first->eta() - myJet.eta()) < .005) && (deltaPhi(jet_it_2->first->phi(),myJet.phi()) < .005)) break;
206     }
207 naodell 1.6
208     TCJet* jetCon = new ((*jet_ak5PF)[PFcount]) TCJet;
209    
210     jetCon->SetP4(myJet.px(), myJet.py(), myJet.pz(), myJet.energy());
211     jetCon->SetVtx(myJet.vx(), myJet.vy(), myJet.vz());
212 naodell 1.4
213 naodell 1.6 jetCon->SetChHadFrac(myJet.chargedHadronEnergyFraction());
214     jetCon->SetNeuHadFrac(myJet.neutralHadronEnergyFraction());
215     jetCon->SetChEmFrac(myJet.chargedEmEnergyFraction());
216     jetCon->SetNeuEmFrac(myJet.neutralEmEnergyFraction());
217 naodell 1.4
218 naodell 1.6 jetCon->SetNumConstit(myJet.chargedMultiplicity() + myJet.neutralMultiplicity());
219     jetCon->SetNumChPart(myJet.chargedMultiplicity());
220 naodell 1.4
221 naodell 1.6 jetCon->SetNumChPart(myJet.chargedMultiplicity());
222 naodell 1.4
223 naodell 1.7 if(jet_it_2 != bTags2.end()) jetCon->SetBDiscrTrkCountHiPure(jet_it_2->second);
224     if(jet_it_1 != bTags1.end()) jetCon->SetBDiscrTrkCountHiEff(jet_it_1->second);
225    
226     float scale2 = correctorL2->correction(myJet);
227     myJet.scaleEnergy(scale2);
228     float scale3 = correctorL3->correction(myJet);
229     myJet.scaleEnergy(scale3);
230    
231     //more corrections?
232    
233 naodell 1.6 jetCon->SetJetCorr(2, scale2);
234     jetCon->SetJetCorr(3, scale3);
235 naodell 1.7
236 naodell 1.1 ++PFcount;
237 naodell 1.4 }
238     }
239 naodell 1.1 }
240    
241 naodell 1.10 if(!isRealData){
242 naodell 1.1
243 naodell 1.10
244 naodell 1.1 Handle<reco::GenJetCollection> GenJets;
245     iEvent.getByLabel(GenJetHandle_, GenJets);
246    
247    
248     for(GenJetCollection::const_iterator jet_iter = GenJets->begin(); jet_iter!= GenJets->end(); ++jet_iter){
249    
250     reco::GenJet myJet = reco::GenJet(*jet_iter);
251    
252 naodell 1.4 if(myJet.pt() > 5){
253 naodell 1.1
254     new ((*jetP4_ak5Gen)[Gencount]) TLorentzVector(myJet.px(), myJet.py(), myJet.pz(), myJet.energy());
255    
256     ++Gencount;
257    
258     }
259 naodell 1.13 }
260 naodell 1.1 }
261 naodell 1.10
262    
263 naodell 1.5 Handle<reco::VertexCollection> primaryVtcs;
264     iEvent.getByLabel(PrimaryVtxHandle_, primaryVtcs);
265 naodell 1.10
266 naodell 1.5 for(VertexCollection::const_iterator vtx_iter = primaryVtcs->begin(); vtx_iter!= primaryVtcs->end(); ++vtx_iter){
267 naodell 1.1
268     reco::Vertex myVtx = reco::Vertex(*vtx_iter);
269 devildog 1.11 if(!myVtx.isValid() || myVtx.isFake()) continue;
270 naodell 1.6 TCPrimaryVtx* vtxCon = new ((*primaryVtx)[Vtxcount]) TCPrimaryVtx;
271 naodell 1.4
272 naodell 1.6 vtxCon->SetPosition(myVtx.x(), myVtx.y(), myVtx.z());
273     vtxCon->SetNDof(myVtx.ndof());
274     vtxCon->SetChi2(myVtx.chi2());
275 devildog 1.12 vtxCon->SetNTrks(myVtx.tracksSize());
276     vtxCon->SetSumPt2Trks(sumPtSquared(myVtx));
277 naodell 1.1 ++Vtxcount;
278    
279     }
280    
281 andrey 1.2
282 naodell 1.7 if(triggerHLT_){
283 andrey 1.2
284 naodell 1.7 //---------- Filling HLT trigger bits! ------------
285 andrey 1.2
286 naodell 1.7 edm::Handle<TriggerResults> hltR;
287     triggerResultsTag_ = InputTag(hlTriggerResults_,"",hltName_);
288     iEvent.getByLabel(triggerResultsTag_,hltR);
289    
290     const TriggerNames & triggerNames = iEvent.triggerNames(*hltR);
291     hlNames=triggerNames.triggerNames();
292    
293 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"};
294 naodell 1.7
295     bool triggerPassed = false;
296     triggerStatus = 0x0;
297 naodell 1.4
298 naodell 1.7 for (uint i=0; i<hlNames.size(); ++i) {
299    
300     triggerPassed = triggerDecision(hltR, i);
301 naodell 1.4
302 naodell 1.7 if(triggerPassed){
303 andrey 1.3
304 naodell 1.7 for (uint j = 0; j != 32; ++j){
305    
306     if (hlNames[i] == MPI_TriggerNames[j])
307     {
308     // cout<<"trigger name: "<<hlNames[i]<<"list: "<<dec<<j+1<<endl;
309     triggerStatus |= 0x01 << j;
310    
311     }
312     }
313 naodell 1.4 }
314 naodell 1.7 }
315     }
316     // cout<< "total status: "<<hex<<triggerStatus<<endl;
317 naodell 1.4
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);