ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/HbbAnalysis/plugins/HbbTreeMaker.cc
Revision: 1.3
Committed: Tue Oct 13 12:14:05 2009 UTC (15 years, 7 months ago) by amagnan
Content type: text/plain
Branch: MAIN
Changes since 1.2: +28 -5 lines
Log Message:
add innerTrack IPd0, IPdz and nHits info for muons.

File Contents

# Content
1 #include "DataFormats/Common/interface/Handle.h"
2 #include "DataFormats/Common/interface/TriggerResults.h"
3 #include "DataFormats/Common/interface/HLTenums.h"
4 #include "DataFormats/Common/interface/ValueMap.h"
5 #include "DataFormats/Candidate/interface/Candidate.h"
6 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
7 #include "DataFormats/MuonReco/interface/Muon.h"
8 #include "DataFormats/MuonReco/interface/MuonFwd.h"
9 #include "DataFormats/VertexReco/interface/Vertex.h"
10 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
11 #include "DataFormats/TrackReco/interface/TrackFwd.h"
12
13 #include "DataFormats/PatCandidates/interface/Lepton.h"
14 #include "DataFormats/PatCandidates/interface/Muon.h"
15 #include "DataFormats/PatCandidates/interface/Electron.h"
16 #include "DataFormats/PatCandidates/interface/Tau.h"
17 #include "DataFormats/PatCandidates/interface/Jet.h"
18 #include "DataFormats/PatCandidates/interface/MET.h"
19
20 #include "FWCore/ServiceRegistry/interface/Service.h"
21 #include "PhysicsTools/UtilAlgos/interface/TFileService.h"
22
23
24 #include "UserCode/HbbAnalysis/interface/Electron.hh"
25 #include "UserCode/HbbAnalysis/interface/Muon.hh"
26 #include "UserCode/HbbAnalysis/interface/Tau.hh"
27 #include "UserCode/HbbAnalysis/interface/Jet.hh"
28 #include "UserCode/HbbAnalysis/interface/Met.hh"
29 #include "UserCode/HbbAnalysis/interface/Trigger.hh"
30
31 #include "UserCode/HbbAnalysis/plugins/HbbTreeMaker.hh"
32
33 using namespace HbbAnalysis;
34
35 HbbTreeMaker::HbbTreeMaker(const edm::ParameterSet & pset):
36 debug_(pset.getParameter<int>("DEBUG")),
37 flavour_(pset.getParameter<unsigned int>("JetFlavour")),
38 doGen_(pset.getParameter<bool>("DOGEN")),
39 genParticleSrc_(pset.getParameter<edm::InputTag>("GenParticles")),
40 electronSrc_(pset.getParameter<edm::InputTag>("Electrons")),
41 muonSrc_(pset.getParameter<edm::InputTag>("Muons")),
42 caloTauSrc_(pset.getParameter<edm::InputTag>("CaloTaus")),
43 pfTauSrc_(pset.getParameter<edm::InputTag>("PFTaus")),
44 caloJetSrc_(pset.getParameter<edm::InputTag>("CaloJets")),
45 jptJetSrc_(pset.getParameter<edm::InputTag>("JPTJets")),
46 pfJetSrc_(pset.getParameter<edm::InputTag>("PFJets")),
47 caloMetSrc_(pset.getParameter<edm::InputTag>("CaloMET")),
48 tcMetSrc_(pset.getParameter<edm::InputTag>("TCMET")),
49 pfMetSrc_(pset.getParameter<edm::InputTag>("PFMET")),
50 //pairSrc_(pset.getParameter<edm::InputTag>("Pair")),
51 //mmPairSrc_(pset.getParameter<edm::InputTag>("MuMuPair")),
52 //etPairSrc_(pset.getParameter<edm::InputTag>("ETauPair")),
53 //mtPairSrc_(pset.getParameter<edm::InputTag>("MuTauPair")),
54 vertexSrc_(pset.getParameter<edm::InputTag>("Vertex")),
55 triggerSrc_(pset.getParameter<edm::InputTag>("Trigger")),
56 hltPaths_(pset.getParameter<std::vector<std::string> >("HLTPaths")),
57 event_(0),
58 tree_(0)
59 //processVec_(pset.getParameter<std::vector<std::string> >("ProcessVec"))
60 {//constructor
61
62
63 }//constructor
64
65 HbbTreeMaker::~HbbTreeMaker(){//destructor
66 }//destructor
67
68
69
70 void HbbTreeMaker::beginJob(const edm::EventSetup&){//beginJob
71
72
73 edm::Service<TFileService> lFileService;
74
75 TFileDirectory lDir = lFileService->mkdir("HbbAnalysis");
76
77 tree_ = lDir.make<TTree>("Tree","Tree");
78 tree_->Branch("HbbEvent","HbbAnalysis::HbbEvent",&event_);
79 event_ = new HbbEvent();
80
81 if (debug_) std::cout << "Initialising JetFlavour : " << std::endl;
82
83 lDir = lFileService->mkdir("JetFlavours");
84
85 jetFlav_.Initialise(lDir, debug_, flavour_);
86
87 }//beginJob
88
89 void HbbTreeMaker::endJob(){//endJob
90 jetFlav_.printSummary();
91
92 //tree_->Write();
93 //delete tree_;
94 //delete event_;
95
96 }//endJob
97
98 void HbbTreeMaker::analyze(const edm::Event& aEvt, const edm::EventSetup& aEvtSetup){//analyze
99
100 //dumpContent(aEvt,"","MuTauPAT");
101 if (debug_) std::cout << "Processing event #" << aEvt.id().event() << std::endl;
102
103
104 static bool firstEvent = true;
105
106 event_->Clear();
107 event_->event(aEvt.id().event());
108
109 edm::Handle<reco::GenParticleCollection> lGenParticles;
110 try {
111 aEvt.getByLabel(genParticleSrc_,lGenParticles);
112 if (debug_) std::cout << "** ngenParticles = " << lGenParticles->size() << std::endl;
113 } catch(cms::Exception& e) {
114 std::cout << "AMM: Collection genParticles not available! Exception : " << e.what() << ". " << std::endl;
115 }
116
117 if (doGen_) HbbParticles(lGenParticles,event_->particles());
118
119 unsigned int lNPartons = jetFlav_.fillPartons(lGenParticles);
120
121 if (debug_) std::cout << "--- Number of partons = " << lNPartons << "." << std::endl;
122
123
124 edm::Handle<std::vector<reco::Vertex> > lRecoVertices;
125 try {
126 aEvt.getByLabel(vertexSrc_, lRecoVertices);
127 if (!lRecoVertices.isValid()){
128 edm::LogInfo("ERROR")<< "Error! Can't get vertex by label. ";
129 }
130 if (debug_) std::cout << "** vertexcollection = " << lRecoVertices->size() << " elements." << std::endl;
131 } catch(cms::Exception& e) {
132 std::cout << "AMM: Collection " << vertexSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
133 }
134
135 edm::Handle<std::vector<pat::Electron> > lElectronCollection;
136
137 try {
138 aEvt.getByLabel(electronSrc_,lElectronCollection);
139 if (!lElectronCollection.isValid()){
140 edm::LogInfo("ERROR")<< "Error! Can't get electron by label. ";
141 }
142 if (debug_) std::cout << "** electroncollection = " << lElectronCollection->size() << " elements." << std::endl;
143 } catch(cms::Exception& e) {
144 std::cout << "AMM: Collection " << electronSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
145 }
146
147 HbbElectrons(lElectronCollection,event_->electrons());
148
149
150 edm::Handle<std::vector<pat::Muon> > lMuonCollection;
151
152 try {
153 aEvt.getByLabel(muonSrc_,lMuonCollection);
154 if (!lMuonCollection.isValid()){
155 edm::LogInfo("ERROR")<< "Error! Can't get muon by label. ";
156 }
157 if (debug_) std::cout << "** muoncollection = " << lMuonCollection->size() << " elements." << std::endl;
158 } catch(cms::Exception& e) {
159 std::cout << "AMM: Collection " << muonSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
160 }
161
162 HbbMuons(lMuonCollection,lRecoVertices,event_->muons());
163
164
165 edm::Handle<std::vector<pat::Tau> > lTauCollection;
166
167 try {
168 aEvt.getByLabel(caloTauSrc_,lTauCollection);
169 if (!lTauCollection.isValid()){
170 edm::LogInfo("ERROR")<< "Error! Can't get caloTau by label. ";
171 }
172 if (debug_) std::cout << "** caloTaucollection = " << lTauCollection->size() << " elements." << std::endl;
173 } catch(cms::Exception& e) {
174 std::cout << "AMM: Collection " << caloTauSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
175 }
176
177 HbbTaus(lTauCollection,lRecoVertices,event_->caloTaus());
178
179 edm::Handle<std::vector<pat::Tau> > lPFTauCollection;
180
181 try {
182 aEvt.getByLabel(pfTauSrc_,lPFTauCollection);
183 if (!lPFTauCollection.isValid()){
184 edm::LogInfo("ERROR")<< "Error! Can't get pfTau by label. ";
185 }
186 if (debug_) std::cout << "** pfTaucollection = " << lPFTauCollection->size() << " elements." << std::endl;
187 } catch(cms::Exception& e) {
188 std::cout << "AMM: Collection " << pfTauSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
189 }
190
191 HbbTaus(lPFTauCollection,lRecoVertices,event_->pfTaus());
192
193 edm::Handle<std::vector<pat::Jet> > lCaloJetCollection;
194
195 try {
196 aEvt.getByLabel(caloJetSrc_,lCaloJetCollection);
197 if (!lCaloJetCollection.isValid()){
198 edm::LogInfo("ERROR")<< "Error! Can't get caloJets by label. ";
199 }
200 if (debug_) std::cout << "** caloJetcollection = " << lCaloJetCollection->size() << " elements." << std::endl;
201 } catch(cms::Exception& e) {
202 std::cout << "AMM: Collection " << caloJetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
203 }
204
205 HbbJets(lCaloJetCollection,jetFlav_,lGenParticles,event_->caloJets());
206
207 edm::Handle<std::vector<pat::Jet> > lJptJetCollection;
208
209 try {
210 aEvt.getByLabel(jptJetSrc_,lJptJetCollection);
211 if (!lJptJetCollection.isValid()){
212 edm::LogInfo("ERROR")<< "Error! Can't get jptJets by label. ";
213 }
214 if (debug_) std::cout << "** jptJetcollection = " << lJptJetCollection->size() << " elements." << std::endl;
215 } catch(cms::Exception& e) {
216 std::cout << "AMM: Collection " << jptJetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
217 }
218
219 HbbJets(lJptJetCollection,jetFlav_,lGenParticles,event_->jptJets());
220
221 edm::Handle<std::vector<pat::Jet> > lPfJetCollection;
222
223 try {
224 aEvt.getByLabel(pfJetSrc_,lPfJetCollection);
225 if (!lPfJetCollection.isValid()){
226 edm::LogInfo("ERROR")<< "Error! Can't get pfJets by label. ";
227 }
228 if (debug_) std::cout << "** pfJetcollection = " << lPfJetCollection->size() << " elements." << std::endl;
229 } catch(cms::Exception& e) {
230 std::cout << "AMM: Collection " << pfJetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
231 }
232
233 HbbJets(lPfJetCollection,jetFlav_,lGenParticles,event_->pfJets());
234
235 edm::Handle<std::vector<pat::MET> > lCaloMetCol;
236
237 try {
238 aEvt.getByLabel(caloMetSrc_,lCaloMetCol);
239 if (!lCaloMetCol.isValid()){
240 edm::LogInfo("ERROR")<< "Error! Can't get caloMet by label. ";
241 }
242 if (debug_) std::cout << "** caloMet = " << lCaloMetCol->size() << " elements." << std::endl;
243 } catch(cms::Exception& e) {
244 std::cout << "AMM: Collection " << caloMetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
245 }
246
247 HbbMet(lCaloMetCol,event_->caloMet());
248
249 edm::Handle<std::vector<pat::MET> > lTcMetCol;
250
251 try {
252 aEvt.getByLabel(tcMetSrc_,lTcMetCol);
253 if (!lTcMetCol.isValid()){
254 edm::LogInfo("ERROR")<< "Error! Can't get tcMet by label. ";
255 }
256 if (debug_) std::cout << "** tcMet = " << lTcMetCol->size() << " elements." << std::endl;
257 } catch(cms::Exception& e) {
258 std::cout << "AMM: Collection " << tcMetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
259 }
260
261 HbbMet(lTcMetCol,event_->tcMet());
262
263 edm::Handle<std::vector<pat::MET> > lPfMetCol;
264
265 try {
266 aEvt.getByLabel(pfMetSrc_,lPfMetCol);
267 if (!lPfMetCol.isValid()){
268 edm::LogInfo("ERROR")<< "Error! Can't get pfMet by label. ";
269 }
270 if (debug_) std::cout << "** pfMet = " << lPfMetCol->size() << " elements." << std::endl;
271 } catch(cms::Exception& e) {
272 std::cout << "AMM: Collection " << pfMetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
273 }
274
275 HbbMet(lPfMetCol,event_->pfMet());
276
277
278
279 //triggers
280 edm::Handle<edm::TriggerResults> lTrigCol;
281 try {
282 aEvt.getByLabel(triggerSrc_,lTrigCol);
283 if (!lTrigCol.isValid()){
284 edm::LogInfo("ERROR")<< "Error! Can't get triggers by label. ";
285 }
286 if (debug_) std::cout << "** triggers = " << lTrigCol->size() << std::endl;
287 } catch(cms::Exception& e) {
288 std::cout << "AMM: Collection " << triggerSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
289 }
290
291 HbbTrigger(lTrigCol,event_->triggers());
292
293 //event_->order();
294 tree_->Fill();
295 firstEvent = false;
296
297 }//analyze
298
299
300 void HbbTreeMaker::HbbElectrons(const edm::Handle<std::vector<pat::Electron> > & aCol,
301 std::vector<HbbAnalysis::Electron> & aVec)
302 {//HbbElectrons
303
304 if (aCol.isValid()){
305 if (aCol->size() > 0) {
306
307 unsigned int iEle = 0;
308 for (std::vector<pat::Electron>::const_iterator iter = aCol->begin();
309 iter != aCol->end();
310 iter++)
311 {
312 //if (debug_ > 1) std::cout << "**** Ele #" << iEle << ", pT,eta,phi = " << (*iter).pt() << " " << (*iter).eta() << " " << (*iter).phi() << std::endl;
313
314 HbbAnalysis::GenVars lGen;
315 if ((*iter).genLepton()){
316 lGen.valid = true;
317 lGen.E = (*iter).genLepton()->energy();
318 lGen.pT = (*iter).genLepton()->pt();
319 lGen.eta = (*iter).genLepton()->eta();
320 lGen.phi = (*iter).genLepton()->phi();
321 lGen.charge = (*iter).genLepton()->charge();
322 lGen.pdgId = (*iter).genLepton()->pdgId();
323 lGen.status = (*iter).genLepton()->status();
324 lGen.mass = (*iter).genLepton()->mass();
325 lGen.vx = (*iter).genLepton()->vx();
326 lGen.vy = (*iter).genLepton()->vy();
327 lGen.vz = (*iter).genLepton()->vz();
328 }
329 else {
330 lGen.valid = false;
331 lGen.E = 0;
332 lGen.pT = 0;
333 lGen.eta = 0;
334 lGen.phi = 0;
335 lGen.charge = 0;
336 lGen.pdgId = 0;
337 lGen.status = 0;
338 lGen.mass = 0;
339 lGen.vx = 0;
340 lGen.vy = 0;
341 lGen.vz = 0;
342 }
343
344 HbbAnalysis::BaseVars lReco;
345 lReco.E = (*iter).energy();
346 lReco.pT = (*iter).pt();
347 lReco.eta = (*iter).eta();
348 lReco.phi = (*iter).phi();
349 lReco.charge = (*iter).charge();
350 lReco.vx = (*iter).vx();
351 lReco.vy = (*iter).vy();
352 lReco.vz = (*iter).vz();
353
354 HbbAnalysis::SCVars lSC;
355 lSC.sigmaEtaEta = (*iter).scSigmaEtaEta();
356 lSC.sigmaIEtaIEta = (*iter).scSigmaIEtaIEta();
357 lSC.e1x5 = (*iter).scE1x5();
358 lSC.e2x5Max = (*iter).scE2x5Max();
359 lSC.e5x5 = (*iter).scE5x5();
360 lSC.eOverP = (*iter).eSuperClusterOverP();
361
362 HbbAnalysis::EleIsoVars lIso;
363 lIso.calo = (*iter).caloIso();
364 lIso.track = (*iter).trackIso();
365 lIso.ecal = (*iter).ecalIso();
366 lIso.hcal = (*iter).hcalIso();
367
368 HbbAnalysis::EleIDVars lID;
369 lID.electronIDs = (*iter).electronIDs();
370 lID.hOverE = (*iter).hadronicOverEm();
371 lID.deltaPhiIn = (*iter).deltaPhiSuperClusterTrackAtVtx();
372 lID.deltaEtaIn = (*iter).deltaEtaSuperClusterTrackAtVtx();
373
374 HbbAnalysis::Electron lObj(lGen,lReco,lSC,lIso,lID);
375 aVec.push_back(lObj);
376 iEle++;
377 }
378 }
379 }
380
381 }//HbbElectrons
382
383
384 void HbbTreeMaker::HbbMuons(const edm::Handle<std::vector<pat::Muon> > & aCol,
385 const edm::Handle<std::vector<reco::Vertex> > & aRecoVertices,
386 std::vector<HbbAnalysis::Muon> & aVec)
387 {//HbbMuons
388
389 if (aCol.isValid()){
390 if (aCol->size() > 0) {
391
392 unsigned int iEle = 0;
393 for (std::vector<pat::Muon>::const_iterator iter = aCol->begin();
394 iter != aCol->end();
395 iter++)
396 {
397 const reco::Muon* recoMuon = dynamic_cast<const reco::Muon*>((*iter).originalObject());
398
399 HbbAnalysis::GenVars lGen;
400 if ((*iter).genLepton()){
401 lGen.valid = true;
402 lGen.E = (*iter).genLepton()->energy();
403 lGen.pT = (*iter).genLepton()->pt();
404 lGen.eta = (*iter).genLepton()->eta();
405 lGen.phi = (*iter).genLepton()->phi();
406 lGen.charge = (*iter).genLepton()->charge();
407 lGen.pdgId = (*iter).genLepton()->pdgId();
408 lGen.status = (*iter).genLepton()->status();
409 lGen.mass = (*iter).genLepton()->mass();
410 lGen.vx = (*iter).genLepton()->vx();
411 lGen.vy = (*iter).genLepton()->vy();
412 lGen.vz = (*iter).genLepton()->vz();
413 }
414 else {
415 lGen.valid = false;
416 lGen.E = 0;
417 lGen.pT = 0;
418 lGen.eta = 0;
419 lGen.phi = 0;
420 lGen.charge = 0;
421 lGen.pdgId = 0;
422 lGen.status = 0;
423 lGen.mass = 0;
424 lGen.vx = 0;
425 lGen.vy = 0;
426 lGen.vz = 0;
427 }
428
429 HbbAnalysis::BaseVars lReco;
430 lReco.E = (*iter).energy();
431 lReco.pT = (*iter).pt();
432 lReco.eta = (*iter).eta();
433 lReco.phi = (*iter).phi();
434 lReco.charge = (*iter).charge();
435 lReco.vx = (*iter).vx();
436 lReco.vy = (*iter).vy();
437 lReco.vz = (*iter).vz();
438
439 HbbAnalysis::MuTrkVars lTrk;
440
441 reco::TrackRef lTrackerTrk = (*iter).innerTrack();
442 if ( lTrackerTrk.isAvailable() && lTrackerTrk.isNonnull() ) {
443 if ( aRecoVertices->size() >= 1 ) {
444 const reco::Vertex& thePrimaryEventVertex = (*aRecoVertices->begin());
445 lTrk.IPd0 = -lTrackerTrk->dxy(thePrimaryEventVertex.position());
446 lTrk.IPdz = lTrackerTrk->dz(thePrimaryEventVertex.position());
447 }
448 else {
449 lTrk.IPd0 = 0;
450 lTrk.IPdz = 0;
451 }
452
453 lTrk.nHits = lTrackerTrk->numberOfValidHits();
454 }
455 else {
456 lTrk.IPd0 = 0;
457 lTrk.IPdz = 0;
458 lTrk.nHits = 0;
459 }
460
461 HbbAnalysis::MuIsoVars lIsoR03;
462 lIsoR03.sumPt = recoMuon->isolationR03().sumPt;
463 lIsoR03.emEt = recoMuon->isolationR03().emEt;
464 lIsoR03.hadEt = recoMuon->isolationR03().hadEt;
465 lIsoR03.nTracks = recoMuon->isolationR03().nTracks;
466 lIsoR03.nJets = recoMuon->isolationR03().nJets;
467
468 HbbAnalysis::MuIsoVars lIsoR05;
469 lIsoR05.sumPt = recoMuon->isolationR05().sumPt;
470 lIsoR05.emEt = recoMuon->isolationR05().emEt;
471 lIsoR05.hadEt = recoMuon->isolationR05().hadEt;
472 lIsoR05.nTracks = recoMuon->isolationR05().nTracks;
473 lIsoR05.nJets = recoMuon->isolationR05().nJets;
474
475 HbbAnalysis::MuIDVars lID;
476 if (recoMuon->isGlobalMuon()) lID.type = 1;
477 else if (recoMuon->isTrackerMuon()) lID.type = 2;
478 else if (recoMuon->isStandAloneMuon()) lID.type = 3;
479 else if (recoMuon->isCaloMuon()) lID.type = 4;
480 else lID.type = 0;
481
482 if (recoMuon->isGood(reco::Muon::AllGlobalMuons)) lID.ids.push_back(1);
483 if (recoMuon->isGood(reco::Muon::AllStandAloneMuons)) lID.ids.push_back(2);
484 if (recoMuon->isGood(reco::Muon::AllTrackerMuons)) lID.ids.push_back(3);
485 if (recoMuon->isGood(reco::Muon::TrackerMuonArbitrated)) lID.ids.push_back(4);
486 if (recoMuon->isGood(reco::Muon::AllArbitrated)) lID.ids.push_back(5);
487 if (recoMuon->isGood(reco::Muon::GlobalMuonPromptTight)) lID.ids.push_back(6);
488 if (recoMuon->isGood(reco::Muon::TMLastStationLoose)) lID.ids.push_back(7);
489 if (recoMuon->isGood(reco::Muon::TMLastStationTight)) lID.ids.push_back(8);
490 if (recoMuon->isGood(reco::Muon::TM2DCompatibilityLoose)) lID.ids.push_back(9);
491 if (recoMuon->isGood(reco::Muon::TM2DCompatibilityTight)) lID.ids.push_back(10);
492 if (recoMuon->isGood(reco::Muon::TMOneStationLoose)) lID.ids.push_back(11);
493 if (recoMuon->isGood(reco::Muon::TMOneStationTight)) lID.ids.push_back(12);
494 if (recoMuon->isGood(reco::Muon::TMLastStationOptimizedLowPtLoose)) lID.ids.push_back(13);
495 if (recoMuon->isGood(reco::Muon::TMLastStationOptimizedLowPtTight)) lID.ids.push_back(14);
496
497 lID.caloCompat = recoMuon->caloCompatibility();
498 lID.segCompat = recoMuon->segmentCompatibility();
499 lID.nChambers = recoMuon->numberOfChambers();
500 lID.nMatchesLoose = recoMuon->numberOfMatches(reco::Muon::NoArbitration);
501 lID.nMatchesMedium = recoMuon->numberOfMatches(reco::Muon::SegmentArbitration);
502 lID.nMatchesTight = recoMuon->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
503
504 HbbAnalysis::Muon lObj(lGen,lReco,lTrk,lIsoR03,lIsoR05,lID);
505 aVec.push_back(lObj);
506 iEle++;
507 }
508 }
509 }
510
511 }//HbbMuons
512
513 void HbbTreeMaker::HbbTaus(const edm::Handle<std::vector<pat::Tau> > & aCol,
514 const edm::Handle<std::vector<reco::Vertex> > & aRecoVertices,
515 std::vector<HbbAnalysis::Tau> & aVec)
516 {//HbbTaus
517
518 if (aCol.isValid()){
519 if (aCol->size() > 0) {
520
521 unsigned int iEle = 0;
522 for (std::vector<pat::Tau>::const_iterator iter = aCol->begin();
523 iter != aCol->end();
524 iter++)
525 {
526 //if (debug_ > 1) std::cout << "**** Ele #" << iEle << ", pT,eta,phi = " << (*iter).pt() << " " << (*iter).eta() << " " << (*iter).phi() << std::endl;
527
528 HbbAnalysis::GenVars lGen;
529 if ((*iter).genLepton()){
530 lGen.valid = true;
531 lGen.E = (*iter).genLepton()->energy();
532 lGen.pT = (*iter).genLepton()->pt();
533 lGen.eta = (*iter).genLepton()->eta();
534 lGen.phi = (*iter).genLepton()->phi();
535 lGen.charge = (*iter).genLepton()->charge();
536 lGen.pdgId = (*iter).genLepton()->pdgId();
537 lGen.status = (*iter).genLepton()->status();
538 lGen.mass = (*iter).genLepton()->mass();
539 lGen.vx = (*iter).genLepton()->vx();
540 lGen.vy = (*iter).genLepton()->vy();
541 lGen.vz = (*iter).genLepton()->vz();
542 }
543 else {
544 lGen.valid = false;
545 lGen.E = 0;
546 lGen.pT = 0;
547 lGen.eta = 0;
548 lGen.phi = 0;
549 lGen.charge = 0;
550 lGen.pdgId = 0;
551 lGen.status = 0;
552 lGen.mass = 0;
553 lGen.vx = 0;
554 lGen.vy = 0;
555 lGen.vz = 0;
556 }
557
558
559 HbbAnalysis::GenVars lGenJet;
560 if ((*iter).genJet()){
561 lGenJet.valid = true;
562 lGenJet.pT = (*iter).genJet()->pt();
563 lGenJet.eta = (*iter).genJet()->eta();
564 lGenJet.phi = (*iter).genJet()->phi();
565 lGenJet.charge = (*iter).genJet()->charge();
566 lGenJet.pdgId = (*iter).genJet()->pdgId();
567 lGenJet.status = (*iter).genJet()->status();
568 lGenJet.mass = (*iter).genJet()->mass();
569 lGenJet.vx = (*iter).genJet()->vx();
570 lGenJet.vy = (*iter).genJet()->vy();
571 lGenJet.vz = (*iter).genJet()->vz();
572 }
573 else {
574 lGenJet.valid = false;
575 lGenJet.pT = 0;
576 lGenJet.eta = 0;
577 lGenJet.phi = 0;
578 lGenJet.charge = 0;
579 lGenJet.pdgId = 0;
580 lGenJet.status = 0;
581 lGenJet.mass = 0;
582 lGenJet.vx = 0;
583 lGenJet.vy = 0;
584 lGenJet.vz = 0;
585 }
586
587 HbbAnalysis::BaseVars lReco;
588 lReco.E = (*iter).energy();
589 lReco.pT = (*iter).pt();
590 lReco.eta = (*iter).eta();
591 lReco.phi = (*iter).phi();
592 lReco.charge = (*iter).charge();
593 lReco.vx = (*iter).vx();
594 lReco.vy = (*iter).vy();
595 lReco.vz = (*iter).vz();
596
597 bool isCalo = (*iter).isCaloTau();
598 bool isPF = (*iter).isPFTau();
599
600 assert (isCalo == !isPF);
601
602 HbbAnalysis::TauLeadTrkVars lLead;
603 if (isCalo) {
604 lLead.signedSipt = (*iter).leadTracksignedSipt();
605
606 reco::TrackRef lCaloTrk = (*iter).leadTrack();
607
608 if ( lCaloTrk.isAvailable() && lCaloTrk.isNonnull() ) {
609 lLead.pT = lCaloTrk->pt();
610 lLead.eta = lCaloTrk->eta();
611 lLead.phi = lCaloTrk->phi();
612
613 lLead.matchDist = reco::deltaR(lCaloTrk->momentum(), (*iter).p4());
614
615 if ( aRecoVertices->size() >= 1 ) {
616 const reco::Vertex& thePrimaryEventVertex = (*aRecoVertices->begin());
617 lLead.IPxy = lCaloTrk->dxy(thePrimaryEventVertex.position());
618 lLead.IPz = lCaloTrk->dz(thePrimaryEventVertex.position());
619 }
620 else {
621 lLead.IPxy = 0;
622 lLead.IPz = 0;
623 }
624 }
625 else {
626 lLead.pT = 0;
627 lLead.eta = 0;
628 lLead.phi = 0;
629 lLead.matchDist = 0;
630 lLead.IPxy = 0;
631 lLead.IPz = 0;
632 lLead.signedSipt = 0;
633 }
634 }//calo
635 else {
636 lLead.signedSipt = (*iter).leadPFChargedHadrCandsignedSipt();
637
638 reco::PFCandidateRef lHadrCand = (*iter).leadPFChargedHadrCand();
639
640 if ( lHadrCand.isAvailable() && lHadrCand.isNonnull() ) {
641 lLead.pT = lHadrCand->pt();
642 lLead.eta = lHadrCand->eta();
643 lLead.phi = lHadrCand->phi();
644
645 lLead.matchDist = reco::deltaR(lHadrCand->p4(), (*iter).p4());
646
647 reco::TrackRef lTrk = lHadrCand->trackRef();
648
649 if ( aRecoVertices->size() >= 1) {
650 const reco::Vertex& thePrimaryEventVertex = (*aRecoVertices->begin());
651 lLead.IPxy = lTrk->dxy(thePrimaryEventVertex.position());
652 lLead.IPz = lTrk->dz(thePrimaryEventVertex.position());
653 }
654 else {
655 lLead.IPxy = 0;
656 lLead.IPz = 0;
657 }
658 }
659 else {
660 lLead.pT = 0;
661 lLead.eta = 0;
662 lLead.phi = 0;
663 lLead.matchDist = 0;
664 lLead.IPxy = 0;
665 lLead.IPz = 0;
666 lLead.signedSipt = 0;
667 }
668
669 }//pf
670
671
672
673 const std::vector<std::pair<std::string,float> > & lIDs = (*iter).tauIDs();
674 bool lPrint = false;
675 //a few warning if additional discriminants are found:
676 if ((isPF && lIDs.size() != 7) ||
677 (isCalo && lIDs.size() != 3))
678 lPrint = true;
679
680 if (lPrint) {
681 std::cout << "!!!!!!! Discriminants changed, please update histograms !!!!!!!" << std::endl;
682 std::cout << "--- isCaloTau = " << isCalo << ", isPFTau = " << isPF << std::endl;
683 std::cout << "------ ID names = " << std::endl;
684 }
685
686
687 HbbAnalysis::PFTauIDVars lpfId;
688 HbbAnalysis::CaloTauIDVars lcaloId;
689
690 for (unsigned int id(0); id<lIDs.size(); id++){
691
692 std::string lName = lIDs.at(id).first;
693 float lDiscri = lIDs.at(id).second;
694
695 if (lPrint) std::cout << "--------- " << lName << " = " << lDiscri << std::endl;
696 //pf
697
698 if (isPF) {
699 if (lName.find("leadingTrackFinding") != lName.npos) lpfId.byLeadingTrackFinding = lDiscri;
700 if (lName.find("leadingTrackPtCut") != lName.npos) lpfId.byLeadingTrackPtCut = lDiscri;
701 if (lName.find("trackIsolation") != lName.npos) lpfId.byTrackIsolation = lDiscri;
702 if (lName.find("ecalIsolation") != lName.npos) lpfId.byECALIsolation = lDiscri;
703 if (lName.find("byIsolation") != lName.npos) lpfId.byIsolation = lDiscri;
704 if (lName.find("againstElectron") != lName.npos) lpfId.againstElectron = lDiscri;
705 if (lName.find("againstMuon") != lName.npos) lpfId.againstMuon = lDiscri;
706 }
707 if (isCalo){
708 if (lName.find("byIsolation") != lName.npos) lcaloId.byIsolation = lDiscri;
709 if (lName.find("leadingTrackFinding") != lName.npos) lcaloId.byLeadingTrackFinding = lDiscri;
710 if (lName.find("leadingTrackPtCut") != lName.npos) lcaloId.byLeadingTrackPtCut = lDiscri;
711 }
712
713
714 }
715
716
717 if (isCalo){
718
719 HbbAnalysis::CaloTauIsoVars lcaloIso;
720 lcaloIso.nIsoTracks = (*iter).isolationTracks().size();
721 lcaloIso.nSigTracks = (*iter).signalTracks().size();
722 lcaloIso.leadTrackHCAL3x3hitsEtSum = (*iter).leadTrackHCAL3x3hitsEtSum();
723 lcaloIso.leadTrackHCAL3x3hottesthitDEta = (*iter).leadTrackHCAL3x3hottesthitDEta();
724 lcaloIso.signalTracksInvariantMass = (*iter).signalTracksInvariantMass();
725 lcaloIso.tracksInvariantMass = (*iter).TracksInvariantMass();
726 lcaloIso.isolationTracksPtSum = (*iter).isolationTracksPtSum();
727 lcaloIso.isolationECALhitsEtSum = (*iter).isolationECALhitsEtSum();
728 lcaloIso.maximumHCALhitEt = (*iter).maximumHCALhitEt();
729
730 HbbAnalysis::Tau lObj(lGen,lGenJet,lReco,lLead,lcaloIso,lcaloId);
731 aVec.push_back(lObj);
732 }
733 else {
734
735 HbbAnalysis::PFTauIsoVars lpfIso;
736 lpfIso.nSigCands = (*iter).signalPFCands().size();
737 lpfIso.nIsoCands = (*iter).isolationPFCands().size();
738 lpfIso.maximumHCALPFClusterEt = (*iter).maximumHCALPFClusterEt();
739 //lpfIso.emFraction = (*iter).emFraction();
740 lpfIso.hcalTotOverPLead = (*iter).hcalTotOverPLead();
741 lpfIso.hcalMaxOverPLead = (*iter).hcalMaxOverPLead();
742 lpfIso.hcal3x3OverPLead = (*iter).hcal3x3OverPLead();
743 lpfIso.ecalStripSumEOverPLead = (*iter).ecalStripSumEOverPLead();
744 lpfIso.bremsRecoveryEOverPLead = (*iter).bremsRecoveryEOverPLead();
745
746 HbbAnalysis::PFTauCandVars lHadr;
747 lHadr.nSigCands = (*iter).signalPFChargedHadrCands().size();
748 lHadr.nIsoCands = (*iter).isolationPFChargedHadrCands().size();
749 lHadr.isolationPtSum = (*iter).isolationPFChargedHadrCandsPtSum();
750
751 HbbAnalysis::PFTauCandVars lNeutr;
752 lNeutr.nSigCands = (*iter).signalPFNeutrHadrCands().size();
753 lNeutr.nIsoCands = (*iter).isolationPFNeutrHadrCands().size();
754 lNeutr.isolationPtSum = 0;
755
756 HbbAnalysis::PFTauCandVars lGamma;
757 lGamma.nSigCands = (*iter).signalPFGammaCands().size();
758 lGamma.nIsoCands = (*iter).isolationPFGammaCands().size();
759 lGamma.isolationPtSum = (*iter).isolationPFGammaCandsEtSum();
760
761 HbbAnalysis::PFTauEleIDVars lEleId;
762 if ( (*iter).electronPreIDTrack().isAvailable() && (*iter).electronPreIDTrack().isNonnull() ) {
763 lEleId.pT = (*iter).electronPreIDTrack()->pt();
764 lEleId.eta = (*iter).electronPreIDTrack()->eta();
765 lEleId.phi = (*iter).electronPreIDTrack()->phi();
766 }
767 else {
768 lEleId.pT = 0;
769 lEleId.eta = 0;
770 lEleId.phi = 0;
771 }
772
773 lEleId.output = (*iter).electronPreIDOutput();
774 lEleId.decision = (*iter).electronPreIDDecision();
775
776 HbbAnalysis::PFTauMuIDVars lMuId;
777 lMuId.caloCompat = (*iter).caloComp();
778 lMuId.segCompat = (*iter).segComp();
779 lMuId.decision = (*iter).muonDecision();
780
781 HbbAnalysis::Tau lObj(lGen,lGenJet,lReco,lLead,lpfIso,lHadr,lNeutr,lGamma,lpfId,lEleId,lMuId);
782 aVec.push_back(lObj);
783 }
784 iEle++;
785 }
786 }
787 }
788
789 }//HbbTaus
790
791 void HbbTreeMaker::HbbJets(const edm::Handle<std::vector<pat::Jet> > & aCol,
792 const HbbAnalysis::JetFlavour & aJetFlav,
793 const edm::Handle<reco::GenParticleCollection> & aGenParticles,
794 std::vector<HbbAnalysis::Jet> & aVec)
795 {//HbbJets
796
797 if (aCol.isValid()){//valid
798 if (aCol->size() > 0) {//non empty
799
800 unsigned int iEle = 0;
801 for (std::vector<pat::Jet>::const_iterator iter = aCol->begin();
802 iter != aCol->end();
803 iter++)
804 {//loop on element
805 //if (debug_ > 1) std::cout << "**** Ele #" << iEle << ", pT,eta,phi = " << (*iter).pt() << " " << (*iter).eta() << " " << (*iter).phi() << std::endl;
806
807 HbbAnalysis::GenVars lGen;
808 if ((*iter).genParton()){
809 lGen.valid = true;
810 lGen.E = (*iter).genParton()->energy();
811 lGen.pT = (*iter).genParton()->pt();
812 lGen.eta = (*iter).genParton()->eta();
813 lGen.phi = (*iter).genParton()->phi();
814 lGen.charge = (*iter).genParton()->charge();
815 lGen.pdgId = (*iter).genParton()->pdgId();
816 lGen.status = (*iter).genParton()->status();
817 lGen.mass = (*iter).genParton()->mass();
818 lGen.vx = (*iter).genParton()->vx();
819 lGen.vy = (*iter).genParton()->vy();
820 lGen.vz = (*iter).genParton()->vz();
821 }
822 else {
823 lGen.valid = false;
824 lGen.E = 0;
825 lGen.pT = 0;
826 lGen.eta = 0;
827 lGen.phi = 0;
828 lGen.charge = 0;
829 lGen.pdgId = 0;
830 lGen.status = 0;
831 lGen.mass = 0;
832 lGen.vx = 0;
833 lGen.vy = 0;
834 lGen.vz = 0;
835 }
836
837
838 HbbAnalysis::GenVars lGenJet;
839 if ((*iter).genJet()){
840 lGenJet.valid = true;
841 lGenJet.pT = (*iter).genJet()->pt();
842 lGenJet.eta = (*iter).genJet()->eta();
843 lGenJet.phi = (*iter).genJet()->phi();
844 lGenJet.charge = (*iter).genJet()->charge();
845 lGenJet.pdgId = (*iter).genJet()->pdgId();
846 lGenJet.status = (*iter).genJet()->status();
847 lGenJet.mass = (*iter).genJet()->mass();
848 lGenJet.vx = (*iter).genJet()->vx();
849 lGenJet.vy = (*iter).genJet()->vy();
850 lGenJet.vz = (*iter).genJet()->vz();
851 }
852 else {
853 lGenJet.valid = false;
854 lGenJet.pT = 0;
855 lGenJet.eta = 0;
856 lGenJet.phi = 0;
857 lGenJet.charge = 0;
858 lGenJet.pdgId = 0;
859 lGenJet.status = 0;
860 lGenJet.mass = 0;
861 lGenJet.vx = 0;
862 lGenJet.vy = 0;
863 lGenJet.vz = 0;
864 }
865
866 HbbAnalysis::BaseVars lReco;
867 lReco.E = (*iter).energy();
868 lReco.pT = (*iter).pt();
869 lReco.eta = (*iter).eta();
870 lReco.phi = (*iter).phi();
871 lReco.charge = (*iter).jetCharge();
872 lReco.vx = (*iter).vx();
873 lReco.vy = (*iter).vy();
874 lReco.vz = (*iter).vz();
875
876 unsigned int lFlav = 0;
877 if (aJetFlav.nPartons()) {
878 lFlav = aJetFlav.partonMatchingGenJet((*iter),aGenParticles,0.4).second;
879 }
880
881 HbbAnalysis::JetVars lCommon;
882 lCommon.flavour = lFlav;
883 lCommon.partonFlavour = (*iter).partonFlavour();
884 lCommon.nAssociatedTracks = ((*iter).associatedTracks()).size();
885 lCommon.rawpT = (*iter).pt();
886 if ((*iter).hasJetCorrFactors())
887 lCommon.rawpT = (*iter).pt()/((*iter).jetCorrFactors().scaleDefault());
888 //lCommon.rawEta = (*iter).;
889 //lCommon.rawPhi = (*iter).;
890 // p_hasJetCorrFactors->Fill(aJet.hasJetCorrFactors());
891
892 HbbAnalysis::JetBtagVars lBtag;
893 lBtag.cSV = (*iter).bDiscriminator("combinedSecondaryVertexBJetTags");
894 lBtag.cSVMVA = (*iter).bDiscriminator("combinedSecondaryVertexMVABJetTags");
895 lBtag.iPMVA = (*iter).bDiscriminator("impactParameterMVABJetTags");
896 lBtag.bProba = (*iter).bDiscriminator("jetBProbabilityBJetTags");
897 lBtag.probability = (*iter).bDiscriminator("jetProbabilityBJetTags");
898 lBtag.sSV = (*iter).bDiscriminator("simpleSecondaryVertexBJetTags");
899 lBtag.softElectron = (*iter).bDiscriminator("softElectronBJetTags");
900 lBtag.softMuon = (*iter).bDiscriminator("softMuonBJetTags");
901 lBtag.softMuonNoIP = (*iter).bDiscriminator("softMuonNoIPBJetTags");
902 lBtag.tCHE = (*iter).bDiscriminator("trackCountingHighEffBJetTags");
903 lBtag.tCHP = (*iter).bDiscriminator("trackCountingHighPurBJetTags");
904
905 bool isCalo = (*iter).isCaloJet();
906 bool isPF = (*iter).isPFJet();
907
908 assert (isCalo == !isPF);
909 if (isCalo) {
910 HbbAnalysis::CaloJetVars lCalo;
911 lCalo.maxEInEmTowers = (*iter).maxEInEmTowers();
912 lCalo.maxEInHadTowers = (*iter).maxEInHadTowers();
913 lCalo.energyFractionHadronic = (*iter).energyFractionHadronic();
914 lCalo.emEnergyFraction = (*iter).emEnergyFraction();
915 lCalo.hadEnergyInHB = (*iter).hadEnergyInHB();
916 lCalo.hadEnergyInHO = (*iter).hadEnergyInHO();
917 lCalo.hadEnergyInHE = (*iter).hadEnergyInHE();
918 lCalo.hadEnergyInHF = (*iter).hadEnergyInHF();
919 lCalo.emEnergyInEB = (*iter).emEnergyInEB();
920 lCalo.emEnergyInEE = (*iter).emEnergyInEE();
921 lCalo.emEnergyInHF = (*iter).emEnergyInHF();
922 lCalo.towersArea = (*iter).towersArea();
923 lCalo.n90 = (*iter).n90();
924 lCalo.n60 = (*iter).n60();
925
926 HbbAnalysis::Jet lObj(lGen,lGenJet,lReco,lCommon,lCalo,lBtag);
927 aVec.push_back(lObj);
928
929 }
930
931 if (isPF){
932 HbbAnalysis::PFJetVars lPF;
933 lPF.chargedHadronEnergy = (*iter).chargedHadronEnergy();
934 lPF.chargedHadronEnergyFraction = (*iter).chargedHadronEnergyFraction();
935 lPF.neutralHadronEnergy = (*iter).neutralHadronEnergy();
936 lPF.neutralHadronEnergyFraction = (*iter).neutralHadronEnergyFraction();
937 lPF.chargedEmEnergy = (*iter).chargedEmEnergy();
938 lPF.chargedEmEnergyFraction = (*iter).chargedEmEnergyFraction();
939 lPF.chargedMuEnergy = (*iter).chargedMuEnergy();
940 lPF.chargedMuEnergyFraction = (*iter).chargedMuEnergyFraction();
941 lPF.neutralEmEnergy = (*iter).neutralEmEnergy();
942 lPF.neutralEmEnergyFraction = (*iter).neutralEmEnergyFraction();
943 lPF.chargedMultiplicity = (*iter).chargedMultiplicity();
944 lPF.neutralMultiplicity = (*iter).neutralMultiplicity();
945 lPF.muonMultiplicity = (*iter).muonMultiplicity();
946
947 HbbAnalysis::Jet lObj(lGen,lGenJet,lReco,lCommon,lPF,lBtag);
948 aVec.push_back(lObj);
949
950 }
951
952 iEle++;
953 }//loop on element
954 }//non empty
955 }//valid
956
957 }//HbbJets
958
959 void HbbTreeMaker::HbbMet(const edm::Handle<std::vector<pat::MET> > & aCol,
960 HbbAnalysis::Met & aMet)
961 {//HbbMet
962 HbbAnalysis::MetVars lGen;
963 HbbAnalysis::MetVars lReco;
964
965 assert(aCol->size() == 1);
966 const pat::MET & lMet = *(aCol->begin());
967
968 lReco.mET = lMet.pt();
969 lReco.mEx = lMet.px();
970 lReco.mEy = lMet.py();
971 lReco.sumET = lMet.sumEt();
972 lReco.phi = lMet.phi();
973 lReco.mEtSig = lMet.mEtSig();
974
975 const reco::GenMET *lGenMET = lMet.genMET();
976
977 if (lGenMET){
978 lGen.mET = lGenMET->pt();
979 lGen.mEx = lGenMET->px();
980 lGen.mEy = lGenMET->py();
981 lGen.sumET = lGenMET->sumEt();
982 lGen.phi = lGenMET->phi();
983 lGen.mEtSig = lGenMET->mEtSig();
984 }
985 else {
986 lGen.mET = 0;
987 lGen.mEx = 0;
988 lGen.mEy = 0;
989 lGen.sumET = 0;
990 lGen.phi = 0;
991 lGen.mEtSig = 0;
992 }
993
994 aMet.genVars(lGen);
995 aMet.recoVars(lReco);
996
997 }//HbbMet
998
999 void HbbTreeMaker::HbbTrigger(const edm::Handle<edm::TriggerResults> & aCol,
1000 std::vector<HbbAnalysis::Trigger> & aVec)
1001 {//HbbTrigger
1002
1003 if (aCol.isValid()){
1004 edm::TriggerNames lNames;
1005 lNames.init(*aCol);
1006
1007 //for (unsigned int j=0; j<hltPaths_.size(); j++) {
1008 //bool valid=false;
1009 for (unsigned int i=0; i<aCol->size(); i++){
1010 std::string trueName = lNames.triggerNames().at(i);
1011
1012 HbbAnalysis::TriggerVars lTrigVars;
1013 lTrigVars.name = trueName;
1014 lTrigVars.index = i;
1015 lTrigVars.accept = aCol->accept(i);
1016
1017 HbbAnalysis::Trigger lObj(lTrigVars);
1018
1019 aVec.push_back(lObj);
1020 }
1021
1022 }
1023
1024 }
1025
1026 void HbbTreeMaker::HbbParticles(const edm::Handle<reco::GenParticleCollection> & aCol,
1027 std::vector<HbbAnalysis::GenParticle> & aVec)
1028 {//genparticles
1029
1030 //add genParticle inside vector
1031 for (unsigned int mccount = 0;mccount<aCol->size();mccount++)
1032 {//loop on particles
1033 const reco::Candidate & p = (*aCol)[mccount];
1034
1035 HbbAnalysis::MCVars lMC;
1036 lMC.index = mccount;
1037 lMC.E = p.energy();
1038 lMC.pT = p.pt();
1039 lMC.eta = p.eta();
1040 lMC.phi = p.phi();
1041 lMC.pdgId = p.pdgId();
1042 lMC.status = p.status();
1043 HbbAnalysis::GenParticle lObj(lMC);
1044 aVec.push_back(lObj);
1045
1046 }//loop on particles
1047
1048 //now need to get the list of parents....
1049
1050 for (unsigned int mccount = 0;mccount<aCol->size();mccount++)
1051 {//loop on particles
1052 const reco::Candidate & p = (*aCol)[mccount];
1053
1054 unsigned int nD = p.numberOfDaughters();
1055 //std::cout << mccount << " has " << nD << " daughter(s) : " ;
1056 for(unsigned int dau=0;dau<nD;dau++){//loop on daughters
1057 const reco::Candidate * pDau = p.daughter(dau);
1058 //std::cout << pDau << " ("<<std::flush;
1059 for (unsigned int mccount2 = 0;mccount2<aCol->size();mccount2++)
1060 {//loop on particles
1061 const reco::Candidate & p2 = (*aCol)[mccount2];
1062 //std::cout<<"DBG: mccount2 = "<<mccount2<<" gen size = "<<data_->gen().size()<<std::endl;
1063
1064 HbbAnalysis::GenParticle & part2 = aVec.at(mccount2);
1065 if (pDau == &p2) {
1066 part2.setParent(mccount);
1067 //std::cout << &p2 << ", index = " << mccount2 << "), "<<std::flush;
1068 break;
1069 }
1070 //else std::cout << std::endl << "****no match " << mccount2 << " " << &p2 << std::endl;
1071 //if (p2.pdgId() == pDau->pdgId() && fabs(p2.energy()-pDau->energy())<0.000001 && fabs(p2.eta()-pDau->eta())<0.000001 && fabs(p2.phi()-pDau->phi())<0.000001){
1072 //part2.parent(mccount);
1073 //}
1074 }//loop on particles
1075 }//loop on daughters
1076 //std::cout << std::endl;
1077
1078 }//loop on particles
1079
1080 if (debug_){
1081 for (unsigned int imc(0); imc < aVec.size(); imc++){
1082 HbbAnalysis::GenParticle & p = aVec.at(imc);
1083 p.print();
1084 }
1085 }
1086
1087
1088 }//genparticles