ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Vuko/WZAnalysis/src/MuonAnalyzer.cc
Revision: 1.3
Committed: Fri Apr 4 08:51:22 2008 UTC (17 years ago) by senka
Content type: text/plain
Branch: MAIN
Changes since 1.2: +53 -1 lines
Log Message:
default values for track characteristics added

File Contents

# User Rev Content
1 senka 1.1 // -*- C++ -*-
2     //
3     // Package: WZAnalyzer
4     // Class: WZAnalyzer
5     //
6     /**\class WZAnalyzer WZAnalyzer.cc Vuko/WZAnalysis/src/WZAnalyzer.cc
7    
8     Description: <one line class summary>
9    
10     Implementation:
11     <Notes on implementation>
12     */
13     //
14     // Original Author: Vuko Brigljevic
15     // Created: Fri Nov 9 11:07:27 CET 2007
16 senka 1.2 // $Id: MuonAnalyzer.cc,v 1.1 2008/03/27 09:35:54 senka Exp $
17 senka 1.1 //
18     //
19    
20    
21    
22     // system include files
23     #include <memory>
24    
25     // user include files
26     #include "FWCore/Framework/interface/Frameworkfwd.h"
27     #include "FWCore/Framework/interface/EDAnalyzer.h"
28    
29     #include "FWCore/Framework/interface/Event.h"
30     #include "FWCore/Framework/interface/MakerMacros.h"
31    
32     #include "FWCore/ParameterSet/interface/ParameterSet.h"
33    
34     #include "Vuko/WZAnalysis/interface/MuonAnalyzer.h"
35     #include "DataFormats/Candidate/interface/OverlapChecker.h"
36    
37     #include "Vuko/WZAnalysis/interface/ElectronProperties.h"
38     #include "Vuko/WZAnalysis/interface/MuonProperties.h"
39    
40     #include "DataFormats/Common/interface/TriggerResults.h"
41    
42     //--- muon AOD:
43     #include "DataFormats/JetReco/interface/CaloJetCollection.h"
44     #include "DataFormats/EgammaCandidates/interface/Electron.h"
45     #include "DataFormats/EgammaCandidates/interface/ElectronFwd.h"
46     #include "DataFormats/MuonReco/interface/MuonFwd.h"
47     #include "DataFormats/MuonReco/interface/Muon.h"
48     #include "DataFormats/MuonReco/interface/MuIsoDeposit.h"
49     #include "DataFormats/EgammaReco/interface/BasicCluster.h"
50     #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
51     #include "DataFormats/EgammaCandidates/interface/PixelMatchGsfElectron.h"
52    
53     #include "DataFormats/METReco/interface/CaloMET.h"
54    
55    
56     #include "TFile.h"
57     #include "TH1F.h"
58     #include "TH2F.h"
59     #include "TTree.h"
60    
61     #include <map>
62    
63     //
64     // constants, enums and typedefs
65     //
66    
67     //
68     // static data member definitions
69     //
70    
71     //
72     // constructors and destructor
73     //
74     MuonAnalyzer::MuonAnalyzer(const edm::ParameterSet& iConfig)
75    
76     {
77    
78     //now do what ever initialization is needed
79    
80     fOutputFileName = iConfig.getUntrackedParameter<string>("HistOutFile");
81     theMCTruthCollection = iConfig.getParameter<edm::InputTag>("genParticles");
82    
83     // theMCTruthCollection = iConfig.getParameter<edm::InputTag>("genParticleCollection");
84    
85     theLooseMuonCollection = iConfig.getParameter<edm::InputTag>("LooseMuons");
86     theGoodMuonCollection = iConfig.getParameter<edm::InputTag>("GoodMuons");
87 senka 1.2 // theLooseElectronCollection = iConfig.getParameter<edm::InputTag>("LooseElectrons");
88     // theGoodElectronCollection = iConfig.getParameter<edm::InputTag>("GoodElectrons");
89     // theTightLeptonCollection = iConfig.getParameter<edm::InputTag>("TightLeptons");
90 senka 1.1 // Z's
91 senka 1.2 // theZeeCollection = iConfig.getParameter<edm::InputTag>("ZtoEE");
92     // theZmumuCollection = iConfig.getParameter<edm::InputTag>("ZtoMuMu");
93     // theZllCollection = iConfig.getParameter<edm::InputTag>("ZtoLL");
94 senka 1.1
95     theJetCollection = iConfig.getParameter<edm::InputTag>("Jets");
96    
97     theWeightCollection = iConfig.getParameter<edm::InputTag>("EventWeight");
98     alpgenIDTag = iConfig.getParameter<edm::InputTag>("AlpgenIDTag");
99    
100     getAlpgenID = false;
101     getAlpgenID = iConfig.getParameter<bool>("getAlpgenID");
102    
103     getEventWeight = false;
104     getEventWeight = iConfig.getParameter<bool>("getEventWeight");
105    
106     }
107    
108    
109     MuonAnalyzer::~MuonAnalyzer()
110     {
111    
112     // do anything here that needs to be done at desctruction time
113     // (e.g. close files, deallocate resources etc.)
114    
115     }
116    
117    
118     //
119     // member functions
120     //
121    
122     // ------------ method called to for each event ------------
123     void
124     MuonAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
125     {
126    
127     cout <<" ulazi u analyze metodu" << endl;
128     using namespace edm;
129     using namespace reco;
130     using namespace std;
131    
132     ///////////////////////////////////////
133     //
134     /// EVENT INITIALISATION
135     ///
136    
137     //read run number and event ID number
138     RunNumber = iEvent.id().run();
139     EventID = iEvent.id().event();
140    
141     ////////////////////////////////////////
142     // Get lists
143    
144     // loose muons
145     Handle<CandidateCollection> looseMuonCands;
146     iEvent.getByLabel( theLooseMuonCollection , looseMuonCands);
147    
148     // good muons
149     Handle<CandidateCollection> goodMuonCands;
150     iEvent.getByLabel( theGoodMuonCollection , goodMuonCands);
151    
152     // tight leptons
153 senka 1.2 // Handle<CandidateCollection> tightLeptonCands;
154     // iEvent.getByLabel( theTightLeptonCollection , tightLeptonCands);
155 senka 1.1
156     // jets
157    
158     Handle<CaloJetCollection> jetCands;
159     iEvent.getByLabel( theJetCollection , jetCands);
160    
161    
162     // event weights (CSA07 soups...)
163     eventWeight = -1;
164    
165     Handle< double> weightHandle;
166     Handle<int> genProcessID;
167     Handle<int> alpgenProcessID;
168     Handle<double> genEventScale;
169    
170     if (getEventWeight) {
171    
172     iEvent.getByLabel (theWeightCollection, weightHandle);
173     if (weightHandle.isValid()) {
174     eventWeight = * weightHandle;
175     }
176    
177     iEvent.getByLabel( "genEventProcID", genProcessID );
178     if (genProcessID.isValid()) {
179     processID = *genProcessID;
180     }
181    
182     if (processID==4 && getAlpgenID) {
183     iEvent.getByLabel( alpgenIDTag , alpgenProcessID );
184     if (alpgenProcessID.isValid()) {
185     alpgenID = *alpgenProcessID;
186     }
187     } else alpgenID=0;
188    
189    
190     iEvent.getByLabel( "genEventScale", genEventScale );
191     if (genEventScale.isValid()) {
192     eventScale = *genEventScale;
193     }
194     }
195    
196     // selected muons
197    
198     // cout << "\t # loose mu : " << looseMuonCands->size()
199     // << "\t # tight mu : " << goodMuonCands->size()
200     // << "\t # loose e : " << looseElectronCands->size()
201     // << "\t # tight e : " << goodElectronCands->size()
202     // << "\t # tight l : " << tightLeptonCands->size()
203     // << "\t # Z->mumu : " << zmumuCands->size()
204     // << "\t # Z->ee : " << zeeCands->size()
205     // << "\t # Z->ll : " << zllCands->size()
206     // << endl;
207    
208     N_looseMuons=0;
209     N_looseElectrons=0;
210     N_goodMuonCands=0;
211     N_looseMuons = looseMuonCands->size();
212     N_goodMuonCands=goodMuonCands->size();
213    
214     // Handle <reco::GenParticleCandidateCollection> mccands;
215     Handle<CandidateCollection> mccands;
216     iEvent.getByLabel( theMCTruthCollection,mccands );
217    
218     collectMCsummary(mccands);
219    
220    
221    
222     // matching:
223    
224     matching(iEvent,looseMuonCands, 13);
225    
226     getMother(&(*mccands)[1]);
227    
228     // dodatak za muon fake rate:
229     matching_JetsWithMuons(iEvent, looseMuonCands, jetCands);
230     matching_JetsWithB(iEvent, jetCands);
231     matching_JetsWithC(iEvent, jetCands);
232    
233     // for(CaloJetCollection::const_iterator jet = jetCands->begin();
234     // jet != jetCands->end(); jet++) {
235    
236     nb_jet=0;
237     nb_jet_eta_cut=0;
238     nb_jet_eta_and_pt_cut=0;
239     for (int i=0;i<jetCands->size();i++){
240     nb_jet++;
241     if (fabs((*jetCands)[i].eta())<2.5 ){
242     nb_jet_eta_cut++;
243     if ((*jetCands)[i].pt()<10){
244     nb_jet_eta_and_pt_cut++;
245     }
246     }
247     }
248    
249     nb_muon=0;
250     for (int i=0;i<looseMuonCands->size();i++){
251     if (fabs((*looseMuonCands)[i].eta())<2.5 ){
252     nb_muon++;
253     }
254     }
255    
256     for (int i=0;i<jetCands->size();i++){
257    
258     // cout <<" # jets ="<<jetCands->size() <<endl;
259     JetmatchedMuon=0;
260     JetmatchedB=0;
261     JetmatchedC=0;
262     MuonmatchedGenMuon=0;
263     dr_muon=0.;
264     dr_B=0.;
265     dr_C=0.;
266     // cout <<" dr na pocetku="<<dr << endl;
267     Jet_Et=-10.;
268     Jet_eta=-10.;
269     Jet_phi=-10.;
270     Jet_Pt=-10.;
271    
272     Jet_maxEInEmTowers=-10.;
273     Jet_maxEInHadTowers=-10.;
274     Jet_energyFractionHadronic=-10.;
275     Jet_emEnergyFraction=-10.;
276     Jet_hadEnergyInHB=-10.;
277     Jet_hadEnergyInHO=-10.;
278     Jet_hadEnergyInHE=-10.;
279     Jet_hadEnergyInHF=-10.;
280     Jet_emEnergyInEB=-10.;
281     Jet_emEnergyInEE=-10.;
282     Jet_emEnergyInHF=-10.;
283     Jet_towersArea=-10.;
284     Jet_n90=-10.;
285     Jet_n60=-10.;
286    
287     muon_Et=-10.;
288     muon_eta=-10.;
289     muon_Pt=-10.;
290     B_Et=-10.;
291     B_eta=-10.;
292     B_Pt=-10.;
293     C_Et=-10.;
294     C_eta=-10.;
295     C_Pt=-10.;
296     gen_muon_Et=-10.;
297     gen_muon_eta=-10.;
298     gen_muon_Pt=-10.;
299     muon_mother=0;
300     muon_mother_id=0;
301 senka 1.2 W_mother_id=0;
302 senka 1.1
303 senka 1.3
304     muon_global_chi2=-10.;
305     muon_global_p=-10.;
306     muon_global_pt=-10.;
307     muon_global_outerP=-10.;
308     muon_global_outerPt=-10.;
309     muon_global_ndof=-10.;
310     muon_global_normalizedChi2=-10.;
311     muon_global_recHitsSize=-10;
312     muon_global_numberOfLostHits=-10;
313     muon_global_numberOfValidHits=-10;
314     muon_global_innerPosition_x=-10.;
315     muon_global_innerPosition_y=-10.;
316     muon_global_innerPosition_z=-10.;
317     muon_global_outerPosition_x=-10.;
318     muon_global_outerPosition_y=-10.;
319     muon_global_outerPosition_z=-10.;
320    
321     muon_STA_chi2=-10.;
322     muon_STA_p=-10.;
323     muon_STA_pt=-10.;
324     muon_STA_outerP=-10.;
325     muon_STA_outerPt=-10.;
326     muon_STA_ndof=-10.;
327     muon_STA_normalizedChi2=-10.;
328     muon_STA_recHitsSize=-10;
329     muon_STA_numberOfLostHits=-10;
330     muon_STA_numberOfValidHits=-10;
331     muon_STA_innerPosition_x=-10.;
332     muon_STA_innerPosition_y=-10.;
333     muon_STA_innerPosition_z=-10.;
334     muon_STA_outerPosition_x=-10.;
335     muon_STA_outerPosition_y=-10.;
336     muon_STA_outerPosition_z=-10.;
337    
338     muon_track_chi2=-10.;
339     muon_track_p=-10.;
340     muon_track_pt=-10.;
341     muon_track_outerP=-10.;
342     muon_track_outerPt=-10.;
343     muon_track_ndof=-10.;
344     muon_track_normalizedChi2=-10.;
345     muon_track_recHitsSize=-10;
346     muon_track_numberOfLostHits=-10;
347     muon_track_numberOfValidHits=-10;
348     muon_track_innerPosition_x=-10.;
349     muon_track_innerPosition_y=-10.;
350     muon_track_innerPosition_z=-10.;
351     muon_track_outerPosition_x=-10.;
352     muon_track_outerPosition_y=-10.;
353     muon_track_outerPosition_z=-10.;
354    
355    
356 senka 1.1 // if (fabs((*jetCands)[i].eta())<2.5) {
357     // cout <<" jet passing eta cut" <<endl;
358     // if jet has matched muon and eta(muon)<2.5
359    
360 senka 1.2 const reco::Candidate * matchedmuon=MatchedMuonWithJet(&(*jetCands)[i]);
361    
362     if (MatchedMuonWithJet(&(*jetCands)[i])) {
363     cout <<" jet IMA matched muon !!!!!!!!!!!!!!" <<endl;
364     JetmatchedMuon=1;
365     dr_muon=dR_Muon_Jet(&(*jetCands)[i]);
366     cout <<" dr="<< dr_muon << " dR_Muon_Jet(&(*jetCands)[i])="<<dR_Muon_Jet(&(*jetCands)[i])<< endl;
367     // muon_Et=MatchedMuonWithJet(&(*jetCands)[i])->et();
368     // muon_eta=MatchedMuonWithJet(&(*jetCands)[i])->eta();
369     // muon_Pt=MatchedMuonWithJet(&(*jetCands)[i])->pt();
370     // muon_STA_chi2=MatchedMuonWithJet(&(*jetCands)[i])->standAloneMuon()->chi2();
371     muon_Et = matchedmuon->et();
372     muon_eta = matchedmuon->eta();
373     muon_Pt = matchedmuon->pt();
374    
375     const reco::Muon* matchedmuon_recoMuon = dynamic_cast<const reco::Muon *> (matchedmuon); // conversion of matchedmuon from reco::Candidate to reco::Muon
376    
377     muon_global_chi2 = matchedmuon_recoMuon->combinedMuon()->chi2();
378     muon_global_p = matchedmuon_recoMuon->combinedMuon()->p();
379     muon_global_pt = matchedmuon_recoMuon->combinedMuon()->pt();
380     muon_global_outerP = matchedmuon_recoMuon->combinedMuon()->outerP();
381     muon_global_outerPt = matchedmuon_recoMuon->combinedMuon()->outerPt();
382     muon_global_ndof = matchedmuon_recoMuon->combinedMuon()->ndof();
383     muon_global_normalizedChi2 = matchedmuon_recoMuon->combinedMuon()->normalizedChi2();
384     muon_global_recHitsSize = matchedmuon_recoMuon->combinedMuon()->recHitsSize();
385     muon_global_numberOfLostHits = matchedmuon_recoMuon->combinedMuon()->numberOfLostHits();
386     muon_global_numberOfValidHits = matchedmuon_recoMuon->combinedMuon()->numberOfValidHits();
387     muon_global_innerPosition_x = matchedmuon_recoMuon->combinedMuon()->innerPosition().x();
388     muon_global_innerPosition_y = matchedmuon_recoMuon->combinedMuon()->innerPosition().y();
389     muon_global_innerPosition_z = matchedmuon_recoMuon->combinedMuon()->innerPosition().z();
390     muon_global_outerPosition_x = matchedmuon_recoMuon->combinedMuon()->outerPosition().x();
391     muon_global_outerPosition_y = matchedmuon_recoMuon->combinedMuon()->outerPosition().y();
392     muon_global_outerPosition_z = matchedmuon_recoMuon->combinedMuon()->outerPosition().z();
393    
394     muon_STA_chi2 = matchedmuon_recoMuon->standAloneMuon()->chi2();
395     muon_STA_p = matchedmuon_recoMuon->standAloneMuon()->p();
396     muon_STA_pt = matchedmuon_recoMuon->standAloneMuon()->pt();
397     muon_STA_outerP = matchedmuon_recoMuon->standAloneMuon()->outerP();
398     muon_STA_outerPt = matchedmuon_recoMuon->standAloneMuon()->outerPt();
399     muon_STA_ndof = matchedmuon_recoMuon->standAloneMuon()->ndof();
400     muon_STA_normalizedChi2 = matchedmuon_recoMuon->standAloneMuon()->normalizedChi2();
401     muon_STA_recHitsSize = matchedmuon_recoMuon->standAloneMuon()->recHitsSize();
402     muon_STA_numberOfLostHits = matchedmuon_recoMuon->standAloneMuon()->numberOfLostHits();
403     muon_STA_numberOfValidHits = matchedmuon_recoMuon->standAloneMuon()->numberOfValidHits();
404     muon_STA_innerPosition_x = matchedmuon_recoMuon->standAloneMuon()->innerPosition().x();
405     muon_STA_innerPosition_y = matchedmuon_recoMuon->standAloneMuon()->innerPosition().y();
406     muon_STA_innerPosition_z = matchedmuon_recoMuon->standAloneMuon()->innerPosition().z();
407     muon_STA_outerPosition_x = matchedmuon_recoMuon->standAloneMuon()->outerPosition().x();
408     muon_STA_outerPosition_y = matchedmuon_recoMuon->standAloneMuon()->outerPosition().y();
409     muon_STA_outerPosition_z = matchedmuon_recoMuon->standAloneMuon()->outerPosition().z();
410    
411     muon_track_chi2 = matchedmuon_recoMuon->track()->chi2();
412     muon_track_p = matchedmuon_recoMuon->track()->p();
413     muon_track_pt = matchedmuon_recoMuon->track()->pt();
414     muon_track_outerP = matchedmuon_recoMuon->track()->outerP();
415     muon_track_outerPt = matchedmuon_recoMuon->track()->outerPt();
416     muon_track_ndof = matchedmuon_recoMuon->track()->ndof();
417     muon_track_normalizedChi2 = matchedmuon_recoMuon->track()->normalizedChi2();
418     muon_track_recHitsSize = matchedmuon_recoMuon->track()->recHitsSize();
419     muon_track_numberOfLostHits = matchedmuon_recoMuon->track()->numberOfLostHits();
420     muon_track_numberOfValidHits = matchedmuon_recoMuon->track()->numberOfValidHits();
421     muon_track_innerPosition_x = matchedmuon_recoMuon->track()->innerPosition().x();
422     muon_track_innerPosition_y = matchedmuon_recoMuon->track()->innerPosition().y();
423     muon_track_innerPosition_z = matchedmuon_recoMuon->track()->innerPosition().z();
424     muon_track_outerPosition_x = matchedmuon_recoMuon->track()->outerPosition().x();
425     muon_track_outerPosition_y = matchedmuon_recoMuon->track()->outerPosition().y();
426     muon_track_outerPosition_z = matchedmuon_recoMuon->track()->outerPosition().z();
427    
428     // if (fabs(MatchedMuonWithJet(&(*jetCands)[i])->eta())<2.5) {
429     // cout <<" muon prolazi eta cut" <<endl;
430    
431     for (map< int,
432 senka 1.1 std::vector< matchedParticles > > ::iterator iter = leptonMatching.begin();
433 senka 1.2 iter!=leptonMatching.end(); iter++) // entire map
434     {
435     // for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
436     // const reco::Candidate * nesto=iter->second[j].pair.second;
437    
438     if ((iter->first)==13)
439     cout <<"IDE PO CIJELOJ MAPI 13" << endl;
440     // }
441     }
442    
443     // ako muon (matched with jet) ima matched gen_muon
444 senka 1.1
445 senka 1.2 if (MatchedGenParticle(matchedmuon)){
446     cout <<" muon has matched gen muon !!!!!!!!!!" <<endl;
447     MuonmatchedGenMuon=1;
448     gen_muon_Et=MatchedGenParticle(matchedmuon)->et();
449     gen_muon_eta=MatchedGenParticle(matchedmuon)->eta();
450     gen_muon_Pt=MatchedGenParticle(matchedmuon)->pt();
451     MuonAnalyzer::LeptonOrigin origin = getMother(MatchedGenParticle(matchedmuon));
452    
453     if (origin == MuonAnalyzer::wboson ){
454     muon_mother=1;
455     cout << "MatchedGenFromW=true" <<endl;
456     W_mother_id = getMother_Id(getMother_particle(MatchedGenParticle(matchedmuon)));
457 senka 1.1 }
458 senka 1.2 else if (origin == MuonAnalyzer::zboson ){
459     muon_mother=2;
460     cout << "MatchedGenFromZ=true" <<endl;
461     }
462     else if (origin == MuonAnalyzer::bdecay ) {
463     muon_mother = 3;
464 senka 1.1 cout << "MatchedGenFromB=true" <<endl;
465 senka 1.2 }
466     else if (origin == MuonAnalyzer::cdecay ) {
467     muon_mother = 4;
468     cout << "MatchedGenFromC=true" <<endl;
469     }
470     else if (origin == MuonAnalyzer::tdecay ) {
471     muon_mother = 5;
472     cout << "MatchedGenFromT=true" <<endl;
473     }
474     muon_mother_id=getMother_Id(MatchedGenParticle(matchedmuon));
475    
476    
477 senka 1.1 // }
478     }
479 senka 1.2 // }
480     cout <<" jet nema matched muon" <<endl;
481     }
482    
483     if (MatchedBhadronWithJet(&(*jetCands)[i])) {
484 senka 1.1 cout <<" ima matched B hadron !!!!!!!!" << endl;
485     /////////////////
486    
487     JetmatchedB=1;
488     dr_B=dR_B_Jet(&(*jetCands)[i]);
489     cout <<" dr="<< dr_B << " dR_B_Jet(&(*jetCands)[i])="<<dR_B_Jet(&(*jetCands)[i])<< endl;
490     B_Et=MatchedBhadronWithJet(&(*jetCands)[i])->et();
491     B_eta=MatchedBhadronWithJet(&(*jetCands)[i])->eta();
492     B_Pt=MatchedBhadronWithJet(&(*jetCands)[i])->pt();
493    
494     }
495    
496     if (MatchedChadronWithJet(&(*jetCands)[i])) {
497     cout <<" ima matched C hadron !!!!!!!!" << endl;
498     /////////////////
499    
500     JetmatchedC=1;
501     dr_C=dR_C_Jet(&(*jetCands)[i]);
502     cout <<" dr="<< dr_C << " dR_C_Jet(&(*jetCands)[i])="<<dR_C_Jet(&(*jetCands)[i])<< endl;
503     C_Et=MatchedChadronWithJet(&(*jetCands)[i])->et();
504     C_eta=MatchedChadronWithJet(&(*jetCands)[i])->eta();
505     C_Pt=MatchedChadronWithJet(&(*jetCands)[i])->pt();
506    
507     }
508    
509     Jet_Et=(*jetCands)[i].et();
510     Jet_eta=(*jetCands)[i].eta();
511     Jet_phi=(*jetCands)[i].phi();
512     Jet_Pt=(*jetCands)[i].pt();
513    
514     Jet_maxEInEmTowers=(*jetCands)[i].maxEInEmTowers();
515     Jet_maxEInHadTowers=(*jetCands)[i].maxEInHadTowers();
516     Jet_energyFractionHadronic=(*jetCands)[i].energyFractionHadronic();
517     Jet_emEnergyFraction=(*jetCands)[i].emEnergyFraction();
518     Jet_hadEnergyInHB=(*jetCands)[i].hadEnergyInHB();
519     Jet_hadEnergyInHO=(*jetCands)[i].hadEnergyInHO();
520     Jet_hadEnergyInHE=(*jetCands)[i].hadEnergyInHE();
521     Jet_hadEnergyInHF=(*jetCands)[i].hadEnergyInHF();
522     Jet_emEnergyInEB=(*jetCands)[i].emEnergyInEB();
523     Jet_emEnergyInEE=(*jetCands)[i].emEnergyInEE();
524     Jet_emEnergyInHF=(*jetCands)[i].emEnergyInHF();
525     Jet_towersArea=(*jetCands)[i].towersArea();
526     Jet_n90=(*jetCands)[i].n90();
527     Jet_n60=(*jetCands)[i].n60();
528    
529     // cout <<" muon mother="<< muon_mother <<endl;
530    
531     muonTree->Fill();
532     }
533    
534     // for (int i=0;i<N_looseMuons; i++){
535     // if (MatchedGenParticle(&(*looseMuonCands)[i])==0) matchedMuon=0;
536     // else matchedMuon=1;
537     // // (&(*looseMuonCands)[0])
538     // muonTree->Fill();
539     // }
540    
541     JetMuTree->Fill();
542    
543    
544     }
545    
546    
547     // ------------ method called once each job just before starting event loop ------------
548     void
549     MuonAnalyzer::beginJob(const edm::EventSetup&)
550     {
551    
552     using namespace wzana;
553    
554     // theFile = new TFile( "wz.root", "RECREATE" ) ;
555     theFile = new TFile( fOutputFileName.c_str(), "RECREATE" ) ;
556    
557     muonTree = new TTree("JetMuon_Tree","muon informations for fake rate");
558     JetMuTree = new TTree("JetMuEvent_Tree","event info");
559    
560     initialiseTree();
561    
562     }
563    
564     // ------------ method called once each job just after ending the event loop ------------
565     void
566     MuonAnalyzer::endJob() {
567    
568     theFile->Write();
569     theFile->Close();
570    
571     }
572    
573    
574     void MuonAnalyzer::initialiseTree() {
575    
576     // muon properties for fake rate
577    
578 senka 1.2 muonTree->Branch("weight", &eventWeight, "weight/F");
579     muonTree->Branch("processID", &processID, "processID/I");
580     muonTree->Branch("alpgenID", &alpgenID, "alpgenID/I");
581     muonTree->Branch("genEventScale", &eventScale, "genEventScale/F");
582     muonTree->Branch("RunNumber", &RunNumber, "RunNumber/I");
583     muonTree->Branch("EventID", &EventID, "EventID/I");
584    
585     muonTree->Branch("N_looseMuons",&N_looseMuons,"N_looseMuons/I");
586     muonTree->Branch("N_goodMuonCands",&N_goodMuonCands,"N_goodMuonCands/I");
587     muonTree->Branch("matchedMuon",&matchedMuon,"matchedMuon/I");
588    
589     muonTree->Branch("nb_jet",&nb_jet,"nb_jet/I");
590     muonTree->Branch("nb_jet_eta_cut",&nb_jet_eta_cut,"nb_jet_eta_cut/I");
591     muonTree->Branch("nb_jet_eta_and_pt_cut",&nb_jet_eta_and_pt_cut,"nb_jet_eta_and_pt_cut/I");
592    
593     muonTree->Branch("nb_muon",&nb_muon,"nb_muon/I");
594 senka 1.1
595     // muonTree->Branch("NbJets",&,"/I");
596     muonTree->Branch("JetmatchedMuon",&JetmatchedMuon,"JetmatchedMuon/I");
597     muonTree->Branch("JetmatchedB",&JetmatchedB,"JetmatchedB/I");
598     muonTree->Branch("JetmatchedC",&JetmatchedC,"JetmatchedC/I");
599     muonTree->Branch("MuonmatchedGenMuon",&MuonmatchedGenMuon,"MuonmatchedGenMuon/I");
600    
601     muonTree->Branch("Jet_Et",&Jet_Et,"Jet_Et/F");
602     muonTree->Branch("Jet_eta",&Jet_eta,"Jet_eta/F");
603     muonTree->Branch("Jet_Pt",&Jet_Pt,"Jet_Pt/F");
604     muonTree->Branch("Jet_maxEInEmTowers",&Jet_maxEInEmTowers,"Jet_maxEInEmTowers/F");
605     muonTree->Branch("Jet_maxEInHadTowers",&Jet_maxEInHadTowers,"Jet_maxEInHadTowers/F");
606     muonTree->Branch("Jet_energyFractionHadronic",&Jet_energyFractionHadronic,"Jet_energyFractionHadronic/F");
607     muonTree->Branch("Jet_emEnergyFraction",&Jet_emEnergyFraction,"Jet_emEnergyFraction/F");
608     muonTree->Branch("Jet_hadEnergyInHB",&Jet_hadEnergyInHB,"Jet_hadEnergyInHB/F");
609     muonTree->Branch("Jet_hadEnergyInHO",&Jet_hadEnergyInHO,"Jet_hadEnergyInHO/F");
610     muonTree->Branch("Jet_hadEnergyInHE",&Jet_hadEnergyInHE,"Jet_hadEnergyInHE/F");
611     muonTree->Branch("Jet_hadEnergyInHF",&Jet_hadEnergyInHF,"Jet_hadEnergyInHF/F");
612     muonTree->Branch("Jet_emEnergyInEB",&Jet_emEnergyInEB,"Jet_emEnergyInEB/F");
613     muonTree->Branch("Jet_emEnergyInEE",&Jet_emEnergyInEE,"Jet_emEnergyInEE/F");
614     muonTree->Branch("Jet_emEnergyInHF",&Jet_emEnergyInHF,"Jet_emEnergyInHF/F");
615     muonTree->Branch("Jet_towersArea",&Jet_towersArea,"Jet_towersArea/F");
616     muonTree->Branch("Jet_n90",&Jet_n90,"Jet_n90/F");
617     muonTree->Branch("Jet_n60",&Jet_n60,"Jet_n60/F");
618    
619     muonTree->Branch("muon_Et",&muon_Et,"muon_Et/F");
620     muonTree->Branch("muon_eta",&muon_eta,"muon_eta/F");
621     muonTree->Branch("muon_Pt",&muon_Pt,"muon_Pt/F");
622    
623     muonTree->Branch("B_Et",&B_Et,"B_Et/F");
624     muonTree->Branch("B_eta",&B_eta,"B_eta/F");
625     muonTree->Branch("B_Pt",&B_Pt,"B_Pt/F");
626    
627     muonTree->Branch("C_Et",&C_Et,"C_Et/F");
628     muonTree->Branch("C_eta",&C_eta,"C_eta/F");
629     muonTree->Branch("C_Pt",&C_Pt,"C_Pt/F");
630    
631     muonTree->Branch("gen_muon_Et",&gen_muon_Et,"gen_muon_Et/F");
632     muonTree->Branch("gen_muon_eta",&gen_muon_eta,"gen_muon_eta/F");
633     muonTree->Branch("gen_muon_Pt",&gen_muon_Pt,"gen_muon_Pt/F");
634    
635     muonTree->Branch("dr_muon",&dr_muon,"dr_muon/F");
636     muonTree->Branch("dr_B",&dr_B,"dr_B/F");
637     muonTree->Branch("dr_C",&dr_C,"dr_C/F");
638    
639     muonTree->Branch("muon_mother",&muon_mother,"muon_mother/I");
640     muonTree->Branch("muon_mother_id",&muon_mother_id,"muon_mother_id/I");
641    
642 senka 1.2 muonTree->Branch("W_mother_id",&W_mother_id,"W_mother_id/I");
643    
644     muonTree->Branch("muon_global_chi2", &muon_global_chi2 , "muon_global_chi2/F");
645     muonTree->Branch("muon_global_p",&muon_global_p," muon_global_p/F");
646     muonTree->Branch("muon_global_pt",&muon_global_pt, "muon_global_pt/F");
647     muonTree->Branch("muon_global_outerP",& muon_global_outerP,"muon_global_outerP/F");
648     muonTree->Branch("muon_global_outerPt",& muon_global_outerPt,"muon_global_outerPt/F");
649     muonTree->Branch("muon_global_ndof",&muon_global_ndof,"muon_global_ndof/F");
650     muonTree->Branch("muon_global_normalizedChi2",& muon_global_normalizedChi2,"muon_global_normalizedChi2/F");
651     muonTree->Branch("muon_global_recHitsSize",&muon_global_recHitsSize,"muon_global_recHitsSize/I");
652     muonTree->Branch("muon_global_numberOfLostHits",& muon_global_numberOfLostHits,"muon_global_numberOfLostHits/I");
653     muonTree->Branch("muon_global_numberOfValidHits",&muon_global_numberOfValidHits,"muon_global_numberOfValidHits/I");
654     muonTree->Branch("muon_global_innerPosition_x",& muon_global_innerPosition_x,"muon_global_innerPosition_x/F");
655     muonTree->Branch("muon_global_innerPosition_y",& muon_global_innerPosition_y,"muon_global_innerPosition_y/F");
656     muonTree->Branch("muon_global_innerPosition_z",& muon_global_innerPosition_z,"muon_global_innerPosition_z/F");
657     muonTree->Branch("muon_global_outerPosition_x",& muon_global_outerPosition_x,"muon_global_outerPosition_x/F");
658     muonTree->Branch("muon_global_outerPosition_y",& muon_global_outerPosition_y,"muon_global_outerPosition_y/F");
659     muonTree->Branch("muon_global_outerPosition_z",& muon_global_outerPosition_z,"muon_global_outerPosition_z/F");
660    
661     muonTree->Branch("muon_STA_chi2", &muon_STA_chi2 , "muon_STA_chi2/F");
662     muonTree->Branch("muon_STA_p",&muon_STA_p," muon_STA_p/F");
663     muonTree->Branch("muon_STA_pt",&muon_STA_pt, "muon_STA_pt/F");
664     muonTree->Branch("muon_STA_outerP",& muon_STA_outerP,"muon_STA_outerP/F");
665     muonTree->Branch("muon_STA_outerPt",& muon_STA_outerPt,"muon_STA_outerPt/F");
666     muonTree->Branch("muon_STA_ndof",&muon_STA_ndof,"muon_STA_ndof/F");
667     muonTree->Branch("muon_STA_normalizedChi2",& muon_STA_normalizedChi2,"muon_STA_normalizedChi2/F");
668     muonTree->Branch("muon_STA_recHitsSize",&muon_STA_recHitsSize,"muon_STA_recHitsSize/I");
669     muonTree->Branch("muon_STA_numberOfLostHits",& muon_STA_numberOfLostHits,"muon_STA_numberOfLostHits/I");
670     muonTree->Branch("muon_STA_numberOfValidHits",&muon_STA_numberOfValidHits,"muon_STA_numberOfValidHits/I");
671     muonTree->Branch("muon_STA_innerPosition_x",& muon_STA_innerPosition_x,"muon_STA_innerPosition_x/F");
672     muonTree->Branch("muon_STA_innerPosition_y",& muon_STA_innerPosition_y,"muon_STA_innerPosition_y/F");
673     muonTree->Branch("muon_STA_innerPosition_z",& muon_STA_innerPosition_z,"muon_STA_innerPosition_z/F");
674     muonTree->Branch("muon_STA_outerPosition_x",& muon_STA_outerPosition_x,"muon_STA_outerPosition_x/F");
675     muonTree->Branch("muon_STA_outerPosition_y",& muon_STA_outerPosition_y,"muon_STA_outerPosition_y/F");
676     muonTree->Branch("muon_STA_outerPosition_z",& muon_STA_outerPosition_z,"muon_STA_outerPosition_z/F");
677    
678     muonTree->Branch("muon_track_chi2", &muon_track_chi2 , "muon_track_chi2/F");
679     muonTree->Branch("muon_track_p",&muon_track_p," muon_track_p/F");
680     muonTree->Branch("muon_track_pt",&muon_track_pt, "muon_track_pt/F");
681     muonTree->Branch("muon_track_outerP",& muon_track_outerP,"muon_track_outerP/F");
682     muonTree->Branch("muon_track_outerPt",& muon_track_outerPt,"muon_track_outerPt/F");
683     muonTree->Branch("muon_track_ndof",&muon_track_ndof,"muon_track_ndof/F");
684     muonTree->Branch("muon_track_normalizedChi2",& muon_track_normalizedChi2,"muon_track_normalizedChi2/F");
685     muonTree->Branch("muon_track_recHitsSize",&muon_track_recHitsSize,"muon_track_recHitsSize/I");
686     muonTree->Branch("muon_track_numberOfLostHits",& muon_track_numberOfLostHits,"muon_track_numberOfLostHits/I");
687     muonTree->Branch("muon_track_numberOfValidHits",&muon_track_numberOfValidHits,"muon_track_numberOfValidHits/I");
688     muonTree->Branch("muon_track_innerPosition_x",& muon_track_innerPosition_x,"muon_track_innerPosition_x/F");
689     muonTree->Branch("muon_track_innerPosition_y",& muon_track_innerPosition_y,"muon_track_innerPosition_y/F");
690     muonTree->Branch("muon_track_innerPosition_z",& muon_track_innerPosition_z,"muon_track_innerPosition_z/F");
691     muonTree->Branch("muon_track_outerPosition_x",& muon_track_outerPosition_x,"muon_track_outerPosition_x/F");
692     muonTree->Branch("muon_track_outerPosition_y",& muon_track_outerPosition_y,"muon_track_outerPosition_y/F");
693     muonTree->Branch("muon_track_outerPosition_z",& muon_track_outerPosition_z,"muon_track_outerPosition_z/F");
694    
695 senka 1.1 JetMuTree->Branch("nb_jet",&nb_jet,"nb_jet/I");
696     JetMuTree->Branch("nb_jet_eta_cut",&nb_jet_eta_cut,"nb_jet_eta_cut/I");
697     JetMuTree->Branch("nb_jet_eta_and_pt_cut",&nb_jet_eta_and_pt_cut,"nb_jet_eta_and_pt_cut/I");
698 senka 1.2
699 senka 1.1 JetMuTree->Branch("nb_muon",&nb_muon,"nb_muon/I");
700    
701 senka 1.2 JetMuTree->Branch("weight", &eventWeight, "weight/F");
702     JetMuTree->Branch("processID", &processID, "processID/I");
703     JetMuTree->Branch("alpgenID", &alpgenID, "alpgenID/I");
704     JetMuTree->Branch("genEventScale", &eventScale, "genEventScale/F");
705     JetMuTree->Branch("RunNumber", &RunNumber, "RunNumber/I");
706     JetMuTree->Branch("EventID", &EventID, "EventID/I");
707    
708 senka 1.1 }
709    
710    
711    
712     ////////////////////////
713     // matching
714     //
715    
716     class MatchingInfo {
717     public:
718    
719     MatchingInfo(float dr, int gen, int reco) : deltaR(dr), genMuon(gen),recoMuon(reco) {}
720    
721     bool operator<(const MatchingInfo m) { return deltaR < m.deltaR;}
722    
723     float deltaR;
724     int genMuon;
725     int recoMuon;
726    
727     };
728    
729    
730     bool betterMatch2(MatchingInfo m1, MatchingInfo m2)
731     { return m1.deltaR < m2.deltaR;}
732    
733    
734     // method that gives the vector<MatchingInfo> matching_Infos containing (dR,i,j); i=genParticle index, j=ParticleCollection index:
735    
736     void MuonAnalyzer::matching(const edm::Event& iEvent, Handle<reco::CandidateCollection> cands, int pdgid){
737    
738     using namespace edm;
739     using namespace reco;
740     using namespace std;
741    
742     //-------------------------
743     Handle<CandidateCollection> mccands;
744     iEvent.getByLabel( theMCTruthCollection,mccands );
745    
746     vector <GenParticleCandidate*> genparticles;
747    
748     for(CandidateCollection::const_iterator p = mccands->begin();
749     p != mccands->end(); p++) {
750    
751     // reco::Candidate* p_tmp = &(*p);
752     reco::GenParticleCandidate* ptr = (reco::GenParticleCandidate*) &(*p);
753    
754     if ( (ptr)->status() == 1
755     /*&& (ptr)->momentum().perp() > 2.*/ ) { // stable particles with pt>2 only
756     if ( abs((ptr)->pdgId()) == pdgid ) {
757     genparticles.push_back((ptr));
758     //cout << "electron MCT\n";
759     }
760     }
761     }
762    
763    
764     int n1=0;
765     vector<MatchingInfo> matching_Infos;
766     for ( unsigned int i=0; i<genparticles.size(); i++ ) // po svim MC muonima
767     {
768    
769     for (unsigned int j = 0; j < cands->size(); j++)// po svim reco muonima
770     {
771    
772     double d_eta2=(*cands)[j].eta()-(genparticles)[i]->eta();
773     double d_phi=fabs((*cands)[j].phi()-(genparticles)[i]->phi());
774     double dR_byhand= sqrt(d_eta2*d_eta2+d_phi*d_phi);
775     double dR=dR_byhand;
776    
777     if (dR<0.15){
778     n1++;
779     matching_Infos.push_back(MatchingInfo(dR,i,j));
780     }
781    
782     }
783     }
784    
785     sort(matching_Infos.begin(),matching_Infos.end(),betterMatch2); // sorting (dR,i,j)
786    
787     if (genparticles.size()!=0 && cands->size()!=0){
788     // Now exclude combinations using same gen or rec muons
789     for (std::vector<MatchingInfo>::iterator it=matching_Infos.begin();
790     it != matching_Infos.end(); it++) {
791     for (std::vector<MatchingInfo>::iterator it2=matching_Infos.begin();
792     it2!=it; it2++) {
793     if ( (it->genMuon == it2->genMuon) || (it->recoMuon == it2->recoMuon) ) {
794     matching_Infos.erase(it);
795     it=matching_Infos.begin();
796     break;
797     }
798     }
799     }
800     }
801    
802     // Now put result in vector of pairs
803    
804     leptonMatching[pdgid].clear();
805    
806     for (vector<MatchingInfo>::iterator match=matching_Infos.begin();
807     match !=matching_Infos.end(); match++) { // po tablici dr,i,j
808    
809     pair<const GenParticleCandidate*, const reco::Candidate *> pair;
810     // vector<pair,float > matchedParticles;
811     pair.first = genparticles[match->genMuon];
812     pair.second = & (*cands)[match->recoMuon];
813    
814     matchedParticles cestice(genparticles[match->genMuon] ,& (*cands)[match->recoMuon] , match->deltaR);
815    
816     leptonMatching[pdgid].push_back(cestice);
817    
818     }
819    
820     }
821    
822 senka 1.2
823 senka 1.1 /////////////////////////////////////////////////////////////////////
824     // made for muon fake rate
825     // method that gives the vector<MatchingInfo> matching_Infos containing (dR,i,j); i=genParticle index, j=ParticleCollection index:
826    
827     void MuonAnalyzer::matching_JetsWithMuons(const edm::Event& iEvent, Handle<reco::CandidateCollection> cands, Handle<reco::CaloJetCollection> jets){
828    
829     using namespace edm;
830     using namespace reco;
831     using namespace std;
832    
833     cout <<" ulazi u matching_JetsWithMuons" <<endl;
834    
835     int n1=0;
836     vector<MatchingInfo> matching_Infos;
837     cout <<"# jets="<<jets->size()<< " # muons="<<cands->size() << endl;
838     for ( unsigned int i=0; i<jets->size(); i++ ) // po svim jetovima
839     {
840     // cout <<" next jet" <<endl;
841     for (unsigned int j = 0; j < cands->size(); j++)// po svim reco muonima
842     {
843     // cout <<" next muon" <<endl;
844     double d_eta2=(*cands)[j].eta()-(*jets)[i].eta();
845     double d_phi=fabs((*cands)[j].phi()-(*jets)[i].phi());
846     double dR_byhand= sqrt(d_eta2*d_eta2+d_phi*d_phi);
847     double dR=dR_byhand;
848     // cout <<"dR="<< dR <<endl;
849     // if (dR<0.15){
850     n1++;
851     matching_Infos.push_back(MatchingInfo(dR,i,j));
852     // }
853    
854     }
855     }
856    
857     sort(matching_Infos.begin(),matching_Infos.end(),betterMatch2); // sorting (dR,i,j)
858    
859     if (jets->size()!=0 && cands->size()!=0){
860     // Now exclude combinations using same gen or rec muons
861     for (std::vector<MatchingInfo>::iterator it=matching_Infos.begin();
862     it != matching_Infos.end(); it++) {
863     for (std::vector<MatchingInfo>::iterator it2=matching_Infos.begin();
864     it2!=it; it2++) {
865     if ( (it->genMuon == it2->genMuon) || (it->recoMuon == it2->recoMuon) ) {
866     matching_Infos.erase(it);
867     it=matching_Infos.begin();
868     break;
869     }
870     }
871     }
872     }
873    
874     // Now t in vector of pairs
875    
876     Muon_Jet_Matching[1].clear();
877    
878     for (vector<MatchingInfo>::iterator match=matching_Infos.begin();
879     match !=matching_Infos.end(); match++) { // po tablici dr,i,j
880    
881     pair<const CaloJet*, const reco::Candidate *> pair;
882     // vector<pair,float > matchedParticles;
883     pair.first = &(*jets)[match->genMuon];
884     pair.second = & (*cands)[match->recoMuon];
885    
886     matchedJets_and_Muons cestice(&(*jets)[match->genMuon] ,& (*cands)[match->recoMuon] , match->deltaR);
887    
888     Muon_Jet_Matching[1].push_back(cestice);
889    
890     }
891    
892     cout <<" nakon pisanja u Muon_Jet_Matching[1]" << endl;
893     for (vector<MatchingInfo>::iterator match=matching_Infos.begin();
894     match !=matching_Infos.end(); match++) { // po tablici dr,i,j
895     // cout <<" nakon pisanja u Muon_Jet_Matching[1]" << endl;
896     for (map< int,
897     std::vector< matchedJets_and_Muons > > ::iterator iter = Muon_Jet_Matching.begin();
898     iter!=Muon_Jet_Matching.end(); iter++) // entire map
899     {
900     for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
901     const reco::Candidate * nesto=iter->second[j].pair.second;
902     // ::OverlapChecker overlap;
903     // if (overlap(*nesto,*daughter)){
904     cout <<"*****dr="<< iter->second[j].delta_R<<endl;
905     // cout <<"*****dr="<< iter->second[j].delta_R<<endl;
906     // }
907     }
908    
909    
910    
911     }
912    
913    
914    
915     }
916    
917     }
918    
919     /////////////////////////////////////////////////////////////////////
920     // made for muon fake rate
921     // method that gives the vector<MatchingInfo> matching_Infos containing (dR,i,j); i=genParticle index, j=ParticleCollection index:
922    
923     void MuonAnalyzer::matching_JetsWithB(const edm::Event& iEvent, Handle<reco::CaloJetCollection> jets){
924    
925     cout <<" ulazi u matching_JetsWithB" <<endl;
926    
927     using namespace edm;
928     using namespace reco;
929     using namespace std;
930    
931     Handle<CandidateCollection> mccands;
932     iEvent.getByLabel( theMCTruthCollection,mccands );
933    
934     vector <GenParticleCandidate*> genBparticles; // vector u kojem su svi B hadroni cija majka nije B hadron (b kvark)
935     // vector <GenParticleCandidate*> genCparticles; // vector u kojem su svi D hadroni cija majka nije D niti B hadron (c kvark)
936    
937    
938     for(CandidateCollection::const_iterator p = mccands->begin();
939     p != mccands->end(); p++) {
940    
941     reco::GenParticleCandidate* ptr = (reco::GenParticleCandidate*) &(*p);
942    
943     ///////////////////////////////////////////////////////
944     if ((((abs(ptr->pdgId()))/100)%10 ==5)||(((abs(ptr->pdgId()))/1000)%10 ==5 )) { // b
945     if ((((abs(ptr->mother()->pdgId()))/100)%10 !=5)&&(((abs(ptr->mother()->pdgId()))/1000)%10 !=5 )) genBparticles.push_back(ptr); // mother is not b
946     }
947    
948     ///////////////////////////////////////////////////////
949    
950    
951     // int pdg_id=ptr->pdgId();
952     // // cout << "particle id : " << pdg_id << endl;
953    
954     // while ( ptr ->mother()!=0 ) { //izvrsava se sve dok ne dodje do nule tj. dok ne trazi majku od protona
955     // // cout << "Going up: ";
956     // // cout << "Mother id : " << genParticle->pdgId() << endl;
957     // ptr= (reco::GenParticleCandidate*) ptr->mother();
958    
959     // if (abs(ptr->pdgId())!=abs(pdg_id)) {
960     // // cout<< " good mother " <<endl;
961    
962     // if ((((abs(ptr->pdgId()))/100)%10 ==5)||(((abs(ptr->pdgId()))/1000)%10 ==5 )) { // mother is b
963     // // while (ptr = const_cast<reco::GenParticleCandidate *>(dynamic_cast<const reco::GenParticleCandidate *>(ptr->mother()))) {
964     // while (ptr ->mother()!=0) {
965    
966     // // ptr= (reco::GenParticleCandidate*) ptr->mother();
967     // // if ((((abs(ptr->pdgId()))/100)%10 !=5)||(((abs(ptr->pdgId()))/1000)%10 !=5 )){ // ako majka nije b
968     // if ((((abs(ptr ->mother()->pdgId()))/100)%10 !=5)||(((abs(ptr ->mother()->pdgId()))/1000)%10 !=5 )){ // ako majka nije b
969     // genBparticles.push_back(ptr);
970     // break;
971     // }
972     // ptr= (reco::GenParticleCandidate*) ptr->mother();
973     // }
974    
975     // }
976    
977     // }
978    
979     // }
980    
981     }
982    
983     int n1=0;
984     vector<MatchingInfo> matching_Infos;
985     cout <<"# jets="<<jets->size()<< " # B hadrons="<<genBparticles.size() << endl;
986     for ( unsigned int i=0; i<genBparticles.size(); i++ ) // po svim MC B hadronima
987     {
988     // cout <<" next B hadron" <<endl;
989     for (unsigned int j = 0; j < jets->size(); j++)// po svim jetovima
990     {
991     // cout <<" next jet" <<endl;
992     double d_eta2=(*jets)[j].eta()-(genBparticles)[i]->eta();
993     double d_phi=fabs((*jets)[j].phi()-(genBparticles)[i]->phi());
994     double dR_byhand= sqrt(d_eta2*d_eta2+d_phi*d_phi);
995     double dR=dR_byhand;
996     // cout <<"dR="<< dR <<endl;
997     // if (dR<0.15){
998     n1++;
999     matching_Infos.push_back(MatchingInfo(dR,i,j));
1000     // }
1001    
1002     }
1003     }
1004    
1005     sort(matching_Infos.begin(),matching_Infos.end(),betterMatch2); // sorting (dR,i,j)
1006    
1007     if (jets->size()!=0 && genBparticles.size()!=0){
1008     // Now exclude combinations using same gen or rec muons
1009     for (std::vector<MatchingInfo>::iterator it=matching_Infos.begin();
1010     it != matching_Infos.end(); it++) {
1011     for (std::vector<MatchingInfo>::iterator it2=matching_Infos.begin();
1012     it2!=it; it2++) {
1013     if ( (it->genMuon == it2->genMuon) || (it->recoMuon == it2->recoMuon) ) {
1014     matching_Infos.erase(it);
1015     it=matching_Infos.begin();
1016     break;
1017     }
1018     }
1019     }
1020     }
1021    
1022     // Now t in vector of pairs
1023    
1024     B_Jet_Matching[1].clear();
1025    
1026     for (vector<MatchingInfo>::iterator match=matching_Infos.begin();
1027     match !=matching_Infos.end(); match++) { // po tablici dr,i,j
1028    
1029     pair<const GenParticleCandidate*, const CaloJet*> pair;
1030     // vector<pair,float > matchedParticles;
1031     pair.first = genBparticles[match->genMuon];
1032     pair.second = &(*jets)[match->recoMuon];
1033    
1034     matchedJets_and_B cestice(genBparticles[match->genMuon],&(*jets)[match->recoMuon] , match->deltaR);
1035    
1036     B_Jet_Matching[1].push_back(cestice);
1037    
1038     }
1039    
1040     cout <<" nakon pisanja u B_Jet_Matching[1]" << endl;
1041     for (vector<MatchingInfo>::iterator match=matching_Infos.begin();
1042     match !=matching_Infos.end(); match++) { // po tablici dr,i,j
1043     // int j=0;
1044     for (map< int,
1045     std::vector< matchedJets_and_B > > ::iterator iter = B_Jet_Matching.begin();
1046     iter!=B_Jet_Matching.end(); iter++) // entire map
1047     {
1048     for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
1049     const reco::Candidate * nesto=iter->second[j].pair.second;
1050     // ::OverlapChecker overlap;
1051     // if (overlap(*nesto,*daughter)){
1052     cout <<"*****dr="<< iter->second[j].delta_R<<endl;
1053     // cout <<"*****dr="<< iter->second[j].delta_R<<endl;
1054     // }
1055     }
1056    
1057     }
1058    
1059     }
1060    
1061     }
1062    
1063     /////////////////////////////////////////////////////////////////////
1064     // made for muon fake rate
1065     // method that gives the vector<MatchingInfo> matching_Infos containing (dR,i,j); i=genParticle index, j=ParticleCollection index:
1066    
1067     void MuonAnalyzer::matching_JetsWithC(const edm::Event& iEvent, Handle<reco::CaloJetCollection> jets){
1068    
1069     cout <<" ulazi u matching_JetsWithC" <<endl;
1070    
1071     using namespace edm;
1072     using namespace reco;
1073     using namespace std;
1074    
1075     Handle<CandidateCollection> mccands;
1076     iEvent.getByLabel( theMCTruthCollection,mccands );
1077    
1078     vector <GenParticleCandidate*> genCparticles; // vector u kojem su svi D hadroni cija majka nije D niti B hadron (c kvark)
1079    
1080     for(CandidateCollection::const_iterator p = mccands->begin();
1081     p != mccands->end(); p++) {
1082    
1083     reco::GenParticleCandidate* ptr = (reco::GenParticleCandidate*) &(*p);
1084    
1085     if ((((abs(ptr->pdgId()))/100)%10 ==4)||(((abs(ptr->pdgId()))/1000)%10 ==4 )) { // c
1086     if ((((abs(ptr->mother()->pdgId()))/100)%10 !=4)&&(((abs(ptr->mother()->pdgId()))/1000)%10 !=4 )
1087     && (((abs(ptr->mother()->pdgId()))/100)%10 !=5)&&(((abs(ptr->mother()->pdgId()))/1000)%10 !=5 ))
1088     genCparticles.push_back(ptr); // mother is not b or c
1089     }
1090    
1091     // int pdg_id=ptr->pdgId();
1092     // // cout << "particle id : " << pdg_id << endl;
1093    
1094     // while (ptr->mother()!=0) { //izvrsava se sve dok ne dodje do nule tj. dok ne trazi majku od protona
1095     // // cout << "Going up: ";
1096     // // cout << "Mother id : " << genParticle->pdgId() << endl;
1097     // ptr= (reco::GenParticleCandidate*) ptr->mother();
1098     // if (abs(ptr->pdgId())!=abs(pdg_id)) {
1099     // // cout<< " good mother " <<endl;
1100    
1101    
1102     // if ((((abs(ptr->pdgId()))/100)%10 ==4)||(((abs(ptr->pdgId()))/1000)%10 ==4 )) { // mother is c
1103     // while (ptr->mother()!=0) {
1104     // // ptr= (reco::GenParticleCandidate*) ptr->mother();
1105     // if (((((abs(ptr->mother()->pdgId()))/100)%10 !=4)||(((abs(ptr->mother()->pdgId()))/1000)%10 !=4 )) &&
1106     // ((((abs(ptr->mother()->pdgId()))/100)%10 !=5)||(((abs(ptr->mother()->pdgId()))/1000)%10 !=5 ))){ // ako majka nije c niti b
1107     // genCparticles.push_back(ptr);
1108     // break;
1109     // }
1110     // ptr= (reco::GenParticleCandidate*) ptr->mother();
1111     // }
1112     // }
1113    
1114    
1115     // }
1116     // }
1117    
1118     }
1119    
1120     int n1=0;
1121     vector<MatchingInfo> matching_Infos;
1122     cout <<"# jets="<<jets->size()<< " # C hadrons="<<genCparticles.size() << endl;
1123     for ( unsigned int i=0; i<genCparticles.size(); i++ ) // po svim MC B hadronima
1124     {
1125     // cout <<" next B hadron" <<endl;
1126     for (unsigned int j = 0; j < jets->size(); j++)// po svim jetovima
1127     {
1128     // cout <<" next jet" <<endl;
1129     double d_eta2=(*jets)[j].eta()-(genCparticles)[i]->eta();
1130     double d_phi=fabs((*jets)[j].phi()-(genCparticles)[i]->phi());
1131     double dR_byhand= sqrt(d_eta2*d_eta2+d_phi*d_phi);
1132     double dR=dR_byhand;
1133     // cout <<"dR="<< dR <<endl;
1134     // if (dR<0.15){
1135     n1++;
1136     matching_Infos.push_back(MatchingInfo(dR,i,j));
1137     // }
1138    
1139     }
1140     }
1141    
1142     sort(matching_Infos.begin(),matching_Infos.end(),betterMatch2); // sorting (dR,i,j)
1143    
1144     if (jets->size()!=0 && genCparticles.size()!=0){
1145     // Now exclude combinations using same gen or rec muons
1146     for (std::vector<MatchingInfo>::iterator it=matching_Infos.begin();
1147     it != matching_Infos.end(); it++) {
1148     for (std::vector<MatchingInfo>::iterator it2=matching_Infos.begin();
1149     it2!=it; it2++) {
1150     if ( (it->genMuon == it2->genMuon) || (it->recoMuon == it2->recoMuon) ) {
1151     matching_Infos.erase(it);
1152     it=matching_Infos.begin();
1153     break;
1154     }
1155     }
1156     }
1157     }
1158    
1159     // Now t in vector of pairs
1160    
1161     C_Jet_Matching[1].clear();
1162    
1163     for (vector<MatchingInfo>::iterator match=matching_Infos.begin();
1164     match !=matching_Infos.end(); match++) { // po tablici dr,i,j
1165    
1166     pair<const GenParticleCandidate*, const CaloJet*> pair;
1167     // vector<pair,float > matchedParticles;
1168     pair.first = genCparticles[match->genMuon];
1169     pair.second = &(*jets)[match->recoMuon];
1170    
1171     matchedJets_and_C cestice(genCparticles[match->genMuon],&(*jets)[match->recoMuon] , match->deltaR);
1172    
1173     C_Jet_Matching[1].push_back(cestice);
1174    
1175     }
1176    
1177     cout <<" nakon pisanja u C_Jet_Matching[1]" << endl;
1178     for (vector<MatchingInfo>::iterator match=matching_Infos.begin();
1179     match !=matching_Infos.end(); match++) { // po tablici dr,i,j
1180     for (map< int,
1181     std::vector< matchedJets_and_B > > ::iterator iter = B_Jet_Matching.begin();
1182     iter!=B_Jet_Matching.end(); iter++) // entire map
1183     {
1184     for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
1185     const reco::Candidate * nesto=iter->second[j].pair.second;
1186     // ::OverlapChecker overlap;
1187     // if (overlap(*nesto,*daughter)){
1188     cout <<"*****dr="<< iter->second[j].delta_R<<endl;
1189     // cout <<"*****dr="<< iter->second[j].delta_R<<endl;
1190     // }
1191     }
1192    
1193     }
1194    
1195    
1196    
1197     }
1198    
1199     }
1200    
1201    
1202     //---------------------------------------------------
1203    
1204    
1205     // get mother of particle
1206     // returns mother = 1 if mother is Z boson
1207     // returns mother = 2 if mother is W boson
1208     // returns mother = 4 if mother is b
1209     // returns mother = 0 else
1210    
1211     MuonAnalyzer::LeptonOrigin MuonAnalyzer::getMother(const reco::Candidate* genParticle){
1212    
1213     int pdg_id=genParticle->pdgId();
1214     // cout << "particle id : " << pdg_id << endl;
1215    
1216     while (genParticle = genParticle->mother() ) { //izvrsava se sve dok ne dodje do nule tj. dok ne trazi majku od protona
1217     // cout << "Going up: ";
1218     // cout << "Mother id : " << genParticle->pdgId() << endl;
1219     if (abs(genParticle->pdgId())!=abs(pdg_id)) {
1220     // cout<< " good mother " <<endl;
1221     if (abs(genParticle->pdgId())==23){ // mother is Z
1222     return zboson;
1223     }
1224     if (abs(genParticle->pdgId())==24) { // mother is W
1225     return wboson;
1226     }
1227     // if (abs(genParticle->pdgId())==23 || abs(genParticle->pdgId())==24) { // mother is W or Z
1228     // WZ_matched=1;
1229     // mother=3;
1230     // }
1231     if ((((abs(genParticle->pdgId()))/100)%10 ==5)||(((abs(genParticle->pdgId()))/1000)%10 ==5 )) { // mother is b
1232     return bdecay;
1233     }
1234     if ((((abs(genParticle->pdgId()))/100)%10 ==4)||(((abs(genParticle->pdgId()))/1000)%10 ==4 )) { // mother is c
1235     return cdecay;
1236     }
1237 senka 1.2 if (abs(genParticle->pdgId()==6)) { // mother is t
1238     return tdecay;
1239     }
1240 senka 1.1 return other;
1241     }
1242     }
1243    
1244     return other;
1245     }
1246    
1247 senka 1.2 const reco::Candidate * MuonAnalyzer::getMother_particle(const reco::Candidate* genParticle){
1248    
1249     int pdg_id=genParticle->pdgId();
1250     // cout << "particle id : " << pdg_id << endl;
1251    
1252     while (genParticle = genParticle->mother() ) { //izvrsava se sve dok ne dodje do nule tj. dok ne trazi majku od protona
1253     // cout << "Going up: ";
1254     // cout << "Mother id : " << genParticle->pdgId() << endl;
1255     if (abs(genParticle->pdgId())!=abs(pdg_id)) {
1256     // cout<< " good mother " <<endl;
1257     return genParticle->mother();
1258     }
1259     }
1260    
1261     // return ;
1262     }
1263    
1264    
1265    
1266 senka 1.1 ///////////////////
1267     int MuonAnalyzer::getMother_Id(const reco::Candidate* genParticle){
1268    
1269     int mother_id=0;
1270    
1271     int pdg_id=genParticle->pdgId();
1272     // cout << "particle id : " << pdg_id << endl;
1273    
1274     while (genParticle = genParticle->mother() ) { //izvrsava se sve dok ne dodje do nule tj. dok ne trazi majku od protona
1275     if (abs(genParticle->pdgId())!=abs(pdg_id)) {
1276     mother_id=genParticle->pdgId();
1277     return mother_id;
1278     }
1279     }
1280    
1281     return mother_id;
1282     }
1283    
1284     /////////////////
1285    
1286    
1287     const reco::Candidate * MuonAnalyzer::MatchedGenParticle(const reco::Candidate * daughter){
1288     ::OverlapChecker overlap;
1289     matched_genParticle=0;
1290    
1291     for (map< int,
1292     std::vector< matchedParticles > > ::iterator iter = leptonMatching.begin();
1293     iter!=leptonMatching.end(); iter++) // entire map
1294     {
1295     for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
1296     const reco::Candidate * nesto=iter->second[j].pair.second;
1297     if (overlap(*nesto,*daughter)){
1298     matched_genParticle=iter->second[j].pair.first;
1299     }
1300     }
1301     }
1302     return matched_genParticle;
1303     }
1304    
1305     /////////////////////////////////////////
1306     // for muon fake rate
1307    
1308     const reco::Candidate * MuonAnalyzer::MatchedMuonWithJet(const reco::CaloJet * daughter){
1309     ::OverlapChecker overlap;
1310     matched_Muon=0;
1311    
1312     for (map< int,
1313     std::vector< matchedJets_and_Muons > > ::iterator iter = Muon_Jet_Matching.begin();
1314     iter!=Muon_Jet_Matching.end(); iter++) // entire map
1315     {
1316     // cout <<" ide preko cijele mape za jet" <<endl;
1317     for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
1318     const reco::CaloJet * nesto=iter->second[j].pair.first;
1319     if (overlap(*nesto,*daughter)){
1320     matched_Muon=iter->second[j].pair.second; /// tu je bila greska
1321     }
1322     }
1323     }
1324     return matched_Muon;
1325     }
1326    
1327    
1328     /////////////////////////////
1329     // Jet and B hadron matching
1330    
1331     const reco::Candidate * MuonAnalyzer::MatchedBhadronWithJet(const reco::CaloJet * daughter){
1332     ::OverlapChecker overlap;
1333     matched_B=0;
1334    
1335     for (map< int,
1336     std::vector< matchedJets_and_B > > ::iterator iter = B_Jet_Matching.begin();
1337     iter!=B_Jet_Matching.end(); iter++) // entire map
1338     {
1339     // cout <<" ide preko cijele mape za jet" <<endl;
1340     for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
1341     const reco::CaloJet * nesto=iter->second[j].pair.second;
1342     if (overlap(*nesto,*daughter)){
1343     matched_B=iter->second[j].pair.first; /// tu je bila greska
1344     }
1345     }
1346     }
1347     return matched_B;
1348     }
1349    
1350     /////////////////////////////
1351     // Jet and C hadron matching
1352    
1353     const reco::Candidate * MuonAnalyzer::MatchedChadronWithJet(const reco::CaloJet * daughter){
1354     ::OverlapChecker overlap;
1355     matched_C=0;
1356    
1357     for (map< int,
1358     std::vector< matchedJets_and_C > > ::iterator iter = C_Jet_Matching.begin();
1359     iter!=C_Jet_Matching.end(); iter++) // entire map
1360     {
1361     // cout <<" ide preko cijele mape za jet" <<endl;
1362     for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
1363     const reco::CaloJet * nesto=iter->second[j].pair.second;
1364     if (overlap(*nesto,*daughter)){
1365     matched_C=iter->second[j].pair.first; /// tu je bila greska
1366     }
1367     }
1368     }
1369     return matched_C;
1370     }
1371    
1372    
1373     float MuonAnalyzer::dR(const reco::Candidate * daughter){
1374    
1375     ::OverlapChecker overlap;
1376    
1377     float delta_r = -10;
1378     for (map< int,
1379     std::vector< matchedParticles > > ::iterator iter = leptonMatching.begin();
1380     iter!=leptonMatching.end(); iter++) // entire map
1381     {
1382     for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
1383     const reco::Candidate * nesto=iter->second[j].pair.second;
1384     if (overlap(*nesto,*daughter)){
1385     delta_r=iter->second[j].delta_R;
1386     }
1387     }
1388     }
1389     return delta_r;
1390     }
1391    
1392     float MuonAnalyzer::dR_Muon_Jet(const reco::CaloJet * daughter){
1393    
1394     ::OverlapChecker overlap;
1395    
1396     float delta_r = -10;
1397     for (map< int,
1398     std::vector<matchedJets_and_Muons > > ::iterator iter = Muon_Jet_Matching.begin();
1399     iter!=Muon_Jet_Matching.end(); iter++) // entire map
1400     {
1401     for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
1402     const reco::CaloJet * nesto=iter->second[j].pair.first;
1403     if (overlap(*nesto,*daughter)){
1404     delta_r=iter->second[j].delta_R;
1405     }
1406     }
1407     }
1408     return delta_r;
1409     }
1410    
1411     float MuonAnalyzer::dR_B_Jet(const reco::CaloJet * daughter){
1412    
1413     ::OverlapChecker overlap;
1414    
1415     float delta_r = -10;
1416     for (map< int,
1417     std::vector<matchedJets_and_B > > ::iterator iter = B_Jet_Matching.begin();
1418     iter!=B_Jet_Matching.end(); iter++) // entire map
1419     {
1420     for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
1421     const reco::CaloJet * nesto=iter->second[j].pair.second;
1422     if (overlap(*nesto,*daughter)){
1423     delta_r=iter->second[j].delta_R;
1424     }
1425     }
1426     }
1427     return delta_r;
1428     }
1429    
1430     float MuonAnalyzer::dR_C_Jet(const reco::CaloJet * daughter){
1431    
1432     ::OverlapChecker overlap;
1433    
1434     float delta_r = -10;
1435     for (map< int,
1436     std::vector<matchedJets_and_C > > ::iterator iter = C_Jet_Matching.begin();
1437     iter!=C_Jet_Matching.end(); iter++) // entire map
1438     {
1439     for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
1440     const reco::CaloJet * nesto=iter->second[j].pair.second;
1441     if (overlap(*nesto,*daughter)){
1442     delta_r=iter->second[j].delta_R;
1443     }
1444     }
1445     }
1446     return delta_r;
1447     }
1448    
1449    
1450    
1451     void MuonAnalyzer::collectMCsummary(Handle <reco::CandidateCollection> mccands) {
1452    
1453     using namespace reco;
1454    
1455     //collections
1456    
1457     MCleptZ1_pdgid = -1;
1458     MCleptZ2_pdgid = -1;
1459     MCleptW_pdgid = -1;
1460    
1461     MCleptZ1_pt = -1;
1462     MCleptZ2_pt = -1;
1463     MCleptW_pt = -1;
1464    
1465     MCleptZ1_eta = -1;
1466     MCleptZ2_eta = -1;
1467     MCleptW_eta = -1;
1468    
1469     MCleptZ1_phi = -1;
1470     MCleptZ2_phi = -1;
1471     MCleptW_phi = -1;
1472    
1473     vector<reco::GenParticleCandidate*> Tau;
1474     vector<reco::GenParticleCandidate*> StableMuons;
1475     vector<reco::GenParticleCandidate*> StableElectrons;
1476    
1477    
1478    
1479     for(CandidateCollection::const_iterator p = mccands->begin();
1480     p != mccands->end(); p++) {
1481    
1482     // reco::Candidate* p_tmp = &(*p);
1483     reco::GenParticleCandidate* ptr = (reco::GenParticleCandidate*) &(*p);
1484    
1485     if ( (ptr)->status() == 1
1486     /*&& (ptr)->momentum().perp() > 2.*/ ) { // stable particles with pt>2 only
1487     if ( abs((ptr)->pdgId()) == 11 ) {
1488    
1489     StableElectrons.push_back((ptr));
1490     //cout << "electron MCT\n";
1491     }
1492     else if ( abs((ptr)->pdgId()) == 13 ) {
1493    
1494     StableMuons.push_back((ptr)) ;
1495     //cout << "muon MCT\n";
1496     }
1497     }
1498     else if ((ptr)->status() == 2) {
1499     if ( abs((ptr)->pdgId()) == 15 ) {//tau
1500     Tau.push_back((ptr));
1501     }
1502     }
1503     }
1504    
1505     // std::cout << "# Electrons : " << StableElectrons.size()
1506     // << "# Muons : " << StableMuons.size() << std::endl
1507     // << "# Tau : " << Tau.size() << std::endl;
1508    
1509    
1510     vector<reco::GenParticleCandidate*> * StableLeptons[3] = {&StableElectrons, &StableMuons, &Tau};
1511    
1512     vector<reco::GenParticleCandidate*> TauChildLeptons;
1513    
1514    
1515     //erase electrons and muons from tau from GenparticleCandidate Vector and put to Tau Vector
1516     for (int i=1; i>=0; i--) {
1517     for ( vector<reco::GenParticleCandidate*>::iterator lepton = StableLeptons[i]->begin();
1518     lepton != StableLeptons[i]->end();lepton++) {
1519    
1520     reco::GenParticleCandidate* mcParticleRef = *lepton;
1521    
1522     int parentItr=0;
1523     while ( mcParticleRef->mother()!=0 && parentItr<2) {
1524    
1525     parentItr++;
1526    
1527     mcParticleRef = (reco::GenParticleCandidate*) mcParticleRef->mother();
1528     if (mcParticleRef==0) break;
1529    
1530     //if tau
1531     if (abs((mcParticleRef)->pdgId())==15 ) {
1532     //put into collection
1533     TauChildLeptons.push_back(*lepton);
1534     StableLeptons[i]->erase(lepton);
1535     lepton--;
1536     break;
1537     }
1538     }
1539     }
1540     }
1541    
1542    
1543     bool firstZlept = true;
1544     int MC_tauDecayTypeZ1 = 0;
1545     int MC_tauDecayTypeZ2 = 0;
1546     int MC_tauDecayTypeW = 0;
1547    
1548    
1549     for (int i=2; i>=0; i--) {
1550     while (StableLeptons[i]->size() > 0) {
1551     float maxPt = 0;
1552     vector<reco::GenParticleCandidate*>::iterator index = StableLeptons[i]->end();
1553    
1554     for ( vector<reco::GenParticleCandidate*>::iterator lepton = StableLeptons[i]->begin();
1555     lepton != StableLeptons[i]->end(); lepton++ ) {
1556    
1557     if ((*lepton)->pt() > maxPt) {
1558     maxPt = (*lepton)->pt();
1559     index = lepton;
1560     }
1561     }
1562     bool Zchild = false;
1563     bool Wchild = false;
1564     bool isTau = false;
1565    
1566    
1567     reco::GenParticleCandidate* mcParticleRef = *index;
1568    
1569     if (abs((*index)->pdgId()) == 15) isTau = true;
1570    
1571     //find W/Z mother
1572     int parentItr=0;
1573     while ( mcParticleRef->mother()!=0 && parentItr<2) {
1574    
1575     if (parentItr>=2) break;
1576     parentItr++;
1577     mcParticleRef = (reco::GenParticleCandidate*) mcParticleRef->mother();
1578     if (mcParticleRef==0) break;
1579    
1580     if (abs((mcParticleRef)->pdgId())==23 ) {
1581     Zchild=true;
1582     break;
1583     }
1584     if (abs((mcParticleRef)->pdgId())==24 ) {
1585     Wchild=true;
1586     break;
1587     }
1588     }
1589    
1590    
1591     if (maxPt >0) {
1592     int * MCtauDecayTypePtr = 0;
1593    
1594     if (Wchild) {
1595     MCleptW_pdgid=(*index)->pdgId();
1596     MCleptW_pt=(float)(*index)->pt();
1597     MCleptW_eta=(float)(*index)->eta();
1598     MCleptW_phi=(float)(*index)->phi();
1599     if (isTau) {
1600     MCtauDecayTypePtr = &MC_tauDecayTypeW;
1601     MC_tauDecayTypeW =3;//default to hadronic decay
1602     }
1603    
1604     }
1605     if (Zchild) {
1606     if (firstZlept) {
1607     firstZlept=false;
1608     MCleptZ1_pdgid=(*index)->pdgId();
1609     MCleptZ1_pt=(float)(*index)->pt();
1610     MCleptZ1_eta=(float)(*index)->eta();
1611     MCleptZ1_phi=(float)(*index)->phi();
1612     if (isTau) {
1613     MCtauDecayTypePtr = &MC_tauDecayTypeZ1;
1614     MC_tauDecayTypeZ1 =3;
1615     }
1616    
1617     }
1618     else {
1619     MCleptZ2_pdgid=(*index)->pdgId();
1620     MCleptZ2_pt=(float)(*index)->pt();
1621     MCleptZ2_eta=(float)(*index)->eta();
1622     MCleptZ2_phi=(float)(*index)->phi();
1623     if (isTau) {
1624     MCtauDecayTypePtr = &MC_tauDecayTypeZ2;
1625     MC_tauDecayTypeZ2 =3;
1626     }
1627    
1628     }
1629     }
1630     //check type of tau decay
1631     if (MCtauDecayTypePtr) {
1632     for( vector<reco::GenParticleCandidate*>::iterator p = TauChildLeptons.begin();p != TauChildLeptons.end(); p++) {
1633     int pitr = 0;
1634     reco::GenParticleCandidate* mcParticleRef = *p;
1635     while (mcParticleRef->mother() && pitr<2) {
1636     pitr++;
1637     mcParticleRef =(reco::GenParticleCandidate*) mcParticleRef->mother();
1638     if (mcParticleRef==0) break;
1639     if (mcParticleRef == *index) {
1640     if (abs((*p)->pdgId()) == 11) *MCtauDecayTypePtr = 1;
1641     if (abs((*p)->pdgId()) == 13) *MCtauDecayTypePtr = 2;
1642     }
1643     }
1644     }
1645     }
1646     }
1647     StableLeptons[i]->erase(index);
1648     }
1649     }
1650     }
1651     //define this as a plug-in
1652     // DEFINE_FWK_MODULE(WZAnalyzer);
1653