ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MPIAnalyzer/src/MPIntuple.cc
Revision: 1.12
Committed: Fri Jul 9 17:30:39 2010 UTC (14 years, 9 months ago) by devildog
Content type: text/plain
Branch: MAIN
Changes since 1.11: +3 -6 lines
Log Message:
Changed accessors to match new set functions in PrimaryVertex

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 devildog 1.12 // $Id: MPIntuple.cc,v 1.11 2010/07/09 17:09:10 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     #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.10 edm::InputTag genParticles_;
92 naodell 1.1
93     // Counting variables //
94    
95 naodell 1.8 int eventNumber, runNumber, lumiSection, bunchCross;
96 naodell 1.1
97 naodell 1.6 TTree* sTree;
98 naodell 1.1 TFile* ntupleFile;
99    
100 naodell 1.4 TClonesArray* jet_ak5PF;
101 naodell 1.1 TClonesArray* jetP4_ak5Gen;
102 naodell 1.10 TClonesArray* genParticleP4;
103 naodell 1.5 TClonesArray* primaryVtx;
104 naodell 1.1
105     bool doGenJets_;
106     bool doPFJets_;
107 naodell 1.7 bool triggerHLT_;
108 naodell 1.8 bool isRealData;
109 naodell 1.1 string rootfilename;
110    
111 andrey 1.2 //Triggers
112 naodell 1.4 string hlTriggerResults_, hltName_, triggerName_;
113     TriggerNames triggerNames;
114 andrey 1.2 HLTConfigProvider hltConfig_;
115 naodell 1.4 vector<string> hlNames;
116 naodell 1.5 unsigned int triggerStatus;
117 andrey 1.2
118    
119 naodell 1.1 };
120    
121     MPIntuple::MPIntuple(const edm::ParameterSet& iConfig):
122    
123     PFJetHandle_(iConfig.getUntrackedParameter<edm::InputTag>("PFJetTag")),
124     GenJetHandle_(iConfig.getUntrackedParameter<edm::InputTag>("GenJetTag")),
125     PrimaryVtxHandle_(iConfig.getUntrackedParameter<edm::InputTag>("PrimaryVtxTag")),
126 naodell 1.10 genParticles_(iConfig.getUntrackedParameter<edm::InputTag>("genParticleTag")),
127 naodell 1.1 doGenJets_(iConfig.getUntrackedParameter<bool>("doGenJets")),
128     doPFJets_(iConfig.getUntrackedParameter<bool>("doPFJets")),
129 naodell 1.7 triggerHLT_(iConfig.getUntrackedParameter<bool>("triggerHLT")),
130 andrey 1.2 rootfilename(iConfig.getUntrackedParameter<string>("rootfilename")),
131     hlTriggerResults_(iConfig.getUntrackedParameter<string>("HLTriggerResults","TriggerResults")),
132     hltName_(iConfig.getUntrackedParameter<string>("hltName","HLT"))
133 naodell 1.1 {
134 naodell 1.10
135 naodell 1.1
136     }
137    
138    
139     MPIntuple::~MPIntuple()
140     {
141    
142     }
143    
144 devildog 1.11 double MPIntuple::sumPtSquared(const Vertex & v)
145     {
146     double sum = 0.;
147     double pT;
148     for (Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
149     pT = (**it).pt();
150     double epT=(**it).ptError(); pT=pT>epT ? pT-epT : 0;
151    
152     sum += pT*pT;
153     }
154     return sum;
155     }
156    
157 naodell 1.1 //
158     // member functions
159     //
160    
161     // ------------ method called to for each event ------------
162     void MPIntuple::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
163     {
164    
165 naodell 1.6 eventNumber = iEvent.id().event();
166     runNumber = iEvent.id().run();
167 naodell 1.1 lumiSection = (unsigned int)iEvent.getLuminosityBlock().luminosityBlock();
168 naodell 1.8 bunchCross = (unsigned int)iEvent.bunchCrossing();
169     isRealData = iEvent.isRealData();
170    
171 naodell 1.10 int PFcount = 0;
172     int Gencount = 0;
173     int GenPcount = 0;
174     int Vtxcount = 0;
175 naodell 1.1
176     if(doPFJets_){
177 naodell 1.4
178 naodell 1.7 edm::Handle<reco::JetTagCollection> bTagHandle1;
179     iEvent.getByLabel("trackCountingHighEffBJetTags", bTagHandle1);
180     const reco::JetTagCollection & bTags1 = *(bTagHandle1.product());
181     reco::JetTagCollection::const_iterator jet_it_1;
182    
183     edm::Handle<reco::JetTagCollection> bTagHandle2;
184     iEvent.getByLabel("trackCountingHighPurBJetTags", bTagHandle2);
185     const reco::JetTagCollection & bTags2 = *(bTagHandle2.product());
186     reco::JetTagCollection::const_iterator jet_it_2;
187    
188 naodell 1.4 const JetCorrector* correctorL2 = JetCorrector::getJetCorrector ("ak5PFL2Relative",iSetup);
189     const JetCorrector* correctorL3 = JetCorrector::getJetCorrector ("ak5PFL3Absolute",iSetup);
190    
191 naodell 1.1 Handle<reco::PFJetCollection> PFJets;
192     iEvent.getByLabel(PFJetHandle_, PFJets);
193    
194     PFcount = 0;
195    
196     for(PFJetCollection::const_iterator jet_iter = PFJets->begin(); jet_iter!= PFJets->end(); ++jet_iter){
197    
198     reco::PFJet myJet = reco::PFJet(*jet_iter);
199 naodell 1.4
200 naodell 1.7 // float scale2 = correctorL2->correction(myJet);
201     // float scale3 = correctorL3->correction(myJet);
202 naodell 1.4
203 naodell 1.1 if(myJet.et() > 5){
204 naodell 1.7
205     for(jet_it_1 = bTags1.begin(); jet_it_1 != bTags1.end(); jet_it_1++) {
206     if( (fabs(jet_it_1->first->eta() - myJet.eta()) < .005) && (deltaPhi(jet_it_1->first->phi(),myJet.phi()) < .005)) break;
207     }
208    
209     for(jet_it_2 = bTags2.begin(); jet_it_2 != bTags2.end(); jet_it_2++) {
210     if( (fabs(jet_it_2->first->eta() - myJet.eta()) < .005) && (deltaPhi(jet_it_2->first->phi(),myJet.phi()) < .005)) break;
211     }
212 naodell 1.6
213     TCJet* jetCon = new ((*jet_ak5PF)[PFcount]) TCJet;
214    
215     jetCon->SetP4(myJet.px(), myJet.py(), myJet.pz(), myJet.energy());
216     jetCon->SetVtx(myJet.vx(), myJet.vy(), myJet.vz());
217 naodell 1.4
218 naodell 1.6 jetCon->SetChHadFrac(myJet.chargedHadronEnergyFraction());
219     jetCon->SetNeuHadFrac(myJet.neutralHadronEnergyFraction());
220     jetCon->SetChEmFrac(myJet.chargedEmEnergyFraction());
221     jetCon->SetNeuEmFrac(myJet.neutralEmEnergyFraction());
222 naodell 1.4
223 naodell 1.6 jetCon->SetNumConstit(myJet.chargedMultiplicity() + myJet.neutralMultiplicity());
224     jetCon->SetNumChPart(myJet.chargedMultiplicity());
225 naodell 1.4
226 naodell 1.6 jetCon->SetNumChPart(myJet.chargedMultiplicity());
227 naodell 1.4
228 naodell 1.7 if(jet_it_2 != bTags2.end()) jetCon->SetBDiscrTrkCountHiPure(jet_it_2->second);
229     if(jet_it_1 != bTags1.end()) jetCon->SetBDiscrTrkCountHiEff(jet_it_1->second);
230    
231     float scale2 = correctorL2->correction(myJet);
232     myJet.scaleEnergy(scale2);
233     float scale3 = correctorL3->correction(myJet);
234     myJet.scaleEnergy(scale3);
235    
236     //more corrections?
237    
238 naodell 1.6 jetCon->SetJetCorr(2, scale2);
239     jetCon->SetJetCorr(3, scale3);
240 naodell 1.7
241 naodell 1.1 ++PFcount;
242 naodell 1.4 }
243     }
244 naodell 1.1 }
245    
246 naodell 1.10 if(!isRealData){
247 naodell 1.1
248 naodell 1.10 Handle<reco::GenParticleCollection > genpcoll;
249     iEvent.getByLabel(genParticles_, genpcoll);
250    
251 naodell 1.1 Handle<reco::GenJetCollection> GenJets;
252     iEvent.getByLabel(GenJetHandle_, GenJets);
253    
254    
255     for(GenJetCollection::const_iterator jet_iter = GenJets->begin(); jet_iter!= GenJets->end(); ++jet_iter){
256    
257     reco::GenJet myJet = reco::GenJet(*jet_iter);
258    
259 naodell 1.4 if(myJet.pt() > 5){
260 naodell 1.1
261     new ((*jetP4_ak5Gen)[Gencount]) TLorentzVector(myJet.px(), myJet.py(), myJet.pz(), myJet.energy());
262    
263     ++Gencount;
264    
265     }
266     }
267 naodell 1.10
268     for(reco::GenParticleCollection::const_iterator iGen = genpcoll->begin(); iGen!=genpcoll->end(); ++iGen) {
269    
270     reco::GenParticle p = reco::GenParticle(*iGen);
271    
272     new ((*genParticleP4)[GenPcount]) TLorentzVector(p.px(), p.py(), p.pz(), p.energy());
273    
274     ++GenPcount;
275    
276     }
277    
278    
279 naodell 1.1 }
280 naodell 1.10
281    
282 naodell 1.5 Handle<reco::VertexCollection> primaryVtcs;
283     iEvent.getByLabel(PrimaryVtxHandle_, primaryVtcs);
284 naodell 1.10
285 naodell 1.5 for(VertexCollection::const_iterator vtx_iter = primaryVtcs->begin(); vtx_iter!= primaryVtcs->end(); ++vtx_iter){
286 naodell 1.1
287     reco::Vertex myVtx = reco::Vertex(*vtx_iter);
288 devildog 1.11 if(!myVtx.isValid() || myVtx.isFake()) continue;
289 naodell 1.6 TCPrimaryVtx* vtxCon = new ((*primaryVtx)[Vtxcount]) TCPrimaryVtx;
290 naodell 1.4
291 naodell 1.6 vtxCon->SetPosition(myVtx.x(), myVtx.y(), myVtx.z());
292     vtxCon->SetNDof(myVtx.ndof());
293     vtxCon->SetChi2(myVtx.chi2());
294 devildog 1.12 vtxCon->SetNTrks(myVtx.tracksSize());
295     vtxCon->SetSumPt2Trks(sumPtSquared(myVtx));
296 naodell 1.1 ++Vtxcount;
297    
298     }
299    
300 andrey 1.2
301 naodell 1.7 if(triggerHLT_){
302 andrey 1.2
303 naodell 1.7 //---------- Filling HLT trigger bits! ------------
304 andrey 1.2
305 naodell 1.7 edm::Handle<TriggerResults> hltR;
306     triggerResultsTag_ = InputTag(hlTriggerResults_,"",hltName_);
307     iEvent.getByLabel(triggerResultsTag_,hltR);
308    
309     const TriggerNames & triggerNames = iEvent.triggerNames(*hltR);
310     hlNames=triggerNames.triggerNames();
311    
312 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"};
313 naodell 1.7
314     bool triggerPassed = false;
315     triggerStatus = 0x0;
316 naodell 1.4
317 naodell 1.7 for (uint i=0; i<hlNames.size(); ++i) {
318    
319     triggerPassed = triggerDecision(hltR, i);
320 naodell 1.4
321 naodell 1.7 if(triggerPassed){
322 andrey 1.3
323 naodell 1.7 for (uint j = 0; j != 32; ++j){
324    
325     if (hlNames[i] == MPI_TriggerNames[j])
326     {
327     // cout<<"trigger name: "<<hlNames[i]<<"list: "<<dec<<j+1<<endl;
328     triggerStatus |= 0x01 << j;
329    
330     }
331     }
332 naodell 1.4 }
333 naodell 1.7 }
334     }
335     // cout<< "total status: "<<hex<<triggerStatus<<endl;
336 naodell 1.4
337     if(PFcount > 0 || Gencount > 0 || Vtxcount > 0); sTree -> Fill();
338    
339     jet_ak5PF->Clear();
340 naodell 1.5 primaryVtx->Clear();
341 naodell 1.4 jetP4_ak5Gen->Clear();
342 andrey 1.2
343 naodell 1.1 }
344    
345    
346     // ------------ method called once each job just before starting event loop ------------
347     void MPIntuple::beginJob()
348     {
349    
350     ntupleFile = new TFile(rootfilename.c_str(), "RECREATE");
351 naodell 1.5 sTree = new TTree("mpiTree", "Tree for Jets");
352 naodell 1.1
353 naodell 1.4 jet_ak5PF = new TClonesArray("TCJet");
354 naodell 1.1 jetP4_ak5Gen = new TClonesArray("TLorentzVector");
355 naodell 1.10 genParticleP4 = new TClonesArray("TLorentzVector");
356 naodell 1.5 primaryVtx = new TClonesArray("TCPrimaryVtx");
357 naodell 1.1
358 naodell 1.4 sTree->Branch("jet_ak5PF",&jet_ak5PF, 6400, 0);
359 naodell 1.1 sTree->Branch("jetP4_ak5Gen",&jetP4_ak5Gen, 6400, 0);
360 naodell 1.10 sTree->Branch("genParticleP4",&genParticleP4, 6400, 0);
361 naodell 1.5 sTree->Branch("primaryVtx",&primaryVtx, 6400, 0);
362 naodell 1.1
363     sTree->Branch("eventNumber",&eventNumber, "eventNumber/I");
364     sTree->Branch("runNumber",&runNumber, "runNumber/I");
365     sTree->Branch("lumiSection",&lumiSection, "lumiSection/I");
366 naodell 1.6 sTree->Branch("triggerStatus",&triggerStatus, "triggerStatus/i");
367 naodell 1.8 sTree->Branch("isRealData",&isRealData, "isRealData/i");
368     sTree->Branch("bunchCross",&bunchCross, "bunchCross/i");
369 naodell 1.1
370     }
371    
372     // ------------ method called once each job just after ending the event loop ------------
373     void MPIntuple::endJob()
374     {
375     ntupleFile->Write();
376     ntupleFile->Close();
377    
378     }
379    
380    
381 andrey 1.2 bool MPIntuple::triggerDecision(edm::Handle<edm::TriggerResults> &hltR, int iTrigger)
382     {
383     bool triggerPassed = false;
384     if(hltR->wasrun(iTrigger) &&
385     hltR->accept(iTrigger) &&
386     !hltR->error(iTrigger) ){
387     triggerPassed = true;
388     }
389     return triggerPassed;
390     }
391    
392    
393 naodell 1.1 //define this as a plug-in
394     DEFINE_FWK_MODULE(MPIntuple);