ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbbAnalysis/HbbAnalyzer/plugins/HbbAnalyzerNew.cc
Revision: 1.80
Committed: Fri Jan 25 22:16:36 2013 UTC (12 years, 3 months ago) by apana
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.79: +79 -10 lines
Log Message: 
Add softlepton info to filterjets; move fetching of collections used for soft lepton finding outside of SimpleJets2 loop

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