ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/HbbAnalysis/plugins/HbbTreeMaker.cc
Revision: 1.41
Committed: Wed Aug 3 17:31:39 2011 UTC (13 years, 9 months ago) by amagnan
Content type: text/plain
Branch: MAIN
Changes since 1.40: +41 -16 lines
Log Message:
add parton jets

File Contents

# Content
1 #include "TLorentzVector.h"
2 #include <limits>
3
4 #include "DataFormats/Common/interface/Handle.h"
5 #include "DataFormats/Common/interface/TriggerResults.h"
6 #include "DataFormats/Common/interface/HLTenums.h"
7 #include "DataFormats/Common/interface/ValueMap.h"
8
9 #include "DataFormats/Candidate/interface/Candidate.h"
10 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
11 #include "DataFormats/Scalers/interface/DcsStatus.h"
12 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
13 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
14
15 #include "DataFormats/MuonReco/interface/Muon.h"
16 #include "DataFormats/MuonReco/interface/MuonFwd.h"
17 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
18 #include "DataFormats/MuonReco/interface/MuonChamberMatch.h"
19 #include "DataFormats/MuonReco/interface/MuonIsolation.h"
20 #include "DataFormats/MuonReco/interface/MuonEnergy.h"
21 #include "DataFormats/MuonReco/interface/MuonTime.h"
22 #include "DataFormats/MuonReco/interface/MuonQuality.h"
23
24 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
25 #include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
26 #include "DataFormats/RecoCandidate/interface/IsoDepositFwd.h"
27
28
29 #include "DataFormats/VertexReco/interface/Vertex.h"
30 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
31 #include "DataFormats/TrackReco/interface/Track.h"
32 #include "DataFormats/TrackReco/interface/TrackFwd.h"
33 #include "DataFormats/JetReco/interface/GenJet.h"
34
35 #include "DataFormats/PatCandidates/interface/Lepton.h"
36 #include "DataFormats/PatCandidates/interface/Muon.h"
37 #include "DataFormats/PatCandidates/interface/Electron.h"
38 #include "DataFormats/PatCandidates/interface/Tau.h"
39 #include "DataFormats/PatCandidates/interface/Jet.h"
40 #include "DataFormats/PatCandidates/interface/JetCorrFactors.h"
41 #include "DataFormats/PatCandidates/interface/MET.h"
42
43 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
44 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
45 #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
46
47 #include "FWCore/ServiceRegistry/interface/Service.h"
48 #include "FWCore/Framework/interface/Run.h"
49 #include "FWCore/Framework/interface/ESHandle.h"
50 #include "CommonTools/UtilAlgos/interface/TFileService.h"
51
52 #include "PhysicsTools/SelectorUtils/interface/SimpleCutBasedElectronIDSelectionFunctor.h"
53 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
54 #include "MagneticField/Engine/interface/MagneticField.h"
55
56 #include "UserCode/HbbAnalysis/interface/Electron.hh"
57 #include "UserCode/HbbAnalysis/interface/Muon.hh"
58 #include "UserCode/HbbAnalysis/interface/Tau.hh"
59 #include "UserCode/HbbAnalysis/interface/Jet.hh"
60 #include "UserCode/HbbAnalysis/interface/Met.hh"
61 #include "UserCode/HbbAnalysis/interface/Trigger.hh"
62 #include "UserCode/HbbAnalysis/interface/Vertex.hh"
63 #include "UserCode/HbbAnalysis/interface/L1Object.hh"
64 #include "UserCode/HbbAnalysis/interface/HLTObject.hh"
65
66 //AJG
67 #include "DataFormats/JetReco/interface/JPTJetCollection.h"
68 #include "DataFormats/JetReco/interface/JPTJet.h"
69 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
70 #include "SimDataFormats/Vertex/interface/SimVertex.h"
71
72 #include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h"
73 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
74
75 #include "UserCode/HbbAnalysis/plugins/HbbTreeMaker.hh"
76
77
78 HbbTreeMaker::HbbTreeMaker(const edm::ParameterSet & pset):
79 debug_(pset.getParameter<int>("DEBUG")),
80 processData_(false),
81 flavour_(pset.getParameter<unsigned int>("JetFlavour")),
82 doGen_(pset.getParameter<bool>("DOGEN")),
83 genParticleSrc_(pset.getParameter<edm::InputTag>("GenParticles")),
84 genJetSrc_(pset.getParameter<edm::InputTag>("GenJets")),
85 partonJetSrc_(pset.getParameter<edm::InputTag>("PartonJets")),
86 trackSrc_(pset.getParameter<edm::InputTag>("Tracks")),
87 dcsSrc_(pset.getParameter<edm::InputTag>("DCSInfo")),
88 electronSrc_(pset.getParameter<edm::InputTag>("Electrons")),
89 muonSrc_(pset.getParameter<edm::InputTag>("Muons")),
90 caloTauSrc_(pset.getParameter<edm::InputTag>("CaloTaus")),
91 pfTauSrc_(pset.getParameter<edm::InputTag>("PFTaus")),
92 caloJetSrc_(pset.getParameter<edm::InputTag>("CaloJets")),
93 jptJetSrc_(pset.getParameter<edm::InputTag>("JPTJets")),
94 ak7JetSrc_(pset.getParameter<edm::InputTag>("AK7Jets")),
95 pfJetSrc_(pset.getParameter<edm::InputTag>("PFJets")),
96 caloMetSrc_(pset.getParameter<edm::InputTag>("CaloMET")),
97 tcMetSrc_(pset.getParameter<edm::InputTag>("TCMET")),
98 pfMetSrc_(pset.getParameter<edm::InputTag>("PFMET")),
99 //pairSrc_(pset.getParameter<edm::InputTag>("Pair")),
100 //mmPairSrc_(pset.getParameter<edm::InputTag>("MuMuPair")),
101 //etPairSrc_(pset.getParameter<edm::InputTag>("ETauPair")),
102 //mtPairSrc_(pset.getParameter<edm::InputTag>("MuTauPair")),
103 vertexSrc_(pset.getParameter<edm::InputTag>("Vertex")),
104 triggerSrc_(pset.getParameter<edm::InputTag>("Trigger")),
105 hltPaths_(pset.getParameter<std::vector<std::string> >("HLTPaths")),
106 l1CenJetSrc_(pset.getParameter<edm::InputTag>("L1CenJet")),
107 l1TauJetSrc_(pset.getParameter<edm::InputTag>("L1TauJet")),
108 l1FwdJetSrc_(pset.getParameter<edm::InputTag>("L1FwdJet")),
109 hltSummarySrc_(pset.getParameter<edm::InputTag>("HLTSummary")),
110 event_(0),
111 tree_(0)
112 //processVec_(pset.getParameter<std::vector<std::string> >("ProcessVec"))
113 {//constructor
114
115 std::string JEC_PATH("CondFormats/JetMETObjects/data/");
116 try {
117 edm::FileInPath fipCalo(JEC_PATH+"Jec10V3_Uncertainty_AK5Calo.txt");
118 edm::FileInPath fipPF(JEC_PATH+"Jec10V3_Uncertainty_AK5PF.txt");
119 edm::FileInPath fipJPT(JEC_PATH+"Jec10V3_Uncertainty_AK5JPT.txt");
120 jecUncCalo = new JetCorrectionUncertainty(fipCalo.fullPath());
121 jecUncPF = new JetCorrectionUncertainty(fipPF.fullPath());
122 jecUncJPT = new JetCorrectionUncertainty(fipJPT.fullPath());
123 } catch(edm::Exception& e) {
124 std::cout << "JEC Uncertainty files unavailable! Exception : " << e.what() << ". " << std::endl;
125 jecUncCalo = 0;
126 jecUncPF = 0;
127 jecUncJPT = 0;
128 }
129
130 }//constructor
131
132 HbbTreeMaker::~HbbTreeMaker()
133 {//destructor
134 delete jecUncCalo;
135 delete jecUncPF;
136 delete jecUncJPT;
137
138 }//destructor
139
140
141
142 void HbbTreeMaker::beginJob(){//beginJob
143
144
145 edm::Service<TFileService> lFileService;
146
147 TFileDirectory lDir = lFileService->mkdir("HbbAnalysis");
148
149 tree_ = lDir.make<TTree>("Tree","Tree");
150 tree_->Branch("HbbEvent","HbbAnalysis::HbbEvent",&event_);
151 event_ = new HbbAnalysis::HbbEvent();
152
153 if (debug_) std::cout << "Initialising JetFlavour : " << std::endl;
154
155 lDir = lFileService->mkdir("JetFlavours");
156 jetFlav_.Initialise(lDir, debug_, flavour_);
157
158
159 }//beginJob
160
161 void HbbTreeMaker::endJob(){//endJob
162 if (!processData_) jetFlav_.printSummary();
163
164 //tree_->Write();
165 //delete tree_;
166 //delete event_;
167
168 }//endJob
169
170 void HbbTreeMaker::analyze(const edm::Event& aEvt, const edm::EventSetup& aEvtSetup){//analyze
171
172 //dumpContent(aEvt,"","MuTauPAT");
173 if (debug_) std::cout << "Processing event #" << aEvt.id().event() << std::endl;
174
175
176 static bool firstEvent = true;
177
178 event_->Clear();
179 event_->event(aEvt.id().event());
180 processData_ = aEvt.isRealData();
181
182 event_->run(aEvt.run());
183 event_->lumiBlock(aEvt.luminosityBlock());
184 event_->isRealData(aEvt.isRealData());
185 event_->bx(aEvt.bunchCrossing());
186
187 edm::Handle<reco::BeamSpot> lBeamSpot;
188 try {
189 aEvt.getByLabel(edm::InputTag("offlineBeamSpot"), lBeamSpot);
190 } catch(cms::Exception& e) {
191 std::cout<< "AMM: No beam spot found. Exception : " << e.what() << "." << std::endl;
192 }
193 HbbBeamSpot(lBeamSpot,event_->beamSpot());
194
195 //AJG - This will only work on RECO sample
196 /*
197 std::cout << "Beamspot = (" << lBeamSpot->x0() << "," << lBeamSpot->y0() << "," << lBeamSpot->z0() << ")" << std::endl;
198
199 edm::Handle<edm::SimVertexContainer> simVertexCollection;
200 aEvt.getByLabel("g4SimHits", simVertexCollection);
201 const edm::SimVertexContainer simVC = *(simVertexCollection.product());
202 std::cout << "SimVertexContainer size = " << simVC.size() << std::endl;
203 std::cout << "SimVertex[0] = (" << simVC.at(0).position().x() << ","
204 << simVC.at(0).position().y() << ","
205 << simVC.at(0).position().z() << ")" << std::endl;
206 */
207
208 //used in jet method....
209 edm::Handle<reco::GenParticleCollection> lGenParticles;
210
211 if (!processData_) {
212 edm::Handle<edm::HepMCProduct> mcProd;
213 edm::Handle<GenRunInfoProduct> lRunProd;
214 try {
215 aEvt.getByLabel("generator", mcProd );
216 } catch(cms::Exception& e) {
217 std::cout << "AMM: Collection HepMCProduct not available! Exception : " << e.what() << ". " << std::endl;
218 }
219 try {
220 aEvt.getRun().getByLabel("generator",lRunProd);
221 } catch(cms::Exception& e) {
222 std::cout << "AMM: Collection GenInfoProduct not available! Exception : " << e.what() << ". " << std::endl;
223 }
224
225 if (mcProd.isValid() && lRunProd.isValid()){
226 HbbGenInfo(mcProd,lRunProd);
227 }
228
229 edm::Handle<reco::GenJetCollection> lGenJets;
230 try {
231 aEvt.getByLabel(genJetSrc_,lGenJets);
232 if (debug_) std::cout << "** ngenJets = " << lGenJets->size() << std::endl;
233 } catch(cms::Exception& e) {
234 std::cout << "AMM: Collection GenJets not available! Exception : " << e.what() << ". " << std::endl;
235 }
236
237 //do not save components
238 if (lGenJets.isValid()) HbbGenJets(lGenJets,event_->genjets(),false);
239
240 edm::Handle<reco::GenJetCollection> lPartonJets;
241 try {
242 aEvt.getByLabel(partonJetSrc_,lPartonJets);
243 if (debug_) std::cout << "** npartonJets = " << lPartonJets->size() << std::endl;
244 } catch(cms::Exception& e) {
245 std::cout << "AMM: Collection PartonJets not available! Exception : " << e.what() << ". " << std::endl;
246 }
247
248 //save components
249 if (lPartonJets.isValid()) HbbGenJets(lPartonJets,event_->partonjets(),true);
250
251 edm::Handle<LHEEventProduct> lLHEEvt;
252 try {
253 aEvt.getByType(lLHEEvt);
254 } catch(cms::Exception& e) {
255 std::cout << "AMM: Collection of LHEEventProduct not available! Exception : " << e.what() << ". " << std::endl;
256 }
257
258 if (lLHEEvt.isValid()){
259 HbbLHEInfo(lLHEEvt,event_->lheParticles());
260 }
261
262 edm::Handle<std::vector<PileupSummaryInfo> > lPUInfo;
263 try {
264 aEvt.getByLabel("addPileupInfo",lPUInfo);
265 } catch(cms::Exception& e) {
266 std::cout << "AMM: Collection of PileupSummaryInfo not available! Exception : " << e.what() << ". " << std::endl;
267 }
268
269 if (lPUInfo.isValid()){
270 HbbPUInfo(lPUInfo,event_->puVars());
271 }
272
273 try {
274 aEvt.getByLabel(genParticleSrc_,lGenParticles);
275 if (debug_) std::cout << "** ngenParticles = " << lGenParticles->size() << std::endl;
276 } catch(cms::Exception& e) {
277 std::cout << "AMM: Collection genParticles not available! Exception : " << e.what() << ". " << std::endl;
278 }
279
280 if (doGen_) HbbParticles(lGenParticles,event_->particles());
281
282 unsigned int lNPartons = 0;
283 lNPartons = jetFlav_.fillPartons(lGenParticles);
284
285 if (debug_) std::cout << "--- Number of partons = " << lNPartons << "." << std::endl;
286 }
287
288 edm::Handle<std::vector<reco::Vertex> > lRecoVertices;
289 try {
290 aEvt.getByLabel(vertexSrc_, lRecoVertices);
291 if (!lRecoVertices.isValid()){
292 edm::LogInfo("ERROR")<< "Error! Can't get vertex by label. ";
293 }
294 if (debug_) std::cout << "** vertexcollection = " << lRecoVertices->size() << " elements." << std::endl;
295 } catch(cms::Exception& e) {
296 std::cout << "AMM: Collection " << vertexSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
297 }
298
299 HbbVertices(lRecoVertices,event_->vertices());
300
301 edm::Handle<std::vector<pat::Electron> > lElectronCollection;
302
303 try {
304 aEvt.getByLabel(electronSrc_,lElectronCollection);
305 if (!lElectronCollection.isValid()){
306 edm::LogInfo("ERROR")<< "Error! Can't get electron by label. ";
307 }
308 if (debug_) std::cout << "** electroncollection = " << lElectronCollection->size() << " elements." << std::endl;
309 } catch(cms::Exception& e) {
310 std::cout << "AMM: Collection " << electronSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
311 }
312
313 HbbElectrons(lElectronCollection,event_->electrons());
314
315 edm::Handle<std::vector<pat::Muon> > lMuonCollection;
316
317 try {
318 aEvt.getByLabel(muonSrc_,lMuonCollection);
319 if (!lMuonCollection.isValid()){
320 edm::LogInfo("ERROR")<< "Error! Can't get muon by label. ";
321 }
322 if (debug_) std::cout << "** muoncollection = " << lMuonCollection->size() << " elements." << std::endl;
323 } catch(cms::Exception& e) {
324 std::cout << "AMM: Collection " << muonSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
325 }
326
327 HbbMuons(lMuonCollection,lRecoVertices,lBeamSpot,event_->muons());
328
329
330 if (!( caloTauSrc_.label()=="" && caloTauSrc_.instance()=="" )) {
331 edm::Handle<std::vector<pat::Tau> > lTauCollection;
332
333 try {
334 aEvt.getByLabel(caloTauSrc_,lTauCollection);
335 if (!lTauCollection.isValid()){
336 edm::LogInfo("ERROR")<< "Error! Can't get caloTau by label. ";
337 }
338 if (debug_) std::cout << "** caloTaucollection = " << lTauCollection->size() << " elements." << std::endl;
339 } catch(cms::Exception& e) {
340 std::cout << "AMM: Collection " << caloTauSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
341 }
342
343 HbbTaus(lTauCollection,lRecoVertices,event_->caloTaus());
344 }
345
346 edm::Handle<std::vector<pat::Tau> > lPFTauCollection;
347
348 try {
349 aEvt.getByLabel(pfTauSrc_,lPFTauCollection);
350 if (!lPFTauCollection.isValid()){
351 edm::LogInfo("ERROR")<< "Error! Can't get pfTau by label. ";
352 }
353 if (debug_) std::cout << "** pfTaucollection = " << lPFTauCollection->size() << " elements." << std::endl;
354 } catch(cms::Exception& e) {
355 std::cout << "AMM: Collection " << pfTauSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
356 }
357
358 HbbTaus(lPFTauCollection,lRecoVertices,event_->pfTaus());
359
360 edm::Handle<std::vector<pat::Jet> > lCaloJetCollection;
361
362 try {
363 aEvt.getByLabel(caloJetSrc_,lCaloJetCollection);
364 if (!lCaloJetCollection.isValid()){
365 edm::LogInfo("ERROR")<< "Error! Can't get caloJets by label. ";
366 }
367 if (debug_) std::cout << "** caloJetcollection = " << lCaloJetCollection->size() << " elements." << std::endl;
368 } catch(cms::Exception& e) {
369 std::cout << "AMM: Collection " << caloJetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
370 }
371
372 HbbJets(lCaloJetCollection,jetFlav_,lGenParticles,event_->caloJets());
373
374 //AJG - Process reco JPT Jet Corrections:
375 //This section should be uncommented for testing purposes
376
377 /*
378 edm::Handle<reco::JPTJetCollection > lJptRL2L3RESJetCollection;
379 edm::Handle<reco::JPTJetCollection > lJptRL1L2L3RESJetCollection;
380
381 edm::Handle<reco::JPTJetCollection > lJptRRAWJetCollection;
382 try {
383 aEvt.getByLabel("JetPlusTrackZSPCorJetAntiKt5",lJptRRAWJetCollection);
384 if (!lJptRRAWJetCollection.isValid()){
385 edm::LogInfo("ERROR")<< "Error! Can't get JetPlusZSPCorJetAntiKt5 by label. ";
386 }
387 if (debug_) std::cout << "** jetPlusTrackZSPCorJetAntiKt5 = " << lJptRRAWJetCollection->size() << " elements." << std::endl;
388 } catch(cms::Exception& e) {
389 std::cout << "AJG: Collection not available! Exception : " << e.what() << ". " << std::endl;
390 }
391
392 try {
393 aEvt.getByLabel("ak5JPTJetsL2L3Residual",lJptRL2L3RESJetCollection);
394 aEvt.getByLabel("JetPlusTrackZSPCorJetAntiKt5",lJptRRAWJetCollection);
395 aEvt.getByLabel("ak5JPTJetsL1L2L3Residual",lJptRL1L2L3RESJetCollection);
396 if (!lJptRRAWJetCollection.isValid()){
397 edm::LogInfo("ERROR")<< "Error! Can't get jptRecoJets by label. ";
398 }
399 if (debug_) std::cout << "** jptRecoJetcollection = " << lJptRRAWJetCollection->size() << " elements." << std::endl;
400 } catch(cms::Exception& e) {
401 std::cout << "AJG: Collection " << jptJetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
402 }
403
404
405 for(reco::JPTJetCollection::const_iterator jptjet = lJptRRAWJetCollection->begin(); jptjet != lJptRRAWJetCollection->end(); ++jptjet ) {
406
407 std::cout << "pT = " << jptjet->pt() << std::endl;
408
409 reco::TrackRefVector tpionsInVInC = jptjet->getPionsInVertexInCalo();
410 std::cout << "Pion Track Count (reco) = " << tpionsInVInC.size() << std::endl;
411 }
412
413 reco::JPTJetCollection::const_iterator jptRAWjet = lJptRRAWJetCollection->begin();
414 reco::JPTJetCollection::const_iterator jptL2L3RESjet = lJptRL2L3RESJetCollection->begin();
415 reco::JPTJetCollection::const_iterator jptL1L2L3RESjet = lJptRL1L2L3RESJetCollection->begin();
416 std::cout << "pT RAW = " << jptRAWjet->pt() << std::endl;
417 std::cout << "pT L2L3RES = " << jptL2L3RESjet->pt() << std::endl;
418 std::cout << "pT L1L2L3RES = " << jptL1L2L3RESjet->pt() << std::endl;
419 */
420
421
422
423 edm::Handle<std::vector<pat::Jet> > lJptJetCollection;
424
425 try {
426 aEvt.getByLabel(jptJetSrc_,lJptJetCollection);
427 if (!lJptJetCollection.isValid()){
428 edm::LogInfo("ERROR")<< "Error! Can't get jptJets by label. ";
429 }
430 if (debug_) std::cout << "** jptJetcollection = " << lJptJetCollection->size() << " elements." << std::endl;
431 } catch(cms::Exception& e) {
432 std::cout << "AMM: Collection " << jptJetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
433 }
434
435 //std::cout << "Processing JPT jets:" << std::endl;
436 HbbJets(lJptJetCollection,jetFlav_,lGenParticles,event_->jptJets());
437
438 edm::Handle<std::vector<pat::Jet> > lAk7JetCollection;
439
440 try {
441 aEvt.getByLabel(ak7JetSrc_,lAk7JetCollection);
442 if (!lAk7JetCollection.isValid()){
443 edm::LogInfo("ERROR")<< "Error! Can't get ak7Jets by label. ";
444 }
445 if (debug_) std::cout << "** ak7Jetcollection = " << lAk7JetCollection->size() << " elements." << std::endl;
446 } catch(cms::Exception& e) {
447 std::cout << "AMM: Collection " << ak7JetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
448 }
449
450 //std::cout << "Processing JPT jets:" << std::endl;
451 HbbJets(lAk7JetCollection,jetFlav_,lGenParticles,event_->ak7Jets());
452
453 edm::Handle<std::vector<pat::Jet> > lPfJetCollection;
454
455 try {
456 aEvt.getByLabel(pfJetSrc_,lPfJetCollection);
457 if (!lPfJetCollection.isValid()){
458 edm::LogInfo("ERROR")<< "Error! Can't get pfJets by label. ";
459 }
460 if (debug_) std::cout << "** pfJetcollection = " << lPfJetCollection->size() << " elements." << std::endl;
461 } catch(cms::Exception& e) {
462 std::cout << "AMM: Collection " << pfJetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
463 }
464
465 //std::cout << "Processing PF jets:" << std::endl;
466 HbbJets(lPfJetCollection,jetFlav_,lGenParticles,event_->pfJets());
467
468 edm::Handle<std::vector<pat::MET> > lCaloMetCol;
469
470 try {
471 aEvt.getByLabel(caloMetSrc_,lCaloMetCol);
472 if (!lCaloMetCol.isValid()){
473 edm::LogInfo("ERROR")<< "Error! Can't get caloMet by label. ";
474 }
475 if (debug_) std::cout << "** caloMet = " << lCaloMetCol->size() << " elements." << std::endl;
476 } catch(cms::Exception& e) {
477 std::cout << "AMM: Collection " << caloMetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
478 }
479
480 HbbMet(lCaloMetCol,event_->caloMet());
481
482 edm::Handle<std::vector<pat::MET> > lTcMetCol;
483
484 try {
485 aEvt.getByLabel(tcMetSrc_,lTcMetCol);
486 if (!lTcMetCol.isValid()){
487 edm::LogInfo("ERROR")<< "Error! Can't get tcMet by label. ";
488 }
489 if (debug_) std::cout << "** tcMet = " << lTcMetCol->size() << " elements." << std::endl;
490 } catch(cms::Exception& e) {
491 std::cout << "AMM: Collection " << tcMetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
492 }
493
494 HbbMet(lTcMetCol,event_->tcMet());
495
496 edm::Handle<std::vector<pat::MET> > lPfMetCol;
497
498 try {
499 aEvt.getByLabel(pfMetSrc_,lPfMetCol);
500 if (!lPfMetCol.isValid()){
501 edm::LogInfo("ERROR")<< "Error! Can't get pfMet by label. ";
502 }
503 if (debug_) std::cout << "** pfMet = " << lPfMetCol->size() << " elements." << std::endl;
504 } catch(cms::Exception& e) {
505 std::cout << "AMM: Collection " << pfMetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
506 }
507
508 HbbMet(lPfMetCol,event_->pfMet());
509
510
511
512 //triggers
513 edm::Handle<edm::TriggerResults> lTrigCol;
514 try {
515 aEvt.getByLabel(triggerSrc_,lTrigCol);
516 if (!lTrigCol.isValid()){
517 edm::LogInfo("ERROR")<< "Error! Can't get triggers by label. ";
518 }
519 if (debug_) std::cout << "** triggers = " << lTrigCol->size() << std::endl;
520 } catch(cms::Exception& e) {
521 std::cout << "AMM: Collection " << triggerSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
522 }
523 if (lTrigCol.isValid()) {
524 const edm::TriggerNames & lNames = aEvt.triggerNames(*lTrigCol);
525 HbbTrigger(lTrigCol,lNames,event_->triggers());
526 }
527
528 //L1 jet objects
529
530 edm::Handle<l1extra::L1JetParticleCollection> lCenJet;
531 edm::Handle<l1extra::L1JetParticleCollection> lFwdJet;
532 edm::Handle<l1extra::L1JetParticleCollection> lTauJet;
533 try {
534 aEvt.getByLabel(l1CenJetSrc_,lCenJet);
535 aEvt.getByLabel(l1TauJetSrc_,lTauJet);
536 aEvt.getByLabel(l1FwdJetSrc_,lFwdJet);
537 } catch(cms::Exception& e) {
538 std::cout << "AMM: L1 Collections " << l1CenJetSrc_ << ", " << l1TauJetSrc_ << ", " << l1FwdJetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
539 }
540 if (lCenJet.isValid()) {
541 HbbL1Objects(lCenJet,event_->l1Jets(),1);
542 }
543 if (lTauJet.isValid()) {
544 HbbL1Objects(lTauJet,event_->l1Jets(),2);
545 }
546 if (lFwdJet.isValid()) {
547 HbbL1Objects(lFwdJet,event_->l1Jets(),3);
548 }
549
550
551 edm::Handle<trigger::TriggerEvent> lHltSummary;
552 try {
553 aEvt.getByLabel(hltSummarySrc_,lHltSummary);
554 } catch(cms::Exception& e) {
555 std::cout << "AMM: collection " << hltSummarySrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
556 }
557 if (lHltSummary.isValid()) {
558 HbbHLTObjects(lHltSummary,event_->hltObjects());
559 }
560
561
562
563 //event_->order();
564 tree_->Fill();
565 firstEvent = false;
566
567 }//analyze
568
569
570 void HbbTreeMaker::HbbBeamSpot(const edm::Handle<reco::BeamSpot> & aBeamSpot,
571 HbbAnalysis::BeamSpot & aBS)
572 {
573
574 if (aBeamSpot.isValid()) {
575 HbbAnalysis::BeamSpotVars lbs;
576 lbs.x0 = aBeamSpot->x0();
577 lbs.y0 = aBeamSpot->y0();
578 lbs.z0 = aBeamSpot->z0();
579 lbs.sigmaZ = aBeamSpot->sigmaZ();
580 lbs.BeamWidthX = aBeamSpot->BeamWidthX();
581 lbs.BeamWidthY = aBeamSpot->BeamWidthY();
582 lbs.x0Error = aBeamSpot->x0Error();
583 lbs.y0Error = aBeamSpot->y0Error();
584 lbs.z0Error = aBeamSpot->z0Error();
585 lbs.sigmaZ0Error = aBeamSpot->sigmaZ0Error();
586 lbs.BeamWidthXError = aBeamSpot->BeamWidthXError();
587 lbs.BeamWidthYError = aBeamSpot->BeamWidthYError();
588
589 aBS.bsVars(lbs);
590 }
591
592 }
593
594
595 void HbbTreeMaker::HbbElectrons(const edm::Handle<std::vector<pat::Electron> > & aCol,
596 std::vector<HbbAnalysis::Electron> & aVec)
597 {//HbbElectrons
598
599 if (aCol.isValid()){
600 if (aCol->size() > 0) {
601
602 unsigned int iEle = 0;
603 for (std::vector<pat::Electron>::const_iterator iter = aCol->begin();
604 iter != aCol->end();
605 ++iter,++iEle)
606 {
607 //if (debug_ > 1) std::cout << "**** Ele #" << iEle << ", pT,eta,phi = " << (*iter).pt() << " " << (*iter).eta() << " " << (*iter).phi() << std::endl;
608 // get the track collection
609
610
611
612 // define your function for WPXX using relative or absolute isolations
613 SimpleCutBasedElectronIDSelectionFunctor patSeleR95(SimpleCutBasedElectronIDSelectionFunctor::relIso95);
614 SimpleCutBasedElectronIDSelectionFunctor patSeleR90(SimpleCutBasedElectronIDSelectionFunctor::relIso90);
615 SimpleCutBasedElectronIDSelectionFunctor patSeleR85(SimpleCutBasedElectronIDSelectionFunctor::relIso85);
616 SimpleCutBasedElectronIDSelectionFunctor patSeleR80(SimpleCutBasedElectronIDSelectionFunctor::relIso80);
617 SimpleCutBasedElectronIDSelectionFunctor patSeleR70(SimpleCutBasedElectronIDSelectionFunctor::relIso70);
618 SimpleCutBasedElectronIDSelectionFunctor patSeleR60(SimpleCutBasedElectronIDSelectionFunctor::relIso60);
619 SimpleCutBasedElectronIDSelectionFunctor patSeleC95(SimpleCutBasedElectronIDSelectionFunctor::cIso95);
620 SimpleCutBasedElectronIDSelectionFunctor patSeleC90(SimpleCutBasedElectronIDSelectionFunctor::cIso90);
621 SimpleCutBasedElectronIDSelectionFunctor patSeleC85(SimpleCutBasedElectronIDSelectionFunctor::cIso85);
622 SimpleCutBasedElectronIDSelectionFunctor patSeleC80(SimpleCutBasedElectronIDSelectionFunctor::cIso80);
623 SimpleCutBasedElectronIDSelectionFunctor patSeleC70(SimpleCutBasedElectronIDSelectionFunctor::cIso70);
624 SimpleCutBasedElectronIDSelectionFunctor patSeleC60(SimpleCutBasedElectronIDSelectionFunctor::cIso60);
625
626
627 HbbAnalysis::GenVars lGen;
628 if (!processData_ && (*iter).genLepton()){
629 lGen.valid = true;
630 lGen.E = (*iter).genLepton()->energy();
631 lGen.pT = (*iter).genLepton()->pt();
632 lGen.eta = (*iter).genLepton()->eta();
633 lGen.y = (*iter).genLepton()->y();
634 lGen.phi = (*iter).genLepton()->phi();
635 lGen.charge = (*iter).genLepton()->charge();
636 lGen.pdgId = (*iter).genLepton()->pdgId();
637 lGen.status = (*iter).genLepton()->status();
638 lGen.mass = (*iter).genLepton()->mass();
639 lGen.vx = (*iter).genLepton()->vx();
640 lGen.vy = (*iter).genLepton()->vy();
641 lGen.vz = (*iter).genLepton()->vz();
642 }
643 else {
644 lGen.valid = false;
645 lGen.E = 0;
646 lGen.pT = 0;
647 lGen.eta = 0;
648 lGen.y = 0;
649 lGen.phi = 0;
650 lGen.charge = 0;
651 lGen.pdgId = 0;
652 lGen.status = 0;
653 lGen.mass = 0;
654 lGen.vx = 0;
655 lGen.vy = 0;
656 lGen.vz = 0;
657 }
658
659 HbbAnalysis::BaseVars lReco;
660 lReco.E = (*iter).energy();
661 lReco.pT = (*iter).pt();
662 lReco.eta = (*iter).eta();
663 lReco.y = (*iter).y();
664 lReco.phi = (*iter).phi();
665 lReco.charge = (*iter).charge();
666 lReco.vx = (*iter).vx();
667 lReco.vy = (*iter).vy();
668 lReco.vz = (*iter).vz();
669 lReco.hltMatch = (*iter).triggerObjectMatches().size();
670
671 HbbAnalysis::SCVars lSC;
672 lSC.sigmaEtaEta = (*iter).scSigmaEtaEta();
673 lSC.sigmaIEtaIEta = (*iter).scSigmaIEtaIEta();
674 lSC.e1x5 = (*iter).scE1x5();
675 lSC.e2x5Max = (*iter).scE2x5Max();
676 lSC.e5x5 = (*iter).scE5x5();
677 lSC.eOverP = (*iter).eSuperClusterOverP();
678
679 HbbAnalysis::EleIsoVars lIso;
680 lIso.calo = (*iter).caloIso();
681 lIso.track = (*iter).trackIso();
682 lIso.ecal = (*iter).ecalIso();
683 lIso.hcal = (*iter).hcalIso();
684
685 HbbAnalysis::EleIDVars lID;
686 lID.idAndIso =
687 (patSeleR95(*iter)<<11) | (patSeleR90(*iter)<<10) |
688 (patSeleR85(*iter)<<9) | (patSeleR80(*iter)<<8) |
689 (patSeleR70(*iter)<<7) | (patSeleR60(*iter)<<6) |
690 (patSeleC95(*iter)<<5) | (patSeleC90(*iter)<<4) |
691 (patSeleC85(*iter)<<3) | (patSeleC80(*iter)<<2) |
692 (patSeleC70(*iter)<<1) | (patSeleC60(*iter));
693 lID.electronIDs = (*iter).electronIDs();
694 lID.hOverE = (*iter).hadronicOverEm();
695 lID.deltaPhiIn = (*iter).deltaPhiSuperClusterTrackAtVtx();
696 lID.deltaEtaIn = (*iter).deltaEtaSuperClusterTrackAtVtx();
697 lID.ecalDrivenSeed = (*iter).ecalDrivenSeed();
698 lID.trackerDrivenSeed = (*iter).trackerDrivenSeed();
699 lID.dB = (*iter).dB();
700
701 HbbAnalysis::Electron lObj(lGen,lReco,lSC,lIso,lID);
702 aVec.push_back(lObj);
703 }
704 }
705 }
706
707 }//HbbElectrons
708
709
710 void HbbTreeMaker::HbbMuons(const edm::Handle<std::vector<pat::Muon> > & aCol,
711 const edm::Handle<std::vector<reco::Vertex> > & aRecoVertices,
712 const edm::Handle<reco::BeamSpot> & aBeamSpot,
713 std::vector<HbbAnalysis::Muon> & aVec)
714 {//HbbMuons
715
716 if (aCol.isValid()){
717 if (aCol->size() > 0) {
718
719 unsigned int iEle = 0;
720 for (std::vector<pat::Muon>::const_iterator iter = aCol->begin();
721 iter != aCol->end();
722 ++iter,++iEle)
723 {
724 const reco::Muon* recoMuon = dynamic_cast<const reco::Muon*>((*iter).originalObject());
725
726 //select only global and tracker muons
727 if (!recoMuon->isGlobalMuon() || !recoMuon->isTrackerMuon()) continue;
728
729 HbbAnalysis::GenVars lGen;
730 if (!processData_ && (*iter).genLepton()){
731 lGen.valid = true;
732 lGen.E = (*iter).genLepton()->energy();
733 lGen.pT = (*iter).genLepton()->pt();
734 lGen.eta = (*iter).genLepton()->eta();
735 lGen.y = (*iter).genLepton()->y();
736 lGen.phi = (*iter).genLepton()->phi();
737 lGen.charge = (*iter).genLepton()->charge();
738 lGen.pdgId = (*iter).genLepton()->pdgId();
739 lGen.status = (*iter).genLepton()->status();
740 lGen.mass = (*iter).genLepton()->mass();
741 lGen.vx = (*iter).genLepton()->vx();
742 lGen.vy = (*iter).genLepton()->vy();
743 lGen.vz = (*iter).genLepton()->vz();
744 }
745 else {
746 lGen.valid = false;
747 lGen.E = 0;
748 lGen.pT = 0;
749 lGen.eta = 0;
750 lGen.y = 0;
751 lGen.phi = 0;
752 lGen.charge = 0;
753 lGen.pdgId = 0;
754 lGen.status = 0;
755 lGen.mass = 0;
756 lGen.vx = 0;
757 lGen.vy = 0;
758 lGen.vz = 0;
759 }
760
761 HbbAnalysis::BaseVars lReco;
762 lReco.E = (*iter).energy();
763 lReco.pT = (*iter).pt();
764 lReco.eta = (*iter).eta();
765 lReco.y = (*iter).y();
766 lReco.phi = (*iter).phi();
767 lReco.charge = (*iter).charge();
768 lReco.vx = (*iter).vx();
769 lReco.vy = (*iter).vy();
770 lReco.vz = (*iter).vz();
771 lReco.hltMatch = (*iter).triggerObjectMatches().size();
772
773 HbbAnalysis::MuTrkVars lTrk;
774 reco::TrackRef glmuon = recoMuon->globalTrack();
775 reco::TrackRef tkmuon = recoMuon->innerTrack();
776
777 if ( glmuon.isAvailable() && glmuon.isNonnull() ) {
778 if (aBeamSpot.isValid()) {
779 lTrk.dxy = glmuon->dxy(aBeamSpot->position());
780 lTrk.dz = glmuon->dz(aBeamSpot->position());
781 }
782 else {
783 lTrk.dxy = 0;
784 lTrk.dz = 0;
785 }
786 lTrk.normalizedChi2 = glmuon->normalizedChi2();
787 lTrk.muonHits = glmuon->hitPattern().numberOfValidMuonHits();
788 lTrk.charge = glmuon->charge();
789 }
790 else {
791 lTrk.dxy = 0;
792 lTrk.dz = 0;
793 lTrk.normalizedChi2 = 0;
794 lTrk.muonHits = 0;
795 lTrk.charge = 0;
796 }
797 if ( tkmuon.isAvailable() && tkmuon.isNonnull() ) {
798 lTrk.trackerHits = tkmuon->hitPattern().numberOfValidTrackerHits();
799 lTrk.pixelHits = tkmuon->hitPattern().numberOfValidPixelHits();
800 }
801 else {
802 lTrk.trackerHits = 0;
803 lTrk.pixelHits = 0;
804 }
805
806 HbbAnalysis::MuIsoVars lIsoR03;
807 lIsoR03.sumPt = recoMuon->isolationR03().sumPt;
808 lIsoR03.emEt = recoMuon->isolationR03().emEt;
809 lIsoR03.hadEt = recoMuon->isolationR03().hadEt;
810 lIsoR03.nTracks = recoMuon->isolationR03().nTracks;
811 lIsoR03.nJets = recoMuon->isolationR03().nJets;
812
813 HbbAnalysis::MuIsoVars lIsoR05;
814 lIsoR05.sumPt = recoMuon->isolationR05().sumPt;
815 lIsoR05.emEt = recoMuon->isolationR05().emEt;
816 lIsoR05.hadEt = recoMuon->isolationR05().hadEt;
817 lIsoR05.nTracks = recoMuon->isolationR05().nTracks;
818 lIsoR05.nJets = recoMuon->isolationR05().nJets;
819
820 HbbAnalysis::MuIDVars lID;
821 //if (recoMuon->isGlobalMuon()) lID.type = 1;
822 //else if (recoMuon->isTrackerMuon()) lID.type = 2;
823 //else if (recoMuon->isStandAloneMuon()) lID.type = 3;
824 //else if (recoMuon->isCaloMuon()) lID.type = 4;
825 //else lID.type = 0;
826 lID.type = (recoMuon->isCaloMuon()<<3) | (recoMuon->isStandAloneMuon()<<2) | (recoMuon->isTrackerMuon()<<1) | recoMuon->isGlobalMuon();
827
828 if (muon::isGoodMuon(*recoMuon,muon::AllGlobalMuons)) lID.ids.push_back(1);
829 if (muon::isGoodMuon(*recoMuon,muon::AllStandAloneMuons)) lID.ids.push_back(2);
830 if (muon::isGoodMuon(*recoMuon,muon::AllTrackerMuons)) lID.ids.push_back(3);
831 if (muon::isGoodMuon(*recoMuon,muon::TrackerMuonArbitrated)) lID.ids.push_back(4);
832 if (muon::isGoodMuon(*recoMuon,muon::AllArbitrated)) lID.ids.push_back(5);
833 if (muon::isGoodMuon(*recoMuon,muon::GlobalMuonPromptTight)) lID.ids.push_back(6);
834 if (muon::isGoodMuon(*recoMuon,muon::TMLastStationLoose)) lID.ids.push_back(7);
835 if (muon::isGoodMuon(*recoMuon,muon::TMLastStationTight)) lID.ids.push_back(8);
836 if (muon::isGoodMuon(*recoMuon,muon::TM2DCompatibilityLoose)) lID.ids.push_back(9);
837 if (muon::isGoodMuon(*recoMuon,muon::TM2DCompatibilityTight)) lID.ids.push_back(10);
838 if (muon::isGoodMuon(*recoMuon,muon::TMOneStationLoose)) lID.ids.push_back(11);
839 if (muon::isGoodMuon(*recoMuon,muon::TMOneStationTight)) lID.ids.push_back(12);
840 if (muon::isGoodMuon(*recoMuon,muon::TMLastStationOptimizedLowPtLoose)) lID.ids.push_back(13);
841 if (muon::isGoodMuon(*recoMuon,muon::TMLastStationOptimizedLowPtTight)) lID.ids.push_back(14);
842 if (muon::isGoodMuon(*recoMuon,muon::GMTkChiCompatibility)) lID.ids.push_back(15);
843 if (muon::isGoodMuon(*recoMuon,muon::GMStaChiCompatibility)) lID.ids.push_back(16);
844 if (muon::isGoodMuon(*recoMuon,muon::GMTkKinkTight)) lID.ids.push_back(17);
845 if (muon::isGoodMuon(*recoMuon,muon::TMLastStationAngLoose)) lID.ids.push_back(18);
846 if (muon::isGoodMuon(*recoMuon,muon::TMLastStationAngTight)) lID.ids.push_back(19);
847 if (muon::isGoodMuon(*recoMuon,muon::TMOneStationAngLoose)) lID.ids.push_back(20);
848 if (muon::isGoodMuon(*recoMuon,muon::TMOneStationAngTight)) lID.ids.push_back(21);
849 if (muon::isGoodMuon(*recoMuon,muon::TMLastStationOptimizedBarrelLowPtLoose)) lID.ids.push_back(22);
850 if (muon::isGoodMuon(*recoMuon,muon::TMLastStationOptimizedBarrelLowPtTight)) lID.ids.push_back(23);
851
852 lID.caloCompat = recoMuon->caloCompatibility();
853 lID.segCompat = muon::segmentCompatibility(*recoMuon);
854 lID.nChambers = recoMuon->numberOfChambers();
855 lID.nMatchesLoose = recoMuon->numberOfMatches(reco::Muon::NoArbitration);
856 lID.nMatchesMedium = recoMuon->numberOfMatches(reco::Muon::SegmentArbitration);
857 lID.nMatchesTight = recoMuon->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
858 lID.dB = (*iter).dB();
859
860 HbbAnalysis::Muon lObj(lGen,lReco,lTrk,lIsoR03,lIsoR05,lID);
861 aVec.push_back(lObj);
862 }
863 }
864 }
865
866 }//HbbMuons
867
868 void HbbTreeMaker::HbbTaus(const edm::Handle<std::vector<pat::Tau> > & aCol,
869 const edm::Handle<std::vector<reco::Vertex> > & aRecoVertices,
870 std::vector<HbbAnalysis::Tau> & aVec)
871 {//HbbTaus
872
873 if (aCol.isValid()){
874 if (aCol->size() > 0) {
875
876 unsigned int iEle = 0;
877 for (std::vector<pat::Tau>::const_iterator iter = aCol->begin();
878 iter != aCol->end();
879 ++iter,++iEle)
880 {
881 //if (debug_ > 1) std::cout << "**** Ele #" << iEle << ", pT,eta,phi = " << (*iter).pt() << " " << (*iter).eta() << " " << (*iter).phi() << std::endl;
882
883 HbbAnalysis::GenVars lGen;
884 if (!processData_ && (*iter).genLepton()){
885 lGen.valid = true;
886 lGen.E = (*iter).genLepton()->energy();
887 lGen.pT = (*iter).genLepton()->pt();
888 lGen.eta = (*iter).genLepton()->eta();
889 lGen.y = (*iter).genLepton()->y();
890 lGen.phi = (*iter).genLepton()->phi();
891 lGen.charge = (*iter).genLepton()->charge();
892 lGen.pdgId = (*iter).genLepton()->pdgId();
893 lGen.status = (*iter).genLepton()->status();
894 lGen.mass = (*iter).genLepton()->mass();
895 lGen.vx = (*iter).genLepton()->vx();
896 lGen.vy = (*iter).genLepton()->vy();
897 lGen.vz = (*iter).genLepton()->vz();
898 }
899 else {
900 lGen.valid = false;
901 lGen.E = 0;
902 lGen.pT = 0;
903 lGen.eta = 0;
904 lGen.y = 0;
905 lGen.phi = 0;
906 lGen.charge = 0;
907 lGen.pdgId = 0;
908 lGen.status = 0;
909 lGen.mass = 0;
910 lGen.vx = 0;
911 lGen.vy = 0;
912 lGen.vz = 0;
913 }
914
915
916 HbbAnalysis::GenVars lGenJet;
917 if (!processData_ && (*iter).genJet()){
918 lGenJet.valid = true;
919 lGenJet.E = (*iter).genJet()->energy();
920 lGenJet.pT = (*iter).genJet()->pt();
921 lGenJet.eta = (*iter).genJet()->eta();
922 lGenJet.y = (*iter).genJet()->y();
923 lGenJet.phi = (*iter).genJet()->phi();
924 lGenJet.charge = (*iter).genJet()->charge();
925 lGenJet.pdgId = (*iter).genJet()->pdgId();
926 lGenJet.status = (*iter).genJet()->status();
927 lGenJet.mass = (*iter).genJet()->mass();
928 lGenJet.vx = (*iter).genJet()->vx();
929 lGenJet.vy = (*iter).genJet()->vy();
930 lGenJet.vz = (*iter).genJet()->vz();
931 }
932 else {
933 lGenJet.valid = false;
934 lGenJet.E = 0;
935 lGenJet.pT = 0;
936 lGenJet.eta = 0;
937 lGenJet.y = 0;
938 lGenJet.phi = 0;
939 lGenJet.charge = 0;
940 lGenJet.pdgId = 0;
941 lGenJet.status = 0;
942 lGenJet.mass = 0;
943 lGenJet.vx = 0;
944 lGenJet.vy = 0;
945 lGenJet.vz = 0;
946 }
947
948 HbbAnalysis::BaseVars lReco;
949 lReco.E = (*iter).energy();
950 lReco.pT = (*iter).pt();
951 lReco.eta = (*iter).eta();
952 lReco.y = (*iter).y();
953 lReco.phi = (*iter).phi();
954 lReco.charge = (*iter).charge();
955 lReco.vx = (*iter).vx();
956 lReco.vy = (*iter).vy();
957 lReco.vz = (*iter).vz();
958
959 bool isCalo = (*iter).isCaloTau();
960 bool isPF = (*iter).isPFTau();
961
962 assert (isCalo == !isPF);
963
964 HbbAnalysis::TrkVars lLead;
965 if (isCalo) {
966 lLead.signedSipt = (*iter).leadTracksignedSipt();
967
968 reco::TrackRef lCaloTrk = (*iter).leadTrack();
969
970 if ( lCaloTrk.isAvailable() && lCaloTrk.isNonnull() ) {
971 lLead.pT = lCaloTrk->pt();
972 lLead.eta = lCaloTrk->eta();
973 lLead.phi = lCaloTrk->phi();
974
975 lLead.matchDist = reco::deltaR(lCaloTrk->momentum(), (*iter).p4());
976
977 if ( aRecoVertices->size() >= 1 ) {
978 const reco::Vertex& thePrimaryEventVertex = (*aRecoVertices->begin());
979 lLead.IPxy = lCaloTrk->dxy(thePrimaryEventVertex.position());
980 lLead.IPz = lCaloTrk->dz(thePrimaryEventVertex.position());
981 }
982 else {
983 lLead.IPxy = 0;
984 lLead.IPz = 0;
985 }
986 }
987 else {
988 lLead.pT = 0;
989 lLead.eta = 0;
990 lLead.phi = 0;
991 lLead.matchDist = 0;
992 lLead.IPxy = 0;
993 lLead.IPz = 0;
994 lLead.signedSipt = 0;
995 }
996 }//calo
997 else {
998 lLead.signedSipt = (*iter).leadPFChargedHadrCandsignedSipt();
999
1000 reco::PFCandidateRef lHadrCand = (*iter).leadPFChargedHadrCand();
1001
1002 if ( lHadrCand.isAvailable() && lHadrCand.isNonnull() ) {
1003 lLead.pT = lHadrCand->pt();
1004 lLead.eta = lHadrCand->eta();
1005 lLead.phi = lHadrCand->phi();
1006
1007 lLead.matchDist = reco::deltaR(lHadrCand->p4(), (*iter).p4());
1008
1009 reco::TrackRef lTrk = lHadrCand->trackRef();
1010
1011 if ( aRecoVertices->size() >= 1) {
1012 const reco::Vertex& thePrimaryEventVertex = (*aRecoVertices->begin());
1013 lLead.IPxy = lTrk->dxy(thePrimaryEventVertex.position());
1014 lLead.IPz = lTrk->dz(thePrimaryEventVertex.position());
1015 }
1016 else {
1017 lLead.IPxy = 0;
1018 lLead.IPz = 0;
1019 }
1020 }
1021 else {
1022 lLead.pT = 0;
1023 lLead.eta = 0;
1024 lLead.phi = 0;
1025 lLead.matchDist = 0;
1026 lLead.IPxy = 0;
1027 lLead.IPz = 0;
1028 lLead.signedSipt = 0;
1029 }
1030
1031 }//pf
1032
1033
1034
1035 const std::vector<std::pair<std::string,float> > & lIDs = (*iter).tauIDs();
1036 bool lPrint = false;
1037 //a few warning if additional discriminants are found:
1038 if ((isPF && lIDs.size() != 16) ||
1039 (isCalo && lIDs.size() != 3))
1040 lPrint = true;
1041
1042 if (lPrint) {
1043 std::cout << "!!!!!!! Discriminants changed, please update tree !!!!!!!" << std::endl;
1044 std::cout << "--- isCaloTau = " << isCalo << ", isPFTau = " << isPF << std::endl;
1045 std::cout << "------ ID names = " << std::endl;
1046 }
1047
1048
1049 HbbAnalysis::CaloTauIDVars lcaloId;
1050 lcaloId.byIsolation = 0;
1051 lcaloId.byLeadingTrackFinding = 0;
1052 lcaloId.byLeadingTrackPtCut = 0;
1053
1054 HbbAnalysis::PFTauIDVars lpfId;
1055 lpfId.byLeadingTrackFinding = 0;
1056 lpfId.byLeadingTrackPtCut = 0;
1057 lpfId.byTrackIsolation = 0;
1058 lpfId.byECALIsolation = 0;
1059 lpfId.byIsolation = 0;
1060 lpfId.againstElectron = 0;
1061 lpfId.againstMuon = 0;
1062 lpfId.byIsolationUsingLeadingPion = 0;
1063 lpfId.byTaNC = 0;
1064 lpfId.byTaNCfrHalfPercent = 0;
1065 lpfId.byTaNCfrOnePercent = 0;
1066 lpfId.byTaNCfrQuarterPercent = 0;
1067 lpfId.byTaNCfrTenthPercent = 0;
1068 lpfId.ecalIsolationUsingLeadingPion = 0;
1069 lpfId.leadingPionPtCut = 0;
1070 lpfId.trackIsolationUsingLeadingPion = 0;
1071
1072
1073
1074 for (unsigned int id(0); id<lIDs.size(); ++id){
1075
1076 std::string lName = lIDs.at(id).first;
1077 float lDiscri = lIDs.at(id).second;
1078
1079 if (lPrint) std::cout << "--------- " << lName << " = " << lDiscri << std::endl;
1080 //pf
1081
1082 if (isPF) {
1083 if (lName.find("leadingTrackFinding") != lName.npos) lpfId.byLeadingTrackFinding = lDiscri;
1084 if (lName.find("leadingTrackPtCut") != lName.npos) lpfId.byLeadingTrackPtCut = lDiscri;
1085 if (lName.find("trackIsolationUsingLeadingPion") != lName.npos) lpfId.trackIsolationUsingLeadingPion = lDiscri;
1086 if (lName.find("trackIsolation") != lName.npos &&
1087 lName.find("trackIsolationUsingLeadingPion") == lName.npos) lpfId.byTrackIsolation = lDiscri;
1088 if (lName.find("ecalIsolationUsingLeadingPion") != lName.npos) lpfId.ecalIsolationUsingLeadingPion = lDiscri;
1089 if (lName.find("ecalIsolation") != lName.npos &&
1090 lName.find("ecalIsolationUsingLeadingPion") == lName.npos) lpfId.byECALIsolation = lDiscri;
1091 if (lName.find("byIsolationUsingLeadingPion") != lName.npos) lpfId.byIsolationUsingLeadingPion = lDiscri;
1092 if (lName.find("byIsolation") != lName.npos &&
1093 lName.find("byIsolationUsingLeadingPion") == lName.npos) lpfId.byIsolation = lDiscri;
1094 if (lName.find("againstElectron") != lName.npos) lpfId.againstElectron = lDiscri;
1095 if (lName.find("againstMuon") != lName.npos) lpfId.againstMuon = lDiscri;
1096 if (lName.find("byTaNCfrHalfPercent") != lName.npos) lpfId.byTaNCfrHalfPercent = lDiscri;
1097 if (lName.find("byTaNCfrOnePercent") != lName.npos) lpfId.byTaNCfrOnePercent = lDiscri;
1098 if (lName.find("byTaNCfrQuarterPercent") != lName.npos) lpfId.byTaNCfrQuarterPercent = lDiscri;
1099 if (lName.find("byTaNCfrTenthPercent") != lName.npos) lpfId.byTaNCfrTenthPercent = lDiscri;
1100 if (lName.find("byTaNC") != lName.npos &&
1101 lName.find("byTaNCfr") == lName.npos) lpfId.byTaNC = lDiscri;
1102 if (lName.find("leadingPionPtCut") != lName.npos) lpfId.leadingPionPtCut = lDiscri;
1103 }
1104 if (isCalo){
1105 if (lName.find("byIsolation") != lName.npos) lcaloId.byIsolation = lDiscri;
1106 if (lName.find("leadingTrackFinding") != lName.npos) lcaloId.byLeadingTrackFinding = lDiscri;
1107 if (lName.find("leadingTrackPtCut") != lName.npos) lcaloId.byLeadingTrackPtCut = lDiscri;
1108 }
1109
1110
1111 }
1112
1113
1114 if (isCalo){
1115
1116 HbbAnalysis::CaloTauIsoVars lcaloIso;
1117 lcaloIso.nIsoTracks = (*iter).isolationTracks().size();
1118 lcaloIso.nSigTracks = (*iter).signalTracks().size();
1119 lcaloIso.leadTrackHCAL3x3hitsEtSum = (*iter).leadTrackHCAL3x3hitsEtSum();
1120 lcaloIso.leadTrackHCAL3x3hottesthitDEta = (*iter).leadTrackHCAL3x3hottesthitDEta();
1121 lcaloIso.signalTracksInvariantMass = (*iter).signalTracksInvariantMass();
1122 lcaloIso.tracksInvariantMass = (*iter).TracksInvariantMass();
1123 lcaloIso.isolationTracksPtSum = (*iter).isolationTracksPtSum();
1124 lcaloIso.isolationECALhitsEtSum = (*iter).isolationECALhitsEtSum();
1125 lcaloIso.maximumHCALhitEt = (*iter).maximumHCALhitEt();
1126
1127 HbbAnalysis::Tau lObj(lGen,lGenJet,lReco,lLead,lcaloIso,lcaloId);
1128 aVec.push_back(lObj);
1129 }
1130 else {
1131
1132 HbbAnalysis::PFTauIsoVars lpfIso;
1133 lpfIso.nSigCands = (*iter).signalPFCands().size();
1134 lpfIso.nIsoCands = (*iter).isolationPFCands().size();
1135 lpfIso.maximumHCALPFClusterEt = (*iter).maximumHCALPFClusterEt();
1136 //lpfIso.emFraction = (*iter).emFraction();
1137 lpfIso.hcalTotOverPLead = (*iter).hcalTotOverPLead();
1138 lpfIso.hcalMaxOverPLead = (*iter).hcalMaxOverPLead();
1139 lpfIso.hcal3x3OverPLead = (*iter).hcal3x3OverPLead();
1140 lpfIso.ecalStripSumEOverPLead = (*iter).ecalStripSumEOverPLead();
1141 lpfIso.bremsRecoveryEOverPLead = (*iter).bremsRecoveryEOverPLead();
1142
1143 HbbAnalysis::PFTauCandVars lHadr;
1144 lHadr.nSigCands = (*iter).signalPFChargedHadrCands().size();
1145 lHadr.nIsoCands = (*iter).isolationPFChargedHadrCands().size();
1146 lHadr.isolationPtSum = (*iter).isolationPFChargedHadrCandsPtSum();
1147
1148 HbbAnalysis::PFTauCandVars lNeutr;
1149 lNeutr.nSigCands = (*iter).signalPFNeutrHadrCands().size();
1150 lNeutr.nIsoCands = (*iter).isolationPFNeutrHadrCands().size();
1151 lNeutr.isolationPtSum = (*iter).neutralHadronIso();
1152
1153 HbbAnalysis::PFTauCandVars lGamma;
1154 lGamma.nSigCands = (*iter).signalPFGammaCands().size();
1155 lGamma.nIsoCands = (*iter).isolationPFGammaCands().size();
1156 lGamma.isolationPtSum = (*iter).isolationPFGammaCandsEtSum();
1157
1158 HbbAnalysis::PFTauEleIDVars lEleId;
1159 if ( (*iter).electronPreIDTrack().isAvailable() && (*iter).electronPreIDTrack().isNonnull() ) {
1160 lEleId.pT = (*iter).electronPreIDTrack()->pt();
1161 lEleId.eta = (*iter).electronPreIDTrack()->eta();
1162 lEleId.phi = (*iter).electronPreIDTrack()->phi();
1163 }
1164 else {
1165 lEleId.pT = 0;
1166 lEleId.eta = 0;
1167 lEleId.phi = 0;
1168 }
1169
1170 lEleId.output = (*iter).electronPreIDOutput();
1171 lEleId.decision = (*iter).electronPreIDDecision();
1172
1173 HbbAnalysis::PFTauMuIDVars lMuId;
1174 lMuId.caloCompat = (*iter).caloComp();
1175 lMuId.segCompat = (*iter).segComp();
1176 lMuId.decision = (*iter).muonDecision();
1177
1178 HbbAnalysis::Tau lObj(lGen,lGenJet,lReco,lLead,lpfIso,lHadr,lNeutr,lGamma,lpfId,lEleId,lMuId);
1179 aVec.push_back(lObj);
1180 }
1181 }
1182 }
1183 }
1184
1185 }//HbbTaus
1186
1187 void HbbTreeMaker::HbbJets(const edm::Handle<std::vector<pat::Jet> > & aCol,
1188 const HbbAnalysis::JetFlavour & aJetFlav,
1189 const edm::Handle<reco::GenParticleCollection> & aGenParticles,
1190 std::vector<HbbAnalysis::Jet> & aVec)
1191 {//HbbJets
1192
1193 if (aCol.isValid()){//valid
1194 if (aCol->size() > 0) {//non empty
1195
1196 unsigned int iEle = 0;
1197 for (std::vector<pat::Jet>::const_iterator iter = aCol->begin();
1198 iter != aCol->end();
1199 ++iter,++iEle)
1200 {//loop on element
1201 //if (debug_ > 1) std::cout << "**** Ele #" << iEle << ", pT,eta,phi = " << (*iter).pt() << " " << (*iter).eta() << " " << (*iter).phi() << std::endl;
1202
1203 HbbAnalysis::GenVars lGen;
1204 if (!processData_ && (*iter).genParton()){
1205 lGen.valid = true;
1206 lGen.E = (*iter).genParton()->energy();
1207 lGen.pT = (*iter).genParton()->pt();
1208 lGen.eta = (*iter).genParton()->eta();
1209 lGen.y = (*iter).genParton()->y();
1210 lGen.phi = (*iter).genParton()->phi();
1211 lGen.charge = (*iter).genParton()->charge();
1212 lGen.pdgId = (*iter).genParton()->pdgId();
1213 lGen.status = (*iter).genParton()->status();
1214 lGen.mass = (*iter).genParton()->mass();
1215 lGen.vx = (*iter).genParton()->vx();
1216 lGen.vy = (*iter).genParton()->vy();
1217 lGen.vz = (*iter).genParton()->vz();
1218 }
1219 else {
1220 lGen.valid = false;
1221 lGen.E = 0;
1222 lGen.pT = 0;
1223 lGen.eta = 0;
1224 lGen.y = 0;
1225 lGen.phi = 0;
1226 lGen.charge = 0;
1227 lGen.pdgId = 0;
1228 lGen.status = 0;
1229 lGen.mass = 0;
1230 lGen.vx = 0;
1231 lGen.vy = 0;
1232 lGen.vz = 0;
1233 }
1234
1235
1236 HbbAnalysis::GenVars lGenJet;
1237 if (!processData_ && (*iter).genJet()){
1238 lGenJet.valid = true;
1239 lGenJet.E = (*iter).genJet()->energy();
1240 lGenJet.pT = (*iter).genJet()->pt();
1241 lGenJet.eta = (*iter).genJet()->eta();
1242 lGenJet.y = (*iter).genJet()->y();
1243 lGenJet.phi = (*iter).genJet()->phi();
1244 lGenJet.charge = (*iter).genJet()->charge();
1245 lGenJet.pdgId = (*iter).genJet()->pdgId();
1246 lGenJet.status = (*iter).genJet()->status();
1247 lGenJet.mass = (*iter).genJet()->mass();
1248 lGenJet.vx = (*iter).genJet()->vx();
1249 lGenJet.vy = (*iter).genJet()->vy();
1250 lGenJet.vz = (*iter).genJet()->vz();
1251 }
1252 else {
1253 lGenJet.valid = false;
1254 lGenJet.E = 0;
1255 lGenJet.pT = 0;
1256 lGenJet.eta = 0;
1257 lGenJet.y = 0;
1258 lGenJet.phi = 0;
1259 lGenJet.charge = 0;
1260 lGenJet.pdgId = 0;
1261 lGenJet.status = 0;
1262 lGenJet.mass = 0;
1263 lGenJet.vx = 0;
1264 lGenJet.vy = 0;
1265 lGenJet.vz = 0;
1266 }
1267
1268 HbbAnalysis::BaseVars lReco;
1269 lReco.E = (*iter).energy();
1270 lReco.pT = (*iter).pt();
1271 lReco.eta = (*iter).eta();
1272 //lReco.etaDet = (*iter).detectorEta((*iter).vz(),(*iter).eta());
1273 lReco.y = (*iter).y();
1274 lReco.phi = (*iter).phi();
1275 lReco.charge = (*iter).jetCharge();
1276 lReco.vx = (*iter).vx();
1277 lReco.vy = (*iter).vy();
1278 lReco.vz = (*iter).vz();
1279 lReco.hltMatch = (*iter).triggerObjectMatches().size();
1280 //lReco.myEtaDet = HbbAnalysis::EtaDetector(lReco);
1281
1282 unsigned int lFlav = 0;
1283 if (!processData_ && aJetFlav.nPartons()) {
1284 lFlav = aJetFlav.partonMatchingGenJet((*iter),aGenParticles,0.4).second;
1285 }
1286
1287 HbbAnalysis::JetVars lCommon;
1288 lCommon.flavour = lFlav;
1289 lCommon.partonFlavour = (*iter).partonFlavour();
1290 lCommon.nAssociatedTracks = ((*iter).associatedTracks()).size();
1291 lCommon.rawpT = (*iter).pt();
1292 //if ((*iter).hasCorrFactors() && (*iter).corrFactor("raw")>0)
1293 //lCommon.rawpT = (*iter).pt()/((*iter).corrFactor("raw"));
1294
1295
1296 // std::cout << " -- Current JEC SET " << (*iter).currentJECSet()
1297 // << " Level " << (*iter).currentJECLevel()
1298 // << std::endl;
1299
1300 HbbAnalysis::JetECorVars lEnergyCorrections;
1301 if ((*iter).jecSetsAvailable()){
1302 float lJecFactor = (*iter).jecFactor("Uncorrected","NONE",(*iter).currentJECSet());
1303 if (lJecFactor > 0) lCommon.rawpT = (*iter).pt()*lJecFactor;
1304 if (debug_){
1305 std::cout << " ---------------------------------------------- " << std::endl
1306 << " - Raw to Cor Jet # " << iEle << " = " << std::endl
1307 << " -- Raw pT = " << lCommon.rawpT << " eta " << lReco.eta << " phi " << lReco.phi << std::endl;
1308 }
1309
1310 std::vector<std::string> lJECnames = (*iter).availableJECSets();
1311 for (unsigned int iSet(0); iSet < lJECnames.size(); ++iSet){
1312 if (debug_) std::cout << " -- JEC SET " << iSet << " " << lJECnames[iSet] << std::endl;
1313 std::vector<std::string> lJEClevels = (*iter).availableJECLevels(iSet);
1314 for (unsigned int iLev(0); iLev<lJEClevels.size(); ++iLev){
1315 if (debug_) std::cout << " ---- JEC level " << iLev << " " << lJEClevels[iLev]
1316 << " , factor = " << (*iter).jecFactor(iLev,pat::JetCorrFactors::NONE,iSet)
1317 << std::endl;
1318 lEnergyCorrections.Levels.push_back(std::pair<std::string,double>(lJEClevels[iLev],(*iter).jecFactor(iLev,pat::JetCorrFactors::NONE,iSet)));
1319 }
1320 }
1321 if (debug_) std::cout << " -- Reco pT = " << lReco.pT
1322 << std::endl;
1323 }
1324
1325 if(jecUncCalo && jecUncPF && jecUncJPT){
1326 if ((*iter).isCaloJet()) {
1327 jecUncCalo->setJetEta((*iter).eta());
1328 jecUncCalo->setJetPt((*iter).pt());
1329 lEnergyCorrections.Levels.push_back(std::pair<std::string,double>("Uncertainty",jecUncCalo->getUncertainty(true)));
1330 } else if ((*iter).isPFJet()) {
1331 jecUncPF->setJetEta((*iter).eta());
1332 jecUncPF->setJetPt((*iter).pt());
1333 lEnergyCorrections.Levels.push_back(std::pair<std::string,double>("Uncertainty",jecUncPF->getUncertainty(true)));
1334 } else if ((*iter).isJPTJet()) {
1335 jecUncJPT->setJetEta((*iter).eta());
1336 jecUncJPT->setJetPt((*iter).pt());
1337 lEnergyCorrections.Levels.push_back(std::pair<std::string,double>("Uncertainty",jecUncJPT->getUncertainty(true)));
1338 }
1339 }
1340
1341 //lEnergyCorrections.L1 = (*iter).jecFactor("L1Offset","NONE",(*iter).currentJECSet());
1342 //lEnergyCorrections.L1JPT = (*iter).jecFactor("L1JPTOffset","NONE",(*iter).currentJECSet());
1343 //lEnergyCorrections.L2 = (*iter).jecFactor("L2Relative","NONE",(*iter).currentJECSet());
1344 //lEnergyCorrections.L3 = (*iter).jecFactor("L3Absolute","NONE",(*iter).currentJECSet());
1345 //if (processData_) lEnergyCorrections.L2L3Residual = (*iter).jecFactor("L2L3Residual","NONE",(*iter).currentJECSet());
1346
1347 lCommon.etaMean = (*iter).etaPhiStatistics().etaMean;
1348 lCommon.phiMean = (*iter).etaPhiStatistics().phiMean;
1349 lCommon.etaEtaMoment = (*iter).etaPhiStatistics().etaEtaMoment;
1350 lCommon.phiPhiMoment = (*iter).etaPhiStatistics().phiPhiMoment;
1351 lCommon.etaPhiMoment = (*iter).etaPhiStatistics().etaPhiMoment;
1352
1353
1354 //lCommon.rawEta = (*iter).;
1355 //lCommon.rawPhi = (*iter).;
1356 // p_hasJetCorrFactors->Fill(aJet.hasJetCorrFactors());
1357
1358 HbbAnalysis::JetBtagVars lBtag;
1359 lBtag.cSV = (*iter).bDiscriminator("combinedSecondaryVertexBJetTags");
1360 lBtag.cSVMVA = (*iter).bDiscriminator("combinedSecondaryVertexMVABJetTags");
1361 lBtag.iPMVA = (*iter).bDiscriminator("impactParameterMVABJetTags");
1362 lBtag.bProba = (*iter).bDiscriminator("jetBProbabilityBJetTags");
1363 lBtag.probability = (*iter).bDiscriminator("jetProbabilityBJetTags");
1364 lBtag.sSVHP = (*iter).bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
1365 lBtag.sSVHE = (*iter).bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
1366 lBtag.softElectronByPt = ((*iter).bDiscriminator("softElectronByPtBJetTags") < -100) ? -100 : (*iter).bDiscriminator("softElectronByPtBJetTags");
1367 lBtag.softElectronByIP3d = ((*iter).bDiscriminator("softElectronByIP3dBJetTags") < -100) ? -100 : (*iter).bDiscriminator("softElectronByIP3dBJetTags");
1368 lBtag.softMuon = ((*iter).bDiscriminator("softMuonBJetTags") < -100) ? -100 : (*iter).bDiscriminator("softMuonBJetTags");
1369 lBtag.softMuonByPt = ((*iter).bDiscriminator("softMuonByPtBJetTags") < -100) ? -100 : (*iter).bDiscriminator("softMuonByPtBJetTags");
1370 lBtag.softMuonByIP3d = ((*iter).bDiscriminator("softMuonByIP3dBJetTags") < -100) ? -100 : (*iter).bDiscriminator("softMuonByIP3dBJetTags");
1371
1372 //std::cout << " -- New values of b-discri for jet : " << iEle << std::endl
1373 // << " ---- softElecs: " << lBtag.softElectronByPt << " " << lBtag.softElectronByIP3d << std::endl
1374 // << " ---- softMus: " << lBtag.softMuon << " " << lBtag.softMuonByPt << " " << lBtag.softMuonByIP3d << std::endl;
1375
1376
1377 lBtag.tCHE = (*iter).bDiscriminator("trackCountingHighEffBJetTags");
1378 lBtag.tCHP = (*iter).bDiscriminator("trackCountingHighPurBJetTags");
1379
1380
1381 bool isCalo = (*iter).isCaloJet();
1382 bool isPF = (*iter).isPFJet();
1383 bool isJPT = (*iter).isJPTJet();
1384
1385 //std::cout << " -- isCalo = " << isCalo << ", isPF = " << isPF << std::endl;
1386
1387 //assert (isCalo == !isPF);
1388 HbbAnalysis::JetIDVars lId;
1389 lId.fHPD = 0;
1390 lId.fRBX = 0;
1391 lId.n90Hits = 0;
1392 lId.hitsInN90 = 0;
1393 lId.restrictedEMF = 0;
1394 lId.nHCALTowers = 0;
1395 lId.nECALTowers = 0;
1396 if (isCalo || isJPT) {
1397 lId.fHPD = (*iter).jetID().fHPD;
1398 lId.fRBX = (*iter).jetID().fRBX;
1399 lId.n90Hits = (*iter).jetID().n90Hits;
1400 lId.hitsInN90 = (*iter).jetID().hitsInN90;
1401 lId.restrictedEMF = (*iter).jetID().restrictedEMF;
1402 lId.nHCALTowers = (*iter).jetID().nHCALTowers;
1403 lId.nECALTowers = (*iter).jetID().nECALTowers;
1404 }
1405 if (isCalo){
1406 HbbAnalysis::CaloJetVars lCalo;
1407 lCalo.maxEInEmTowers = (*iter).maxEInEmTowers();
1408 lCalo.maxEInHadTowers = (*iter).maxEInHadTowers();
1409 lCalo.energyFractionHadronic = (*iter).energyFractionHadronic();
1410 lCalo.emEnergyFraction = (*iter).emEnergyFraction();
1411 lCalo.hadEnergyInHB = (*iter).hadEnergyInHB();
1412 lCalo.hadEnergyInHO = (*iter).hadEnergyInHO();
1413 lCalo.hadEnergyInHE = (*iter).hadEnergyInHE();
1414 lCalo.hadEnergyInHF = (*iter).hadEnergyInHF();
1415 lCalo.emEnergyInEB = (*iter).emEnergyInEB();
1416 lCalo.emEnergyInEE = (*iter).emEnergyInEE();
1417 lCalo.emEnergyInHF = (*iter).emEnergyInHF();
1418 lCalo.towersArea = (*iter).towersArea();
1419 lCalo.n90 = (*iter).n90();
1420 lCalo.n60 = (*iter).n60();
1421
1422 HbbAnalysis::Jet lObj(lGen,lGenJet,lReco,lCommon,lCalo,lBtag,lId,lEnergyCorrections);
1423 aVec.push_back(lObj);
1424 }
1425 HbbAnalysis::JPTPFJetVars lJPTPF;
1426 lJPTPF.chargedHadronEnergy = 0;
1427 lJPTPF.chargedHadronEnergyFraction = 0;
1428 lJPTPF.neutralHadronEnergy = 0;
1429 lJPTPF.neutralHadronEnergyFraction = 0;
1430 lJPTPF.chargedEmEnergy = 0;
1431 lJPTPF.chargedEmEnergyFraction = 0;
1432 lJPTPF.neutralEmEnergy = 0;
1433 lJPTPF.neutralEmEnergyFraction = 0;
1434 lJPTPF.chargedMultiplicity = 0;
1435 lJPTPF.muonMultiplicity = 0;
1436
1437 if (isPF || isJPT){
1438 lJPTPF.chargedEmEnergy = (*iter).chargedEmEnergy();
1439 lJPTPF.neutralEmEnergy = (*iter).neutralEmEnergy();
1440 lJPTPF.muonMultiplicity = (*iter).muonMultiplicity();
1441 lJPTPF.chargedHadronEnergy = (*iter).chargedHadronEnergy();
1442 lJPTPF.neutralHadronEnergy = (*iter).neutralHadronEnergy();
1443 lJPTPF.chargedMultiplicity = (*iter).chargedMultiplicity();
1444
1445 lJPTPF.chargedHadronEnergyFraction = (*iter).chargedHadronEnergyFraction();
1446 lJPTPF.neutralHadronEnergyFraction = (*iter).neutralHadronEnergyFraction();
1447 lJPTPF.chargedEmEnergyFraction = (*iter).chargedEmEnergyFraction();
1448 lJPTPF.neutralEmEnergyFraction = (*iter).neutralEmEnergyFraction();
1449 }
1450 if (isJPT) {
1451 HbbAnalysis::JPTJetVars lJPT;
1452 reco::TrackRefVector pionsInVInC = (*iter).pionsInVertexInCalo();
1453 reco::TrackRefVector pionsInVOuC = (*iter).pionsInVertexOutCalo();
1454 reco::TrackRefVector pionsOuVInC = (*iter).pionsOutVertexInCalo();
1455
1456 reco::TrackRefVector muonsInVInC = (*iter).muonsInVertexInCalo();
1457 reco::TrackRefVector muonsInVOuC = (*iter).muonsInVertexOutCalo();
1458 reco::TrackRefVector muonsOuVInC = (*iter).muonsOutVertexInCalo();
1459
1460 reco::TrackRefVector elecsInVInC = (*iter).elecsInVertexInCalo();
1461 reco::TrackRefVector elecsInVOuC = (*iter).elecsInVertexOutCalo();
1462 reco::TrackRefVector elecsOuVInC = (*iter).elecsOutVertexInCalo();
1463
1464 /*
1465 unsigned test_track_count = pionsInVInC.size() + pionsInVOuC.size() + pionsOuVInC.size();
1466 std::cout << "Total Pion Track Count = " << test_track_count << std::endl;
1467 */
1468
1469 for (reco::TrackRefVector::const_iterator trk_iter = pionsInVInC.begin();
1470 trk_iter != pionsInVInC.end();
1471 ++trk_iter) {
1472 lJPT.particleStatusPt.push_back(std::pair<unsigned int,double>(1,(*trk_iter)->pt()));
1473 }
1474 for (reco::TrackRefVector::const_iterator trk_iter = pionsInVOuC.begin();
1475 trk_iter != pionsInVOuC.end();
1476 ++trk_iter) {
1477 lJPT.particleStatusPt.push_back(std::pair<unsigned int,double>(2,(*trk_iter)->pt()));
1478 }
1479 for (reco::TrackRefVector::const_iterator trk_iter = pionsOuVInC.begin();
1480 trk_iter != pionsOuVInC.end();
1481 ++trk_iter) {
1482 lJPT.particleStatusPt.push_back(std::pair<unsigned int,double>(3,(*trk_iter)->pt()));
1483 }
1484
1485 for (reco::TrackRefVector::const_iterator trk_iter = muonsInVInC.begin();
1486 trk_iter != muonsInVInC.end();
1487 ++trk_iter) {
1488 lJPT.particleStatusPt.push_back(std::pair<unsigned int,double>(4,(*trk_iter)->pt()));
1489 }
1490 for (reco::TrackRefVector::const_iterator trk_iter = muonsInVOuC.begin();
1491 trk_iter != muonsInVOuC.end();
1492 ++trk_iter) {
1493 lJPT.particleStatusPt.push_back(std::pair<unsigned int,double>(5,(*trk_iter)->pt()));
1494 }
1495 for (reco::TrackRefVector::const_iterator trk_iter = muonsOuVInC.begin();
1496 trk_iter != muonsOuVInC.end();
1497 ++trk_iter) {
1498 lJPT.particleStatusPt.push_back(std::pair<unsigned int,double>(6,(*trk_iter)->pt()));
1499 }
1500
1501 for (reco::TrackRefVector::const_iterator trk_iter = elecsInVInC.begin();
1502 trk_iter != elecsInVInC.end();
1503 ++trk_iter) {
1504 lJPT.particleStatusPt.push_back(std::pair<unsigned int,double>(7,(*trk_iter)->pt()));
1505 }
1506 for (reco::TrackRefVector::const_iterator trk_iter = elecsInVOuC.begin();
1507 trk_iter != elecsInVOuC.end();
1508 ++trk_iter) {
1509 lJPT.particleStatusPt.push_back(std::pair<unsigned int,double>(8,(*trk_iter)->pt()));
1510 }
1511 for (reco::TrackRefVector::const_iterator trk_iter = elecsOuVInC.begin();
1512 trk_iter != elecsOuVInC.end();
1513 ++trk_iter) {
1514 lJPT.particleStatusPt.push_back(std::pair<unsigned int,double>(9,(*trk_iter)->pt()));
1515 }
1516
1517 // const reco::TrackRefVector& pionsInVertexInCalo () const{return jptSpecific().pionsInVertexInCalo; }
1518 // const reco::TrackRefVector& pionsInVertexOutCalo() const{return jptSpecific().pionsInVertexOutCalo;}
1519 // const reco::TrackRefVector& pionsOutVertexInCalo() const{return jptSpecific().pionsOutVertexInCalo;}
1520 // const reco::TrackRefVector& muonsInVertexInCalo () const{return jptSpecific().muonsInVertexInCalo; }
1521 // const reco::TrackRefVector& muonsInVertexOutCalo() const{return jptSpecific().muonsInVertexOutCalo;}
1522 // const reco::TrackRefVector& muonsOutVertexInCalo() const{return jptSpecific().muonsOutVertexInCalo;}
1523 // const reco::TrackRefVector& elecsInVertexInCalo () const{return jptSpecific().elecsInVertexInCalo; }
1524 // const reco::TrackRefVector& elecsInVertexOutCalo() const{return jptSpecific().elecsInVertexOutCalo;}
1525 // const reco::TrackRefVector& elecsOutVertexInCalo() const{return jptSpecific().elecsOutVertexInCalo;}
1526 lJPT.zspCorrection = (*iter).zspCorrection();
1527 lJPT.elecMultiplicity = (*iter).elecMultiplicity();
1528
1529 HbbAnalysis::CaloJetVars lCalo;
1530 const reco::JPTJet * lJetRef = dynamic_cast<const reco::JPTJet *>((*iter).originalObject());
1531 edm::RefToBase<reco::Jet> jptjetRef = lJetRef->getCaloJetRef();
1532
1533 // Dynamic_cast to CaloJet object
1534 reco::CaloJet const * rawcalojet =
1535 dynamic_cast<reco::CaloJet const *>( &* jptjetRef);
1536 lJPT.calopT = rawcalojet->pt();
1537 lJPT.caloEta = rawcalojet->eta();
1538 if (debug_) {
1539 std::cout << " -- Raw CaloJet pT = " << rawcalojet->pt() << std::endl;
1540 std::cout << " -- Raw CaloJet Eta = " << rawcalojet->eta() << std::endl;
1541 }
1542 lCalo.maxEInEmTowers = rawcalojet->maxEInEmTowers();
1543 lCalo.maxEInHadTowers = rawcalojet->maxEInHadTowers();
1544 lCalo.energyFractionHadronic = rawcalojet->energyFractionHadronic();
1545 lCalo.emEnergyFraction = rawcalojet->emEnergyFraction();
1546 lCalo.hadEnergyInHB = rawcalojet->hadEnergyInHB();
1547 lCalo.hadEnergyInHO = rawcalojet->hadEnergyInHO();
1548 lCalo.hadEnergyInHE = rawcalojet->hadEnergyInHE();
1549 lCalo.hadEnergyInHF = rawcalojet->hadEnergyInHF();
1550 lCalo.emEnergyInEB = rawcalojet->emEnergyInEB();
1551 lCalo.emEnergyInEE = rawcalojet->emEnergyInEE();
1552 lCalo.emEnergyInHF = rawcalojet->emEnergyInHF();
1553 lCalo.towersArea = rawcalojet->towersArea();
1554 lCalo.n90 = rawcalojet->n90();
1555 lCalo.n60 = rawcalojet->n60();
1556
1557 lCommon.etaMean = rawcalojet->etaPhiStatistics().etaMean;
1558 lCommon.phiMean = rawcalojet->etaPhiStatistics().phiMean;
1559 lCommon.etaEtaMoment = rawcalojet->etaPhiStatistics().etaEtaMoment;
1560 lCommon.phiPhiMoment = rawcalojet->etaPhiStatistics().phiPhiMoment;
1561 lCommon.etaPhiMoment = rawcalojet->etaPhiStatistics().etaPhiMoment;
1562
1563 HbbAnalysis::Jet lObj(lGen,lGenJet,lReco,lCommon,lCalo,lJPT,lJPTPF,lBtag,lId,lEnergyCorrections);
1564 aVec.push_back(lObj);
1565 }
1566 if (isPF) {
1567 HbbAnalysis::PFJetVars lPF;
1568 lPF.chargedMuEnergy = (*iter).chargedMuEnergy();
1569 lPF.chargedMuEnergyFraction = (*iter).chargedMuEnergyFraction();
1570 lPF.neutralMultiplicity = (*iter).neutralMultiplicity();
1571 lPF.numberOfDaughters = (*iter).numberOfDaughters();
1572 lPF.HFHadronEnergy = (*iter).HFHadronEnergy();
1573
1574 HbbAnalysis::Jet lObj(lGen,lGenJet,lReco,lCommon,lPF,lJPTPF,lBtag,lEnergyCorrections);
1575 aVec.push_back(lObj);
1576 }
1577
1578 }//loop on element
1579 }//non empty
1580 }//valid
1581
1582 }//HbbJets
1583
1584 void HbbTreeMaker::HbbMet(const edm::Handle<std::vector<pat::MET> > & aCol,
1585 HbbAnalysis::Met & aMet)
1586 {//HbbMet
1587 if (aCol.isValid()) {
1588 HbbAnalysis::MetVars lGen;
1589 HbbAnalysis::MetVars lReco;
1590
1591 assert(aCol->size() == 1);
1592 const pat::MET & lMet = *(aCol->begin());
1593
1594 lReco.mET = lMet.pt();
1595 lReco.mEx = lMet.px();
1596 lReco.mEy = lMet.py();
1597 lReco.sumET = lMet.sumEt();
1598 lReco.phi = lMet.phi();
1599 lReco.mEtSig = lMet.mEtSig();
1600
1601 const reco::GenMET *lGenMET = lMet.genMET();
1602
1603 if (!processData_ && lGenMET){
1604 lGen.mET = lGenMET->pt();
1605 lGen.mEx = lGenMET->px();
1606 lGen.mEy = lGenMET->py();
1607 lGen.sumET = lGenMET->sumEt();
1608 lGen.phi = lGenMET->phi();
1609 lGen.mEtSig = lGenMET->mEtSig();
1610 }
1611 else {
1612 lGen.mET = 0;
1613 lGen.mEx = 0;
1614 lGen.mEy = 0;
1615 lGen.sumET = 0;
1616 lGen.phi = 0;
1617 lGen.mEtSig = 0;
1618 }
1619
1620 aMet.genVars(lGen);
1621 aMet.recoVars(lReco);
1622 }
1623 }//HbbMet
1624
1625 void HbbTreeMaker::HbbTrigger(const edm::Handle<edm::TriggerResults> & aCol,
1626 const edm::TriggerNames & aNames,
1627 std::vector<HbbAnalysis::Trigger> & aVec)
1628 {//HbbTrigger
1629
1630 if (aCol.isValid()){
1631 //edm::TriggerNames lNames;
1632 //lNames.init(*aCol);
1633
1634 //for (unsigned int j=0; j<hltPaths_.size(); j++) {
1635 //bool valid=false;
1636 for (unsigned int i=0; i<aCol->size(); ++i){
1637 std::string trueName = aNames.triggerName(i);//.at(i);
1638
1639 HbbAnalysis::TriggerVars lTrigVars;
1640 lTrigVars.name = trueName;
1641 lTrigVars.index = i;
1642 lTrigVars.accept = aCol->accept(i);
1643
1644 HbbAnalysis::Trigger lObj(lTrigVars);
1645
1646 aVec.push_back(lObj);
1647 }
1648
1649 }
1650
1651 }
1652
1653 void HbbTreeMaker::HbbParticles(const edm::Handle<reco::GenParticleCollection> & aCol,
1654 std::vector<HbbAnalysis::GenParticle> & aVec)
1655 {//genparticles
1656
1657 //add genParticle inside vector
1658 for (unsigned int mccount = 0;mccount<aCol->size();++mccount)
1659 {//loop on particles
1660 const reco::Candidate & p = (*aCol)[mccount];
1661
1662 HbbAnalysis::MCVars lMC;
1663 lMC.index = mccount;
1664 lMC.E = p.energy();
1665 lMC.pT = p.pt();
1666 lMC.eta = p.eta();
1667 lMC.y = p.y();
1668 lMC.phi = p.phi();
1669 lMC.pdgId = p.pdgId();
1670 lMC.status = p.status();
1671 if (p.status()==3 ||
1672 (p.status()==2 && (abs(p.pdgId())==4 || abs(p.pdgId())==5)) ||
1673 (p.status()==1 && ( (p.pt()>5 && abs(p.pdgId()) > 16) || abs(p.pdgId()) <= 16) )
1674 ){
1675 HbbAnalysis::GenParticle lObj(lMC);
1676 aVec.push_back(lObj);
1677 }
1678 }//loop on particles
1679
1680 //now need to get the list of parents....
1681
1682 for (unsigned int mccount = 0;mccount<aCol->size();++mccount)
1683 {//loop on particles
1684 const reco::Candidate & p = (*aCol)[mccount];
1685
1686 unsigned int nD = p.numberOfDaughters();
1687 //std::cout << mccount << " has " << nD << " daughter(s) : " ;
1688 for(unsigned int dau=0;dau<nD;++dau){//loop on daughters
1689 const reco::Candidate * pDau = p.daughter(dau);
1690 //std::cout << pDau << " ("<<std::flush;
1691 for (unsigned int mccount2 = 0;mccount2<aCol->size();++mccount2)
1692 {//loop on particles
1693 const reco::Candidate & p2 = (*aCol)[mccount2];
1694 //std::cout<<"DBG: mccount2 = "<<mccount2<<" gen size = "<<data_->gen().size()<<std::endl;
1695 for (unsigned int imc(0); imc < aVec.size(); ++imc){
1696 HbbAnalysis::GenParticle & part2 = aVec.at(imc);
1697 //part2 = aVec.at(imc);
1698 if (part2.partVars().index == mccount2) {
1699 if (pDau == &p2) {
1700 part2.setParent(mccount);
1701 //std::cout << &p2 << ", index = " << mccount2 << "), "<<std::flush;
1702 break;
1703 }
1704 }
1705 }
1706 }//loop on particles
1707 //else std::cout << std::endl << "****no match " << mccount2 << " " << &p2 << std::endl;
1708 //if (p2.pdgId() == pDau->pdgId() && fabs(p2.energy()-pDau->energy())<0.000001 && fabs(p2.eta()-pDau->eta())<0.000001 && fabs(p2.phi()-pDau->phi())<0.000001){
1709 //part2.parent(mccount);
1710 //}
1711 }//loop on daughters
1712 //std::cout << std::endl;
1713
1714 }//loop on particles
1715
1716 if (debug_>1){
1717 for (unsigned int imc(0); imc < aVec.size(); ++imc){
1718 HbbAnalysis::GenParticle & p = aVec.at(imc);
1719 p.print();
1720 }
1721 }
1722
1723
1724 }//genparticles
1725
1726 void HbbTreeMaker::HbbGenJets(const edm::Handle<reco::GenJetCollection> & aCol,
1727 std::vector<HbbAnalysis::GenJet> & aVec,
1728 const bool aSaveConstituants)
1729 {
1730
1731 if (aCol.isValid()){//valid
1732 if (aCol->size() > 0) {//non empty
1733
1734 unsigned int iEle = 0;
1735 for (std::vector<reco::GenJet>::const_iterator iter = aCol->begin();
1736 iter != aCol->end();
1737 ++iter,++iEle)
1738 {//loop on element
1739 //if (debug_ > 1) std::cout << "**** Ele #" << iEle << ", pT,eta,phi = " << (*iter).pt() << " " << (*iter).eta() << " " << (*iter).phi() << std::endl;
1740
1741 HbbAnalysis::GenVars lGenJet;
1742 if (!processData_ ){
1743 lGenJet.valid = true;
1744 lGenJet.E = (*iter).energy();
1745 lGenJet.pT = (*iter).pt();
1746 lGenJet.eta = (*iter).eta();
1747 lGenJet.y = (*iter).y();
1748 lGenJet.phi = (*iter).phi();
1749 lGenJet.charge = (*iter).charge();
1750
1751 //only status 1 saved !!!!
1752 //get flavour from list of constituants.
1753 //std::vector <const reco::GenParticle*> lConst = (*iter).getGenConstituents();
1754 //unsigned int nb=0, nc=0;
1755 // for (unsigned int iMC(0); iMC<lConst.size(); ++iMC){
1756 // const reco::GenParticle* lPart = lConst[iMC];
1757 // //std::cout << lPart->pdgId() << " " << fabs(lPart->pdgId()) << " " << abs(lPart->pdgId()) << std::endl;
1758 // if (fabs(lPart->pdgId()) == 5){
1759 // ++nb;
1760 // break;
1761 // }
1762 // if (fabs(lPart->pdgId()) == 4) ++nc;
1763 // }
1764 // if (nb > 0) lGenJet.pdgId = 5;
1765 // else if (nc > 0) lGenJet.pdgId = 4;
1766 // else lGenJet.pdgId = 3;
1767 lGenJet.pdgId = (*iter).pdgId();
1768 lGenJet.status = (*iter).status();
1769 lGenJet.mass = (*iter).mass();
1770 lGenJet.vx = (*iter).vx();
1771 lGenJet.vy = (*iter).vy();
1772 lGenJet.vz = (*iter).vz();
1773 }
1774 else {
1775 lGenJet.valid = false;
1776 lGenJet.E = 0;
1777 lGenJet.pT = 0;
1778 lGenJet.eta = 0;
1779 lGenJet.y = 0;
1780 lGenJet.phi = 0;
1781 lGenJet.charge = 0;
1782 lGenJet.pdgId = 0;
1783 lGenJet.status = 0;
1784 lGenJet.mass = 0;
1785 lGenJet.vx = 0;
1786 lGenJet.vy = 0;
1787 lGenJet.vz = 0;
1788 }
1789
1790 HbbAnalysis::GenJet lObj(lGenJet);
1791 if (aSaveConstituants){
1792 for (unsigned int iC(0); iC<(*iter).getGenConstituents().size(); ++iC){
1793 const reco::GenParticle* lPart = (*iter).getGenConstituent(iC);
1794
1795 HbbAnalysis::MCVars lMC;
1796 lMC.pT = lPart->pt();
1797 lMC.eta = lPart->eta();
1798 lMC.phi = lPart->phi();
1799 lMC.pdgId = lPart->pdgId();
1800 lObj.addConstituent(lMC);
1801 }
1802 }
1803 else lObj.nConstituents((*iter).getGenConstituents().size());
1804 aVec.push_back(lObj);
1805 }
1806 }
1807 }
1808 }
1809
1810 void HbbTreeMaker::HbbVertices(const edm::Handle<std::vector<reco::Vertex> > & aCol,
1811 std::vector<HbbAnalysis::Vertex> & aVec)
1812 {
1813
1814 if (aCol.isValid()) {
1815 for (std::vector<reco::Vertex>::const_iterator iter = aCol->begin();
1816 iter != aCol->end();
1817 ++iter)
1818 {
1819
1820 if (!((*iter).isValid()) || (*iter).isFake()) continue;
1821
1822 HbbAnalysis::VertexVars lVtx;
1823
1824 if ((*iter).tracksSize() > 0) {
1825 for (std::vector<reco::TrackBaseRef>::const_iterator lTrk = (*iter).tracks_begin();
1826 lTrk != (*iter).tracks_end();
1827 ++lTrk) {
1828 lVtx.trackWeights.push_back((*iter).trackWeight(*lTrk));
1829 lVtx.trackpT.push_back((*(*lTrk)).pt());
1830 }
1831 }
1832
1833 if (lVtx.trackWeights.size() != (*iter).tracksSize()) {
1834 std::cout<< " -- Problem with tracks, size is not as expected ! "
1835 << std::endl
1836 << " --- Size of recoVertex tracks = " << (*iter).tracksSize()
1837 << std::endl
1838 << " --- Size of trackWeights = " << lVtx.trackWeights.size()
1839 << std::endl;
1840 }
1841 if (lVtx.trackpT.size() != (*iter).tracksSize()) {
1842 std::cout<< " -- Problem with tracks, size is not as expected ! "
1843 << std::endl
1844 << " --- Size of recoVertex tracks = " << (*iter).tracksSize()
1845 << std::endl
1846 << " --- Size of trackpT = " << lVtx.trackpT.size()
1847 << std::endl;
1848 }
1849
1850 lVtx.chi2 = (*iter).chi2();
1851 lVtx.ndof = (*iter).ndof();
1852 lVtx.x = (*iter).x();
1853 lVtx.y = (*iter).y();
1854 lVtx.z = (*iter).z();
1855 lVtx.rho = (*iter).position().Rho();
1856 lVtx.xError = (*iter).xError();
1857 lVtx.yError = (*iter).yError();
1858 lVtx.zError = (*iter).zError();
1859 lVtx.cov01 = (*iter).covariance(0,1);
1860 lVtx.cov02 = (*iter).covariance(0,2);
1861 lVtx.cov12 = (*iter).covariance(1,2);
1862
1863 //AJG
1864 /*
1865 std::cout << "Vertex (" << (*iter).x() << ","
1866 << (*iter).y() << ","
1867 << (*iter).z() << ")" << std::endl;
1868 */
1869
1870 HbbAnalysis::Vertex lObj(lVtx);
1871 aVec.push_back(lObj);
1872 }
1873 }
1874
1875 }
1876
1877 void HbbTreeMaker::HbbL1Objects(const edm::Handle<l1extra::L1JetParticleCollection> & aCol,
1878 std::vector<HbbAnalysis::L1Object> & aVec,
1879 const unsigned int aType)
1880 {
1881
1882 for (l1extra::L1JetParticleCollection::const_iterator iter = aCol->begin();
1883 iter != aCol->end();
1884 ++iter)
1885 {
1886
1887 HbbAnalysis::L1Vars lVars;
1888 lVars.pT = (*iter).pt();
1889 lVars.ET = (*iter).et();
1890 lVars.eta = (*iter).eta();
1891 lVars.phi = (*iter).phi();
1892 lVars.bx = (*iter).bx();
1893 lVars.type = aType;
1894
1895 HbbAnalysis::L1Object lObj(lVars);
1896 aVec.push_back(lObj);
1897
1898 }
1899
1900 }
1901
1902
1903 void HbbTreeMaker::HbbHLTObjects(const edm::Handle<trigger::TriggerEvent> & aCol,
1904 std::vector<HbbAnalysis::HLTObject> & aVec)
1905 {
1906
1907
1908 //cout
1909 //const trigger::TriggerObjectCollection & toc(aCol->getObjects());
1910 //for (unsigned i=0; i<aCol->sizeFilters(); i++){
1911 //const edm::InputTag & lTag = aCol->filterTag(i);
1912 //std::cout << " -- ele " << i << ": " << lTag << ", key ids = " ;
1913 ////const int index = aCol->filterIndex(lTag);
1914 ////if (index >= aCol->sizeFilters()) {
1915 ////std::cout << " -- Filter " << aTag << " not found !" << std::endl;
1916 ////continue;
1917 ////}
1918 //const trigger::Keys & k = aCol->filterKeys(i);
1919 //for (trigger::Keys::const_iterator ki = k.begin(); ki !=k.end(); ++ki ) {
1920 //std::cout << toc[*ki].id() << ", " ;
1921 //}
1922 //std::cout << std::endl;
1923 //}
1924
1925 }
1926
1927 void HbbTreeMaker::fillHLTVector(const edm::InputTag & aTag,
1928 const edm::Handle<trigger::TriggerEvent> & aCol,
1929 std::vector<HbbAnalysis::HLTObject> & aVec,
1930 const unsigned int aType){
1931
1932 const trigger::TriggerObjectCollection & toc(aCol->getObjects());
1933 const int index = aCol->filterIndex(aTag);
1934 if (index >= aCol->sizeFilters()) {
1935 //std::cout << " -- Filter " << aTag << " not found !" << std::endl;
1936 return;
1937 }
1938
1939 const trigger::Keys & k = aCol->filterKeys(index);
1940
1941 for (trigger::Keys::const_iterator ki = k.begin(); ki !=k.end(); ++ki ) {
1942 HbbAnalysis::HLTVars lVars;
1943 lVars.pT = toc[*ki].pt();
1944 lVars.eta = toc[*ki].eta();
1945 lVars.phi = toc[*ki].phi();
1946 lVars.id = toc[*ki].id();
1947 lVars.mass = toc[*ki].mass();
1948 lVars.type = aType;
1949 HbbAnalysis::HLTObject lObj(lVars);
1950 aVec.push_back(lObj);
1951 }
1952 }
1953
1954 void HbbTreeMaker::HbbGenInfo(const edm::Handle<edm::HepMCProduct> & amcProd,
1955 const edm::Handle<GenRunInfoProduct> & aRunProd)
1956 {
1957 const HepMC::GenEvent * genEvt = amcProd->GetEvent();
1958 HbbAnalysis::GenInfo lInfo;
1959
1960 lInfo.processID(genEvt->signal_process_id());
1961 lInfo.ptHat(genEvt->event_scale());
1962
1963 lInfo.internalXS(aRunProd->internalXSec().value(),
1964 aRunProd->internalXSec().error());
1965
1966 lInfo.externalXSLO(aRunProd->externalXSecLO().value(),
1967 aRunProd->externalXSecLO().error());
1968
1969 lInfo.externalXSNLO(aRunProd->externalXSecNLO().value(),
1970 aRunProd->externalXSecNLO().error());
1971
1972 lInfo.filterEff(aRunProd->filterEfficiency());
1973
1974 event_->generatorInfo(lInfo);
1975
1976 }
1977
1978 void HbbTreeMaker::HbbPUInfo(const edm::Handle<std::vector<PileupSummaryInfo> > & aPUInfo,
1979 std::vector<HbbAnalysis::PUVars> & aVec) {
1980 std::vector<PileupSummaryInfo>::const_iterator PVI;
1981 for (PVI = aPUInfo->begin(); PVI != aPUInfo->end(); ++PVI) {
1982 HbbAnalysis::PUVars aVars;
1983 aVars.bunchCrossing = PVI->getBunchCrossing();
1984 aVars.numInteractions = PVI->getPU_NumInteractions();
1985 aVec.push_back(aVars);
1986 }
1987 }
1988
1989 void HbbTreeMaker::HbbLHEInfo(const edm::Handle<LHEEventProduct> & aLHEEvt,
1990 std::vector<HbbAnalysis::GenParticle> & aVec) {
1991
1992 const lhef::HEPEUP lhepeup = aLHEEvt->hepeup();
1993
1994 //4-momentum (px,py,pz,E)+mass
1995 const std::vector<lhef::HEPEUP::FiveVector> lpup = lhepeup.PUP;
1996 //pdgid int
1997 const std::vector<int> lidup = lhepeup.IDUP;
1998 //status int
1999 const std::vector<int> listup = lhepeup.ISTUP;
2000 //mother: pair<int,int>
2001 const std::vector<std::pair<int,int> > lmothup = lhepeup.MOTHUP;
2002
2003
2004 //add genParticle inside vector
2005 for (unsigned int imc = 0;imc<lpup.size();++imc)
2006 {//loop on particles
2007
2008 HbbAnalysis::MCVars lMC;
2009 lMC.index = imc;
2010 lMC.E = (lpup[imc])[3];
2011
2012 TLorentzVector lVec((lpup[imc])[0],(lpup[imc])[1],(lpup[imc])[2],(lpup[imc])[3]);
2013
2014 lMC.pT = lVec.Pt();
2015 // Special case if pT = 0 (otherwise ROOT complains)
2016 if (fabs(0.0 - lVec.Pt()) < std::numeric_limits<double>::epsilon()){
2017 if (lVec.Pz() >= 0) {
2018 lMC.eta = 1E10;
2019 lMC.y = 1E10;
2020 } else {
2021 lMC.eta = -1E10;
2022 lMC.y = -1E10;
2023 }
2024 } else {
2025 lMC.eta = lVec.Eta();
2026 lMC.y = lVec.Rapidity();
2027 }
2028 lMC.phi = lVec.Phi();
2029 lMC.pdgId = lidup[imc];
2030 lMC.status = listup[imc];
2031
2032 HbbAnalysis::GenParticle lObj(lMC);
2033 lObj.setParent(lmothup[imc].first);
2034 lObj.setParent(lmothup[imc].second);
2035 aVec.push_back(lObj);
2036
2037 }//loop on particles
2038
2039 if (debug_>1){
2040 for (unsigned int imc(0); imc < aVec.size(); ++imc){
2041 HbbAnalysis::GenParticle & p = aVec.at(imc);
2042 p.print();
2043 }
2044 }
2045 }
2046
2047 #include "FWCore/Framework/interface/MakerMacros.h"
2048 DEFINE_FWK_MODULE(HbbTreeMaker);
2049
2050