ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbbAnalysis/HbbAnalyzer/plugins/HbbAnalyzerNew.cc
Revision: 1.75
Committed: Thu May 24 15:55:16 2012 UTC (12 years, 11 months ago) by arizzi
Content type: text/plain
Branch: MAIN
CVS Tags: Step2ForV33_v2, Step2ForV33_v1, EdmV33Jun12v0, Step2ForV32_v2, Step2ForV32_v0, EdmV32May24v0
Changes since 1.74: +7 -6 lines
Log Message:
store PU true number of PU and change pfMET input tag to be really pfMET

File Contents

# Content
1 // -*- C++ -*-
2 //
3 // Package: HbbAnalyzerNew
4 // Class: HbbAnalyzerNew
5 //
6 /**\class HbbAnalyzerNew HbbAnalyzerNew.cc Analysis/HbbAnalyzer/src/HbbAnalyzerNew.cc
7
8 Description: <one line class summary>
9
10 Implementation:
11 <Notes on implementation>
12 */
13 //
14 // Original Author: David Lopes Pegna,Address unknown,NONE,
15 // Created: Thu Mar 5 13:51:28 EST 2009
16 // $Id: HbbAnalyzerNew.cc,v 1.74 2012/05/17 10:58:34 arizzi Exp $
17 //
18 //
19
20
21 //uncomment to save also jet collections 1 and 4
22 //#define ENABLE_SIMPLEJETS1
23 //#define ENABLE_SIMPLEJETS4
24
25 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
26 #include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h"
27 #include "JetMETCorrections/Objects/interface/JetCorrector.h"
28 #include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"
29 #include "DataFormats/TrackReco/interface/TrackFwd.h"
30
31 #include "VHbbAnalysis/HbbAnalyzer/interface/HbbAnalyzerNew.h"
32 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
33
34 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
35 #include "DataFormats/Math/interface/deltaR.h"
36 #include "DataFormats/Math/interface/LorentzVector.h"
37 #include "DataFormats/Math/interface/Vector3D.h"
38 #include "Math/GenVector/PxPyPzM4D.h"
39
40 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
41 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
42 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
43 #include "TrackingTools/IPTools/interface/IPTools.h"
44
45 #include "CMGTools/External/interface/PileupJetIdentifier.h"
46 #include "CMGTools/External/interface/PileupJetIdAlgo.h"
47 //#include "CMGTools/External/interface/PileupJetIdProducer.h"
48
49 #include <cmath>
50
51
52
53
54 #define GENPTOLOR(a) TLorentzVector((a).px(), (a).py(), (a).pz(), (a).energy())
55 #define GENPTOLORP(a) TLorentzVector((a)->px(), (a)->py(), (a)->pz(), (a)->energy())
56
57
58
59 struct CompareJetPtMuons {
60 bool operator()( const VHbbEvent::MuonInfo& j1, const VHbbEvent::MuonInfo& j2 ) const {
61 return j1.p4.Pt() > j2.p4.Pt();
62 }
63 };
64 struct CompareJetPtElectrons {
65 bool operator()( const VHbbEvent::ElectronInfo& j1, const VHbbEvent::ElectronInfo& j2 ) const {
66 return j1.p4.Pt() > j2.p4.Pt();
67 }
68 };
69 struct CompareJetPtTaus {
70 bool operator()( const VHbbEvent::TauInfo& j1, const VHbbEvent::TauInfo& j2 ) const {
71 return j1.p4.Pt() > j2.p4.Pt();
72 }
73 };
74
75
76
77 HbbAnalyzerNew::HbbAnalyzerNew(const edm::ParameterSet& iConfig):
78 eleLabel_(iConfig.getParameter<edm::InputTag>("electronTag")),
79 muoLabel_(iConfig.getParameter<edm::InputTag>("muonTag")),
80 lep_ptCutForBjets_(iConfig.getParameter<double>("lep_ptCutForBjets")),
81 elenoCutsLabel_(iConfig.getParameter<edm::InputTag>("electronNoCutsTag")),
82 muonoCutsLabel_(iConfig.getParameter<edm::InputTag>("muonNoCutsTag")),
83 jetLabel_(iConfig.getParameter<edm::InputTag>("jetTag")),
84 subjetLabel_(iConfig.getParameter<edm::InputTag>("subjetTag")),
85 filterjetLabel_(iConfig.getParameter<edm::InputTag>("filterjetTag")),
86 simplejet1Label_(iConfig.getParameter<edm::InputTag>("simplejet1Tag")),
87 simplejet2Label_(iConfig.getParameter<edm::InputTag>("simplejet2Tag")),
88 simplejet3Label_(iConfig.getParameter<edm::InputTag>("simplejet3Tag")),
89 simplejet4Label_(iConfig.getParameter<edm::InputTag>("simplejet4Tag")),
90 tauLabel_(iConfig.getParameter<edm::InputTag>("tauTag")),
91 metLabel_(iConfig.getParameter<edm::InputTag>("metTag")),
92 phoLabel_(iConfig.getParameter<edm::InputTag>("photonTag")),
93 hltResults_(iConfig.getParameter<edm::InputTag>("hltResultsTag")),
94 runOnMC_(iConfig.getParameter<bool>("runOnMC")), verbose_(iConfig.getUntrackedParameter<bool>("verbose")) {
95
96 //
97 // put the setwhatproduced etc etc
98
99 produces<VHbbEvent>();
100 produces<VHbbEventAuxInfo>();
101
102
103 }
104
105
106 HbbAnalyzerNew::~HbbAnalyzerNew(){
107
108 // do anything here that needs to be done at desctruction time
109 // (e.g. close files, deallocate resources etc.)
110
111 }
112
113
114 //
115 // member functions
116 //
117
118 // ------------ method called to for each event ------------
119 void
120 HbbAnalyzerNew::produce(edm::Event& iEvent, const edm::EventSetup& iSetup){
121 using namespace edm;
122 using namespace reco;
123
124
125 // JEC Uncertainty
126
127 // JetCorrectionUncertainty *jecUnc=0;
128 edm::ESHandle<JetCorrectorParametersCollection> JetCorParColl;
129 iSetup.get<JetCorrectionsRecord>().get("AK5PFchs",JetCorParColl);
130 JetCorrectionUncertainty *jecUnc=0;
131 // if (!runOnMC_){
132 JetCorrectorParameters const & JetCorPar = (*JetCorParColl)["Uncertainty"];
133 jecUnc = new JetCorrectionUncertainty(JetCorPar);
134 // }
135
136 std::auto_ptr<VHbbEvent> hbbInfo( new VHbbEvent() );
137 std::auto_ptr<VHbbEventAuxInfo> auxInfo( new VHbbEventAuxInfo() );
138
139
140 if (runOnMC_){
141 Handle<GenEventInfoProduct> evt_info;
142 iEvent.getByType(evt_info);
143 auxInfo->weightMCProd = evt_info->weight();
144 }
145 else
146 { auxInfo->weightMCProd =1.;}
147 //
148 // ??
149
150 // trigger
151
152 // trigger
153 edm::Handle<edm::TriggerResults> hltresults;
154 //iEvent.getByLabel("TriggerResults", hltresults);
155
156 //edm::InputTag tag("TriggerResults::HLT");
157 // edm::InputTag tag("TriggerResults::HLT0");
158 iEvent.getByLabel(hltResults_, hltresults);
159
160 const edm::TriggerNames & triggerNames_ = iEvent.triggerNames(*hltresults);
161
162 int ntrigs = hltresults->size();
163 if (ntrigs==0){std::cerr << "%HLTInfo -- No trigger name given in TriggerResults of the input " << std::endl;}
164
165 BeamSpot vertexBeamSpot;
166 edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
167 iEvent.getByLabel("offlineBeamSpot",recoBeamSpotHandle);
168 vertexBeamSpot = *recoBeamSpotHandle;
169 /*
170 double BSx=vertexBeamSpot.x0();
171 double BSy=vertexBeamSpot.y0();
172 double BSz=vertexBeamSpot.z0();
173 */
174
175 double MinVtxProb=-999.;
176 int VtxIn=-99;
177
178 Handle<reco::VertexCollection> recVtxs;
179 iEvent.getByLabel("offlinePrimaryVertices", recVtxs);
180
181 auxInfo->pvInfo.nVertices = recVtxs->size();
182
183 for(size_t i = 0; i < recVtxs->size(); ++ i) {
184 const Vertex &vtx = (*recVtxs)[i];
185 double RecVtxProb=TMath::Prob(vtx.chi2(),vtx.ndof());
186 if(RecVtxProb>MinVtxProb){
187 VtxIn=i;
188 MinVtxProb=RecVtxProb;
189 }
190 }
191
192 const Vertex &RecVtx = (*recVtxs)[VtxIn];
193 const Vertex &RecVtxFirst = (*recVtxs)[0];
194 const Vertex &vertex = RecVtxFirst; //used in ele id 2012
195
196 auxInfo->pvInfo.firstPVInPT2 = TVector3(RecVtxFirst.x(), RecVtxFirst.y(), RecVtxFirst.z());
197 auxInfo->pvInfo.firstPVInProb = TVector3(RecVtx.x(), RecVtx.y(), RecVtx.z());
198
199 (auxInfo->pvInfo).efirstPVInPT2 = (RecVtxFirst.error());
200 (auxInfo->pvInfo).efirstPVInProb = RecVtx.error();
201
202 edm::Handle<double> rhoHandle;
203 iEvent.getByLabel(edm::InputTag("kt6PFJets", "rho"),rhoHandle);
204 auxInfo->puInfo.rho = *rhoHandle;
205
206 edm::Handle<double> rho25Handle;
207 iEvent.getByLabel(edm::InputTag("kt6PFJets25", "rho"),rho25Handle);
208 auxInfo->puInfo.rho25 = *rho25Handle;
209
210 edm::Handle<double> rhoNeutralHandle;
211 iEvent.getByLabel(edm::InputTag("kt6PFJetsCentralNeutral", "rho"),rhoNeutralHandle);
212 auxInfo->puInfo.rhoNeutral = *rhoNeutralHandle;
213
214
215 edm::Handle<std::vector< PileupSummaryInfo> > puHandle;
216
217 if (runOnMC_){
218 iEvent.getByType(puHandle);
219 if (puHandle.isValid()){
220
221 std::vector< PileupSummaryInfo> pu = (*puHandle);
222 for (std::vector<PileupSummaryInfo>::const_iterator it= pu.begin(); it!=pu.end(); ++it){
223 int bx = (*it).getBunchCrossing();
224 if(bx == 0) { auxInfo->puInfo.truePU = (*it).getTrueNumInteractions();}
225 unsigned int num = (*it).getPU_NumInteractions();
226 // std::cout <<" PU PUSHING "<<bx<<" " <<num<<std::endl;
227 auxInfo->puInfo.pus[bx] =num;
228 }
229 }
230 }
231
232 //// real start
233
234
235 Handle<GenParticleCollection> genParticles;
236
237 bool printJet=0;
238
239
240 if(runOnMC_){
241
242 iEvent.getByLabel("genParticles", genParticles);
243
244 for(size_t i = 0; i < genParticles->size(); ++ i) {
245
246 const GenParticle & p = (*genParticles)[i];
247 int id = p.pdgId();
248 int st = p.status();
249
250 if(id==25){
251
252 VHbbEventAuxInfo::ParticleMCInfo htemp;
253 htemp.status=st;
254 htemp.charge=p.charge();
255 if(p.mother(0)!=0) htemp.momid=p.mother(0)->pdgId();
256 if(p.mother(0)!=0 && p.mother(0)->mother(0)!=0) htemp.gmomid=p.mother(0)->mother(0)->pdgId();
257 htemp.p4 = GENPTOLOR(p);
258
259 int ndau = p.numberOfDaughters();
260 for(int j = 0; j < ndau; ++ j) {
261 const Candidate * Hdau = p.daughter( j );
262 htemp.dauid.push_back(Hdau->pdgId());
263 htemp.dauFourMomentum.push_back(GENPTOLORP(Hdau));
264 }
265 (auxInfo->mcH).push_back(htemp);
266 }
267
268
269 if(abs(id)==24){
270
271 VHbbEventAuxInfo::ParticleMCInfo wtemp;
272 wtemp.status=st;
273 wtemp.charge=p.charge();
274 if(p.mother(0)!=0) wtemp.momid=p.mother(0)->pdgId();
275 if(p.mother(0)!=0 && p.mother(0)->mother(0)!=0) wtemp.gmomid=p.mother(0)->mother(0)->pdgId();
276 wtemp.p4=GENPTOLOR(p);
277
278 int ndau = p.numberOfDaughters();
279 for(int j = 0; j < ndau; ++ j) {
280 const Candidate * Wdau = p.daughter( j );
281 wtemp.dauid.push_back(Wdau->pdgId());
282 wtemp.dauFourMomentum.push_back(GENPTOLORP(Wdau));
283 }
284 auxInfo->mcW.push_back(wtemp);
285 }
286
287 if(abs(id)==15) {
288 VHbbEventAuxInfo::ParticleMCInfo tautemp;
289 tautemp.status=st;
290 tautemp.charge=p.charge();
291 if(p.mother(0)!=0) tautemp.momid=p.mother(0)->pdgId();
292 if(p.mother(0)!=0 && p.mother(0)->mother(0)!=0) tautemp.gmomid=p.mother(0)->mother(0)->pdgId();
293 tautemp.p4=GENPTOLOR(p);
294
295 int ndau = p.numberOfDaughters();
296 for(int j = 0; j < ndau; ++ j) {
297 const Candidate * Taudau = p.daughter( j );
298 tautemp.dauid.push_back(Taudau->pdgId());
299 tautemp.dauFourMomentum.push_back(GENPTOLORP(Taudau));
300 }
301 auxInfo->mcTau.push_back(tautemp);
302 }
303
304 if(abs(id)==23){
305
306
307 VHbbEventAuxInfo::ParticleMCInfo ztemp;
308 ztemp.status=st;
309 ztemp.charge=p.charge();
310 if(p.mother(0)!=0) ztemp.momid=p.mother(0)->pdgId();
311 if(p.mother(0)!=0 && p.mother(0)->mother(0)!=0) ztemp.gmomid=p.mother(0)->mother(0)->pdgId();
312 ztemp.p4=GENPTOLOR(p);
313
314 int ndau = p.numberOfDaughters();
315 for(int j = 0; j < ndau; ++ j) {
316 const Candidate * Zdau = p.daughter( j );
317 ztemp.dauid.push_back(Zdau->pdgId());
318 ztemp.dauFourMomentum.push_back(GENPTOLORP(Zdau));
319 }
320 auxInfo->mcZ.push_back(ztemp);
321 }
322 //
323 // binfo
324 //
325
326
327 if(id==5){
328
329 VHbbEventAuxInfo::ParticleMCInfo btemp;
330 btemp.status=st;
331 btemp.charge=p.charge();
332 if(p.mother(0)!=0) btemp.momid=p.mother(0)->pdgId();
333 if(p.mother(0)!=0 && p.mother(0)->mother(0)!=0) btemp.gmomid=p.mother(0)->mother(0)->pdgId();
334
335 btemp.p4=GENPTOLOR(p);
336
337 int nHDaubdau = p.numberOfDaughters();
338 for(int j = 0; j < nHDaubdau; ++ j) {
339 const Candidate * Bdau = p.daughter( j );
340 btemp.dauid.push_back(Bdau->pdgId());
341 }
342 auxInfo->mcB.push_back(btemp);
343 }
344
345 if(id==-5){
346
347 VHbbEventAuxInfo::ParticleMCInfo bbtemp;
348
349 bbtemp.status=st;
350 bbtemp.charge=p.charge();
351 if(p.mother(0)!=0) bbtemp.momid=p.mother(0)->pdgId();
352 if(p.mother(0)!=0 && p.mother(0)->mother(0)!=0) bbtemp.gmomid=p.mother(0)->mother(0)->pdgId();
353
354 bbtemp.p4=GENPTOLOR(p);
355
356 int nHDaubdau = p.numberOfDaughters();
357 for(int j = 0; j < nHDaubdau; ++ j) {
358 const Candidate * Bdau = p.daughter( j );
359 bbtemp.dauid.push_back(Bdau->pdgId());
360 }
361
362
363 auxInfo->mcBbar.push_back(bbtemp);
364 }
365
366 if(abs(id)==4){
367 VHbbEventAuxInfo::ParticleMCInfo ctemp;
368 ctemp.status=st;
369 ctemp.charge=p.charge();
370 if(p.mother(0)!=0) ctemp.momid=p.mother(0)->pdgId();
371 if(p.mother(0)!=0 && p.mother(0)->mother(0)!=0) ctemp.gmomid=p.mother(0)->mother(0)->pdgId();
372
373 ctemp.p4=GENPTOLOR(p);
374
375 int nHDaubdau = p.numberOfDaughters();
376 for(int j = 0; j < nHDaubdau; ++ j) {
377 const Candidate * Bdau = p.daughter( j );
378 ctemp.dauid.push_back(Bdau->pdgId());
379 }
380
381 auxInfo->mcC.push_back(ctemp);
382
383 }
384
385 }
386
387 } // isMC
388
389 /////// end generator block
390
391 /// photon used in isolation
392 edm::Handle<edm::View<reco::PFCandidate> > photonIsoH;
393 iEvent.getByLabel("pfAllPhotons",photonIsoH);
394 edm::View<reco::PFCandidate> photonsForIso = *photonIsoH;
395
396 edm::Handle<edm::View<pat::Muon> > muonHandle;
397 iEvent.getByLabel(muoLabel_,muonHandle);
398 edm::View<pat::Muon> muons = *muonHandle;
399
400 // hard jet
401 edm::Handle<edm::View<pat::Jet> > jetHandle;
402 iEvent.getByLabel(jetLabel_,jetHandle);
403 edm::View<pat::Jet> jets = *jetHandle;
404
405 // sub jet
406 edm::Handle<edm::View<pat::Jet> > subjetHandle;
407 iEvent.getByLabel(subjetLabel_,subjetHandle);
408 edm::View<pat::Jet> subjets = *subjetHandle;
409
410 // filter jet
411 edm::Handle<edm::View<pat::Jet> > filterjetHandle;
412 iEvent.getByLabel(filterjetLabel_,filterjetHandle);
413 edm::View<pat::Jet> filterjets = *filterjetHandle;
414
415 // standard jets
416
417
418 edm::Handle<edm::View<pat::Jet> > simplejet2Handle;
419 iEvent.getByLabel(simplejet2Label_,simplejet2Handle);
420 edm::View<pat::Jet> simplejets2 = *simplejet2Handle;
421
422 edm::Handle<edm::View<pat::Jet> > simplejet3Handle;
423 iEvent.getByLabel(simplejet3Label_,simplejet3Handle);
424 edm::View<pat::Jet> simplejets3 = *simplejet3Handle;
425
426
427
428 edm::Handle<edm::View<pat::Electron> > electronHandle;
429 iEvent.getByLabel(eleLabel_,electronHandle);
430 edm::View<pat::Electron> electrons = *electronHandle;
431
432
433 // edm::Handle<edm::View<pat::Photon> > phoHandle;
434 // iEvent.getByLabel(phoLabel_,phoHandle);
435 // edm::View<pat::Photon> photons = *phoHandle;
436
437 edm::Handle<edm::View<pat::Tau> > tauHandle;
438 iEvent.getByLabel(tauLabel_,tauHandle);
439 edm::View<pat::Tau> taus = *tauHandle;
440
441 //Get the computer for the CSV
442 ESHandle<JetTagComputer> handle;
443 iSetup.get<JetTagComputerRecord>().get("combinedSecondaryVertex", handle);
444 computer = dynamic_cast<const GenericMVAJetTagComputer*>(handle.product());
445
446 //BTAGGING SCALE FACTOR FROM DATABASE
447 //Combined Secondary Vertex Loose
448 edm::ESHandle<BtagPerformance> bTagSF_CSVL_;
449 iSetup.get<BTagPerformanceRecord>().get("BTAGCSVL",bTagSF_CSVL_);
450 //Combined Secondary Vertex Medium
451 edm::ESHandle<BtagPerformance> bTagSF_CSVM_;
452 iSetup.get<BTagPerformanceRecord>().get("BTAGCSVM",bTagSF_CSVM_);
453 //Combined Secondary Vertex Tight
454 edm::ESHandle<BtagPerformance> bTagSF_CSVT_;
455 iSetup.get<BTagPerformanceRecord>().get("BTAGCSVT",bTagSF_CSVT_);
456
457 edm::ESHandle<BtagPerformance> mistagSF_CSVL_;
458 iSetup.get<BTagPerformanceRecord>().get("MISTAGCSVL",mistagSF_CSVL_);
459 //Combined Secondary Vertex Medium
460 edm::ESHandle<BtagPerformance> mistagSF_CSVM_;
461 iSetup.get<BTagPerformanceRecord>().get("MISTAGCSVM",mistagSF_CSVM_);
462 //Combined Secondary Vertex Tight
463 edm::ESHandle<BtagPerformance> mistagSF_CSVT_;
464 iSetup.get<BTagPerformanceRecord>().get("MISTAGCSVT",mistagSF_CSVT_);
465
466 BTagSFContainer btagSFs;
467 btagSFs.BTAGSF_CSVL = (bTagSF_CSVL_.product());
468 btagSFs.BTAGSF_CSVM = (bTagSF_CSVM_.product());
469 btagSFs.BTAGSF_CSVT = (bTagSF_CSVT_.product());
470 btagSFs.MISTAGSF_CSVL = (mistagSF_CSVL_.product());
471 btagSFs.MISTAGSF_CSVM = (mistagSF_CSVM_.product());
472 btagSFs.MISTAGSF_CSVT = (mistagSF_CSVT_.product());
473
474 #ifdef ENABLE_SIMPLEJETS1
475 edm::Handle<edm::View<pat::Jet> > simplejet1Handle;
476 iEvent.getByLabel(simplejet1Label_,simplejet1Handle);
477 edm::View<pat::Jet> simplejets1 = *simplejet1Handle;
478 for(edm::View<pat::Jet>::const_iterator jet_iter = simplejets1.begin(); jet_iter!=simplejets1.end(); ++jet_iter){
479 // if(jet_iter->pt()>50)
480 // njetscounter++;
481 VHbbEvent::SimpleJet sj;
482 // std::cout <<" sj1"<<std::endl;
483 fillSimpleJet(sj,jet_iter);
484 // if(!runOnMC_)
485
486 setJecUnc(sj,jecUnc);
487
488 Particle::LorentzVector p4Jet = jet_iter->p4();
489
490 if(runOnMC_){
491
492 fillScaleFactors(sj, btagSFs);
493
494 //PAT genJet matching
495 //genJet
496 const reco::GenJet *gJ = jet_iter->genJet();
497 //physical parton for mother info ONLY
498 if( (jet_iter->genParton()) ){
499 sj.bestMCid = jet_iter->genParton()->pdgId();
500 if( (jet_iter->genParton()->mother()) )
501 sj.bestMCmomid=jet_iter->genParton()->mother()->pdgId();
502 }
503 TLorentzVector gJp4;
504 if(gJ){
505 gJp4.SetPtEtaPhiE(gJ->pt(),gJ->eta(),gJ->phi(),gJ->energy());
506 sj.bestMCp4 = gJp4;
507 if(verbose_){
508 std::clog << "genJet matched Pt = " << gJp4.Pt() << std::endl;
509 std::clog << "genJet matched eta = " << gJp4.Eta() << std::endl;
510 std::clog << "genJet matched deltaR = " <<gJp4.DeltaR(sj.p4) << std::endl;
511 std::clog << "genJet matched mother id = " << sj.bestMCmomid << std::endl;
512 }
513 }
514
515 } //isMC
516 hbbInfo->simpleJets.push_back(sj);
517
518 }
519 #endif //ENABLE_SIMPLEJETS1
520
521 for(edm::View<pat::Jet>::const_iterator jet_iter = simplejets3.begin(); jet_iter!=simplejets3.end(); ++jet_iter){
522 // if(jet_iter->pt()>50)
523 // njetscounter++;
524 VHbbEvent::SimpleJet sj;
525 // std::cout <<" sj3"<<std::endl;
526 fillSimpleJet(sj,jet_iter);
527 // if(!runOnMC_)
528 setJecUnc(sj,jecUnc);
529
530 Particle::LorentzVector p4Jet = jet_iter->p4();
531
532 if(runOnMC_){
533
534 fillScaleFactors(sj, btagSFs);
535
536 //PAT genJet matching
537 //genJet
538 const reco::GenJet *gJ = jet_iter->genJet();
539 //physical parton for mother info ONLY
540 if( (jet_iter->genParton()) ){
541 sj.bestMCid = jet_iter->genParton()->pdgId();
542 if( (jet_iter->genParton()->mother()) )
543 sj.bestMCmomid=jet_iter->genParton()->mother()->pdgId();
544 }
545 TLorentzVector gJp4;
546 if(gJ){
547 gJp4.SetPtEtaPhiE(gJ->pt(),gJ->eta(),gJ->phi(),gJ->energy());
548 sj.bestMCp4 = gJp4;
549 if(verbose_){
550 std::clog << "genJet matched Pt = " << gJp4.Pt() << std::endl;
551 std::clog << "genJet matched eta = " << gJp4.Eta() << std::endl;
552 std::clog << "genJet matched deltaR = " <<gJp4.DeltaR(sj.p4) << std::endl;
553 std::clog << "genJet matched mother id = " << sj.bestMCmomid << std::endl;
554 }
555 }
556
557 } //isMC
558 //
559
560
561 hbbInfo->simpleJets3.push_back(sj);
562
563 }
564
565 #ifdef ENABLE_SIMPLEJETS4
566 edm::Handle<edm::View<pat::Jet> > simplejet4Handle;
567 iEvent.getByLabel(simplejet4Label_,simplejet4Handle);
568 edm::View<pat::Jet> simplejets4 = *simplejet4Handle;
569 for(edm::View<pat::Jet>::const_iterator jet_iter = simplejets4.begin(); jet_iter!=simplejets4.end(); ++jet_iter){
570 // if(jet_iter->pt()>50)
571 // njetscounter++;
572
573 VHbbEvent::SimpleJet sj;
574 // std::cout <<" sj4"<<std::endl;
575 fillSimpleJet(sj,jet_iter);
576 // if(!runOnMC_)
577 setJecUnc(sj,jecUnc);
578
579
580
581 Particle::LorentzVector p4Jet = jet_iter->p4();
582
583 if(runOnMC_){
584
585 fillScaleFactors(sj, btagSFs);
586
587 //PAT genJet matching
588 //genJet
589 const reco::GenJet *gJ = jet_iter->genJet();
590 //physical parton for mother info ONLY
591 if( (jet_iter->genParton()) ){
592 sj.bestMCid = jet_iter->genParton()->pdgId();
593 if( (jet_iter->genParton()->mother()) )
594 sj.bestMCmomid=jet_iter->genParton()->mother()->pdgId();
595 }
596 TLorentzVector gJp4;
597 if(gJ){
598 gJp4.SetPtEtaPhiE(gJ->pt(),gJ->eta(),gJ->phi(),gJ->energy());
599 sj.bestMCp4 = gJp4;
600 if(verbose_){
601 std::clog << "genJet matched Pt = " << gJp4.Pt() << std::endl;
602 std::clog << "genJet matched eta = " << gJp4.Eta() << std::endl;
603 std::clog << "genJet matched deltaR = " <<gJp4.DeltaR(sj.p4) << std::endl;
604 std::clog << "genJet matched mother id = " << sj.bestMCmomid << std::endl;
605 }
606 }
607
608 } //isMC
609 hbbInfo->simpleJets4.push_back(sj);
610
611 }
612 #endif //ENABLE SIMPLEJETS4
613
614
615 for(edm::View<pat::Jet>::const_iterator jet_iter = simplejets2.begin(); jet_iter!=simplejets2.end(); ++jet_iter){
616
617
618
619 VHbbEvent::SimpleJet sj;
620 // std::cout <<" sj2"<<std::endl;
621 fillSimpleJet(sj,jet_iter);
622
623 ///########### PU JET ID #################
624 // add puId...
625 edm::Handle<edm::ValueMap<float> > puJetIdMVA;
626 iEvent.getByLabel("puJetMva","fullDiscriminant", puJetIdMVA);
627
628 edm::Handle<edm::ValueMap<int> > puJetIdFlag;
629 iEvent.getByLabel("puJetMva", "fullId", puJetIdFlag);
630
631 // cout << " pt " << jet_iter->pt() << " eta " << jet_iter->eta() << std::endl;
632 unsigned int idx = jet_iter - simplejets2.begin();
633
634
635
636 sj.puJetIdMva = (*puJetIdMVA)[simplejets2.refAt(idx)];
637 int idflag = (*puJetIdFlag)[simplejets2.refAt(idx)];
638
639
640 // cout << " PU JetID MVA " << mva;
641 if( PileupJetIdentifier::passJetId( idflag, PileupJetIdentifier::kLoose )) {
642 //cout << " pass loose wp";
643 sj.puJetIdL =1;
644 }
645 if( PileupJetIdentifier::passJetId( idflag, PileupJetIdentifier::kMedium )) {
646 // cout << " pass medium wp";
647 sj.puJetIdM =1;
648 }
649 if( PileupJetIdentifier::passJetId( idflag, PileupJetIdentifier::kTight )) {
650 // cout << " pass tight wp";
651 sj.puJetIdT =1;
652 }
653 // cout << endl;
654 // ############# END OF PU JET ID ######################
655
656
657 // if(!runOnMC_)
658 setJecUnc(sj,jecUnc);
659 /* sj.flavour = jet_iter->partonFlavour();
660
661
662 sj.tche=jet_iter->bDiscriminator("trackCountingHighEffBJetTags");
663 sj.tchp=jet_iter->bDiscriminator("trackCountingHighPurBJetTags");
664 sj.jp=jet_iter->bDiscriminator("jetProbabilityBJetTags");
665 sj.jpb=jet_iter->bDiscriminator("jetBProbabilityBJetTags");
666 sj.ssvhe=jet_iter->bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
667 sj.csv=jet_iter->bDiscriminator("combinedSecondaryVertexBJetTags");
668 sj.csvmva=jet_iter->bDiscriminator("combinedSecondaryVertexMVABJetTags");
669 sj.charge=jet_iter->jetCharge();
670 sj.ntracks=jet_iter->associatedTracks().size();
671 sj.p4=GENPTOLORP(jet_iter);
672 sj.chargedTracksFourMomentum=(getChargedTracksMomentum(&*(jet_iter)));
673 sj.SF_CSVL=1;
674 sj.SF_CSVM=1;
675 sj.SF_CSVT=1;
676 sj.SF_CSVLerr=0;
677 sj.SF_CSVMerr=0;
678 sj.SF_CSVTerr=0;
679
680
681 //
682 // addtaginfo for csv
683 //
684
685 if (jet_iter->hasTagInfo("SimpleSecondaryVertex")) {
686
687 const reco::SecondaryVertexTagInfo * tf = jet_iter->tagInfoSecondaryVertex();
688 sj.vtxMass = tf->secondaryVertex(0).p4().mass();
689 sj.vtxNTracks = tf->secondaryVertex(0).nTracks();
690 Measurement1D m = tf->flightDistance(0);
691 sj.vtx3dL = m.value();
692 sj.vtx3deL = m.error();
693 }
694
695
696 //
697 // add tVector
698 //
699 sj.tVector = getTvect(&(*jet_iter));
700 */
701 Particle::LorentzVector p4Jet = jet_iter->p4();
702
703 if(runOnMC_){
704
705 //BTV scale factors
706 fillScaleFactors(sj, btagSFs);
707
708 //PAT genJet matching
709 //genJet
710 const reco::GenJet *gJ = jet_iter->genJet();
711 //physical parton for mother info ONLY
712 if( (jet_iter->genParton()) ){
713 sj.bestMCid = jet_iter->genParton()->pdgId();
714 if( (jet_iter->genParton()->mother()) )
715 sj.bestMCmomid=jet_iter->genParton()->mother()->pdgId();
716 }
717 TLorentzVector gJp4;
718 if(gJ){
719 gJp4.SetPtEtaPhiE(gJ->pt(),gJ->eta(),gJ->phi(),gJ->energy());
720 sj.bestMCp4 = gJp4;
721 if(verbose_){
722 std::clog << "genJet matched Pt = " << gJp4.Pt() << std::endl;
723 std::clog << "genJet matched eta = " << gJp4.Eta() << std::endl;
724 std::clog << "genJet matched deltaR = " << gJp4.DeltaR(sj.p4) << std::endl;
725 std::clog << "genJet matched mother id = " << sj.bestMCmomid << std::endl;
726 }
727 }
728
729 // add flag if a mc lepton is find inside a cone around the jets...
730 iEvent.getByLabel("genParticles", genParticles);
731
732 for(size_t i = 0; i < genParticles->size(); ++ i) {
733
734 const GenParticle & p = (*genParticles)[i];
735 int id = 0;
736 p.pt()> lep_ptCutForBjets_ ? id= p.pdgId(): 0;
737
738 // std::cout<< "found a muon with pt " << mu->pt() << std::endl;
739 if ((abs(id)==13 || abs(id)==11) && deltaR(p.eta(), p.phi(), sj.p4.Eta(), sj.p4.Phi() ) <0.5) sj.isSemiLeptMCtruth=1;
740 }
741
742 } //isMC
743
744 // add flag if a reco lepton is find inside a cone around the jets...
745 edm::Handle<edm::View<reco::Candidate> > muonNoCutsHandle;
746 iEvent.getByLabel(muonoCutsLabel_,muonNoCutsHandle);
747 edm::View<reco::Candidate> muonsNoCuts = *muonNoCutsHandle;
748
749
750
751 for(edm::View<reco::Candidate>::const_iterator mu = muonsNoCuts.begin(); mu!=muonsNoCuts.end() && sj.isSemiLept!=1; ++mu){
752 // std::cout<< "found a muon with pt " << mu->pt() << std::endl;
753 const pat::Muon& m = static_cast <const pat::Muon&> (*mu);
754 float Smpt = m.pt();
755 float Smeta = m.eta();
756 float Smphi = m.phi();
757
758 float SmJdR = deltaR(Smeta, Smphi, sj.p4.Eta(), sj.p4.Phi());
759
760 if ( Smpt> lep_ptCutForBjets_ && SmJdR <0.5) {
761 sj.isSemiLept=1;
762 //isSemiLept(-99), isSemiLeptMCtruth(-99), SoftLeptPt(-99), SoftLeptdR(-99), SoftLeptptRel(-99), SoftLeptpdgId(-99), SoftLeptIdlooseMu(-99), SoftLeptId95(-99), SoftLeptRelCombIso(-99),
763 sj.SoftLeptpdgId =13;
764 sj.SoftLeptdR= SmJdR;
765 sj.SoftLeptPt=Smpt;
766 TVector3 mvec ( m.p4().Vect().X(), m.p4().Vect().Y(), m.p4().Vect().Z() );
767 sj.SoftLeptptRel= sj.p4.Perp( mvec );
768 sj.SoftLeptRelCombIso = (m.trackIso() + m.ecalIso() + m.hcalIso() ) / Smpt ;
769 sj.SoftLeptIdlooseMu=m.muonID("TMLastStationLoose");
770 }
771 }
772
773
774 edm::Handle<edm::View<reco::Candidate> > eleNoCutsHandle;
775 iEvent.getByLabel(elenoCutsLabel_,eleNoCutsHandle);
776 edm::View<reco::Candidate> elesNoCuts = *eleNoCutsHandle;
777
778
779
780 for(edm::View<reco::Candidate>::const_iterator ele = elesNoCuts.begin(); ele!=elesNoCuts.end() && sj.isSemiLept!=1; ++ele){
781
782 const pat::Electron& e = static_cast <const pat::Electron&> (*ele);
783 float Smpt = e.pt();
784 float Smeta = e.eta();
785 float Smphi = e.phi();
786
787 float SmJdR = deltaR(Smeta, Smphi, sj.p4.Eta(), sj.p4.Phi());
788 if ( Smpt> lep_ptCutForBjets_ && SmJdR <0.5) {
789 sj.isSemiLept=1;
790 sj.SoftLeptpdgId =11;
791 sj.SoftLeptdR= SmJdR;
792 sj.SoftLeptPt=Smpt;
793 TVector3 mvec ( e.p4().Vect().X(), e.p4().Vect().Y(), e.p4().Vect().Z() );
794 sj.SoftLeptptRel= sj.p4.Perp( mvec );
795 sj.SoftLeptRelCombIso = (e.trackIso() + e.ecalIso() + e.hcalIso() ) / Smpt ;
796 // sj.SoftLeptId95=e.electronID("eidVBTFCom95");
797 //std::cout << "before ele id " << std::endl;
798 // std::cout << " e.e.sigmaIetaIeta " << e.sigmaIetaIeta() << std::endl;
799 //std::cout << " e.isEB() " << e.isEB() << std::endl;
800 if (
801 ( fabs(Smeta)<2.5 && !( abs(Smeta)>1.4442 && abs(Smeta)<1.566)) &&
802
803 (( abs(Smeta)>1.566 && (e.sigmaIetaIeta()<0.01) && ( e.deltaPhiSuperClusterTrackAtVtx()<0.8 && e.deltaPhiSuperClusterTrackAtVtx()>-0.8) && ( e.deltaEtaSuperClusterTrackAtVtx()<0.007 && e.deltaEtaSuperClusterTrackAtVtx()>-0.007 ) )
804 || ( abs(Smeta)<1.4442 && (e.sigmaIetaIeta()<0.03) && ( e.deltaPhiSuperClusterTrackAtVtx()<0.7 && e.deltaPhiSuperClusterTrackAtVtx()>-0.7 ) && ( e.deltaEtaSuperClusterTrackAtVtx()<0.01 && e.deltaEtaSuperClusterTrackAtVtx()>-0.01 ) ))
805 )
806 sj.SoftLeptId95=1;
807 }
808 }
809
810
811
812
813
814 hbbInfo->simpleJets2.push_back(sj);
815
816 }
817
818
819 /* const GenJet* jet1Mc = bjet1.genJet();
820 const GenJet* jet2Mc = bjet2.genJet();
821 if(jet1Mc!=0){
822 MCbJet1MomId=jet1Mc->mother()->pdgId();
823 MCbJet1GMomId=jet1Mc->mother()->mother()->pdgId();
824 }
825
826 if(jet2Mc!=0){
827 MCbJet2MomId=jet2Mc->mother()->pdgId();
828 MCbJet2GMomId=jet2Mc->mother()->mother()->pdgId();
829 }
830 */
831
832
833
834 /////// hard jet
835
836
837 double matEta[1000*30],matPhi[1000*30];
838 for(int i=0;i<1000;i++){for(int j=0;j<30;j++){matEta[i*j]=-99.;matPhi[i*j]=-99.;}}
839
840 for(edm::View<pat::Jet>::const_iterator jet_iter = jets.begin(); jet_iter!=jets.end(); ++jet_iter){
841
842 if(printJet) {std::cout << "Jet Pt: " << jet_iter->pt() << " E,M: " << jet_iter->p4().E() << " " << jet_iter->p4().M() << "\n";}
843
844 reco::Jet::Constituents constituents = jet_iter->getJetConstituents();
845
846 // if(printJet) {std::cout << "NsubJets: " << constituents.size() << "\n";}
847 VHbbEvent::HardJet hj;
848 hj.constituents=constituents.size();
849 hj.p4 =GENPTOLORP(jet_iter);
850
851 for (unsigned int iJC(0); iJC<constituents.size(); ++iJC ){
852 Jet::Constituent icandJet = constituents[iJC];
853
854 if(printJet) {std::cout << "subJet Pt: " << icandJet->pt() << " subJet E,M,eta,phi: " << icandJet->p4().E() << ","
855 << icandJet->p4().M() << "," << icandJet->eta() << "," << icandJet->phi() << "\n"; }
856
857
858 hj.subFourMomentum.push_back(GENPTOLORP(icandJet));
859 hj.etaSub.push_back(icandJet->eta());
860 hj.phiSub.push_back(icandJet->phi());
861
862 }
863 hbbInfo->hardJets.push_back(hj);
864 }
865
866 // HardJetSubEta2.SetMatrixArray(matEta);
867 // HardJetSubPhi2.SetMatrixArray(matPhi);
868 // TMatrixDRow a1(HardJetSubEta2,0);
869 // for(int i=0;i<30;i++){
870 // std::cout << "test: " << a1[i] << "\n";
871 // }
872
873 // pat subJets with Btag
874
875
876 for(edm::View<pat::Jet>::const_iterator subjet_iter = subjets.begin(); subjet_iter!=subjets.end(); ++subjet_iter){
877
878 if(printJet) {std::cout << "SubJetTagged Pt: " << subjet_iter->pt() << " E,M,eta,phi,Btag: " << subjet_iter->p4().E()
879 << "," << subjet_iter->p4().M() << "," << subjet_iter->eta() << "," << subjet_iter->phi()
880 << "," << subjet_iter->bDiscriminator("combinedSecondaryVertexBJetTags") << "\n";}
881
882 VHbbEvent::SimpleJet sj;
883 // std::cout <<" sub jet "<<std::endl;
884 fillSimpleJet(sj,subjet_iter);
885 // if(!runOnMC_)
886 setJecUnc(sj,jecUnc);
887 /* sj.flavour = subjet_iter->partonFlavour();
888 sj.tVector = getTvect(&(*subjet_iter));
889 sj.tche=subjet_iter->bDiscriminator("trackCountingHighEffBJetTags");
890 sj.tchp=subjet_iter->bDiscriminator("trackCountingHighPurBJetTags");
891 sj.jp=subjet_iter->bDiscriminator("jetProbabilityBJetTags");
892 sj.jpb=subjet_iter->bDiscriminator("jetBProbabilityBJetTags");
893 sj.ssvhe=subjet_iter->bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
894 sj.csv=subjet_iter->bDiscriminator("combinedSecondaryVertexBJetTags");
895 sj.csvmva=subjet_iter->bDiscriminator("combinedSecondaryVertexMVABJetTags");
896 sj.charge=subjet_iter->jetCharge();
897 sj.ntracks=subjet_iter->associatedTracks().size();
898 sj.p4=GENPTOLORP(subjet_iter);
899 sj.p4=(getChargedTracksMomentum(&*(subjet_iter)));
900
901 //
902 // addtaginfo for csv
903 //
904
905 if (subjet_iter->hasTagInfo("SimpleSecondaryVertex")) {
906
907 const reco::SecondaryVertexTagInfo * tf = subjet_iter->tagInfoSecondaryVertex();
908 sj.vtxMass = tf->secondaryVertex(0).p4().mass();
909 sj.vtxNTracks = tf->secondaryVertex(0).nTracks();
910 Measurement1D m = tf->flightDistance(0);
911 sj.vtx3dL = m.value();
912 sj.vtx3deL = m.error();
913 }
914 */
915 hbbInfo->subJets.push_back(sj);
916
917 }
918
919 for(edm::View<pat::Jet>::const_iterator filterjet_iter = filterjets.begin(); filterjet_iter!=filterjets.end(); ++filterjet_iter){
920
921 if(printJet) {std::cout << "FilterjetTagged Pt: " << filterjet_iter->pt() << " E,M,eta,phi,Btag: " << filterjet_iter->p4().E() << "," << filterjet_iter->p4().M() << "," << filterjet_iter->eta() << "," << filterjet_iter->phi() << "," << filterjet_iter->bDiscriminator("combinedSecondaryVertexBJetTags") << "\n";}
922
923 VHbbEvent::SimpleJet fj;
924 // std::cout <<" sub jet "<<std::endl;
925 fillSimpleJet(fj,filterjet_iter);
926 // if(!runOnMC_)
927 setJecUnc(fj,jecUnc);
928
929 hbbInfo->filterJets.push_back(fj);
930
931
932 }
933
934 //
935 // add charged met
936 //
937
938 edm::Handle<edm::View<reco::MET> > metChargedHandle;
939 iEvent.getByLabel("pfMETNoPUCharge",metChargedHandle);
940 edm::View<reco::MET> metsCh = *metChargedHandle;
941 if(metsCh.size()){
942 hbbInfo->metCh.sumEt=(metsCh[0]).sumEt();
943 hbbInfo->metCh.metSig=metSignificance(& (metsCh[0]));
944 hbbInfo->metCh.eLong=(metsCh[0]).e_longitudinal();
945 hbbInfo->metCh.p4=GENPTOLOR((metsCh[0]));
946 if (verbose_) std::cout <<" METCharged "<< hbbInfo->metCh.metSig <<" " << hbbInfo->metCh.sumEt<<std::endl;
947 }
948
949 // type 1 corr met
950 edm::Handle<edm::View<reco::MET> > pfmetType1corrHandle;
951 iEvent.getByLabel("patType1CorrectedPFMet",pfmetType1corrHandle);
952 edm::View<reco::MET> pfmetsType1corr = *pfmetType1corrHandle;
953 if(pfmetsType1corr.size()){
954 hbbInfo->pfmetType1corr.sumEt=(pfmetsType1corr[0]).sumEt();
955 hbbInfo->pfmetType1corr.metSig=metSignificance(& (pfmetsType1corr[0]));
956 hbbInfo->pfmetType1corr.eLong=(pfmetsType1corr[0]).e_longitudinal();
957 hbbInfo->pfmetType1corr.p4=GENPTOLOR((pfmetsType1corr[0]));
958 if (verbose_) std::cout <<" type 1 corrected pfMET "<< hbbInfo->pfmetType1corr.metSig <<" " << hbbInfo->pfmetType1corr.sumEt<<std::endl;
959 }
960
961
962 // type 1 + 2 corr met
963 edm::Handle<edm::View<reco::MET> > pfmetType1p2corrHandle;
964 iEvent.getByLabel("patType1p2CorrectedPFMet",pfmetType1p2corrHandle);
965 edm::View<reco::MET> pfmetsType1p2corr = *pfmetType1p2corrHandle;
966 if(pfmetsType1p2corr.size()){
967 hbbInfo->pfmetType1p2corr.sumEt=(pfmetsType1p2corr[0]).sumEt();
968 hbbInfo->pfmetType1p2corr.metSig=metSignificance(& (pfmetsType1p2corr[0]));
969 hbbInfo->pfmetType1p2corr.eLong=(pfmetsType1p2corr[0]).e_longitudinal();
970 hbbInfo->pfmetType1p2corr.p4=GENPTOLOR((pfmetsType1p2corr[0]));
971 if (verbose_) std::cout <<" type 1 +2 corrected pfMET "<< hbbInfo->pfmetType1p2corr.metSig <<" " << hbbInfo->pfmetType1p2corr.sumEt<<std::endl;
972 }
973
974 // type 1 corr met NoPU
975 /* edm::Handle<edm::View<reco::MET> > pfmetNoPUType1corrHandle;
976 iEvent.getByLabel("patType1CorrectedPFMetNoPU",pfmetNoPUType1corrHandle);
977 edm::View<reco::MET> pfmetsNoPUType1corr = *pfmetNoPUType1corrHandle;
978 if(pfmetsNoPUType1corr.size()){
979 hbbInfo->pfmetNoPUType1corr.sumEt=(pfmetsNoPUType1corr[0]).sumEt();
980 hbbInfo->pfmetNoPUType1corr.metSig=metSignificance(& (pfmetsNoPUType1corr[0]));
981 hbbInfo->pfmetNoPUType1corr.eLong=(pfmetsNoPUType1corr[0]).e_longitudinal();
982 hbbInfo->pfmetNoPUType1corr.p4=GENPTOLOR((pfmetsNoPUType1corr[0]));
983 if (verbose_) std::cout <<" type 1 corrected pfMET NoPU"<< hbbInfo->pfmetNoPUType1corr.metSig <<" " << hbbInfo->pfmetNoPUType1corr.sumEt<<std::endl;
984 }
985
986
987 // type 1 + 2 corr met
988 edm::Handle<edm::View<reco::MET> > pfmetNoPUType1p2corrHandle;
989 iEvent.getByLabel("patType1p2CorrectedPFMetNoPU",pfmetNoPUType1p2corrHandle);
990 edm::View<reco::MET> pfmetsNoPUType1p2corr = *pfmetNoPUType1p2corrHandle;
991 if(pfmetsNoPUType1p2corr.size()){
992 hbbInfo->pfmetNoPUType1p2corr.sumEt=(pfmetsNoPUType1p2corr[0]).sumEt();
993 hbbInfo->pfmetNoPUType1p2corr.metSig=metSignificance(& (pfmetsNoPUType1p2corr[0]));
994 hbbInfo->pfmetNoPUType1p2corr.eLong=(pfmetsNoPUType1p2corr[0]).e_longitudinal();
995 hbbInfo->pfmetNoPUType1p2corr.p4=GENPTOLOR((pfmetsNoPUType1p2corr[0]));
996 if (verbose_) std::cout <<" type 1 +2 corrected pfMET "<< hbbInfo->pfmetNoPUType1p2corr.metSig <<" " << hbbInfo->pfmetNoPUType1p2corr.sumEt<<std::endl;
997 }
998
999 */
1000
1001 /*
1002 // MET uncertainty vector
1003 vector<pat::MET> "patType1CorrectedPFMet" "" "VH"
1004 vector<pat::MET> "patType1CorrectedPFMetElectronEnDown" "" "VH"
1005 vector<pat::MET> "patType1CorrectedPFMetElectronEnUp" "" "VH"
1006 vector<pat::MET> "patType1CorrectedPFMetJetEnDown" "" "VH"
1007 vector<pat::MET> "patType1CorrectedPFMetJetEnUp" "" "VH"
1008 vector<pat::MET> "patType1CorrectedPFMetJetResDown" "" "VH"
1009 vector<pat::MET> "patType1CorrectedPFMetJetResUp" "" "VH"
1010 vector<pat::MET> "patType1CorrectedPFMetMuonEnDown" "" "VH"
1011 vector<pat::MET> "patType1CorrectedPFMetMuonEnUp" "" "VH"
1012 vector<pat::MET> "patType1CorrectedPFMetNoPU" "" "VH"
1013 vector<pat::MET> "patType1CorrectedPFMetTauEnDown" "" "VH"
1014 vector<pat::MET> "patType1CorrectedPFMetTauEnUp" "" "VH"
1015 vector<pat::MET> "patType1CorrectedPFMetUnclusteredEnDown" "" "VH"
1016 vector<pat::MET> "patType1CorrectedPFMetUnclusteredEnUp" "" "VH"
1017 vector<pat::MET> "patType1p2CorrectedPFMet" "" "VH"
1018 vector<pat::MET> "patType1p2CorrectedPFMetElectronEnDown" "" "VH"
1019 vector<pat::MET> "patType1p2CorrectedPFMetElectronEnUp" "" "VH"
1020 vector<pat::MET> "patType1p2CorrectedPFMetJetEnDown" "" "VH"
1021 vector<pat::MET> "patType1p2CorrectedPFMetJetEnUp" "" "VH"
1022 vector<pat::MET> "patType1p2CorrectedPFMetJetResDown" "" "VH"
1023 vector<pat::MET> "patType1p2CorrectedPFMetJetResUp" "" "VH"
1024 vector<pat::MET> "patType1p2CorrectedPFMetMuonEnDown" "" "VH"
1025 vector<pat::MET> "patType1p2CorrectedPFMetMuonEnUp" "" "VH"
1026 vector<pat::MET> "patType1p2CorrectedPFMetNoPU" "" "VH"
1027 vector<pat::MET> "patType1p2CorrectedPFMetTauEnDown" "" "VH"
1028 vector<pat::MET> "patType1p2CorrectedPFMetTauEnUp" "" "VH"
1029 vector<pat::MET> "patType1p2CorrectedPFMetUnclusteredEnDown" "" "VH"
1030 vector<pat::MET> "patType1p2CorrectedPFMetUnclusteredEnUp" "" "VH"
1031 */
1032
1033 VHbbEvent::METInfo metunc;
1034 edm::Handle<edm::View<reco::MET> > patType1CorrectedPFMetElectronEnDownHandle;
1035 iEvent.getByLabel("patType1CorrectedPFMetElectronEnDown",patType1CorrectedPFMetElectronEnDownHandle);
1036 edm::View<reco::MET> patType1CorrectedPFMetsElectronEnDown = *patType1CorrectedPFMetElectronEnDownHandle;
1037 if(patType1CorrectedPFMetsElectronEnDown.size()){
1038 metunc.sumEt =(patType1CorrectedPFMetsElectronEnDown[0]).sumEt();
1039 metunc.metSig=metSignificance(& (patType1CorrectedPFMetsElectronEnDown[0]));
1040 metunc.eLong=(patType1CorrectedPFMetsElectronEnDown[0]).e_longitudinal();
1041 metunc.p4=GENPTOLOR((patType1CorrectedPFMetsElectronEnDown[0]));
1042 hbbInfo->metUncInfo.push_back(metunc);
1043 }
1044
1045 edm::Handle<edm::View<reco::MET> > patType1CorrectedPFMetElectronEnUpHandle;
1046 iEvent.getByLabel("patType1CorrectedPFMetElectronEnUp",patType1CorrectedPFMetElectronEnUpHandle);
1047 edm::View<reco::MET> patType1CorrectedPFMetsElectronEnUp = *patType1CorrectedPFMetElectronEnUpHandle;
1048 if(patType1CorrectedPFMetsElectronEnUp.size()){
1049 metunc.sumEt =(patType1CorrectedPFMetsElectronEnUp[0]).sumEt();
1050 metunc.metSig=metSignificance(& (patType1CorrectedPFMetsElectronEnUp[0]));
1051 metunc.eLong=(patType1CorrectedPFMetsElectronEnUp[0]).e_longitudinal();
1052 metunc.p4=GENPTOLOR((patType1CorrectedPFMetsElectronEnUp[0]));
1053 hbbInfo->metUncInfo.push_back(metunc);
1054 }
1055
1056
1057
1058 edm::Handle<edm::View<reco::MET> > patType1CorrectedPFMetMuonEnDownHandle;
1059 iEvent.getByLabel("patType1CorrectedPFMetMuonEnDown",patType1CorrectedPFMetMuonEnDownHandle);
1060 edm::View<reco::MET> patType1CorrectedPFMetsMuonEnDown = *patType1CorrectedPFMetMuonEnDownHandle;
1061 if(patType1CorrectedPFMetsMuonEnDown.size()){
1062 metunc.sumEt =(patType1CorrectedPFMetsMuonEnDown[0]).sumEt();
1063 metunc.metSig=metSignificance(& (patType1CorrectedPFMetsMuonEnDown[0]));
1064 metunc.eLong=(patType1CorrectedPFMetsMuonEnDown[0]).e_longitudinal();
1065 metunc.p4=GENPTOLOR((patType1CorrectedPFMetsMuonEnDown[0]));
1066 hbbInfo->metUncInfo.push_back(metunc);
1067 }
1068
1069 edm::Handle<edm::View<reco::MET> > patType1CorrectedPFMetMuonEnUpHandle;
1070 iEvent.getByLabel("patType1CorrectedPFMetMuonEnUp",patType1CorrectedPFMetMuonEnUpHandle);
1071 edm::View<reco::MET> patType1CorrectedPFMetsMuonEnUp = *patType1CorrectedPFMetMuonEnUpHandle;
1072 if(patType1CorrectedPFMetsMuonEnUp.size()){
1073 metunc.sumEt =(patType1CorrectedPFMetsMuonEnUp[0]).sumEt();
1074 metunc.metSig=metSignificance(& (patType1CorrectedPFMetsMuonEnUp[0]));
1075 metunc.eLong=(patType1CorrectedPFMetsMuonEnUp[0]).e_longitudinal();
1076 metunc.p4=GENPTOLOR((patType1CorrectedPFMetsMuonEnUp[0]));
1077 hbbInfo->metUncInfo.push_back(metunc);
1078 }
1079
1080
1081
1082 edm::Handle<edm::View<reco::MET> > patType1CorrectedPFMetTauEnDownHandle;
1083 iEvent.getByLabel("patType1CorrectedPFMetTauEnDown",patType1CorrectedPFMetTauEnDownHandle);
1084 edm::View<reco::MET> patType1CorrectedPFMetsTauEnDown = *patType1CorrectedPFMetTauEnDownHandle;
1085 if(patType1CorrectedPFMetsTauEnDown.size()){
1086 metunc.sumEt =(patType1CorrectedPFMetsTauEnDown[0]).sumEt();
1087 metunc.metSig=metSignificance(& (patType1CorrectedPFMetsTauEnDown[0]));
1088 metunc.eLong=(patType1CorrectedPFMetsTauEnDown[0]).e_longitudinal();
1089 metunc.p4=GENPTOLOR((patType1CorrectedPFMetsTauEnDown[0]));
1090 hbbInfo->metUncInfo.push_back(metunc);
1091 }
1092
1093 edm::Handle<edm::View<reco::MET> > patType1CorrectedPFMetTauEnUpHandle;
1094 iEvent.getByLabel("patType1CorrectedPFMetTauEnUp",patType1CorrectedPFMetTauEnUpHandle);
1095 edm::View<reco::MET> patType1CorrectedPFMetsTauEnUp = *patType1CorrectedPFMetTauEnUpHandle;
1096 if(patType1CorrectedPFMetsTauEnUp.size()){
1097 metunc.sumEt =(patType1CorrectedPFMetsTauEnUp[0]).sumEt();
1098 metunc.metSig=metSignificance(& (patType1CorrectedPFMetsTauEnUp[0]));
1099 metunc.eLong=(patType1CorrectedPFMetsTauEnUp[0]).e_longitudinal();
1100 metunc.p4=GENPTOLOR((patType1CorrectedPFMetsTauEnUp[0]));
1101 hbbInfo->metUncInfo.push_back(metunc);
1102 }
1103
1104
1105 edm::Handle<edm::View<reco::MET> > patType1CorrectedPFMetJetEnDownHandle;
1106 iEvent.getByLabel("patType1CorrectedPFMetJetEnDown",patType1CorrectedPFMetJetEnDownHandle);
1107 edm::View<reco::MET> patType1CorrectedPFMetsJetEnDown = *patType1CorrectedPFMetJetEnDownHandle;
1108 if(patType1CorrectedPFMetsJetEnDown.size()){
1109 metunc.sumEt =(patType1CorrectedPFMetsJetEnDown[0]).sumEt();
1110 metunc.metSig=metSignificance(& (patType1CorrectedPFMetsJetEnDown[0]));
1111 metunc.eLong=(patType1CorrectedPFMetsJetEnDown[0]).e_longitudinal();
1112 metunc.p4=GENPTOLOR((patType1CorrectedPFMetsJetEnDown[0]));
1113 hbbInfo->metUncInfo.push_back(metunc);
1114 }
1115
1116 edm::Handle<edm::View<reco::MET> > patType1CorrectedPFMetJetEnUpHandle;
1117 iEvent.getByLabel("patType1CorrectedPFMetJetEnUp",patType1CorrectedPFMetJetEnUpHandle);
1118 edm::View<reco::MET> patType1CorrectedPFMetsJetEnUp = *patType1CorrectedPFMetJetEnUpHandle;
1119 if(patType1CorrectedPFMetsJetEnUp.size()){
1120 metunc.sumEt =(patType1CorrectedPFMetsJetEnUp[0]).sumEt();
1121 metunc.metSig=metSignificance(& (patType1CorrectedPFMetsJetEnUp[0]));
1122 metunc.eLong=(patType1CorrectedPFMetsJetEnUp[0]).e_longitudinal();
1123 metunc.p4=GENPTOLOR((patType1CorrectedPFMetsJetEnUp[0]));
1124 hbbInfo->metUncInfo.push_back(metunc);
1125 }
1126
1127
1128 edm::Handle<edm::View<reco::MET> > patType1CorrectedPFMetJetResDownHandle;
1129 iEvent.getByLabel("patType1CorrectedPFMetJetResDown",patType1CorrectedPFMetJetResDownHandle);
1130 edm::View<reco::MET> patType1CorrectedPFMetsJetResDown = *patType1CorrectedPFMetJetResDownHandle;
1131 if(patType1CorrectedPFMetsJetResDown.size()){
1132 metunc.sumEt =(patType1CorrectedPFMetsJetResDown[0]).sumEt();
1133 metunc.metSig=metSignificance(& (patType1CorrectedPFMetsJetResDown[0]));
1134 metunc.eLong=(patType1CorrectedPFMetsJetResDown[0]).e_longitudinal();
1135 metunc.p4=GENPTOLOR((patType1CorrectedPFMetsJetResDown[0]));
1136 hbbInfo->metUncInfo.push_back(metunc);
1137 }
1138
1139 edm::Handle<edm::View<reco::MET> > patType1CorrectedPFMetJetResUpHandle;
1140 iEvent.getByLabel("patType1CorrectedPFMetJetResUp",patType1CorrectedPFMetJetResUpHandle);
1141 edm::View<reco::MET> patType1CorrectedPFMetsJetResUp = *patType1CorrectedPFMetJetResUpHandle;
1142 if(patType1CorrectedPFMetsJetResUp.size()){
1143 metunc.sumEt =(patType1CorrectedPFMetsJetResUp[0]).sumEt();
1144 metunc.metSig=metSignificance(& (patType1CorrectedPFMetsJetResUp[0]));
1145 metunc.eLong=(patType1CorrectedPFMetsJetResUp[0]).e_longitudinal();
1146 metunc.p4=GENPTOLOR((patType1CorrectedPFMetsJetResUp[0]));
1147 hbbInfo->metUncInfo.push_back(metunc);
1148 }
1149
1150
1151 edm::Handle<edm::View<reco::MET> > patType1CorrectedPFMetUnclusteredEnDownHandle;
1152 iEvent.getByLabel("patType1CorrectedPFMetUnclusteredEnDown",patType1CorrectedPFMetUnclusteredEnDownHandle);
1153 edm::View<reco::MET> patType1CorrectedPFMetsUnclusteredEnDown = *patType1CorrectedPFMetUnclusteredEnDownHandle;
1154 if(patType1CorrectedPFMetsUnclusteredEnDown.size()){
1155 metunc.sumEt =(patType1CorrectedPFMetsUnclusteredEnDown[0]).sumEt();
1156 metunc.metSig=metSignificance(& (patType1CorrectedPFMetsUnclusteredEnDown[0]));
1157 metunc.eLong=(patType1CorrectedPFMetsUnclusteredEnDown[0]).e_longitudinal();
1158 metunc.p4=GENPTOLOR((patType1CorrectedPFMetsUnclusteredEnDown[0]));
1159 hbbInfo->metUncInfo.push_back(metunc);
1160 }
1161
1162 edm::Handle<edm::View<reco::MET> > patType1CorrectedPFMetUnclusteredEnUpHandle;
1163 iEvent.getByLabel("patType1CorrectedPFMetUnclusteredEnUp",patType1CorrectedPFMetUnclusteredEnUpHandle);
1164 edm::View<reco::MET> patType1CorrectedPFMetsUnclusteredEnUp = *patType1CorrectedPFMetUnclusteredEnUpHandle;
1165 if(patType1CorrectedPFMetsUnclusteredEnUp.size()){
1166 metunc.sumEt =(patType1CorrectedPFMetsUnclusteredEnUp[0]).sumEt();
1167 metunc.metSig=metSignificance(& (patType1CorrectedPFMetsUnclusteredEnUp[0]));
1168 metunc.eLong=(patType1CorrectedPFMetsUnclusteredEnUp[0]).e_longitudinal();
1169 metunc.p4=GENPTOLOR((patType1CorrectedPFMetsUnclusteredEnUp[0]));
1170 hbbInfo->metUncInfo.push_back(metunc);
1171 }
1172
1173
1174 edm::Handle<edm::View<reco::MET> > patType1p2CorrectedPFMetElectronEnDownHandle;
1175 iEvent.getByLabel("patType1p2CorrectedPFMetElectronEnDown",patType1p2CorrectedPFMetElectronEnDownHandle);
1176 edm::View<reco::MET> patType1p2CorrectedPFMetsElectronEnDown = *patType1p2CorrectedPFMetElectronEnDownHandle;
1177 if(patType1p2CorrectedPFMetsElectronEnDown.size()){
1178 metunc.sumEt =(patType1p2CorrectedPFMetsElectronEnDown[0]).sumEt();
1179 metunc.metSig=metSignificance(& (patType1p2CorrectedPFMetsElectronEnDown[0]));
1180 metunc.eLong=(patType1p2CorrectedPFMetsElectronEnDown[0]).e_longitudinal();
1181 metunc.p4=GENPTOLOR((patType1p2CorrectedPFMetsElectronEnDown[0]));
1182 hbbInfo->metUncInfo.push_back(metunc);
1183 }
1184
1185 edm::Handle<edm::View<reco::MET> > patType1p2CorrectedPFMetElectronEnUpHandle;
1186 iEvent.getByLabel("patType1p2CorrectedPFMetElectronEnUp",patType1p2CorrectedPFMetElectronEnUpHandle);
1187 edm::View<reco::MET> patType1p2CorrectedPFMetsElectronEnUp = *patType1p2CorrectedPFMetElectronEnUpHandle;
1188 if(patType1p2CorrectedPFMetsElectronEnUp.size()){
1189 metunc.sumEt =(patType1p2CorrectedPFMetsElectronEnUp[0]).sumEt();
1190 metunc.metSig=metSignificance(& (patType1p2CorrectedPFMetsElectronEnUp[0]));
1191 metunc.eLong=(patType1p2CorrectedPFMetsElectronEnUp[0]).e_longitudinal();
1192 metunc.p4=GENPTOLOR((patType1p2CorrectedPFMetsElectronEnUp[0]));
1193 hbbInfo->metUncInfo.push_back(metunc);
1194 }
1195
1196
1197
1198 edm::Handle<edm::View<reco::MET> > patType1p2CorrectedPFMetMuonEnDownHandle;
1199 iEvent.getByLabel("patType1p2CorrectedPFMetMuonEnDown",patType1p2CorrectedPFMetMuonEnDownHandle);
1200 edm::View<reco::MET> patType1p2CorrectedPFMetsMuonEnDown = *patType1p2CorrectedPFMetMuonEnDownHandle;
1201 if(patType1p2CorrectedPFMetsMuonEnDown.size()){
1202 metunc.sumEt =(patType1p2CorrectedPFMetsMuonEnDown[0]).sumEt();
1203 metunc.metSig=metSignificance(& (patType1p2CorrectedPFMetsMuonEnDown[0]));
1204 metunc.eLong=(patType1p2CorrectedPFMetsMuonEnDown[0]).e_longitudinal();
1205 metunc.p4=GENPTOLOR((patType1p2CorrectedPFMetsMuonEnDown[0]));
1206 hbbInfo->metUncInfo.push_back(metunc);
1207 }
1208
1209 edm::Handle<edm::View<reco::MET> > patType1p2CorrectedPFMetMuonEnUpHandle;
1210 iEvent.getByLabel("patType1p2CorrectedPFMetMuonEnUp",patType1p2CorrectedPFMetMuonEnUpHandle);
1211 edm::View<reco::MET> patType1p2CorrectedPFMetsMuonEnUp = *patType1p2CorrectedPFMetMuonEnUpHandle;
1212 if(patType1p2CorrectedPFMetsMuonEnUp.size()){
1213 metunc.sumEt =(patType1p2CorrectedPFMetsMuonEnUp[0]).sumEt();
1214 metunc.metSig=metSignificance(& (patType1p2CorrectedPFMetsMuonEnUp[0]));
1215 metunc.eLong=(patType1p2CorrectedPFMetsMuonEnUp[0]).e_longitudinal();
1216 metunc.p4=GENPTOLOR((patType1p2CorrectedPFMetsMuonEnUp[0]));
1217 hbbInfo->metUncInfo.push_back(metunc);
1218 }
1219
1220
1221
1222 edm::Handle<edm::View<reco::MET> > patType1p2CorrectedPFMetTauEnDownHandle;
1223 iEvent.getByLabel("patType1p2CorrectedPFMetTauEnDown",patType1p2CorrectedPFMetTauEnDownHandle);
1224 edm::View<reco::MET> patType1p2CorrectedPFMetsTauEnDown = *patType1p2CorrectedPFMetTauEnDownHandle;
1225 if(patType1p2CorrectedPFMetsTauEnDown.size()){
1226 metunc.sumEt =(patType1p2CorrectedPFMetsTauEnDown[0]).sumEt();
1227 metunc.metSig=metSignificance(& (patType1p2CorrectedPFMetsTauEnDown[0]));
1228 metunc.eLong=(patType1p2CorrectedPFMetsTauEnDown[0]).e_longitudinal();
1229 metunc.p4=GENPTOLOR((patType1p2CorrectedPFMetsTauEnDown[0]));
1230 hbbInfo->metUncInfo.push_back(metunc);
1231 }
1232
1233 edm::Handle<edm::View<reco::MET> > patType1p2CorrectedPFMetTauEnUpHandle;
1234 iEvent.getByLabel("patType1p2CorrectedPFMetTauEnUp",patType1p2CorrectedPFMetTauEnUpHandle);
1235 edm::View<reco::MET> patType1p2CorrectedPFMetsTauEnUp = *patType1p2CorrectedPFMetTauEnUpHandle;
1236 if(patType1p2CorrectedPFMetsTauEnUp.size()){
1237 metunc.sumEt =(patType1p2CorrectedPFMetsTauEnUp[0]).sumEt();
1238 metunc.metSig=metSignificance(& (patType1p2CorrectedPFMetsTauEnUp[0]));
1239 metunc.eLong=(patType1p2CorrectedPFMetsTauEnUp[0]).e_longitudinal();
1240 metunc.p4=GENPTOLOR((patType1p2CorrectedPFMetsTauEnUp[0]));
1241 hbbInfo->metUncInfo.push_back(metunc);
1242 }
1243
1244
1245 edm::Handle<edm::View<reco::MET> > patType1p2CorrectedPFMetJetEnDownHandle;
1246 iEvent.getByLabel("patType1p2CorrectedPFMetJetEnDown",patType1p2CorrectedPFMetJetEnDownHandle);
1247 edm::View<reco::MET> patType1p2CorrectedPFMetsJetEnDown = *patType1p2CorrectedPFMetJetEnDownHandle;
1248 if(patType1p2CorrectedPFMetsJetEnDown.size()){
1249 metunc.sumEt =(patType1p2CorrectedPFMetsJetEnDown[0]).sumEt();
1250 metunc.metSig=metSignificance(& (patType1p2CorrectedPFMetsJetEnDown[0]));
1251 metunc.eLong=(patType1p2CorrectedPFMetsJetEnDown[0]).e_longitudinal();
1252 metunc.p4=GENPTOLOR((patType1p2CorrectedPFMetsJetEnDown[0]));
1253 hbbInfo->metUncInfo.push_back(metunc);
1254 }
1255
1256 edm::Handle<edm::View<reco::MET> > patType1p2CorrectedPFMetJetEnUpHandle;
1257 iEvent.getByLabel("patType1p2CorrectedPFMetJetEnUp",patType1p2CorrectedPFMetJetEnUpHandle);
1258 edm::View<reco::MET> patType1p2CorrectedPFMetsJetEnUp = *patType1p2CorrectedPFMetJetEnUpHandle;
1259 if(patType1p2CorrectedPFMetsJetEnUp.size()){
1260 metunc.sumEt =(patType1p2CorrectedPFMetsJetEnUp[0]).sumEt();
1261 metunc.metSig=metSignificance(& (patType1p2CorrectedPFMetsJetEnUp[0]));
1262 metunc.eLong=(patType1p2CorrectedPFMetsJetEnUp[0]).e_longitudinal();
1263 metunc.p4=GENPTOLOR((patType1p2CorrectedPFMetsJetEnUp[0]));
1264 hbbInfo->metUncInfo.push_back(metunc);
1265 }
1266
1267
1268 edm::Handle<edm::View<reco::MET> > patType1p2CorrectedPFMetJetResDownHandle;
1269 iEvent.getByLabel("patType1p2CorrectedPFMetJetResDown",patType1p2CorrectedPFMetJetResDownHandle);
1270 edm::View<reco::MET> patType1p2CorrectedPFMetsJetResDown = *patType1p2CorrectedPFMetJetResDownHandle;
1271 if(patType1p2CorrectedPFMetsJetResDown.size()){
1272 metunc.sumEt =(patType1p2CorrectedPFMetsJetResDown[0]).sumEt();
1273 metunc.metSig=metSignificance(& (patType1p2CorrectedPFMetsJetResDown[0]));
1274 metunc.eLong=(patType1p2CorrectedPFMetsJetResDown[0]).e_longitudinal();
1275 metunc.p4=GENPTOLOR((patType1p2CorrectedPFMetsJetResDown[0]));
1276 hbbInfo->metUncInfo.push_back(metunc);
1277 }
1278
1279 edm::Handle<edm::View<reco::MET> > patType1p2CorrectedPFMetJetResUpHandle;
1280 iEvent.getByLabel("patType1p2CorrectedPFMetJetResUp",patType1p2CorrectedPFMetJetResUpHandle);
1281 edm::View<reco::MET> patType1p2CorrectedPFMetsJetResUp = *patType1p2CorrectedPFMetJetResUpHandle;
1282 if(patType1p2CorrectedPFMetsJetResUp.size()){
1283 metunc.sumEt =(patType1p2CorrectedPFMetsJetResUp[0]).sumEt();
1284 metunc.metSig=metSignificance(& (patType1p2CorrectedPFMetsJetResUp[0]));
1285 metunc.eLong=(patType1p2CorrectedPFMetsJetResUp[0]).e_longitudinal();
1286 metunc.p4=GENPTOLOR((patType1p2CorrectedPFMetsJetResUp[0]));
1287 hbbInfo->metUncInfo.push_back(metunc);
1288 }
1289
1290
1291 edm::Handle<edm::View<reco::MET> > patType1p2CorrectedPFMetUnclusteredEnDownHandle;
1292 iEvent.getByLabel("patType1p2CorrectedPFMetUnclusteredEnDown",patType1p2CorrectedPFMetUnclusteredEnDownHandle);
1293 edm::View<reco::MET> patType1p2CorrectedPFMetsUnclusteredEnDown = *patType1p2CorrectedPFMetUnclusteredEnDownHandle;
1294 if(patType1p2CorrectedPFMetsUnclusteredEnDown.size()){
1295 metunc.sumEt =(patType1p2CorrectedPFMetsUnclusteredEnDown[0]).sumEt();
1296 metunc.metSig=metSignificance(& (patType1p2CorrectedPFMetsUnclusteredEnDown[0]));
1297 metunc.eLong=(patType1p2CorrectedPFMetsUnclusteredEnDown[0]).e_longitudinal();
1298 metunc.p4=GENPTOLOR((patType1p2CorrectedPFMetsUnclusteredEnDown[0]));
1299 hbbInfo->metUncInfo.push_back(metunc);
1300 }
1301
1302 edm::Handle<edm::View<reco::MET> > patType1p2CorrectedPFMetUnclusteredEnUpHandle;
1303 iEvent.getByLabel("patType1p2CorrectedPFMetUnclusteredEnUp",patType1p2CorrectedPFMetUnclusteredEnUpHandle);
1304 edm::View<reco::MET> patType1p2CorrectedPFMetsUnclusteredEnUp = *patType1p2CorrectedPFMetUnclusteredEnUpHandle;
1305 if(patType1p2CorrectedPFMetsUnclusteredEnUp.size()){
1306 metunc.sumEt =(patType1p2CorrectedPFMetsUnclusteredEnUp[0]).sumEt();
1307 metunc.metSig=metSignificance(& (patType1p2CorrectedPFMetsUnclusteredEnUp[0]));
1308 metunc.eLong=(patType1p2CorrectedPFMetsUnclusteredEnUp[0]).e_longitudinal();
1309 metunc.p4=GENPTOLOR((patType1p2CorrectedPFMetsUnclusteredEnUp[0]));
1310 hbbInfo->metUncInfo.push_back(metunc);
1311 }
1312
1313
1314
1315 //
1316 // met is calomet
1317 //
1318
1319 edm::Handle<edm::View<pat::MET> > metTCHandle;
1320 iEvent.getByLabel("patMETsTC",metTCHandle);
1321 edm::View<pat::MET> metsTC = *metTCHandle;
1322 if(metsTC.size()){
1323 hbbInfo->tcmet.sumEt=(metsTC[0]).sumEt();
1324 hbbInfo->tcmet.metSig=metSignificance(&(metsTC[0]));
1325 hbbInfo->tcmet.eLong=(metsTC[0]).e_longitudinal();
1326 hbbInfo->tcmet.p4=GENPTOLOR((metsTC[0]));
1327 if (verbose_) std::cout <<" METTC "<< hbbInfo->tcmet.metSig <<" " << hbbInfo->tcmet.sumEt<<std::endl;
1328 }
1329
1330 edm::Handle<edm::View<reco::MET> > pfMETNoPUHandle;
1331 iEvent.getByLabel("pfMETNoPU",pfMETNoPUHandle);
1332 edm::View<reco::MET> metspfMETNoPU = *pfMETNoPUHandle;
1333 if(metspfMETNoPU.size()){
1334 hbbInfo->metNoPU.sumEt=(metspfMETNoPU[0]).sumEt();
1335 hbbInfo->metNoPU.metSig=metSignificance(&(metspfMETNoPU[0]));
1336 hbbInfo->metNoPU.eLong=(metspfMETNoPU[0]).e_longitudinal();
1337 hbbInfo->metNoPU.p4=GENPTOLOR((metspfMETNoPU[0]));
1338 if (verbose_) std::cout <<" pfMETNoPU "<< hbbInfo->metNoPU.metSig <<" " << hbbInfo->metNoPU.sumEt<<std::endl;
1339 }
1340
1341 edm::Handle<edm::View<reco::MET> > mHTHandle;
1342 iEvent.getByLabel("patMETsHT",mHTHandle);
1343 edm::View<reco::MET> metsHT = *mHTHandle;
1344 if(metsHT.size()){
1345 hbbInfo->mht.sumEt=(metsHT[0]).sumEt();
1346 hbbInfo->mht.metSig=metSignificance(&(metsHT[0]));
1347 hbbInfo->mht.eLong=(metsHT[0]).e_longitudinal();
1348 hbbInfo->mht.p4=GENPTOLOR((metsHT[0]));
1349 if (verbose_) std::cout <<" METHT "<< hbbInfo->mht.metSig <<" " << hbbInfo->mht.sumEt<<std::endl;
1350 }
1351
1352 edm::Handle<edm::View<reco::MET> > metHandle;
1353 iEvent.getByLabel(metLabel_,metHandle);
1354 edm::View<reco::MET> mets = *metHandle;
1355
1356 if(mets.size()){
1357 hbbInfo->calomet.sumEt=(mets[0]).sumEt();
1358 hbbInfo->calomet.metSig=metSignificance(&(mets[0]));
1359 hbbInfo->calomet.eLong=(mets[0]).e_longitudinal();
1360 hbbInfo->calomet.p4=GENPTOLOR((mets[0]));
1361 if (verbose_) std::cout <<" METCALO "<< hbbInfo->calomet.metSig <<" " << hbbInfo->calomet.sumEt<<std::endl;
1362 }
1363
1364 edm::Handle<edm::View<pat::MET> > metPFHandle;
1365 iEvent.getByLabel("patMETsPFlow",metPFHandle);
1366 edm::View<pat::MET> metsPF = *metPFHandle;
1367
1368 if(metsPF.size()){
1369 hbbInfo->pfmet.sumEt=(metsPF[0]).sumEt();
1370 hbbInfo->pfmet.metSig=metSignificance(&(metsPF[0]));
1371 hbbInfo->pfmet.eLong=(metsPF[0]).e_longitudinal();
1372 hbbInfo->pfmet.p4=GENPTOLOR((metsPF[0]));
1373 if (verbose_) std::cout <<" METPF "<< hbbInfo->pfmet.metSig <<" " << hbbInfo->pfmet.sumEt<<std::endl;
1374 }
1375
1376
1377 if(verbose_){
1378 std::cout << "METs: calomet "<<mets.size()<<" tcmet"<<metsTC.size()<<" pfmet "<<metsPF.size()<<" MHT" <<metsHT.size()<<std::endl;
1379 }
1380
1381 if(verbose_)
1382 std::cout << " INPUT MUONS "<<muons.size()<<std::endl;
1383
1384 for(edm::View<pat::Muon>::const_iterator mu = muons.begin(); mu!=muons.end(); ++mu){
1385 VHbbEvent::MuonInfo mf;
1386 mf.p4 =GENPTOLORP( mu);
1387 mf.charge=mu->charge();
1388 mf.tIso=mu->trackIso();
1389 mf.eIso=mu->ecalIso();
1390 mf.hIso=mu->hcalIso();
1391 mf.pfChaIso=mu->chargedHadronIso();
1392 mf.pfChaPUIso=mu->puChargedHadronIso(); //userIso(5);
1393 mf.pfPhoIso=mu->photonIso();
1394 mf.pfNeuIso=mu->neutralHadronIso();
1395 Geom::Phi<double> deltaphi(mu->phi()-atan2(mf.p4.Px(), mf.p4.Py()));
1396 double acop = deltaphi.value();
1397 mf.acop=acop;
1398
1399 mf.emEnergy = mu->calEnergy().em;
1400 mf.hadEnergy = mu->calEnergy().had;
1401
1402 mf.nMatches = mu->numberOfMatches();
1403
1404 mf.ipDb=mu->dB();
1405 mf.ipErrDb=mu->edB();
1406 mf.cat=0;
1407
1408 if(mu->isGlobalMuon()) mf.cat|=1;
1409 if(mu->isTrackerMuon()) mf.cat|=2;
1410 if(mu->isStandAloneMuon()) mf.cat|=4;
1411 TrackRef trkMu1Ref = mu->get<TrackRef>();
1412 if(trkMu1Ref.isNonnull()){
1413 const Track* MuTrk1 = mu->get<TrackRef>().get();
1414 mf.zPVPt=MuTrk1->dz(RecVtxFirst.position());
1415 mf.zPVProb=MuTrk1->dz(RecVtx.position());
1416 mf.nHits=MuTrk1->numberOfValidHits();
1417 mf.chi2=MuTrk1->normalizedChi2();
1418 TrackRef iTrack1 = mu->innerTrack();
1419 const reco::HitPattern& p1 = iTrack1->hitPattern();
1420 mf.nPixelHits=p1.pixelLayersWithMeasurement();
1421
1422 mf.nValidTracker = p1.numberOfValidTrackerHits();
1423 mf.nValidPixel = p1.numberOfValidPixelHits();
1424 mf.nValidLayers = p1.trackerLayersWithMeasurement();
1425 mf.isPF = mu->isPFMuon();
1426
1427
1428
1429 }
1430 if(mu->isGlobalMuon()){
1431 TrackRef gTrack = mu->globalTrack();
1432 const reco::HitPattern& q = gTrack->hitPattern();
1433 mf.globChi2=gTrack.get()->normalizedChi2();
1434 mf.globNHits=q.numberOfValidMuonHits();
1435 mf.validMuStations = q.muonStationsWithValidHits();
1436 }else{
1437 mf.globChi2=-99;
1438 mf.globNHits=-99;
1439 }
1440
1441 //Muon trigger matching
1442 for (int itrig = 0; itrig != ntrigs; ++itrig){
1443 std::string trigName=triggerNames_.triggerName(itrig);
1444 if( (mu->triggerObjectMatchesByPath(trigName,false,false).size() != 0) ){
1445 mf.hltMatchedBits.push_back(itrig);
1446 if(verbose_){
1447 std::clog << "Trigger Matching box" << std::endl;
1448 std::clog << "+++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
1449 std::clog << "Matching parameters are defined in the cfg" << std::endl;
1450 std::clog << "Trigger bit = " << itrig << std::endl;
1451 std::clog << "Trigger name = " << trigName << std::endl;
1452 std::clog << "Trigger object matched collection size = " << mu->triggerObjectMatchesByPath(trigName,false,false).size() << std::endl;
1453 std::clog << "Pat Muon pt = " << mf.p4.Pt() << " HLT object matched = " << mu->triggerObjectMatch(0)->p4().Pt() << std::endl;
1454 std::clog << "+++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
1455 }
1456 }
1457 }
1458 //
1459
1460 // add stamuon
1461
1462 // if (mu->isStandAloneMuon()) {
1463 // reco::TrackRef sta = mu->standAloneMuon();
1464 //
1465 // }
1466
1467
1468 // int muInfo[12];
1469 // fillMuBlock(mu, muInfo);
1470 if(runOnMC_){
1471 const GenParticle* muMc = mu->genLepton();
1472 if(muMc!=0){
1473 mf.mcId=muMc->pdgId();
1474 mf.mcFourMomentum=GENPTOLORP(muMc);
1475 if(muMc->mother()!=0) mf.mcMomId=muMc->mother()->pdgId();
1476 if(muMc->mother()!=0 && muMc->mother()->mother()!=0) mf.mcgMomId=muMc->mother()->mother()->pdgId();
1477 } }
1478 hbbInfo->muInfo.push_back(mf);
1479 }
1480
1481 if(verbose_)
1482 std::cout << " INPUT electrons "<<electrons.size()<<std::endl;
1483 InputTag reducedEBRecHitCollection(string("reducedEcalRecHitsEB"));
1484 InputTag reducedEERecHitCollection(string("reducedEcalRecHitsEE"));
1485 EcalClusterLazyTools lazyTools(iEvent, iSetup, reducedEBRecHitCollection, reducedEERecHitCollection);
1486 edm::ESHandle<TransientTrackBuilder> builder;
1487 iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder", builder);
1488 const TransientTrackBuilder & transientTrackBuilder= *(builder.product());
1489
1490 for(edm::View<pat::Electron>::const_iterator elec = electrons.begin(); elec!=electrons.end(); ++elec){
1491 VHbbEvent::ElectronInfo ef;
1492 ef.p4=GENPTOLORP(elec);
1493 ef.scEta =elec->superCluster()->eta();
1494 ef.scPhi =elec->superCluster()->phi();
1495 // if(ElecEta[eleccont]!=0) ElecEt[eleccont]=elec->superCluster()->energy()/cosh(elec->superCluster()->eta());
1496 ef.charge=elec->charge();
1497 ef.tIso=elec->trackIso();
1498 ef.eIso=elec->ecalIso();
1499 ef.hIso=elec->hcalIso();
1500 ef.pfChaIso=elec->chargedHadronIso();
1501 ef.pfChaPUIso=elec->puChargedHadronIso();//userIso(5);
1502 ef.pfPhoIso=elec->photonIso();
1503 ef.pfPhoIsoDoubleCounted=0;
1504
1505 /* Check if there are photons sharing the super cluster*/
1506 for(size_t k=0;k<photonsForIso.size();k++) {
1507 if(deltaR(elec->eta(),elec->phi(),photonsForIso[k].eta(),photonsForIso[k].phi()) < 0.05 && abs(photonsForIso[k].pt()-elec->pt())/(photonsForIso[k].pt()+elec->pt()) < 0.05 ) {
1508 std::cout << "Double counting of supercluster!" << std::endl;
1509 ef.pfPhoIsoDoubleCounted+=photonsForIso[k].pt();
1510 }
1511 }
1512 ef.pfNeuIso=elec->neutralHadronIso();
1513
1514 //
1515 // ip info
1516 //
1517
1518 ef.ipDb=elec->dB();
1519 ef.ipErrDb=elec->edB();
1520
1521
1522
1523 Geom::Phi<double> deltaphi(elec->superCluster()->phi()-atan2(hbbInfo->pfmet.p4.Py(),hbbInfo->pfmet.p4.Px()));
1524 ef.acop = deltaphi.value();
1525 //
1526 ef.sihih = elec->sigmaIetaIeta();
1527 ef.Dphi = elec->deltaPhiSuperClusterTrackAtVtx();
1528 ef.Deta = elec->deltaEtaSuperClusterTrackAtVtx();
1529 ef.HoE = elec->hadronicOverEm();
1530 ef.convDist = elec->convDist();
1531 ef.convDcot = elec->convDcot();
1532 if(elec->gsfTrack().isNonnull())
1533 {
1534 ef.innerHits = elec->gsfTrack()->trackerExpectedHitsInner().numberOfHits();
1535 }
1536 ef.isEB = elec->isEB();
1537 ef.isEE = elec->isEE();
1538 /* 2012 ELEID*/
1539
1540 const pat::Electron & ele = *elec;
1541 bool validKF= false;
1542 reco::TrackRef myTrackRef = ele.closestCtfTrackRef();
1543 validKF = (myTrackRef.isAvailable());
1544 validKF = (myTrackRef.isNonnull());
1545
1546 // Pure tracking variables
1547 ef.fMVAVar_fbrem = ele.fbrem();
1548 ef.fMVAVar_kfchi2 = (validKF) ? myTrackRef->normalizedChi2() : 0 ;
1549 ef.fMVAVar_kfhits = (validKF) ? myTrackRef->hitPattern().trackerLayersWithMeasurement() : -1. ;
1550 // fMVAVar_kfhitsall = (validKF) ? myTrackRef->numberOfValidHits() : -1. ; // save also this in your ntuple as possible alternative
1551 ef.fMVAVar_gsfchi2 = ele.gsfTrack()->normalizedChi2();
1552
1553
1554 // Geometrical matchings
1555 ef.fMVAVar_deta = ele.deltaEtaSuperClusterTrackAtVtx();
1556 ef.fMVAVar_dphi = ele.deltaPhiSuperClusterTrackAtVtx();
1557 ef.fMVAVar_detacalo = ele.deltaEtaSeedClusterTrackAtCalo();
1558 // fMVAVar_dphicalo = ele.deltaPhiSeedClusterTrackAtCalo(); // save also this in your ntuple
1559
1560
1561 // Pure ECAL -> shower shapes
1562 ef.fMVAVar_see = ele.sigmaIetaIeta(); //EleSigmaIEtaIEta
1563 std::vector<float> vCov = lazyTools.localCovariances(*(ele.superCluster()->seed())) ;
1564 if (!isnan(vCov[2])) ef.fMVAVar_spp = sqrt (vCov[2]); //EleSigmaIPhiIPhi
1565 else ef.fMVAVar_spp = 0.;
1566 // fMVAVar_sigmaIEtaIPhi = vCov[1]; // save also this in your ntuple
1567
1568 ef.fMVAVar_etawidth = ele.superCluster()->etaWidth();
1569 ef.fMVAVar_phiwidth = ele.superCluster()->phiWidth();
1570 ef.fMVAVar_e1x5e5x5 = (ele.e5x5()) !=0. ? 1.-(ele.e1x5()/ele.e5x5()) : -1. ;
1571 ef.fMVAVar_R9 = lazyTools.e3x3(*(ele.superCluster()->seed())) / ele.superCluster()->rawEnergy();
1572 //fMVAVar_nbrems = fabs(ele.numberOfBrems()); // save also this in your ntuple
1573
1574 // Energy matching
1575 ef.fMVAVar_HoE = ele.hadronicOverEm();
1576 ef.fMVAVar_EoP = ele.eSuperClusterOverP();
1577 // fMVAVar_IoEmIoP = (1.0/(ele.superCluster()->energy())) - (1.0 / ele.p()); // in the future to be changed with ele.gsfTrack()->p()
1578 ef.fMVAVar_IoEmIoP = (1.0/ele.ecalEnergy()) - (1.0 / ele.p()); // in the future to be changed with ele.gsfTrack()->p() // 24/04/2012 changed to correctly access the corrected supercluster energy from CMSSW_52X
1579
1580 ef.fMVAVar_eleEoPout = ele.eEleClusterOverPout();
1581 ef.fMVAVar_PreShowerOverRaw= ele.superCluster()->preshowerEnergy() / ele.superCluster()->rawEnergy();
1582 // fMVAVar_EoPout = ele.eSeedClusterOverPout(); // save also this in your ntuple
1583
1584
1585 // Spectators
1586 ef.fMVAVar_eta = ele.superCluster()->eta();
1587 ef.fMVAVar_pt = ele.pt();
1588
1589 //additional for cut based
1590 ef.dxy = elec->gsfTrack()->dxy(vertex.position());
1591 ef.dz = elec->gsfTrack()->dz(vertex.position());
1592
1593
1594 //d0
1595 if (ele.gsfTrack().isNonnull()) {
1596 ef.fMVAVar_d0 = (-1.0)*ele.gsfTrack()->dxy(vertex.position());
1597 } else if (ele.closestCtfTrackRef().isNonnull()) {
1598 ef.fMVAVar_d0 = (-1.0)*ele.closestCtfTrackRef()->dxy(vertex.position());
1599 } else {
1600 ef.fMVAVar_d0 = -9999.0;
1601
1602 //default values for IP3D
1603 ef.fMVAVar_ip3d = -999.0;
1604 // fMVAVar_ip3dSig = 0.0;
1605 if (ele.gsfTrack().isNonnull()) {
1606 const double gsfsign = ( (-ele.gsfTrack()->dxy(vertex.position())) >=0 ) ? 1. : -1.;
1607
1608 const reco::TransientTrack &tt = transientTrackBuilder.build(ele.gsfTrack());
1609 const std::pair<bool,Measurement1D> &ip3dpv = IPTools::absoluteImpactParameter3D(tt,vertex);
1610 if (ip3dpv.first) {
1611 double ip3d = gsfsign*ip3dpv.second.value();
1612 //double ip3derr = ip3dpv.second.error();
1613 ef.fMVAVar_ip3d = ip3d;
1614 // fMVAVar_ip3dSig = ip3d/ip3derr;
1615 }
1616 }
1617 }
1618
1619
1620 /* end of 2012 ELEID*/
1621
1622 //
1623 // fill eleids
1624 //
1625 /* ef.id95 = elec->electronID("simpleEleId95cIso");
1626 ef.id85 = elec->electronID("simpleEleId85cIso");
1627 ef.id70 = elec->electronID("simpleEleId70cIso");
1628 ef.id95r = elec->electronID("simpleEleId95relIso");
1629 ef.id70r = elec->electronID("simpleEleId70relIso");
1630 ef.id85r = elec->electronID("simpleEleId85relIso");
1631 */
1632 ef.id95 =elec->electronID("eidVBTFCom95");
1633 ef.id95r=elec->electronID("eidVBTFRel95");
1634 ef.id85 =elec->electronID("eidVBTFCom85");
1635 ef.id85r=elec->electronID("eidVBTFRel85");
1636 ef.id80 =elec->electronID("eidVBTFCom80");
1637 ef.id80r=elec->electronID("eidVBTFRel80");
1638 ef.id70 =elec->electronID("eidVBTFCom70");
1639 ef.id70r=elec->electronID("eidVBTFRel70");
1640 ef.mvaOut=elec->electronID("mvaNonTrigV0");
1641 ef.mvaOutTrig=elec->electronID("mvaTrigV0");
1642
1643 //Electron trigger matching
1644 for (int itrig = 0; itrig != ntrigs; ++itrig){
1645 std::string trigName=triggerNames_.triggerName(itrig);
1646 if( (elec->triggerObjectMatchesByPath(trigName).size() != 0) ){
1647 ef.hltMatchedBits.push_back(itrig);
1648 if(verbose_){
1649 std::clog << "Trigger Matching box" << std::endl;
1650 std::clog << "+++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
1651 std::clog << "Matching parameters are defined in the cfg" << std::endl;
1652 std::clog << "Trigger bit = " << itrig << std::endl;
1653 std::clog << "Trigger name = " << trigName << std::endl;
1654 std::clog << "Trigger object matched collection size = " << elec->triggerObjectMatchesByPath(trigName).size() << std::endl;
1655 std::clog << "Pat Electron pt = " << ef.p4.Pt() << " HLT object matched = " << elec->triggerObjectMatch(0)->p4().Pt() << std::endl;
1656 std::clog << "+++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
1657 }
1658 }
1659 }
1660
1661 if(runOnMC_){
1662 const GenParticle* elecMc = elec->genLepton();
1663 if(elecMc!=0){
1664 ef.mcId=elecMc->pdgId();
1665 ef.mcFourMomentum=GENPTOLORP(elecMc);
1666 if(elecMc->mother()!=0) ef.mcMomId=elecMc->mother()->pdgId();
1667 if(elecMc->mother()!=0 && elecMc->mother()->mother()!=0) ef.mcgMomId=elecMc->mother()->mother()->pdgId();
1668 }}
1669 hbbInfo->eleInfo.push_back(ef);
1670 }
1671
1672 if(verbose_)
1673 std::cout << " INPUT taus "<<taus.size()<<std::endl;
1674 for(edm::View<pat::Tau>::const_iterator tau = taus.begin(); tau!=taus.end(); ++tau){
1675 VHbbEvent::TauInfo tf;
1676 tf.p4=GENPTOLORP(tau);
1677 tf.charge=tau->charge();
1678 tf.tIso=tau->trackIso();
1679 tf.eIso=tau->ecalIso();
1680 tf.hIso=tau->hcalIso();
1681 Geom::Phi<double> deltaphi(tau->phi()-atan2(hbbInfo->pfmet.p4.Py(),hbbInfo->pfmet.p4.Px()));
1682 double acop = deltaphi.value();
1683 tf.acop=acop;
1684 if (tau->isTauIDAvailable("againstElectronLoose")) tf.againstElectronLoose=tau->tauID("againstElectronLoose");
1685 if (tau->isTauIDAvailable("againstElectronMedium")) tf.againstElectronMedium=tau->tauID("againstElectronMedium");
1686 if (tau->isTauIDAvailable("againstElectronTight")) tf.againstElectronTight=tau->tauID("againstElectronTight");
1687 if (tau->isTauIDAvailable("againstMuonLoose")) tf.againstMuonLoose=tau->tauID("againstMuonLoose");
1688 if (tau->isTauIDAvailable("againstMuonTight")) tf.againstMuonTight=tau->tauID("againstMuonTight");
1689 if (tau->isTauIDAvailable("byLooseIsolation")) tf.byLooseIsolation=tau->tauID("byLooseIsolation");
1690 if (tau->isTauIDAvailable("byMediumIsolation")) tf.byMediumIsolation=tau->tauID("byMediumIsolation");
1691 if (tau->isTauIDAvailable("byTightIsolation")) tf.byTightIsolation=tau->tauID("byTightIsolation");
1692 if (tau->isTauIDAvailable("byVLooseIsolation")) tf.byVLooseIsolation=tau->tauID("byVLooseIsolation");
1693 if (tau->isTauIDAvailable("decayModeFinding")) tf.decayModeFinding=tau->tauID("decayModeFinding");
1694 if (tau->isTauIDAvailable("byIsolation")) tf.byIsolation=tau->tauID("byIsolation");
1695 if (tau->isTauIDAvailable("trackIsolation")) tf.trackIsolation=tau->tauID("trackIsolation");
1696 if (tau->isTauIDAvailable("byTaNCfrOnePercent")) tf.byTaNCfrOnePercent=tau->tauID("byTaNCfrOnePercent");
1697 if (tau->isTauIDAvailable("byTaNCfrHalfPercent")) tf.byTaNCfrHalfPercent=tau->tauID("byTaNCfrHalfPercent");
1698 if (tau->isTauIDAvailable("byTaNCfrQuarterPercent")) tf.byTaNCfrQuarterPercent=tau->tauID("byTaNCfrQuarterPercent");
1699 if (tau->isTauIDAvailable("byTaNCfrTenthPercent")) tf.byTaNCfrTenthPercent=tau->tauID("byTaNCfrTenthPercent");
1700 if (tau->isTauIDAvailable("byTaNC")) tf.byTaNC=tau->tauID("byTaNC");
1701 if (tau->isTauIDAvailable("byLooseCombinedIsolationDeltaBetaCorr")) tf.byLooseCombinedIsolationDeltaBetaCorr=tau->tauID("byLooseCombinedIsolationDeltaBetaCorr");
1702 if (tau->isTauIDAvailable("againstElectronMVA")) tf.againstElectronMVA=tau->tauID("againstElectronMVA");
1703 if (tau->isPFTau()) {
1704 tf.isolationPFChargedHadrCandsPtSum = tau->isolationPFChargedHadrCandsPtSum();
1705 tf.isolationPFGammaCandsEtSum = tau->isolationPFGammaCandsEtSum();
1706 if (tau->leadPFChargedHadrCand().isAvailable()) tf.leadPFChargedHadrCandPt = tau->leadPFChargedHadrCand()->pt();
1707 tf.NsignalPFChargedHadrCands = tau->signalPFChargedHadrCands().size();
1708 tf.NsignalPFGammaCands = tau->signalPFGammaCands().size();
1709 }
1710 hbbInfo->tauInfo.push_back(tf);
1711 if (verbose_) {
1712 std::cout << "SCZ DEBUG: againstElectronLoose is " << tf.againstElectronLoose << std::endl;
1713 std::cout << "SCZ DEBUG: againstElectronMedium is " << tf.againstElectronMedium << std::endl;
1714 std::cout << "SCZ DEBUG: againstElectronTight is " << tf.againstElectronTight << std::endl;
1715 std::cout << "SCZ DEBUG: againstMuonLoose is " << tf.againstMuonLoose << std::endl;
1716 std::cout << "SCZ DEBUG: againstMuonTight is " << tf.againstMuonTight << std::endl;
1717 std::cout << "SCZ DEBUG: byLooseIsolation is " << tf.byLooseIsolation << std::endl;
1718 std::cout << "SCZ DEBUG: byMediumIsolation is " << tf.byMediumIsolation << std::endl;
1719 std::cout << "SCZ DEBUG: byTightIsolation is " << tf.byTightIsolation << std::endl;
1720 std::cout << "SCZ DEBUG: byVLooseIsolation is " << tf.byVLooseIsolation << std::endl;
1721 std::cout << "SCZ DEBUG: decayModeFinding is " << tf.decayModeFinding << std::endl;
1722 std::cout << "SCZ DEBUG: byIsolation is " << tf.byIsolation<< std::endl;
1723 std::cout << "SCZ DEBUG: trackIsolation is " << tf.trackIsolation << std::endl;
1724 std::cout << "SCZ DEBUG: byTaNCfrOnePercent is " << tf.byTaNCfrOnePercent << std::endl;
1725 std::cout << "SCZ DEBUG: byTaNCfrHalfPercent is " << tf.byTaNCfrHalfPercent << std::endl;
1726 std::cout << "SCZ DEBUG: byTaNCfrQuarterPercent is " << tf.byTaNCfrQuarterPercent << std::endl;
1727 std::cout << "SCZ DEBUG: byTaNCfrTenthPercent is " << tf.byTaNCfrTenthPercent << std::endl;
1728 std::cout << "SCZ DEBUG: byTaNC is " << tf.byTaNC << std::endl;
1729 std::cout << "SCZ DEBUG: isolationPFChargedHadrCandsPtSum is " << tf.isolationPFChargedHadrCandsPtSum << std::endl;
1730 std::cout << "SCZ DEBUG: isolationPFGammaCandsEtSum is " << tf.isolationPFGammaCandsEtSum << std::endl;
1731 std::cout << "SCZ DEBUG: isolationPFGammaCandsEtSum is " << tf.leadPFChargedHadrCandPt << std::endl;
1732 std::cout << "SCZ DEBUG: NsignalPFChargedHadrCands is " << tf.NsignalPFChargedHadrCands << std::endl;
1733 std::cout << "SCZ DEBUG: NsignalPFGammaCands is " << tf.NsignalPFGammaCands << std::endl;
1734 std::cout << "SCZ DEBUG: byLooseCombinedIsolationDeltaBetaCorr is " << tf.byLooseCombinedIsolationDeltaBetaCorr << std::endl;
1735 std::cout << "SCZ DEBUG: againstElectronMVA is " << tf.againstElectronMVA << std::endl;
1736 }
1737 }
1738
1739 CompareJetPtMuons ptComparatorMu;
1740 CompareJetPtElectrons ptComparatorE;
1741 CompareJetPtTaus ptComparatorTau;
1742
1743 std::sort(hbbInfo->muInfo.begin(), hbbInfo->muInfo.end(), ptComparatorMu);
1744 std::sort(hbbInfo->eleInfo.begin(), hbbInfo->eleInfo.end(), ptComparatorE);
1745 std::sort(hbbInfo->tauInfo.begin(), hbbInfo->tauInfo.end(), ptComparatorTau);
1746
1747
1748
1749
1750 if (verbose_){
1751 std::cout <<" Pushing hbbInfo "<<std::endl;
1752 std::cout <<" SimpleJets1 = "<<hbbInfo->simpleJets.size()<<std::endl<<
1753 " SimpleJets2 = "<<hbbInfo->simpleJets2.size()<<std::endl<<
1754 " SubJets = "<<hbbInfo->subJets.size()<<std::endl<<
1755 " HardJets = "<<hbbInfo->hardJets.size()<<std::endl<<
1756 " FilterJets = "<<hbbInfo->filterJets.size()<<std::endl<<
1757 " Muons = "<<hbbInfo->muInfo.size()<<std::endl<<
1758 " Electrons = "<<hbbInfo->eleInfo.size()<<std::endl<<
1759 " Taus = "<<hbbInfo->tauInfo.size()<<std::endl<<
1760 " Electrons = "<<hbbInfo->eleInfo.size()<<std::endl<<
1761 "--------------------- "<<std::endl;
1762 }
1763
1764
1765 iEvent.put(hbbInfo);
1766 iEvent.put(auxInfo);
1767
1768
1769 delete jecUnc;
1770
1771 }
1772
1773 void
1774 HbbAnalyzerNew::fillMuBlock(edm::View<pat::Muon>::const_iterator mu, int muInfo[15])
1775 {
1776 if(muon::isGoodMuon(*mu,muon::TMLastStationLoose)) muInfo[0]=1;
1777 if(muon::isGoodMuon(*mu,muon::TMLastStationTight)) muInfo[1]=1;
1778 if(muon::isGoodMuon(*mu,muon::TM2DCompatibilityLoose)) muInfo[2]=1;
1779 if(muon::isGoodMuon(*mu,muon::TM2DCompatibilityTight)) muInfo[3]=1;
1780 if(muon::isGoodMuon(*mu,muon::TMOneStationLoose)) muInfo[4]=1;
1781 if(muon::isGoodMuon(*mu,muon::TMOneStationTight)) muInfo[5]=1;
1782 if(muon::isGoodMuon(*mu,muon::TMLastStationOptimizedLowPtLoose)) muInfo[6]=1;
1783 if(muon::isGoodMuon(*mu,muon::TMLastStationOptimizedLowPtTight))muInfo[7]=1;
1784 if(muon::isGoodMuon(*mu,muon::TMOneStationAngLoose)) muInfo[8]=1;
1785 if(muon::isGoodMuon(*mu,muon::TMOneStationAngTight)) muInfo[9]=1;
1786 if(muon::isGoodMuon(*mu,muon::TMLastStationAngLoose)) muInfo[10]=1;
1787 if(muon::isGoodMuon(*mu,muon::TMLastStationAngTight)) muInfo[11]=1;
1788 if(muon::isGoodMuon(*mu,muon::GMTkChiCompatibility)) muInfo[12]=1;
1789 if(muon::isGoodMuon(*mu,muon::GMStaChiCompatibility)) muInfo[13]=1;
1790 if(muon::isGoodMuon(*mu,muon::GMTkKinkTight)) muInfo[14]=1;
1791 }
1792
1793 // ------------ method called once each job just before starting event loop ------------
1794 void
1795 HbbAnalyzerNew::beginJob(){
1796 }
1797
1798
1799 // ------------ method called once each job just after ending the event loop ------------
1800 void
1801 HbbAnalyzerNew::endJob() {
1802 }
1803
1804 TVector2 HbbAnalyzerNew::getTvect( const pat::Jet* patJet ){
1805
1806 TVector2 t_Vect(0,0);
1807 TVector2 null(0,0);
1808 TVector2 ci(0,0);
1809 TLorentzVector pi(0,0,0,0);
1810 TLorentzVector J(0,0,0,0);
1811 TVector2 r(0,0);
1812 double patJetpfcPt = 1e10;
1813 double r_mag = 1e10;
1814 unsigned int nOfconst = 0;
1815
1816
1817 if (patJet->isPFJet() == false) {
1818 return t_Vect;
1819 }
1820
1821
1822 //re-reconstruct the jet direction with the charged tracks
1823 std::vector<reco::PFCandidatePtr>
1824 patJetpfc = patJet->getPFConstituents();
1825 for(size_t idx = 0; idx < patJetpfc.size(); idx++){
1826 if( patJetpfc.at(idx)->charge() != 0 ){
1827 pi.SetPtEtaPhiE( patJetpfc.at(idx)->pt(), patJetpfc.at(idx)->eta(), patJetpfc.at(idx)->phi(), patJetpfc.at(idx)->energy() );
1828 J += pi;
1829 nOfconst++;
1830 }
1831 }
1832 // if there are less than two charged tracks do not calculate the pull (there is not enough info). It returns a null vector
1833
1834 if( nOfconst < 2 )
1835 return null;
1836
1837
1838
1839 TVector2 v_J( J.Rapidity(), J.Phi() );
1840 //calculate TVector using only charged tracks
1841 for(size_t idx = 0; idx < patJetpfc.size(); idx++){
1842 if( patJetpfc.at(idx)->charge() != 0 ){
1843 patJetpfcPt = patJetpfc.at(idx)->pt();
1844 pi.SetPtEtaPhiE( patJetpfc.at(idx)->pt(), patJetpfc.at(idx)->eta(), patJetpfc.at(idx)->phi(), patJetpfc.at(idx)->energy() );
1845 r.Set( pi.Rapidity() - J.Rapidity(), Geom::deltaPhi( patJetpfc.at(idx)->phi(), J.Phi() ) );
1846 r_mag = r.Mod();
1847 t_Vect += ( patJetpfcPt / J.Pt() ) * r_mag * r;
1848 }
1849 }
1850
1851
1852 return t_Vect;
1853
1854 }
1855
1856 TLorentzVector HbbAnalyzerNew::getChargedTracksMomentum(const pat::Jet* patJet ){
1857 // return TLorentzVector();
1858 TLorentzVector pi(0,0,0,0);
1859 TLorentzVector v_j1(0,0,0,0);
1860
1861
1862 // std::cout <<"fff ECCCCCCOOOOO "<<patJet->isPFJet()<<std::endl;
1863
1864 if (patJet->isPFJet() == false ){
1865 v_j1 = GENPTOLORP(patJet);
1866 return v_j1;
1867 }
1868 std::vector<reco::PFCandidatePtr>
1869 j1pfc = patJet->getPFConstituents();
1870 for(size_t idx = 0; idx < j1pfc.size(); idx++){
1871 if( j1pfc.at(idx)->charge() != 0 ){
1872 pi.SetPtEtaPhiE( j1pfc.at(idx)->pt(), j1pfc.at(idx)->eta(), j1pfc.at(idx)->phi(), j1pfc.at(idx)->energy() );
1873 v_j1 += pi;
1874 }
1875 }
1876 return v_j1;
1877 //re-
1878 }
1879
1880
1881 //Btagging scale factors
1882 void HbbAnalyzerNew::fillScaleFactors(VHbbEvent::SimpleJet& sj, BTagSFContainer iSF){
1883
1884
1885 BinningPointByMap measurePoint;
1886 //for a USDG
1887 //for CB jets
1888 //scale factor 1 for CB jets over 240GeV/c
1889 if( TMath::Abs(sj.flavour) == 4 or TMath::Abs(sj.flavour) == 5 ){
1890 measurePoint.insert( BinningVariables::JetEt, sj.p4.Et() );
1891 measurePoint.insert( BinningVariables::JetAbsEta, fabs(sj.p4.Eta()) );
1892 if( iSF.BTAGSF_CSVL->isResultOk(PerformanceResult::BTAGBEFFCORR , measurePoint) ){
1893 sj.SF_CSVL = iSF.BTAGSF_CSVL->getResult(PerformanceResult::BTAGBEFFCORR , measurePoint);
1894 sj.SF_CSVLerr = iSF.BTAGSF_CSVL->getResult(PerformanceResult::BTAGBERRCORR , measurePoint);
1895 if(verbose_){
1896 std::clog << "C/B Jet flavour = " << sj.flavour << std::endl;
1897 std::clog << "C/B Jet Et = " << sj.p4.Et() << std::endl;
1898 std::clog << "C/B Jet eta = " << sj.p4.Eta() << std::endl;
1899 std::clog << "C/B CSVL Scale Factor = " << sj.SF_CSVL << std::endl;
1900 std::clog << "C/B CSVL Scale Factor error = " << sj.SF_CSVLerr << std::endl;
1901 }
1902 }
1903 if( iSF.BTAGSF_CSVM->isResultOk(PerformanceResult::BTAGBEFFCORR , measurePoint) ){
1904 sj.SF_CSVM = iSF.BTAGSF_CSVM->getResult(PerformanceResult::BTAGBEFFCORR , measurePoint);
1905 sj.SF_CSVMerr = iSF.BTAGSF_CSVM->getResult(PerformanceResult::BTAGBERRCORR , measurePoint);
1906 }
1907 if( iSF.BTAGSF_CSVT->isResultOk(PerformanceResult::BTAGBEFFCORR , measurePoint) ){
1908 sj.SF_CSVT = iSF.BTAGSF_CSVT->getResult(PerformanceResult::BTAGBEFFCORR , measurePoint);
1909 sj.SF_CSVTerr = iSF.BTAGSF_CSVT->getResult(PerformanceResult::BTAGBERRCORR , measurePoint);
1910 }
1911 else{
1912 if(verbose_){
1913 std::cerr << "No SF found in the database for this jet" << std::endl;
1914 std::clog << "No SF found: Jet flavour = " << sj.flavour << std::endl;
1915 std::clog << "No SF found: Jet Et = " << sj.p4.Et() << std::endl;
1916 std::clog << "No SF found: Jet eta = " << sj.p4.Eta() << std::endl;
1917 }
1918 }
1919 }
1920 else {
1921 measurePoint.insert( BinningVariables::JetEt, sj.p4.Et() );
1922 measurePoint.insert( BinningVariables::JetAbsEta, fabs(sj.p4.Eta()) );
1923 if( iSF.MISTAGSF_CSVL->isResultOk(PerformanceResult::BTAGLEFFCORR , measurePoint) ){
1924 sj.SF_CSVL = iSF.MISTAGSF_CSVL->getResult(PerformanceResult::BTAGLEFFCORR , measurePoint);
1925 sj.SF_CSVLerr = iSF.MISTAGSF_CSVL->getResult(PerformanceResult::BTAGLERRCORR , measurePoint);
1926 if(verbose_){
1927 std::clog << "Light Jet flavour = " << sj.flavour << std::endl;
1928 std::clog << "Light Jet Et = " << sj.p4.Et() << std::endl;
1929 std::clog << "Light Jet eta = " << sj.p4.Eta() << std::endl;
1930 std::clog << "Light CSVL Scale Factor = " << sj.SF_CSVL << std::endl;
1931 std::clog << "Light CSVL Scale Factor error = " << sj.SF_CSVLerr << std::endl;
1932 }
1933 }
1934 if( iSF.MISTAGSF_CSVM->isResultOk(PerformanceResult::BTAGLEFFCORR , measurePoint) ){
1935 sj.SF_CSVM = iSF.MISTAGSF_CSVM->getResult(PerformanceResult::BTAGLEFFCORR , measurePoint);
1936 sj.SF_CSVMerr = iSF.MISTAGSF_CSVM->getResult(PerformanceResult::BTAGLERRCORR , measurePoint);
1937 }
1938 if( iSF.MISTAGSF_CSVT->isResultOk(PerformanceResult::BTAGLEFFCORR , measurePoint) ){
1939 sj.SF_CSVT = iSF.MISTAGSF_CSVT->getResult(PerformanceResult::BTAGLEFFCORR , measurePoint);
1940 sj.SF_CSVTerr = iSF.MISTAGSF_CSVT->getResult(PerformanceResult::BTAGLERRCORR , measurePoint);
1941 }
1942 else{
1943 if(verbose_){
1944 std::cerr << "No SF found in the database for this jet" << std::endl;
1945 std::clog << "No SF found: Jet flavour = " << sj.flavour << std::endl;
1946 std::clog << "No SF found: Jet Et = " << sj.p4.Et() << std::endl;
1947 std::clog << "No SF found: Jet eta = " << sj.p4.Eta() << std::endl;
1948 }
1949 }
1950 }
1951
1952 }
1953
1954 void HbbAnalyzerNew::setJecUnc(VHbbEvent::SimpleJet& sj,JetCorrectionUncertainty* jecunc){
1955 //
1956 // test
1957 //
1958
1959 // return;
1960 double eta = sj.p4.Eta();
1961 double pt = sj.p4.Pt();
1962
1963 jecunc->setJetEta(eta);
1964 jecunc->setJetPt(pt); // here you must use the CORRECTED jet pt
1965 double unc = jecunc->getUncertainty(true);
1966 sj.jecunc= unc;
1967 }
1968
1969
1970 void HbbAnalyzerNew ::fillSimpleJet (VHbbEvent::SimpleJet& sj, edm::View<pat::Jet>::const_iterator jet_iter){
1971 sj.flavour = jet_iter->partonFlavour();
1972
1973 sj.tche=jet_iter->bDiscriminator("trackCountingHighEffBJetTags");
1974 sj.tchp=jet_iter->bDiscriminator("trackCountingHighPurBJetTags");
1975 sj.jp=jet_iter->bDiscriminator("jetProbabilityBJetTags");
1976 sj.jpb=jet_iter->bDiscriminator("jetBProbabilityBJetTags");
1977 sj.ssvhe=jet_iter->bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
1978 sj.csv=jet_iter->bDiscriminator("combinedSecondaryVertexBJetTags");
1979 sj.csvmva=jet_iter->bDiscriminator("combinedSecondaryVertexMVABJetTags");
1980 sj.csvivf=jet_iter->bDiscriminator("combinedInclusiveSecondaryVertexBJetTags");
1981 sj.cmva=jet_iter->bDiscriminator("combinedMVABJetTags");
1982 sj.charge=jet_iter->jetCharge();
1983 sj.ntracks=jet_iter->associatedTracks().size();
1984 sj.p4=GENPTOLORP(jet_iter);
1985 // std::cout << " ECCO "<<sj.csv<< " "<< sj.p4.Pt()<<std::endl;
1986 sj.chargedTracksFourMomentum=(getChargedTracksMomentum(&*(jet_iter)));
1987 sj.SF_CSVL=1;
1988 sj.SF_CSVM=1;
1989 sj.SF_CSVT=1;
1990 sj.SF_CSVLerr=0;
1991 sj.SF_CSVMerr=0;
1992 sj.SF_CSVTerr=0;
1993
1994
1995
1996
1997
1998 if (jet_iter->isPFJet() == true) {
1999
2000 sj.chargedHadronEFraction = jet_iter-> chargedHadronEnergyFraction();
2001 sj.neutralHadronEFraction = jet_iter-> neutralHadronEnergyFraction ();
2002 sj.chargedEmEFraction = jet_iter-> chargedEmEnergyFraction ();
2003 sj.neutralEmEFraction = jet_iter-> neutralEmEnergyFraction ();
2004 sj.nConstituents = jet_iter->getPFConstituents().size();
2005
2006 }
2007 sj.jetArea = jet_iter->jetArea();
2008 //
2009 // addtaginfo for csv
2010 //
2011
2012 // if (jet_iter->hasTagInfo("SimpleSecondaryVertex")) {
2013
2014 const reco::SecondaryVertexTagInfo * tf = jet_iter->tagInfoSecondaryVertex();
2015 if (tf){
2016 math::XYZTLorentzVectorD vertexSum;
2017 for(size_t vi=0;vi< tf->nVertices();vi++)
2018 {
2019 vertexSum+=tf->secondaryVertex(vi).p4();
2020 }
2021 sj.vtxP4 = GENPTOLOR(vertexSum);
2022
2023 if (tf->nVertices() >0){
2024 sj.vtxPosition = TVector3(tf->secondaryVertex(0).position().x(),tf->secondaryVertex(0).position().y(),tf->secondaryVertex(0).position().z());
2025 sj.vtxMass = tf->secondaryVertex(0).p4().mass();
2026 sj.vtxNTracks = tf->secondaryVertex(0).nTracks();
2027 std::vector<reco::TrackBaseRef >::const_iterator tit = tf->secondaryVertex(0).tracks_begin();
2028 for (; tit< tf->secondaryVertex(0).tracks_end(); ++tit){
2029 sj.vtxTrackIds.push_back(tit->key());
2030 }
2031 Measurement1D m = tf->flightDistance(0);
2032 sj.vtx3dL = m.value();
2033 sj.vtx3deL = m.error();
2034 }
2035 }
2036
2037 // CSV track info
2038 const reco::SecondaryVertexTagInfo * svTagInfos = jet_iter->tagInfoSecondaryVertex();
2039 const reco::TrackIPTagInfo * ipTagInfos = jet_iter->tagInfoTrackIP();
2040 for (edm::RefVector<reco::TrackCollection>::const_iterator t = ipTagInfos->selectedTracks().begin(); t != ipTagInfos->selectedTracks().end(); t++){
2041 sj.btagTrackIds.push_back(t->key());
2042 }// all btag IP selected tracks
2043 std::vector<const reco::BaseTagInfo*> tagInfos;
2044 tagInfos.push_back(dynamic_cast<const reco::BaseTagInfo*>(ipTagInfos));
2045 tagInfos.push_back(dynamic_cast<const reco::BaseTagInfo*>(svTagInfos));
2046 JetTagComputer::TagInfoHelper helper(tagInfos);
2047 reco::TaggingVariableList varList = computer->taggingVariables(helper); // computer for getting CSV variables
2048
2049 for(reco::TaggingVariableList::const_iterator iter = varList.begin(); iter != varList.end(); ++iter)
2050 {
2051 //std::cout << reco::TaggingVariableTokens[iter->first] << " = " << iter->second << std::endl;
2052 for (edm::RefVector<reco::TrackCollection>::const_iterator t = ipTagInfos->selectedTracks().begin(); t != ipTagInfos->selectedTracks().end(); t++){
2053
2054 if (strcmp(reco::TaggingVariableTokens[iter->first], "trackMomentum") == 0 && (fabs((float)iter->second - (float)(*t)->p()) < 0.0001) ){
2055 sj.csvTrackIds.push_back(t->key());
2056 }// if tagged track
2057 }// loop on IPtracks
2058 }// loop on CSV variables
2059
2060
2061 sj.btagNTracks= ipTagInfos->selectedTracks().size();
2062 sj.csvNTracks = sj.csvTrackIds.size();
2063
2064 //
2065 // add tVector
2066 //
2067 sj.tVector = getTvect(&(*jet_iter));
2068
2069 sj.ptRaw = jet_iter->correctedJet(0).pt();
2070
2071 sj.ptLeadTrack =-9999.;
2072 if (jet_iter->isPFJet() == true) {
2073 std::vector <reco::PFCandidatePtr> constituents = jet_iter->getPFConstituents ();
2074 for (unsigned ic = 0; ic < constituents.size (); ++ic) {
2075 if ( constituents[ic]->particleId() > 3 ) continue;
2076 reco::TrackRef trackRef = constituents[ic]->trackRef();
2077 if ( trackRef.isNonnull() ) { if(trackRef->pt() > sj.ptLeadTrack) sj.ptLeadTrack=trackRef->pt(); }
2078 }
2079 }
2080
2081
2082 }
2083
2084 float HbbAnalyzerNew::metSignificance(const reco::MET * met)
2085 {
2086 double sigmaX2= met->getSignificanceMatrix()(0,0);
2087 double sigmaY2= met->getSignificanceMatrix()(1,1);
2088 double significance = 0;
2089 try {
2090 if(sigmaX2<1.e10 && sigmaY2<1.e10) significance = met->significance();
2091 }
2092 catch(...)
2093 {
2094 std::cout << "PROBLEM WITH MET SIGNIFICANCE sigma X2 and Y2 are: " << sigmaX2 << " " << sigmaY2 << std::endl;
2095 }
2096 return significance;
2097 }
2098
2099
2100 //define this as a plug-in
2101 DEFINE_FWK_MODULE(HbbAnalyzerNew);