ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbbAnalysis/HbbAnalyzer/plugins/HbbAnalyzerNew.cc
Revision: 1.77
Committed: Mon Aug 27 21:34:22 2012 UTC (12 years, 8 months ago) by dlopes
Content type: text/plain
Branch: MAIN
CVS Tags: EDMV42_Step2_V6, EDMV42_Step2_V5a, EDMV42_Step2_V5, tauCandV42, hbbsubstructDev, EDMV42_Step2_V4a, EDMV42_Step2_V4, EDMV42_Step2_V3, EDMV42_Step2_V2, EDMV42_Step2_V1, EdmV42
Branch point for: V42TauCandidate, hbbsubstructDevPostHCP
Changes since 1.76: +28 -1 lines
Log Message:
add back gen jet info on filter jet

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