ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MPIAnalyzer/src/MPIntuple.cc
Revision: 1.15
Committed: Tue Oct 12 19:03:10 2010 UTC (14 years, 6 months ago) by naodell
Content type: text/plain
Branch: MAIN
Changes since 1.14: +56 -6 lines
Log Message:
Included jet-track associator.

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