ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MPIAnalyzer/src/MPIntuple.cc
Revision: 1.11
Committed: Fri Jul 9 17:09:10 2010 UTC (14 years, 9 months ago) by devildog
Content type: text/plain
Branch: MAIN
Changes since 1.10: +22 -3 lines
Log Message:
Added some new Vertex variables, sumPtSquared function

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.11 // $Id: MPIntuple.cc,v 1.10 2010/07/08 09:36:20 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     #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.11 //placeholders
295     // vtxCon->SetNTracks(myVtx.tracksSize());
296     // vtxCon->SetTracksSumPtSquared(sumPtSquared(*myVtx));
297     int test1 = myVtx.tracksSize();
298     double test2 = sumPtSquared(myVtx);
299 naodell 1.1 ++Vtxcount;
300    
301     }
302    
303 andrey 1.2
304 naodell 1.7 if(triggerHLT_){
305 andrey 1.2
306 naodell 1.7 //---------- Filling HLT trigger bits! ------------
307 andrey 1.2
308 naodell 1.7 edm::Handle<TriggerResults> hltR;
309     triggerResultsTag_ = InputTag(hlTriggerResults_,"",hltName_);
310     iEvent.getByLabel(triggerResultsTag_,hltR);
311    
312     const TriggerNames & triggerNames = iEvent.triggerNames(*hltR);
313     hlNames=triggerNames.triggerNames();
314    
315 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"};
316 naodell 1.7
317     bool triggerPassed = false;
318     triggerStatus = 0x0;
319 naodell 1.4
320 naodell 1.7 for (uint i=0; i<hlNames.size(); ++i) {
321    
322     triggerPassed = triggerDecision(hltR, i);
323 naodell 1.4
324 naodell 1.7 if(triggerPassed){
325 andrey 1.3
326 naodell 1.7 for (uint j = 0; j != 32; ++j){
327    
328     if (hlNames[i] == MPI_TriggerNames[j])
329     {
330     // cout<<"trigger name: "<<hlNames[i]<<"list: "<<dec<<j+1<<endl;
331     triggerStatus |= 0x01 << j;
332    
333     }
334     }
335 naodell 1.4 }
336 naodell 1.7 }
337     }
338     // cout<< "total status: "<<hex<<triggerStatus<<endl;
339 naodell 1.4
340     if(PFcount > 0 || Gencount > 0 || Vtxcount > 0); sTree -> Fill();
341    
342     jet_ak5PF->Clear();
343 naodell 1.5 primaryVtx->Clear();
344 naodell 1.4 jetP4_ak5Gen->Clear();
345 andrey 1.2
346 naodell 1.1 }
347    
348    
349     // ------------ method called once each job just before starting event loop ------------
350     void MPIntuple::beginJob()
351     {
352    
353     ntupleFile = new TFile(rootfilename.c_str(), "RECREATE");
354 naodell 1.5 sTree = new TTree("mpiTree", "Tree for Jets");
355 naodell 1.1
356 naodell 1.4 jet_ak5PF = new TClonesArray("TCJet");
357 naodell 1.1 jetP4_ak5Gen = new TClonesArray("TLorentzVector");
358 naodell 1.10 genParticleP4 = new TClonesArray("TLorentzVector");
359 naodell 1.5 primaryVtx = new TClonesArray("TCPrimaryVtx");
360 naodell 1.1
361 naodell 1.4 sTree->Branch("jet_ak5PF",&jet_ak5PF, 6400, 0);
362 naodell 1.1 sTree->Branch("jetP4_ak5Gen",&jetP4_ak5Gen, 6400, 0);
363 naodell 1.10 sTree->Branch("genParticleP4",&genParticleP4, 6400, 0);
364 naodell 1.5 sTree->Branch("primaryVtx",&primaryVtx, 6400, 0);
365 naodell 1.1
366     sTree->Branch("eventNumber",&eventNumber, "eventNumber/I");
367     sTree->Branch("runNumber",&runNumber, "runNumber/I");
368     sTree->Branch("lumiSection",&lumiSection, "lumiSection/I");
369 naodell 1.6 sTree->Branch("triggerStatus",&triggerStatus, "triggerStatus/i");
370 naodell 1.8 sTree->Branch("isRealData",&isRealData, "isRealData/i");
371     sTree->Branch("bunchCross",&bunchCross, "bunchCross/i");
372 naodell 1.1
373     }
374    
375     // ------------ method called once each job just after ending the event loop ------------
376     void MPIntuple::endJob()
377     {
378     ntupleFile->Write();
379     ntupleFile->Close();
380    
381     }
382    
383    
384 andrey 1.2 bool MPIntuple::triggerDecision(edm::Handle<edm::TriggerResults> &hltR, int iTrigger)
385     {
386     bool triggerPassed = false;
387     if(hltR->wasrun(iTrigger) &&
388     hltR->accept(iTrigger) &&
389     !hltR->error(iTrigger) ){
390     triggerPassed = true;
391     }
392     return triggerPassed;
393     }
394    
395    
396 naodell 1.1 //define this as a plug-in
397     DEFINE_FWK_MODULE(MPIntuple);