ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/HbbAnalysis/plugins/HbbTreeMaker.cc
Revision: 1.48
Committed: Wed Jan 25 18:45:49 2012 UTC (13 years, 3 months ago) by agilbert
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.47: +53 -18 lines
Log Message:
Lots of updates

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