ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/HbbAnalysis/plugins/HbbTreeMaker.cc
Revision: 1.1
Committed: Thu Sep 17 19:39:51 2009 UTC (15 years, 7 months ago) by amagnan
Content type: text/plain
Branch: MAIN
CVS Tags: v00-03-00
Log Message:
add Tree

File Contents

# User Rev Content
1 amagnan 1.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,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.pT = (*iter).genLepton()->pt();
318     lGen.eta = (*iter).genLepton()->eta();
319     lGen.phi = (*iter).genLepton()->phi();
320     lGen.charge = (*iter).genLepton()->charge();
321     lGen.pdgId = (*iter).genLepton()->pdgId();
322     lGen.status = (*iter).genLepton()->status();
323     lGen.mass = (*iter).genLepton()->mass();
324     lGen.vx = (*iter).genLepton()->vx();
325     lGen.vy = (*iter).genLepton()->vy();
326     lGen.vz = (*iter).genLepton()->vz();
327     }
328     else {
329     lGen.valid = false;
330     lGen.pT = 0;
331     lGen.eta = 0;
332     lGen.phi = 0;
333     lGen.charge = 0;
334     lGen.pdgId = 0;
335     lGen.status = 0;
336     lGen.mass = 0;
337     lGen.vx = 0;
338     lGen.vy = 0;
339     lGen.vz = 0;
340     }
341    
342     HbbAnalysis::BaseVars lReco;
343     lReco.pT = (*iter).pt();
344     lReco.eta = (*iter).eta();
345     lReco.phi = (*iter).phi();
346     lReco.charge = (*iter).charge();
347     lReco.vx = (*iter).vx();
348     lReco.vy = (*iter).vy();
349     lReco.vz = (*iter).vz();
350    
351     HbbAnalysis::SCVars lSC;
352     lSC.sigmaEtaEta = (*iter).scSigmaEtaEta();
353     lSC.sigmaIEtaIEta = (*iter).scSigmaIEtaIEta();
354     lSC.e1x5 = (*iter).scE1x5();
355     lSC.e2x5Max = (*iter).scE2x5Max();
356     lSC.e5x5 = (*iter).scE5x5();
357     lSC.eOverP = (*iter).eSuperClusterOverP();
358    
359     HbbAnalysis::EleIsoVars lIso;
360     lIso.calo = (*iter).caloIso();
361     lIso.track = (*iter).trackIso();
362     lIso.ecal = (*iter).ecalIso();
363     lIso.hcal = (*iter).hcalIso();
364    
365     HbbAnalysis::EleIDVars lID;
366     lID.electronIDs = (*iter).electronIDs();
367     lID.hOverE = (*iter).hadronicOverEm();
368     lID.deltaPhiIn = (*iter).deltaPhiSuperClusterTrackAtVtx();
369     lID.deltaEtaIn = (*iter).deltaEtaSuperClusterTrackAtVtx();
370    
371     HbbAnalysis::Electron lObj(lGen,lReco,lSC,lIso,lID);
372     aVec.push_back(lObj);
373     iEle++;
374     }
375     }
376     }
377    
378     }//HbbElectrons
379    
380    
381     void HbbTreeMaker::HbbMuons(const edm::Handle<std::vector<pat::Muon> > & aCol,
382     std::vector<HbbAnalysis::Muon> & aVec)
383     {//HbbMuons
384    
385     if (aCol.isValid()){
386     if (aCol->size() > 0) {
387    
388     unsigned int iEle = 0;
389     for (std::vector<pat::Muon>::const_iterator iter = aCol->begin();
390     iter != aCol->end();
391     iter++)
392     {
393     const reco::Muon* recoMuon = dynamic_cast<const reco::Muon*>((*iter).originalObject());
394    
395     HbbAnalysis::GenVars lGen;
396     if ((*iter).genLepton()){
397     lGen.valid = true;
398     lGen.pT = (*iter).genLepton()->pt();
399     lGen.eta = (*iter).genLepton()->eta();
400     lGen.phi = (*iter).genLepton()->phi();
401     lGen.charge = (*iter).genLepton()->charge();
402     lGen.pdgId = (*iter).genLepton()->pdgId();
403     lGen.status = (*iter).genLepton()->status();
404     lGen.mass = (*iter).genLepton()->mass();
405     lGen.vx = (*iter).genLepton()->vx();
406     lGen.vy = (*iter).genLepton()->vy();
407     lGen.vz = (*iter).genLepton()->vz();
408     }
409     else {
410     lGen.valid = false;
411     lGen.pT = 0;
412     lGen.eta = 0;
413     lGen.phi = 0;
414     lGen.charge = 0;
415     lGen.pdgId = 0;
416     lGen.status = 0;
417     lGen.mass = 0;
418     lGen.vx = 0;
419     lGen.vy = 0;
420     lGen.vz = 0;
421     }
422    
423     HbbAnalysis::BaseVars lReco;
424     lReco.pT = (*iter).pt();
425     lReco.eta = (*iter).eta();
426     lReco.phi = (*iter).phi();
427     lReco.charge = (*iter).charge();
428     lReco.vx = (*iter).vx();
429     lReco.vy = (*iter).vy();
430     lReco.vz = (*iter).vz();
431    
432     HbbAnalysis::MuIsoVars lIsoR03;
433     lIsoR03.sumPt = recoMuon->isolationR03().sumPt;
434     lIsoR03.emEt = recoMuon->isolationR03().emEt;
435     lIsoR03.hadEt = recoMuon->isolationR03().hadEt;
436     lIsoR03.nTracks = recoMuon->isolationR03().nTracks;
437     lIsoR03.nJets = recoMuon->isolationR03().nJets;
438    
439     HbbAnalysis::MuIsoVars lIsoR05;
440     lIsoR05.sumPt = recoMuon->isolationR05().sumPt;
441     lIsoR05.emEt = recoMuon->isolationR05().emEt;
442     lIsoR05.hadEt = recoMuon->isolationR05().hadEt;
443     lIsoR05.nTracks = recoMuon->isolationR05().nTracks;
444     lIsoR05.nJets = recoMuon->isolationR05().nJets;
445    
446     HbbAnalysis::MuIDVars lID;
447     if (recoMuon->isGlobalMuon()) lID.type = 1;
448     else if (recoMuon->isTrackerMuon()) lID.type = 2;
449     else if (recoMuon->isStandAloneMuon()) lID.type = 3;
450     else if (recoMuon->isCaloMuon()) lID.type = 4;
451     else lID.type = 0;
452    
453     if (recoMuon->isGood(reco::Muon::AllGlobalMuons)) lID.ids.push_back(1);
454     if (recoMuon->isGood(reco::Muon::AllStandAloneMuons)) lID.ids.push_back(2);
455     if (recoMuon->isGood(reco::Muon::AllTrackerMuons)) lID.ids.push_back(3);
456     if (recoMuon->isGood(reco::Muon::TrackerMuonArbitrated)) lID.ids.push_back(4);
457     if (recoMuon->isGood(reco::Muon::AllArbitrated)) lID.ids.push_back(5);
458     if (recoMuon->isGood(reco::Muon::GlobalMuonPromptTight)) lID.ids.push_back(6);
459     if (recoMuon->isGood(reco::Muon::TMLastStationLoose)) lID.ids.push_back(7);
460     if (recoMuon->isGood(reco::Muon::TMLastStationTight)) lID.ids.push_back(8);
461     if (recoMuon->isGood(reco::Muon::TM2DCompatibilityLoose)) lID.ids.push_back(9);
462     if (recoMuon->isGood(reco::Muon::TM2DCompatibilityTight)) lID.ids.push_back(10);
463     if (recoMuon->isGood(reco::Muon::TMOneStationLoose)) lID.ids.push_back(11);
464     if (recoMuon->isGood(reco::Muon::TMOneStationTight)) lID.ids.push_back(12);
465     if (recoMuon->isGood(reco::Muon::TMLastStationOptimizedLowPtLoose)) lID.ids.push_back(13);
466     if (recoMuon->isGood(reco::Muon::TMLastStationOptimizedLowPtTight)) lID.ids.push_back(14);
467    
468     lID.caloCompat = recoMuon->caloCompatibility();
469     lID.segCompat = recoMuon->segmentCompatibility();
470     lID.nChambers = recoMuon->numberOfChambers();
471     lID.nMatchesLoose = recoMuon->numberOfMatches(reco::Muon::NoArbitration);
472     lID.nMatchesMedium = recoMuon->numberOfMatches(reco::Muon::SegmentArbitration);
473     lID.nMatchesTight = recoMuon->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
474    
475     HbbAnalysis::Muon lObj(lGen,lReco,lIsoR03,lIsoR05,lID);
476     aVec.push_back(lObj);
477     iEle++;
478     }
479     }
480     }
481    
482     }//HbbMuons
483    
484     void HbbTreeMaker::HbbTaus(const edm::Handle<std::vector<pat::Tau> > & aCol,
485     const edm::Handle<std::vector<reco::Vertex> > & aRecoVertices,
486     std::vector<HbbAnalysis::Tau> & aVec)
487     {//HbbTaus
488    
489     if (aCol.isValid()){
490     if (aCol->size() > 0) {
491    
492     unsigned int iEle = 0;
493     for (std::vector<pat::Tau>::const_iterator iter = aCol->begin();
494     iter != aCol->end();
495     iter++)
496     {
497     //if (debug_ > 1) std::cout << "**** Ele #" << iEle << ", pT,eta,phi = " << (*iter).pt() << " " << (*iter).eta() << " " << (*iter).phi() << std::endl;
498    
499     HbbAnalysis::GenVars lGen;
500     if ((*iter).genLepton()){
501     lGen.valid = true;
502     lGen.pT = (*iter).genLepton()->pt();
503     lGen.eta = (*iter).genLepton()->eta();
504     lGen.phi = (*iter).genLepton()->phi();
505     lGen.charge = (*iter).genLepton()->charge();
506     lGen.pdgId = (*iter).genLepton()->pdgId();
507     lGen.status = (*iter).genLepton()->status();
508     lGen.mass = (*iter).genLepton()->mass();
509     lGen.vx = (*iter).genLepton()->vx();
510     lGen.vy = (*iter).genLepton()->vy();
511     lGen.vz = (*iter).genLepton()->vz();
512     }
513     else {
514     lGen.valid = false;
515     lGen.pT = 0;
516     lGen.eta = 0;
517     lGen.phi = 0;
518     lGen.charge = 0;
519     lGen.pdgId = 0;
520     lGen.status = 0;
521     lGen.mass = 0;
522     lGen.vx = 0;
523     lGen.vy = 0;
524     lGen.vz = 0;
525     }
526    
527    
528     HbbAnalysis::GenVars lGenJet;
529     if ((*iter).genJet()){
530     lGenJet.valid = true;
531     lGenJet.pT = (*iter).genJet()->pt();
532     lGenJet.eta = (*iter).genJet()->eta();
533     lGenJet.phi = (*iter).genJet()->phi();
534     lGenJet.charge = (*iter).genJet()->charge();
535     lGenJet.pdgId = (*iter).genJet()->pdgId();
536     lGenJet.status = (*iter).genJet()->status();
537     lGenJet.mass = (*iter).genJet()->mass();
538     lGenJet.vx = (*iter).genJet()->vx();
539     lGenJet.vy = (*iter).genJet()->vy();
540     lGenJet.vz = (*iter).genJet()->vz();
541     }
542     else {
543     lGenJet.valid = false;
544     lGenJet.pT = 0;
545     lGenJet.eta = 0;
546     lGenJet.phi = 0;
547     lGenJet.charge = 0;
548     lGenJet.pdgId = 0;
549     lGenJet.status = 0;
550     lGenJet.mass = 0;
551     lGenJet.vx = 0;
552     lGenJet.vy = 0;
553     lGenJet.vz = 0;
554     }
555    
556     HbbAnalysis::BaseVars lReco;
557     lReco.pT = (*iter).pt();
558     lReco.eta = (*iter).eta();
559     lReco.phi = (*iter).phi();
560     lReco.charge = (*iter).charge();
561     lReco.vx = (*iter).vx();
562     lReco.vy = (*iter).vy();
563     lReco.vz = (*iter).vz();
564    
565     bool isCalo = (*iter).isCaloTau();
566     bool isPF = (*iter).isPFTau();
567    
568     assert (isCalo == !isPF);
569    
570     HbbAnalysis::TauLeadTrkVars lLead;
571     if (isCalo) {
572     lLead.signedSipt = (*iter).leadTracksignedSipt();
573    
574     reco::TrackRef lCaloTrk = (*iter).leadTrack();
575    
576     if ( lCaloTrk.isAvailable() && lCaloTrk.isNonnull() ) {
577     lLead.pT = lCaloTrk->pt();
578     lLead.eta = lCaloTrk->eta();
579     lLead.phi = lCaloTrk->phi();
580    
581     lLead.matchDist = reco::deltaR(lCaloTrk->momentum(), (*iter).p4());
582    
583     if ( aRecoVertices->size() >= 1 ) {
584     const reco::Vertex& thePrimaryEventVertex = (*aRecoVertices->begin());
585     lLead.IPxy = lCaloTrk->dxy(thePrimaryEventVertex.position());
586     lLead.IPz = lCaloTrk->dz(thePrimaryEventVertex.position());
587     }
588     else {
589     lLead.IPxy = 0;
590     lLead.IPz = 0;
591     }
592     }
593     else {
594     lLead.pT = 0;
595     lLead.eta = 0;
596     lLead.phi = 0;
597     lLead.matchDist = 0;
598     lLead.IPxy = 0;
599     lLead.IPz = 0;
600     lLead.signedSipt = 0;
601     }
602     }//calo
603     else {
604     lLead.signedSipt = (*iter).leadPFChargedHadrCandsignedSipt();
605    
606     reco::PFCandidateRef lHadrCand = (*iter).leadPFChargedHadrCand();
607    
608     if ( lHadrCand.isAvailable() && lHadrCand.isNonnull() ) {
609     lLead.pT = lHadrCand->pt();
610     lLead.eta = lHadrCand->eta();
611     lLead.phi = lHadrCand->phi();
612    
613     lLead.matchDist = reco::deltaR(lHadrCand->p4(), (*iter).p4());
614    
615     reco::TrackRef lTrk = lHadrCand->trackRef();
616    
617     if ( aRecoVertices->size() >= 1) {
618     const reco::Vertex& thePrimaryEventVertex = (*aRecoVertices->begin());
619     lLead.IPxy = lTrk->dxy(thePrimaryEventVertex.position());
620     lLead.IPz = lTrk->dz(thePrimaryEventVertex.position());
621     }
622     else {
623     lLead.IPxy = 0;
624     lLead.IPz = 0;
625     }
626     }
627     else {
628     lLead.pT = 0;
629     lLead.eta = 0;
630     lLead.phi = 0;
631     lLead.matchDist = 0;
632     lLead.IPxy = 0;
633     lLead.IPz = 0;
634     lLead.signedSipt = 0;
635     }
636    
637     }//pf
638    
639    
640    
641     const std::vector<std::pair<std::string,float> > & lIDs = (*iter).tauIDs();
642     bool lPrint = false;
643     //a few warning if additional discriminants are found:
644     if ((isPF && lIDs.size() != 7) ||
645     (isCalo && lIDs.size() != 3))
646     lPrint = true;
647    
648     if (lPrint) {
649     std::cout << "!!!!!!! Discriminants changed, please update histograms !!!!!!!" << std::endl;
650     std::cout << "--- isCaloTau = " << isCalo << ", isPFTau = " << isPF << std::endl;
651     std::cout << "------ ID names = " << std::endl;
652     }
653    
654    
655     HbbAnalysis::PFTauIDVars lpfId;
656     HbbAnalysis::CaloTauIDVars lcaloId;
657    
658     for (unsigned int id(0); id<lIDs.size(); id++){
659    
660     std::string lName = lIDs.at(id).first;
661     float lDiscri = lIDs.at(id).second;
662    
663     if (lPrint) std::cout << "--------- " << lName << " = " << lDiscri << std::endl;
664     //pf
665    
666     if (isPF) {
667     if (lName.find("leadingTrackFinding") != lName.npos) lpfId.byLeadingTrackFinding = lDiscri;
668     if (lName.find("leadingTrackPtCut") != lName.npos) lpfId.byLeadingTrackPtCut = lDiscri;
669     if (lName.find("trackIsolation") != lName.npos) lpfId.byTrackIsolation = lDiscri;
670     if (lName.find("ecalIsolation") != lName.npos) lpfId.byECALIsolation = lDiscri;
671     if (lName.find("byIsolation") != lName.npos) lpfId.byIsolation = lDiscri;
672     if (lName.find("againstElectron") != lName.npos) lpfId.againstElectron = lDiscri;
673     if (lName.find("againstMuon") != lName.npos) lpfId.againstMuon = lDiscri;
674     }
675     if (isCalo){
676     if (lName.find("byIsolation") != lName.npos) lcaloId.byIsolation = lDiscri;
677     if (lName.find("leadingTrackFinding") != lName.npos) lcaloId.byLeadingTrackFinding = lDiscri;
678     if (lName.find("leadingTrackPtCut") != lName.npos) lcaloId.byLeadingTrackPtCut = lDiscri;
679     }
680    
681    
682     }
683    
684    
685     if (isCalo){
686    
687     HbbAnalysis::CaloTauIsoVars lcaloIso;
688     lcaloIso.nIsoTracks = (*iter).isolationTracks().size();
689     lcaloIso.nSigTracks = (*iter).signalTracks().size();
690     lcaloIso.leadTrackHCAL3x3hitsEtSum = (*iter).leadTrackHCAL3x3hitsEtSum();
691     lcaloIso.leadTrackHCAL3x3hottesthitDEta = (*iter).leadTrackHCAL3x3hottesthitDEta();
692     lcaloIso.signalTracksInvariantMass = (*iter).signalTracksInvariantMass();
693     lcaloIso.tracksInvariantMass = (*iter).TracksInvariantMass();
694     lcaloIso.isolationTracksPtSum = (*iter).isolationTracksPtSum();
695     lcaloIso.isolationECALhitsEtSum = (*iter).isolationECALhitsEtSum();
696     lcaloIso.maximumHCALhitEt = (*iter).maximumHCALhitEt();
697    
698     HbbAnalysis::Tau lObj(lGen,lGenJet,lReco,lLead,lcaloIso,lcaloId);
699     aVec.push_back(lObj);
700     }
701     else {
702    
703     HbbAnalysis::PFTauIsoVars lpfIso;
704     lpfIso.nSigCands = (*iter).signalPFCands().size();
705     lpfIso.nIsoCands = (*iter).isolationPFCands().size();
706     lpfIso.maximumHCALPFClusterEt = (*iter).maximumHCALPFClusterEt();
707     //lpfIso.emFraction = (*iter).emFraction();
708     lpfIso.hcalTotOverPLead = (*iter).hcalTotOverPLead();
709     lpfIso.hcalMaxOverPLead = (*iter).hcalMaxOverPLead();
710     lpfIso.hcal3x3OverPLead = (*iter).hcal3x3OverPLead();
711     lpfIso.ecalStripSumEOverPLead = (*iter).ecalStripSumEOverPLead();
712     lpfIso.bremsRecoveryEOverPLead = (*iter).bremsRecoveryEOverPLead();
713    
714     HbbAnalysis::PFTauCandVars lHadr;
715     lHadr.nSigCands = (*iter).signalPFChargedHadrCands().size();
716     lHadr.nIsoCands = (*iter).isolationPFChargedHadrCands().size();
717     lHadr.isolationPtSum = (*iter).isolationPFChargedHadrCandsPtSum();
718    
719     HbbAnalysis::PFTauCandVars lNeutr;
720     lNeutr.nSigCands = (*iter).signalPFNeutrHadrCands().size();
721     lNeutr.nIsoCands = (*iter).isolationPFNeutrHadrCands().size();
722     lNeutr.isolationPtSum = 0;
723    
724     HbbAnalysis::PFTauCandVars lGamma;
725     lGamma.nSigCands = (*iter).signalPFGammaCands().size();
726     lGamma.nIsoCands = (*iter).isolationPFGammaCands().size();
727     lGamma.isolationPtSum = (*iter).isolationPFGammaCandsEtSum();
728    
729     HbbAnalysis::PFTauEleIDVars lEleId;
730     if ( (*iter).electronPreIDTrack().isAvailable() && (*iter).electronPreIDTrack().isNonnull() ) {
731     lEleId.pT = (*iter).electronPreIDTrack()->pt();
732     lEleId.eta = (*iter).electronPreIDTrack()->eta();
733     lEleId.phi = (*iter).electronPreIDTrack()->phi();
734     }
735     else {
736     lEleId.pT = 0;
737     lEleId.eta = 0;
738     lEleId.phi = 0;
739     }
740    
741     lEleId.output = (*iter).electronPreIDOutput();
742     lEleId.decision = (*iter).electronPreIDDecision();
743    
744     HbbAnalysis::PFTauMuIDVars lMuId;
745     lMuId.caloCompat = (*iter).caloComp();
746     lMuId.segCompat = (*iter).segComp();
747     lMuId.decision = (*iter).muonDecision();
748    
749     HbbAnalysis::Tau lObj(lGen,lGenJet,lReco,lLead,lpfIso,lHadr,lNeutr,lGamma,lpfId,lEleId,lMuId);
750     aVec.push_back(lObj);
751     }
752     iEle++;
753     }
754     }
755     }
756    
757     }//HbbTaus
758    
759     void HbbTreeMaker::HbbJets(const edm::Handle<std::vector<pat::Jet> > & aCol,
760     const HbbAnalysis::JetFlavour & aJetFlav,
761     const edm::Handle<reco::GenParticleCollection> & aGenParticles,
762     std::vector<HbbAnalysis::Jet> & aVec)
763     {//HbbJets
764    
765     if (aCol.isValid()){//valid
766     if (aCol->size() > 0) {//non empty
767    
768     unsigned int iEle = 0;
769     for (std::vector<pat::Jet>::const_iterator iter = aCol->begin();
770     iter != aCol->end();
771     iter++)
772     {//loop on element
773     //if (debug_ > 1) std::cout << "**** Ele #" << iEle << ", pT,eta,phi = " << (*iter).pt() << " " << (*iter).eta() << " " << (*iter).phi() << std::endl;
774    
775     HbbAnalysis::GenVars lGen;
776     if ((*iter).genParton()){
777     lGen.valid = true;
778     lGen.pT = (*iter).genParton()->pt();
779     lGen.eta = (*iter).genParton()->eta();
780     lGen.phi = (*iter).genParton()->phi();
781     lGen.charge = (*iter).genParton()->charge();
782     lGen.pdgId = (*iter).genParton()->pdgId();
783     lGen.status = (*iter).genParton()->status();
784     lGen.mass = (*iter).genParton()->mass();
785     lGen.vx = (*iter).genParton()->vx();
786     lGen.vy = (*iter).genParton()->vy();
787     lGen.vz = (*iter).genParton()->vz();
788     }
789     else {
790     lGen.valid = false;
791     lGen.pT = 0;
792     lGen.eta = 0;
793     lGen.phi = 0;
794     lGen.charge = 0;
795     lGen.pdgId = 0;
796     lGen.status = 0;
797     lGen.mass = 0;
798     lGen.vx = 0;
799     lGen.vy = 0;
800     lGen.vz = 0;
801     }
802    
803    
804     HbbAnalysis::GenVars lGenJet;
805     if ((*iter).genJet()){
806     lGenJet.valid = true;
807     lGenJet.pT = (*iter).genJet()->pt();
808     lGenJet.eta = (*iter).genJet()->eta();
809     lGenJet.phi = (*iter).genJet()->phi();
810     lGenJet.charge = (*iter).genJet()->charge();
811     lGenJet.pdgId = (*iter).genJet()->pdgId();
812     lGenJet.status = (*iter).genJet()->status();
813     lGenJet.mass = (*iter).genJet()->mass();
814     lGenJet.vx = (*iter).genJet()->vx();
815     lGenJet.vy = (*iter).genJet()->vy();
816     lGenJet.vz = (*iter).genJet()->vz();
817     }
818     else {
819     lGenJet.valid = false;
820     lGenJet.pT = 0;
821     lGenJet.eta = 0;
822     lGenJet.phi = 0;
823     lGenJet.charge = 0;
824     lGenJet.pdgId = 0;
825     lGenJet.status = 0;
826     lGenJet.mass = 0;
827     lGenJet.vx = 0;
828     lGenJet.vy = 0;
829     lGenJet.vz = 0;
830     }
831    
832     HbbAnalysis::BaseVars lReco;
833     lReco.pT = (*iter).pt();
834     lReco.eta = (*iter).eta();
835     lReco.phi = (*iter).phi();
836     lReco.charge = (*iter).jetCharge();
837     lReco.vx = (*iter).vx();
838     lReco.vy = (*iter).vy();
839     lReco.vz = (*iter).vz();
840    
841     unsigned int lFlav = 0;
842     if (aJetFlav.nPartons()) {
843     lFlav = aJetFlav.partonMatchingGenJet((*iter),aGenParticles,0.4).second;
844     }
845    
846     HbbAnalysis::JetVars lCommon;
847     lCommon.flavour = lFlav;
848     lCommon.partonFlavour = (*iter).partonFlavour();
849     lCommon.nAssociatedTracks = ((*iter).associatedTracks()).size();
850     lCommon.rawpT = (*iter).pt();
851     if ((*iter).hasJetCorrFactors())
852     lCommon.rawpT = (*iter).pt()/((*iter).jetCorrFactors().scaleDefault());
853     //lCommon.rawEta = (*iter).;
854     //lCommon.rawPhi = (*iter).;
855     // p_hasJetCorrFactors->Fill(aJet.hasJetCorrFactors());
856    
857     HbbAnalysis::JetBtagVars lBtag;
858     lBtag.cSV = (*iter).bDiscriminator("combinedSecondaryVertexBJetTags");
859     lBtag.cSVMVA = (*iter).bDiscriminator("combinedSecondaryVertexMVABJetTags");
860     lBtag.iPMVA = (*iter).bDiscriminator("impactParameterMVABJetTags");
861     lBtag.bProba = (*iter).bDiscriminator("jetBProbabilityBJetTags");
862     lBtag.probability = (*iter).bDiscriminator("jetProbabilityBJetTags");
863     lBtag.sSV = (*iter).bDiscriminator("simpleSecondaryVertexBJetTags");
864     lBtag.softElectron = (*iter).bDiscriminator("softElectronBJetTags");
865     lBtag.softMuon = (*iter).bDiscriminator("softMuonBJetTags");
866     lBtag.softMuonNoIP = (*iter).bDiscriminator("softMuonNoIPBJetTags");
867     lBtag.tCHE = (*iter).bDiscriminator("trackCountingHighEffBJetTags");
868     lBtag.tCHP = (*iter).bDiscriminator("trackCountingHighPurBJetTags");
869    
870     bool isCalo = (*iter).isCaloJet();
871     bool isPF = (*iter).isPFJet();
872    
873     assert (isCalo == !isPF);
874     if (isCalo) {
875     HbbAnalysis::CaloJetVars lCalo;
876     lCalo.maxEInEmTowers = (*iter).maxEInEmTowers();
877     lCalo.maxEInHadTowers = (*iter).maxEInHadTowers();
878     lCalo.energyFractionHadronic = (*iter).energyFractionHadronic();
879     lCalo.emEnergyFraction = (*iter).emEnergyFraction();
880     lCalo.hadEnergyInHB = (*iter).hadEnergyInHB();
881     lCalo.hadEnergyInHO = (*iter).hadEnergyInHO();
882     lCalo.hadEnergyInHE = (*iter).hadEnergyInHE();
883     lCalo.hadEnergyInHF = (*iter).hadEnergyInHF();
884     lCalo.emEnergyInEB = (*iter).emEnergyInEB();
885     lCalo.emEnergyInEE = (*iter).emEnergyInEE();
886     lCalo.emEnergyInHF = (*iter).emEnergyInHF();
887     lCalo.towersArea = (*iter).towersArea();
888     lCalo.n90 = (*iter).n90();
889     lCalo.n60 = (*iter).n60();
890    
891     HbbAnalysis::Jet lObj(lGen,lGenJet,lReco,lCommon,lCalo,lBtag);
892     aVec.push_back(lObj);
893    
894     }
895    
896     if (isPF){
897     HbbAnalysis::PFJetVars lPF;
898     lPF.chargedHadronEnergy = (*iter).chargedHadronEnergy();
899     lPF.chargedHadronEnergyFraction = (*iter).chargedHadronEnergyFraction();
900     lPF.neutralHadronEnergy = (*iter).neutralHadronEnergy();
901     lPF.neutralHadronEnergyFraction = (*iter).neutralHadronEnergyFraction();
902     lPF.chargedEmEnergy = (*iter).chargedEmEnergy();
903     lPF.chargedEmEnergyFraction = (*iter).chargedEmEnergyFraction();
904     lPF.chargedMuEnergy = (*iter).chargedMuEnergy();
905     lPF.chargedMuEnergyFraction = (*iter).chargedMuEnergyFraction();
906     lPF.neutralEmEnergy = (*iter).neutralEmEnergy();
907     lPF.neutralEmEnergyFraction = (*iter).neutralEmEnergyFraction();
908     lPF.chargedMultiplicity = (*iter).chargedMultiplicity();
909     lPF.neutralMultiplicity = (*iter).neutralMultiplicity();
910     lPF.muonMultiplicity = (*iter).muonMultiplicity();
911    
912     HbbAnalysis::Jet lObj(lGen,lGenJet,lReco,lCommon,lPF,lBtag);
913     aVec.push_back(lObj);
914    
915     }
916    
917     iEle++;
918     }//loop on element
919     }//non empty
920     }//valid
921    
922     }//HbbJets
923    
924     void HbbTreeMaker::HbbMet(const edm::Handle<std::vector<pat::MET> > & aCol,
925     HbbAnalysis::Met & aMet)
926     {//HbbMet
927     HbbAnalysis::MetVars lGen;
928     HbbAnalysis::MetVars lReco;
929    
930     assert(aCol->size() == 1);
931     const pat::MET & lMet = *(aCol->begin());
932    
933     lReco.mET = lMet.pt();
934     lReco.mEx = lMet.px();
935     lReco.mEy = lMet.py();
936     lReco.sumET = lMet.sumEt();
937     lReco.phi = lMet.phi();
938     lReco.mEtSig = lMet.mEtSig();
939    
940     const reco::GenMET *lGenMET = lMet.genMET();
941    
942     if (lGenMET){
943     lGen.mET = lGenMET->pt();
944     lGen.mEx = lGenMET->px();
945     lGen.mEy = lGenMET->py();
946     lGen.sumET = lGenMET->sumEt();
947     lGen.phi = lGenMET->phi();
948     lGen.mEtSig = lGenMET->mEtSig();
949     }
950     else {
951     lGen.mET = 0;
952     lGen.mEx = 0;
953     lGen.mEy = 0;
954     lGen.sumET = 0;
955     lGen.phi = 0;
956     lGen.mEtSig = 0;
957     }
958    
959     aMet.genVars(lGen);
960     aMet.recoVars(lReco);
961    
962     }//HbbMet
963    
964     void HbbTreeMaker::HbbTrigger(const edm::Handle<edm::TriggerResults> & aCol,
965     std::vector<HbbAnalysis::Trigger> & aVec)
966     {//HbbTrigger
967    
968     if (aCol.isValid()){
969     edm::TriggerNames lNames;
970     lNames.init(*aCol);
971    
972     //for (unsigned int j=0; j<hltPaths_.size(); j++) {
973     //bool valid=false;
974     for (unsigned int i=0; i<aCol->size(); i++){
975     std::string trueName = lNames.triggerNames().at(i);
976    
977     HbbAnalysis::TriggerVars lTrigVars;
978     lTrigVars.name = trueName;
979     lTrigVars.index = i;
980     lTrigVars.accept = aCol->accept(i);
981    
982     HbbAnalysis::Trigger lObj(lTrigVars);
983    
984     aVec.push_back(lObj);
985     }
986    
987     }
988    
989     }
990    
991     void HbbTreeMaker::HbbParticles(const edm::Handle<reco::GenParticleCollection> & aCol,
992     std::vector<HbbAnalysis::GenParticle> & aVec)
993     {//genparticles
994    
995     //add genParticle inside vector
996     for (unsigned int mccount = 0;mccount<aCol->size();mccount++)
997     {//loop on particles
998     const reco::Candidate & p = (*aCol)[mccount];
999    
1000     HbbAnalysis::MCVars lMC;
1001     lMC.index = mccount;
1002     lMC.pT = p.pt();
1003     lMC.eta = p.eta();
1004     lMC.phi = p.phi();
1005     lMC.pdgId = p.pdgId();
1006     lMC.status = p.status();
1007     HbbAnalysis::GenParticle lObj(lMC);
1008     aVec.push_back(lObj);
1009    
1010     }//loop on particles
1011    
1012     //now need to get the list of parents....
1013    
1014     for (unsigned int mccount = 0;mccount<aCol->size();mccount++)
1015     {//loop on particles
1016     const reco::Candidate & p = (*aCol)[mccount];
1017    
1018     unsigned int nD = p.numberOfDaughters();
1019     //std::cout << mccount << " has " << nD << " daughter(s) : " ;
1020     for(unsigned int dau=0;dau<nD;dau++){//loop on daughters
1021     const reco::Candidate * pDau = p.daughter(dau);
1022     //std::cout << pDau << " ("<<std::flush;
1023     for (unsigned int mccount2 = 0;mccount2<aCol->size();mccount2++)
1024     {//loop on particles
1025     const reco::Candidate & p2 = (*aCol)[mccount2];
1026     //std::cout<<"DBG: mccount2 = "<<mccount2<<" gen size = "<<data_->gen().size()<<std::endl;
1027    
1028     HbbAnalysis::GenParticle & part2 = aVec.at(mccount2);
1029     if (pDau == &p2) {
1030     part2.setParent(mccount);
1031     //std::cout << &p2 << ", index = " << mccount2 << "), "<<std::flush;
1032     break;
1033     }
1034     //else std::cout << std::endl << "****no match " << mccount2 << " " << &p2 << std::endl;
1035     //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){
1036     //part2.parent(mccount);
1037     //}
1038     }//loop on particles
1039     }//loop on daughters
1040     //std::cout << std::endl;
1041    
1042     }//loop on particles
1043    
1044     if (debug_){
1045     for (unsigned int imc(0); imc < aVec.size(); imc++){
1046     HbbAnalysis::GenParticle & p = aVec.at(imc);
1047     p.print();
1048     }
1049     }
1050    
1051    
1052     }//genparticles