ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Vuko/WZAnalysis/src/MuonAnalyzer.cc
Revision: 1.2
Committed: Thu Apr 3 14:39:59 2008 UTC (17 years, 1 month ago) by senka
Content type: text/plain
Branch: MAIN
Changes since 1.1: +229 -747 lines
Log Message:
infos about hits and global, STA, tracker tracks 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 senka 1.2 muon_STA_chi2=-10.;
291 senka 1.1 B_Et=-10.;
292     B_eta=-10.;
293     B_Pt=-10.;
294     C_Et=-10.;
295     C_eta=-10.;
296     C_Pt=-10.;
297     gen_muon_Et=-10.;
298     gen_muon_eta=-10.;
299     gen_muon_Pt=-10.;
300     muon_mother=0;
301     muon_mother_id=0;
302 senka 1.2 W_mother_id=0;
303 senka 1.1
304     // if (fabs((*jetCands)[i].eta())<2.5) {
305     // cout <<" jet passing eta cut" <<endl;
306     // if jet has matched muon and eta(muon)<2.5
307    
308 senka 1.2 const reco::Candidate * matchedmuon=MatchedMuonWithJet(&(*jetCands)[i]);
309    
310     if (MatchedMuonWithJet(&(*jetCands)[i])) {
311     cout <<" jet IMA matched muon !!!!!!!!!!!!!!" <<endl;
312     JetmatchedMuon=1;
313     dr_muon=dR_Muon_Jet(&(*jetCands)[i]);
314     cout <<" dr="<< dr_muon << " dR_Muon_Jet(&(*jetCands)[i])="<<dR_Muon_Jet(&(*jetCands)[i])<< endl;
315     // muon_Et=MatchedMuonWithJet(&(*jetCands)[i])->et();
316     // muon_eta=MatchedMuonWithJet(&(*jetCands)[i])->eta();
317     // muon_Pt=MatchedMuonWithJet(&(*jetCands)[i])->pt();
318     // muon_STA_chi2=MatchedMuonWithJet(&(*jetCands)[i])->standAloneMuon()->chi2();
319     muon_Et = matchedmuon->et();
320     muon_eta = matchedmuon->eta();
321     muon_Pt = matchedmuon->pt();
322    
323     const reco::Muon* matchedmuon_recoMuon = dynamic_cast<const reco::Muon *> (matchedmuon); // conversion of matchedmuon from reco::Candidate to reco::Muon
324    
325     muon_global_chi2 = matchedmuon_recoMuon->combinedMuon()->chi2();
326     muon_global_p = matchedmuon_recoMuon->combinedMuon()->p();
327     muon_global_pt = matchedmuon_recoMuon->combinedMuon()->pt();
328     muon_global_outerP = matchedmuon_recoMuon->combinedMuon()->outerP();
329     muon_global_outerPt = matchedmuon_recoMuon->combinedMuon()->outerPt();
330     muon_global_ndof = matchedmuon_recoMuon->combinedMuon()->ndof();
331     muon_global_normalizedChi2 = matchedmuon_recoMuon->combinedMuon()->normalizedChi2();
332     muon_global_recHitsSize = matchedmuon_recoMuon->combinedMuon()->recHitsSize();
333     muon_global_numberOfLostHits = matchedmuon_recoMuon->combinedMuon()->numberOfLostHits();
334     muon_global_numberOfValidHits = matchedmuon_recoMuon->combinedMuon()->numberOfValidHits();
335     muon_global_innerPosition_x = matchedmuon_recoMuon->combinedMuon()->innerPosition().x();
336     muon_global_innerPosition_y = matchedmuon_recoMuon->combinedMuon()->innerPosition().y();
337     muon_global_innerPosition_z = matchedmuon_recoMuon->combinedMuon()->innerPosition().z();
338     muon_global_outerPosition_x = matchedmuon_recoMuon->combinedMuon()->outerPosition().x();
339     muon_global_outerPosition_y = matchedmuon_recoMuon->combinedMuon()->outerPosition().y();
340     muon_global_outerPosition_z = matchedmuon_recoMuon->combinedMuon()->outerPosition().z();
341    
342     muon_STA_chi2 = matchedmuon_recoMuon->standAloneMuon()->chi2();
343     muon_STA_p = matchedmuon_recoMuon->standAloneMuon()->p();
344     muon_STA_pt = matchedmuon_recoMuon->standAloneMuon()->pt();
345     muon_STA_outerP = matchedmuon_recoMuon->standAloneMuon()->outerP();
346     muon_STA_outerPt = matchedmuon_recoMuon->standAloneMuon()->outerPt();
347     muon_STA_ndof = matchedmuon_recoMuon->standAloneMuon()->ndof();
348     muon_STA_normalizedChi2 = matchedmuon_recoMuon->standAloneMuon()->normalizedChi2();
349     muon_STA_recHitsSize = matchedmuon_recoMuon->standAloneMuon()->recHitsSize();
350     muon_STA_numberOfLostHits = matchedmuon_recoMuon->standAloneMuon()->numberOfLostHits();
351     muon_STA_numberOfValidHits = matchedmuon_recoMuon->standAloneMuon()->numberOfValidHits();
352     muon_STA_innerPosition_x = matchedmuon_recoMuon->standAloneMuon()->innerPosition().x();
353     muon_STA_innerPosition_y = matchedmuon_recoMuon->standAloneMuon()->innerPosition().y();
354     muon_STA_innerPosition_z = matchedmuon_recoMuon->standAloneMuon()->innerPosition().z();
355     muon_STA_outerPosition_x = matchedmuon_recoMuon->standAloneMuon()->outerPosition().x();
356     muon_STA_outerPosition_y = matchedmuon_recoMuon->standAloneMuon()->outerPosition().y();
357     muon_STA_outerPosition_z = matchedmuon_recoMuon->standAloneMuon()->outerPosition().z();
358    
359     muon_track_chi2 = matchedmuon_recoMuon->track()->chi2();
360     muon_track_p = matchedmuon_recoMuon->track()->p();
361     muon_track_pt = matchedmuon_recoMuon->track()->pt();
362     muon_track_outerP = matchedmuon_recoMuon->track()->outerP();
363     muon_track_outerPt = matchedmuon_recoMuon->track()->outerPt();
364     muon_track_ndof = matchedmuon_recoMuon->track()->ndof();
365     muon_track_normalizedChi2 = matchedmuon_recoMuon->track()->normalizedChi2();
366     muon_track_recHitsSize = matchedmuon_recoMuon->track()->recHitsSize();
367     muon_track_numberOfLostHits = matchedmuon_recoMuon->track()->numberOfLostHits();
368     muon_track_numberOfValidHits = matchedmuon_recoMuon->track()->numberOfValidHits();
369     muon_track_innerPosition_x = matchedmuon_recoMuon->track()->innerPosition().x();
370     muon_track_innerPosition_y = matchedmuon_recoMuon->track()->innerPosition().y();
371     muon_track_innerPosition_z = matchedmuon_recoMuon->track()->innerPosition().z();
372     muon_track_outerPosition_x = matchedmuon_recoMuon->track()->outerPosition().x();
373     muon_track_outerPosition_y = matchedmuon_recoMuon->track()->outerPosition().y();
374     muon_track_outerPosition_z = matchedmuon_recoMuon->track()->outerPosition().z();
375    
376     // if (fabs(MatchedMuonWithJet(&(*jetCands)[i])->eta())<2.5) {
377     // cout <<" muon prolazi eta cut" <<endl;
378    
379     for (map< int,
380 senka 1.1 std::vector< matchedParticles > > ::iterator iter = leptonMatching.begin();
381 senka 1.2 iter!=leptonMatching.end(); iter++) // entire map
382     {
383     // for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
384     // const reco::Candidate * nesto=iter->second[j].pair.second;
385    
386     if ((iter->first)==13)
387     cout <<"IDE PO CIJELOJ MAPI 13" << endl;
388     // }
389     }
390    
391     // ako muon (matched with jet) ima matched gen_muon
392 senka 1.1
393 senka 1.2 if (MatchedGenParticle(matchedmuon)){
394     cout <<" muon has matched gen muon !!!!!!!!!!" <<endl;
395     MuonmatchedGenMuon=1;
396     gen_muon_Et=MatchedGenParticle(matchedmuon)->et();
397     gen_muon_eta=MatchedGenParticle(matchedmuon)->eta();
398     gen_muon_Pt=MatchedGenParticle(matchedmuon)->pt();
399     MuonAnalyzer::LeptonOrigin origin = getMother(MatchedGenParticle(matchedmuon));
400    
401     if (origin == MuonAnalyzer::wboson ){
402     muon_mother=1;
403     cout << "MatchedGenFromW=true" <<endl;
404     W_mother_id = getMother_Id(getMother_particle(MatchedGenParticle(matchedmuon)));
405 senka 1.1 }
406 senka 1.2 else if (origin == MuonAnalyzer::zboson ){
407     muon_mother=2;
408     cout << "MatchedGenFromZ=true" <<endl;
409     }
410     else if (origin == MuonAnalyzer::bdecay ) {
411     muon_mother = 3;
412 senka 1.1 cout << "MatchedGenFromB=true" <<endl;
413 senka 1.2 }
414     else if (origin == MuonAnalyzer::cdecay ) {
415     muon_mother = 4;
416     cout << "MatchedGenFromC=true" <<endl;
417     }
418     else if (origin == MuonAnalyzer::tdecay ) {
419     muon_mother = 5;
420     cout << "MatchedGenFromT=true" <<endl;
421     }
422     muon_mother_id=getMother_Id(MatchedGenParticle(matchedmuon));
423    
424    
425 senka 1.1 // }
426     }
427 senka 1.2 // }
428     cout <<" jet nema matched muon" <<endl;
429     }
430    
431     if (MatchedBhadronWithJet(&(*jetCands)[i])) {
432 senka 1.1 cout <<" ima matched B hadron !!!!!!!!" << endl;
433     /////////////////
434    
435     JetmatchedB=1;
436     dr_B=dR_B_Jet(&(*jetCands)[i]);
437     cout <<" dr="<< dr_B << " dR_B_Jet(&(*jetCands)[i])="<<dR_B_Jet(&(*jetCands)[i])<< endl;
438     B_Et=MatchedBhadronWithJet(&(*jetCands)[i])->et();
439     B_eta=MatchedBhadronWithJet(&(*jetCands)[i])->eta();
440     B_Pt=MatchedBhadronWithJet(&(*jetCands)[i])->pt();
441    
442     }
443    
444     if (MatchedChadronWithJet(&(*jetCands)[i])) {
445     cout <<" ima matched C hadron !!!!!!!!" << endl;
446     /////////////////
447    
448     JetmatchedC=1;
449     dr_C=dR_C_Jet(&(*jetCands)[i]);
450     cout <<" dr="<< dr_C << " dR_C_Jet(&(*jetCands)[i])="<<dR_C_Jet(&(*jetCands)[i])<< endl;
451     C_Et=MatchedChadronWithJet(&(*jetCands)[i])->et();
452     C_eta=MatchedChadronWithJet(&(*jetCands)[i])->eta();
453     C_Pt=MatchedChadronWithJet(&(*jetCands)[i])->pt();
454    
455     }
456    
457     Jet_Et=(*jetCands)[i].et();
458     Jet_eta=(*jetCands)[i].eta();
459     Jet_phi=(*jetCands)[i].phi();
460     Jet_Pt=(*jetCands)[i].pt();
461    
462     Jet_maxEInEmTowers=(*jetCands)[i].maxEInEmTowers();
463     Jet_maxEInHadTowers=(*jetCands)[i].maxEInHadTowers();
464     Jet_energyFractionHadronic=(*jetCands)[i].energyFractionHadronic();
465     Jet_emEnergyFraction=(*jetCands)[i].emEnergyFraction();
466     Jet_hadEnergyInHB=(*jetCands)[i].hadEnergyInHB();
467     Jet_hadEnergyInHO=(*jetCands)[i].hadEnergyInHO();
468     Jet_hadEnergyInHE=(*jetCands)[i].hadEnergyInHE();
469     Jet_hadEnergyInHF=(*jetCands)[i].hadEnergyInHF();
470     Jet_emEnergyInEB=(*jetCands)[i].emEnergyInEB();
471     Jet_emEnergyInEE=(*jetCands)[i].emEnergyInEE();
472     Jet_emEnergyInHF=(*jetCands)[i].emEnergyInHF();
473     Jet_towersArea=(*jetCands)[i].towersArea();
474     Jet_n90=(*jetCands)[i].n90();
475     Jet_n60=(*jetCands)[i].n60();
476    
477     // cout <<" muon mother="<< muon_mother <<endl;
478    
479     muonTree->Fill();
480     }
481    
482     // for (int i=0;i<N_looseMuons; i++){
483     // if (MatchedGenParticle(&(*looseMuonCands)[i])==0) matchedMuon=0;
484     // else matchedMuon=1;
485     // // (&(*looseMuonCands)[0])
486     // muonTree->Fill();
487     // }
488    
489     JetMuTree->Fill();
490    
491    
492     }
493    
494    
495     // ------------ method called once each job just before starting event loop ------------
496     void
497     MuonAnalyzer::beginJob(const edm::EventSetup&)
498     {
499    
500     using namespace wzana;
501    
502     // theFile = new TFile( "wz.root", "RECREATE" ) ;
503     theFile = new TFile( fOutputFileName.c_str(), "RECREATE" ) ;
504    
505     muonTree = new TTree("JetMuon_Tree","muon informations for fake rate");
506     JetMuTree = new TTree("JetMuEvent_Tree","event info");
507    
508     initialiseTree();
509    
510     }
511    
512     // ------------ method called once each job just after ending the event loop ------------
513     void
514     MuonAnalyzer::endJob() {
515    
516     theFile->Write();
517     theFile->Close();
518    
519     }
520    
521    
522     void MuonAnalyzer::initialiseTree() {
523    
524     // muon properties for fake rate
525    
526 senka 1.2 muonTree->Branch("weight", &eventWeight, "weight/F");
527     muonTree->Branch("processID", &processID, "processID/I");
528     muonTree->Branch("alpgenID", &alpgenID, "alpgenID/I");
529     muonTree->Branch("genEventScale", &eventScale, "genEventScale/F");
530     muonTree->Branch("RunNumber", &RunNumber, "RunNumber/I");
531     muonTree->Branch("EventID", &EventID, "EventID/I");
532    
533     muonTree->Branch("N_looseMuons",&N_looseMuons,"N_looseMuons/I");
534     muonTree->Branch("N_goodMuonCands",&N_goodMuonCands,"N_goodMuonCands/I");
535     muonTree->Branch("matchedMuon",&matchedMuon,"matchedMuon/I");
536    
537     muonTree->Branch("nb_jet",&nb_jet,"nb_jet/I");
538     muonTree->Branch("nb_jet_eta_cut",&nb_jet_eta_cut,"nb_jet_eta_cut/I");
539     muonTree->Branch("nb_jet_eta_and_pt_cut",&nb_jet_eta_and_pt_cut,"nb_jet_eta_and_pt_cut/I");
540    
541     muonTree->Branch("nb_muon",&nb_muon,"nb_muon/I");
542 senka 1.1
543     // muonTree->Branch("NbJets",&,"/I");
544     muonTree->Branch("JetmatchedMuon",&JetmatchedMuon,"JetmatchedMuon/I");
545     muonTree->Branch("JetmatchedB",&JetmatchedB,"JetmatchedB/I");
546     muonTree->Branch("JetmatchedC",&JetmatchedC,"JetmatchedC/I");
547     muonTree->Branch("MuonmatchedGenMuon",&MuonmatchedGenMuon,"MuonmatchedGenMuon/I");
548    
549     muonTree->Branch("Jet_Et",&Jet_Et,"Jet_Et/F");
550     muonTree->Branch("Jet_eta",&Jet_eta,"Jet_eta/F");
551     muonTree->Branch("Jet_Pt",&Jet_Pt,"Jet_Pt/F");
552     muonTree->Branch("Jet_maxEInEmTowers",&Jet_maxEInEmTowers,"Jet_maxEInEmTowers/F");
553     muonTree->Branch("Jet_maxEInHadTowers",&Jet_maxEInHadTowers,"Jet_maxEInHadTowers/F");
554     muonTree->Branch("Jet_energyFractionHadronic",&Jet_energyFractionHadronic,"Jet_energyFractionHadronic/F");
555     muonTree->Branch("Jet_emEnergyFraction",&Jet_emEnergyFraction,"Jet_emEnergyFraction/F");
556     muonTree->Branch("Jet_hadEnergyInHB",&Jet_hadEnergyInHB,"Jet_hadEnergyInHB/F");
557     muonTree->Branch("Jet_hadEnergyInHO",&Jet_hadEnergyInHO,"Jet_hadEnergyInHO/F");
558     muonTree->Branch("Jet_hadEnergyInHE",&Jet_hadEnergyInHE,"Jet_hadEnergyInHE/F");
559     muonTree->Branch("Jet_hadEnergyInHF",&Jet_hadEnergyInHF,"Jet_hadEnergyInHF/F");
560     muonTree->Branch("Jet_emEnergyInEB",&Jet_emEnergyInEB,"Jet_emEnergyInEB/F");
561     muonTree->Branch("Jet_emEnergyInEE",&Jet_emEnergyInEE,"Jet_emEnergyInEE/F");
562     muonTree->Branch("Jet_emEnergyInHF",&Jet_emEnergyInHF,"Jet_emEnergyInHF/F");
563     muonTree->Branch("Jet_towersArea",&Jet_towersArea,"Jet_towersArea/F");
564     muonTree->Branch("Jet_n90",&Jet_n90,"Jet_n90/F");
565     muonTree->Branch("Jet_n60",&Jet_n60,"Jet_n60/F");
566    
567     muonTree->Branch("muon_Et",&muon_Et,"muon_Et/F");
568     muonTree->Branch("muon_eta",&muon_eta,"muon_eta/F");
569     muonTree->Branch("muon_Pt",&muon_Pt,"muon_Pt/F");
570    
571     muonTree->Branch("B_Et",&B_Et,"B_Et/F");
572     muonTree->Branch("B_eta",&B_eta,"B_eta/F");
573     muonTree->Branch("B_Pt",&B_Pt,"B_Pt/F");
574    
575     muonTree->Branch("C_Et",&C_Et,"C_Et/F");
576     muonTree->Branch("C_eta",&C_eta,"C_eta/F");
577     muonTree->Branch("C_Pt",&C_Pt,"C_Pt/F");
578    
579     muonTree->Branch("gen_muon_Et",&gen_muon_Et,"gen_muon_Et/F");
580     muonTree->Branch("gen_muon_eta",&gen_muon_eta,"gen_muon_eta/F");
581     muonTree->Branch("gen_muon_Pt",&gen_muon_Pt,"gen_muon_Pt/F");
582    
583     muonTree->Branch("dr_muon",&dr_muon,"dr_muon/F");
584     muonTree->Branch("dr_B",&dr_B,"dr_B/F");
585     muonTree->Branch("dr_C",&dr_C,"dr_C/F");
586    
587     muonTree->Branch("muon_mother",&muon_mother,"muon_mother/I");
588     muonTree->Branch("muon_mother_id",&muon_mother_id,"muon_mother_id/I");
589    
590 senka 1.2 muonTree->Branch("W_mother_id",&W_mother_id,"W_mother_id/I");
591    
592     muonTree->Branch("muon_global_chi2", &muon_global_chi2 , "muon_global_chi2/F");
593     muonTree->Branch("muon_global_p",&muon_global_p," muon_global_p/F");
594     muonTree->Branch("muon_global_pt",&muon_global_pt, "muon_global_pt/F");
595     muonTree->Branch("muon_global_outerP",& muon_global_outerP,"muon_global_outerP/F");
596     muonTree->Branch("muon_global_outerPt",& muon_global_outerPt,"muon_global_outerPt/F");
597     muonTree->Branch("muon_global_ndof",&muon_global_ndof,"muon_global_ndof/F");
598     muonTree->Branch("muon_global_normalizedChi2",& muon_global_normalizedChi2,"muon_global_normalizedChi2/F");
599     muonTree->Branch("muon_global_recHitsSize",&muon_global_recHitsSize,"muon_global_recHitsSize/I");
600     muonTree->Branch("muon_global_numberOfLostHits",& muon_global_numberOfLostHits,"muon_global_numberOfLostHits/I");
601     muonTree->Branch("muon_global_numberOfValidHits",&muon_global_numberOfValidHits,"muon_global_numberOfValidHits/I");
602     muonTree->Branch("muon_global_innerPosition_x",& muon_global_innerPosition_x,"muon_global_innerPosition_x/F");
603     muonTree->Branch("muon_global_innerPosition_y",& muon_global_innerPosition_y,"muon_global_innerPosition_y/F");
604     muonTree->Branch("muon_global_innerPosition_z",& muon_global_innerPosition_z,"muon_global_innerPosition_z/F");
605     muonTree->Branch("muon_global_outerPosition_x",& muon_global_outerPosition_x,"muon_global_outerPosition_x/F");
606     muonTree->Branch("muon_global_outerPosition_y",& muon_global_outerPosition_y,"muon_global_outerPosition_y/F");
607     muonTree->Branch("muon_global_outerPosition_z",& muon_global_outerPosition_z,"muon_global_outerPosition_z/F");
608    
609     muonTree->Branch("muon_STA_chi2", &muon_STA_chi2 , "muon_STA_chi2/F");
610     muonTree->Branch("muon_STA_p",&muon_STA_p," muon_STA_p/F");
611     muonTree->Branch("muon_STA_pt",&muon_STA_pt, "muon_STA_pt/F");
612     muonTree->Branch("muon_STA_outerP",& muon_STA_outerP,"muon_STA_outerP/F");
613     muonTree->Branch("muon_STA_outerPt",& muon_STA_outerPt,"muon_STA_outerPt/F");
614     muonTree->Branch("muon_STA_ndof",&muon_STA_ndof,"muon_STA_ndof/F");
615     muonTree->Branch("muon_STA_normalizedChi2",& muon_STA_normalizedChi2,"muon_STA_normalizedChi2/F");
616     muonTree->Branch("muon_STA_recHitsSize",&muon_STA_recHitsSize,"muon_STA_recHitsSize/I");
617     muonTree->Branch("muon_STA_numberOfLostHits",& muon_STA_numberOfLostHits,"muon_STA_numberOfLostHits/I");
618     muonTree->Branch("muon_STA_numberOfValidHits",&muon_STA_numberOfValidHits,"muon_STA_numberOfValidHits/I");
619     muonTree->Branch("muon_STA_innerPosition_x",& muon_STA_innerPosition_x,"muon_STA_innerPosition_x/F");
620     muonTree->Branch("muon_STA_innerPosition_y",& muon_STA_innerPosition_y,"muon_STA_innerPosition_y/F");
621     muonTree->Branch("muon_STA_innerPosition_z",& muon_STA_innerPosition_z,"muon_STA_innerPosition_z/F");
622     muonTree->Branch("muon_STA_outerPosition_x",& muon_STA_outerPosition_x,"muon_STA_outerPosition_x/F");
623     muonTree->Branch("muon_STA_outerPosition_y",& muon_STA_outerPosition_y,"muon_STA_outerPosition_y/F");
624     muonTree->Branch("muon_STA_outerPosition_z",& muon_STA_outerPosition_z,"muon_STA_outerPosition_z/F");
625    
626     muonTree->Branch("muon_track_chi2", &muon_track_chi2 , "muon_track_chi2/F");
627     muonTree->Branch("muon_track_p",&muon_track_p," muon_track_p/F");
628     muonTree->Branch("muon_track_pt",&muon_track_pt, "muon_track_pt/F");
629     muonTree->Branch("muon_track_outerP",& muon_track_outerP,"muon_track_outerP/F");
630     muonTree->Branch("muon_track_outerPt",& muon_track_outerPt,"muon_track_outerPt/F");
631     muonTree->Branch("muon_track_ndof",&muon_track_ndof,"muon_track_ndof/F");
632     muonTree->Branch("muon_track_normalizedChi2",& muon_track_normalizedChi2,"muon_track_normalizedChi2/F");
633     muonTree->Branch("muon_track_recHitsSize",&muon_track_recHitsSize,"muon_track_recHitsSize/I");
634     muonTree->Branch("muon_track_numberOfLostHits",& muon_track_numberOfLostHits,"muon_track_numberOfLostHits/I");
635     muonTree->Branch("muon_track_numberOfValidHits",&muon_track_numberOfValidHits,"muon_track_numberOfValidHits/I");
636     muonTree->Branch("muon_track_innerPosition_x",& muon_track_innerPosition_x,"muon_track_innerPosition_x/F");
637     muonTree->Branch("muon_track_innerPosition_y",& muon_track_innerPosition_y,"muon_track_innerPosition_y/F");
638     muonTree->Branch("muon_track_innerPosition_z",& muon_track_innerPosition_z,"muon_track_innerPosition_z/F");
639     muonTree->Branch("muon_track_outerPosition_x",& muon_track_outerPosition_x,"muon_track_outerPosition_x/F");
640     muonTree->Branch("muon_track_outerPosition_y",& muon_track_outerPosition_y,"muon_track_outerPosition_y/F");
641     muonTree->Branch("muon_track_outerPosition_z",& muon_track_outerPosition_z,"muon_track_outerPosition_z/F");
642    
643 senka 1.1 JetMuTree->Branch("nb_jet",&nb_jet,"nb_jet/I");
644     JetMuTree->Branch("nb_jet_eta_cut",&nb_jet_eta_cut,"nb_jet_eta_cut/I");
645     JetMuTree->Branch("nb_jet_eta_and_pt_cut",&nb_jet_eta_and_pt_cut,"nb_jet_eta_and_pt_cut/I");
646 senka 1.2
647 senka 1.1 JetMuTree->Branch("nb_muon",&nb_muon,"nb_muon/I");
648    
649 senka 1.2 JetMuTree->Branch("weight", &eventWeight, "weight/F");
650     JetMuTree->Branch("processID", &processID, "processID/I");
651     JetMuTree->Branch("alpgenID", &alpgenID, "alpgenID/I");
652     JetMuTree->Branch("genEventScale", &eventScale, "genEventScale/F");
653     JetMuTree->Branch("RunNumber", &RunNumber, "RunNumber/I");
654     JetMuTree->Branch("EventID", &EventID, "EventID/I");
655    
656 senka 1.1 }
657    
658    
659    
660     ////////////////////////
661     // matching
662     //
663    
664     class MatchingInfo {
665     public:
666    
667     MatchingInfo(float dr, int gen, int reco) : deltaR(dr), genMuon(gen),recoMuon(reco) {}
668    
669     bool operator<(const MatchingInfo m) { return deltaR < m.deltaR;}
670    
671     float deltaR;
672     int genMuon;
673     int recoMuon;
674    
675     };
676    
677    
678     bool betterMatch2(MatchingInfo m1, MatchingInfo m2)
679     { return m1.deltaR < m2.deltaR;}
680    
681    
682     // method that gives the vector<MatchingInfo> matching_Infos containing (dR,i,j); i=genParticle index, j=ParticleCollection index:
683    
684     void MuonAnalyzer::matching(const edm::Event& iEvent, Handle<reco::CandidateCollection> cands, int pdgid){
685    
686     using namespace edm;
687     using namespace reco;
688     using namespace std;
689    
690     //-------------------------
691     Handle<CandidateCollection> mccands;
692     iEvent.getByLabel( theMCTruthCollection,mccands );
693    
694     vector <GenParticleCandidate*> genparticles;
695    
696     for(CandidateCollection::const_iterator p = mccands->begin();
697     p != mccands->end(); p++) {
698    
699     // reco::Candidate* p_tmp = &(*p);
700     reco::GenParticleCandidate* ptr = (reco::GenParticleCandidate*) &(*p);
701    
702     if ( (ptr)->status() == 1
703     /*&& (ptr)->momentum().perp() > 2.*/ ) { // stable particles with pt>2 only
704     if ( abs((ptr)->pdgId()) == pdgid ) {
705     genparticles.push_back((ptr));
706     //cout << "electron MCT\n";
707     }
708     }
709     }
710    
711    
712     int n1=0;
713     vector<MatchingInfo> matching_Infos;
714     for ( unsigned int i=0; i<genparticles.size(); i++ ) // po svim MC muonima
715     {
716    
717     for (unsigned int j = 0; j < cands->size(); j++)// po svim reco muonima
718     {
719    
720     double d_eta2=(*cands)[j].eta()-(genparticles)[i]->eta();
721     double d_phi=fabs((*cands)[j].phi()-(genparticles)[i]->phi());
722     double dR_byhand= sqrt(d_eta2*d_eta2+d_phi*d_phi);
723     double dR=dR_byhand;
724    
725     if (dR<0.15){
726     n1++;
727     matching_Infos.push_back(MatchingInfo(dR,i,j));
728     }
729    
730     }
731     }
732    
733     sort(matching_Infos.begin(),matching_Infos.end(),betterMatch2); // sorting (dR,i,j)
734    
735     if (genparticles.size()!=0 && cands->size()!=0){
736     // Now exclude combinations using same gen or rec muons
737     for (std::vector<MatchingInfo>::iterator it=matching_Infos.begin();
738     it != matching_Infos.end(); it++) {
739     for (std::vector<MatchingInfo>::iterator it2=matching_Infos.begin();
740     it2!=it; it2++) {
741     if ( (it->genMuon == it2->genMuon) || (it->recoMuon == it2->recoMuon) ) {
742     matching_Infos.erase(it);
743     it=matching_Infos.begin();
744     break;
745     }
746     }
747     }
748     }
749    
750     // Now put result in vector of pairs
751    
752     leptonMatching[pdgid].clear();
753    
754     for (vector<MatchingInfo>::iterator match=matching_Infos.begin();
755     match !=matching_Infos.end(); match++) { // po tablici dr,i,j
756    
757     pair<const GenParticleCandidate*, const reco::Candidate *> pair;
758     // vector<pair,float > matchedParticles;
759     pair.first = genparticles[match->genMuon];
760     pair.second = & (*cands)[match->recoMuon];
761    
762     matchedParticles cestice(genparticles[match->genMuon] ,& (*cands)[match->recoMuon] , match->deltaR);
763    
764     leptonMatching[pdgid].push_back(cestice);
765    
766     }
767    
768     }
769    
770 senka 1.2
771 senka 1.1 /////////////////////////////////////////////////////////////////////
772     // made for muon fake rate
773     // method that gives the vector<MatchingInfo> matching_Infos containing (dR,i,j); i=genParticle index, j=ParticleCollection index:
774    
775     void MuonAnalyzer::matching_JetsWithMuons(const edm::Event& iEvent, Handle<reco::CandidateCollection> cands, Handle<reco::CaloJetCollection> jets){
776    
777     using namespace edm;
778     using namespace reco;
779     using namespace std;
780    
781     cout <<" ulazi u matching_JetsWithMuons" <<endl;
782    
783     int n1=0;
784     vector<MatchingInfo> matching_Infos;
785     cout <<"# jets="<<jets->size()<< " # muons="<<cands->size() << endl;
786     for ( unsigned int i=0; i<jets->size(); i++ ) // po svim jetovima
787     {
788     // cout <<" next jet" <<endl;
789     for (unsigned int j = 0; j < cands->size(); j++)// po svim reco muonima
790     {
791     // cout <<" next muon" <<endl;
792     double d_eta2=(*cands)[j].eta()-(*jets)[i].eta();
793     double d_phi=fabs((*cands)[j].phi()-(*jets)[i].phi());
794     double dR_byhand= sqrt(d_eta2*d_eta2+d_phi*d_phi);
795     double dR=dR_byhand;
796     // cout <<"dR="<< dR <<endl;
797     // if (dR<0.15){
798     n1++;
799     matching_Infos.push_back(MatchingInfo(dR,i,j));
800     // }
801    
802     }
803     }
804    
805     sort(matching_Infos.begin(),matching_Infos.end(),betterMatch2); // sorting (dR,i,j)
806    
807     if (jets->size()!=0 && cands->size()!=0){
808     // Now exclude combinations using same gen or rec muons
809     for (std::vector<MatchingInfo>::iterator it=matching_Infos.begin();
810     it != matching_Infos.end(); it++) {
811     for (std::vector<MatchingInfo>::iterator it2=matching_Infos.begin();
812     it2!=it; it2++) {
813     if ( (it->genMuon == it2->genMuon) || (it->recoMuon == it2->recoMuon) ) {
814     matching_Infos.erase(it);
815     it=matching_Infos.begin();
816     break;
817     }
818     }
819     }
820     }
821    
822     // Now t in vector of pairs
823    
824     Muon_Jet_Matching[1].clear();
825    
826     for (vector<MatchingInfo>::iterator match=matching_Infos.begin();
827     match !=matching_Infos.end(); match++) { // po tablici dr,i,j
828    
829     pair<const CaloJet*, const reco::Candidate *> pair;
830     // vector<pair,float > matchedParticles;
831     pair.first = &(*jets)[match->genMuon];
832     pair.second = & (*cands)[match->recoMuon];
833    
834     matchedJets_and_Muons cestice(&(*jets)[match->genMuon] ,& (*cands)[match->recoMuon] , match->deltaR);
835    
836     Muon_Jet_Matching[1].push_back(cestice);
837    
838     }
839    
840     cout <<" nakon pisanja u Muon_Jet_Matching[1]" << endl;
841     for (vector<MatchingInfo>::iterator match=matching_Infos.begin();
842     match !=matching_Infos.end(); match++) { // po tablici dr,i,j
843     // cout <<" nakon pisanja u Muon_Jet_Matching[1]" << endl;
844     for (map< int,
845     std::vector< matchedJets_and_Muons > > ::iterator iter = Muon_Jet_Matching.begin();
846     iter!=Muon_Jet_Matching.end(); iter++) // entire map
847     {
848     for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
849     const reco::Candidate * nesto=iter->second[j].pair.second;
850     // ::OverlapChecker overlap;
851     // if (overlap(*nesto,*daughter)){
852     cout <<"*****dr="<< iter->second[j].delta_R<<endl;
853     // cout <<"*****dr="<< iter->second[j].delta_R<<endl;
854     // }
855     }
856    
857    
858    
859     }
860    
861    
862    
863     }
864    
865     }
866    
867     /////////////////////////////////////////////////////////////////////
868     // made for muon fake rate
869     // method that gives the vector<MatchingInfo> matching_Infos containing (dR,i,j); i=genParticle index, j=ParticleCollection index:
870    
871     void MuonAnalyzer::matching_JetsWithB(const edm::Event& iEvent, Handle<reco::CaloJetCollection> jets){
872    
873     cout <<" ulazi u matching_JetsWithB" <<endl;
874    
875     using namespace edm;
876     using namespace reco;
877     using namespace std;
878    
879     Handle<CandidateCollection> mccands;
880     iEvent.getByLabel( theMCTruthCollection,mccands );
881    
882     vector <GenParticleCandidate*> genBparticles; // vector u kojem su svi B hadroni cija majka nije B hadron (b kvark)
883     // vector <GenParticleCandidate*> genCparticles; // vector u kojem su svi D hadroni cija majka nije D niti B hadron (c kvark)
884    
885    
886     for(CandidateCollection::const_iterator p = mccands->begin();
887     p != mccands->end(); p++) {
888    
889     reco::GenParticleCandidate* ptr = (reco::GenParticleCandidate*) &(*p);
890    
891     ///////////////////////////////////////////////////////
892     if ((((abs(ptr->pdgId()))/100)%10 ==5)||(((abs(ptr->pdgId()))/1000)%10 ==5 )) { // b
893     if ((((abs(ptr->mother()->pdgId()))/100)%10 !=5)&&(((abs(ptr->mother()->pdgId()))/1000)%10 !=5 )) genBparticles.push_back(ptr); // mother is not b
894     }
895    
896     ///////////////////////////////////////////////////////
897    
898    
899     // int pdg_id=ptr->pdgId();
900     // // cout << "particle id : " << pdg_id << endl;
901    
902     // while ( ptr ->mother()!=0 ) { //izvrsava se sve dok ne dodje do nule tj. dok ne trazi majku od protona
903     // // cout << "Going up: ";
904     // // cout << "Mother id : " << genParticle->pdgId() << endl;
905     // ptr= (reco::GenParticleCandidate*) ptr->mother();
906    
907     // if (abs(ptr->pdgId())!=abs(pdg_id)) {
908     // // cout<< " good mother " <<endl;
909    
910     // if ((((abs(ptr->pdgId()))/100)%10 ==5)||(((abs(ptr->pdgId()))/1000)%10 ==5 )) { // mother is b
911     // // while (ptr = const_cast<reco::GenParticleCandidate *>(dynamic_cast<const reco::GenParticleCandidate *>(ptr->mother()))) {
912     // while (ptr ->mother()!=0) {
913    
914     // // ptr= (reco::GenParticleCandidate*) ptr->mother();
915     // // if ((((abs(ptr->pdgId()))/100)%10 !=5)||(((abs(ptr->pdgId()))/1000)%10 !=5 )){ // ako majka nije b
916     // if ((((abs(ptr ->mother()->pdgId()))/100)%10 !=5)||(((abs(ptr ->mother()->pdgId()))/1000)%10 !=5 )){ // ako majka nije b
917     // genBparticles.push_back(ptr);
918     // break;
919     // }
920     // ptr= (reco::GenParticleCandidate*) ptr->mother();
921     // }
922    
923     // }
924    
925     // }
926    
927     // }
928    
929     }
930    
931     int n1=0;
932     vector<MatchingInfo> matching_Infos;
933     cout <<"# jets="<<jets->size()<< " # B hadrons="<<genBparticles.size() << endl;
934     for ( unsigned int i=0; i<genBparticles.size(); i++ ) // po svim MC B hadronima
935     {
936     // cout <<" next B hadron" <<endl;
937     for (unsigned int j = 0; j < jets->size(); j++)// po svim jetovima
938     {
939     // cout <<" next jet" <<endl;
940     double d_eta2=(*jets)[j].eta()-(genBparticles)[i]->eta();
941     double d_phi=fabs((*jets)[j].phi()-(genBparticles)[i]->phi());
942     double dR_byhand= sqrt(d_eta2*d_eta2+d_phi*d_phi);
943     double dR=dR_byhand;
944     // cout <<"dR="<< dR <<endl;
945     // if (dR<0.15){
946     n1++;
947     matching_Infos.push_back(MatchingInfo(dR,i,j));
948     // }
949    
950     }
951     }
952    
953     sort(matching_Infos.begin(),matching_Infos.end(),betterMatch2); // sorting (dR,i,j)
954    
955     if (jets->size()!=0 && genBparticles.size()!=0){
956     // Now exclude combinations using same gen or rec muons
957     for (std::vector<MatchingInfo>::iterator it=matching_Infos.begin();
958     it != matching_Infos.end(); it++) {
959     for (std::vector<MatchingInfo>::iterator it2=matching_Infos.begin();
960     it2!=it; it2++) {
961     if ( (it->genMuon == it2->genMuon) || (it->recoMuon == it2->recoMuon) ) {
962     matching_Infos.erase(it);
963     it=matching_Infos.begin();
964     break;
965     }
966     }
967     }
968     }
969    
970     // Now t in vector of pairs
971    
972     B_Jet_Matching[1].clear();
973    
974     for (vector<MatchingInfo>::iterator match=matching_Infos.begin();
975     match !=matching_Infos.end(); match++) { // po tablici dr,i,j
976    
977     pair<const GenParticleCandidate*, const CaloJet*> pair;
978     // vector<pair,float > matchedParticles;
979     pair.first = genBparticles[match->genMuon];
980     pair.second = &(*jets)[match->recoMuon];
981    
982     matchedJets_and_B cestice(genBparticles[match->genMuon],&(*jets)[match->recoMuon] , match->deltaR);
983    
984     B_Jet_Matching[1].push_back(cestice);
985    
986     }
987    
988     cout <<" nakon pisanja u B_Jet_Matching[1]" << endl;
989     for (vector<MatchingInfo>::iterator match=matching_Infos.begin();
990     match !=matching_Infos.end(); match++) { // po tablici dr,i,j
991     // int j=0;
992     for (map< int,
993     std::vector< matchedJets_and_B > > ::iterator iter = B_Jet_Matching.begin();
994     iter!=B_Jet_Matching.end(); iter++) // entire map
995     {
996     for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
997     const reco::Candidate * nesto=iter->second[j].pair.second;
998     // ::OverlapChecker overlap;
999     // if (overlap(*nesto,*daughter)){
1000     cout <<"*****dr="<< iter->second[j].delta_R<<endl;
1001     // cout <<"*****dr="<< iter->second[j].delta_R<<endl;
1002     // }
1003     }
1004    
1005     }
1006    
1007     }
1008    
1009     }
1010    
1011     /////////////////////////////////////////////////////////////////////
1012     // made for muon fake rate
1013     // method that gives the vector<MatchingInfo> matching_Infos containing (dR,i,j); i=genParticle index, j=ParticleCollection index:
1014    
1015     void MuonAnalyzer::matching_JetsWithC(const edm::Event& iEvent, Handle<reco::CaloJetCollection> jets){
1016    
1017     cout <<" ulazi u matching_JetsWithC" <<endl;
1018    
1019     using namespace edm;
1020     using namespace reco;
1021     using namespace std;
1022    
1023     Handle<CandidateCollection> mccands;
1024     iEvent.getByLabel( theMCTruthCollection,mccands );
1025    
1026     vector <GenParticleCandidate*> genCparticles; // vector u kojem su svi D hadroni cija majka nije D niti B hadron (c kvark)
1027    
1028     for(CandidateCollection::const_iterator p = mccands->begin();
1029     p != mccands->end(); p++) {
1030    
1031     reco::GenParticleCandidate* ptr = (reco::GenParticleCandidate*) &(*p);
1032    
1033     if ((((abs(ptr->pdgId()))/100)%10 ==4)||(((abs(ptr->pdgId()))/1000)%10 ==4 )) { // c
1034     if ((((abs(ptr->mother()->pdgId()))/100)%10 !=4)&&(((abs(ptr->mother()->pdgId()))/1000)%10 !=4 )
1035     && (((abs(ptr->mother()->pdgId()))/100)%10 !=5)&&(((abs(ptr->mother()->pdgId()))/1000)%10 !=5 ))
1036     genCparticles.push_back(ptr); // mother is not b or c
1037     }
1038    
1039     // int pdg_id=ptr->pdgId();
1040     // // cout << "particle id : " << pdg_id << endl;
1041    
1042     // while (ptr->mother()!=0) { //izvrsava se sve dok ne dodje do nule tj. dok ne trazi majku od protona
1043     // // cout << "Going up: ";
1044     // // cout << "Mother id : " << genParticle->pdgId() << endl;
1045     // ptr= (reco::GenParticleCandidate*) ptr->mother();
1046     // if (abs(ptr->pdgId())!=abs(pdg_id)) {
1047     // // cout<< " good mother " <<endl;
1048    
1049    
1050     // if ((((abs(ptr->pdgId()))/100)%10 ==4)||(((abs(ptr->pdgId()))/1000)%10 ==4 )) { // mother is c
1051     // while (ptr->mother()!=0) {
1052     // // ptr= (reco::GenParticleCandidate*) ptr->mother();
1053     // if (((((abs(ptr->mother()->pdgId()))/100)%10 !=4)||(((abs(ptr->mother()->pdgId()))/1000)%10 !=4 )) &&
1054     // ((((abs(ptr->mother()->pdgId()))/100)%10 !=5)||(((abs(ptr->mother()->pdgId()))/1000)%10 !=5 ))){ // ako majka nije c niti b
1055     // genCparticles.push_back(ptr);
1056     // break;
1057     // }
1058     // ptr= (reco::GenParticleCandidate*) ptr->mother();
1059     // }
1060     // }
1061    
1062    
1063     // }
1064     // }
1065    
1066     }
1067    
1068     int n1=0;
1069     vector<MatchingInfo> matching_Infos;
1070     cout <<"# jets="<<jets->size()<< " # C hadrons="<<genCparticles.size() << endl;
1071     for ( unsigned int i=0; i<genCparticles.size(); i++ ) // po svim MC B hadronima
1072     {
1073     // cout <<" next B hadron" <<endl;
1074     for (unsigned int j = 0; j < jets->size(); j++)// po svim jetovima
1075     {
1076     // cout <<" next jet" <<endl;
1077     double d_eta2=(*jets)[j].eta()-(genCparticles)[i]->eta();
1078     double d_phi=fabs((*jets)[j].phi()-(genCparticles)[i]->phi());
1079     double dR_byhand= sqrt(d_eta2*d_eta2+d_phi*d_phi);
1080     double dR=dR_byhand;
1081     // cout <<"dR="<< dR <<endl;
1082     // if (dR<0.15){
1083     n1++;
1084     matching_Infos.push_back(MatchingInfo(dR,i,j));
1085     // }
1086    
1087     }
1088     }
1089    
1090     sort(matching_Infos.begin(),matching_Infos.end(),betterMatch2); // sorting (dR,i,j)
1091    
1092     if (jets->size()!=0 && genCparticles.size()!=0){
1093     // Now exclude combinations using same gen or rec muons
1094     for (std::vector<MatchingInfo>::iterator it=matching_Infos.begin();
1095     it != matching_Infos.end(); it++) {
1096     for (std::vector<MatchingInfo>::iterator it2=matching_Infos.begin();
1097     it2!=it; it2++) {
1098     if ( (it->genMuon == it2->genMuon) || (it->recoMuon == it2->recoMuon) ) {
1099     matching_Infos.erase(it);
1100     it=matching_Infos.begin();
1101     break;
1102     }
1103     }
1104     }
1105     }
1106    
1107     // Now t in vector of pairs
1108    
1109     C_Jet_Matching[1].clear();
1110    
1111     for (vector<MatchingInfo>::iterator match=matching_Infos.begin();
1112     match !=matching_Infos.end(); match++) { // po tablici dr,i,j
1113    
1114     pair<const GenParticleCandidate*, const CaloJet*> pair;
1115     // vector<pair,float > matchedParticles;
1116     pair.first = genCparticles[match->genMuon];
1117     pair.second = &(*jets)[match->recoMuon];
1118    
1119     matchedJets_and_C cestice(genCparticles[match->genMuon],&(*jets)[match->recoMuon] , match->deltaR);
1120    
1121     C_Jet_Matching[1].push_back(cestice);
1122    
1123     }
1124    
1125     cout <<" nakon pisanja u C_Jet_Matching[1]" << endl;
1126     for (vector<MatchingInfo>::iterator match=matching_Infos.begin();
1127     match !=matching_Infos.end(); match++) { // po tablici dr,i,j
1128     for (map< int,
1129     std::vector< matchedJets_and_B > > ::iterator iter = B_Jet_Matching.begin();
1130     iter!=B_Jet_Matching.end(); iter++) // entire map
1131     {
1132     for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
1133     const reco::Candidate * nesto=iter->second[j].pair.second;
1134     // ::OverlapChecker overlap;
1135     // if (overlap(*nesto,*daughter)){
1136     cout <<"*****dr="<< iter->second[j].delta_R<<endl;
1137     // cout <<"*****dr="<< iter->second[j].delta_R<<endl;
1138     // }
1139     }
1140    
1141     }
1142    
1143    
1144    
1145     }
1146    
1147     }
1148    
1149    
1150     //---------------------------------------------------
1151    
1152    
1153     // get mother of particle
1154     // returns mother = 1 if mother is Z boson
1155     // returns mother = 2 if mother is W boson
1156     // returns mother = 4 if mother is b
1157     // returns mother = 0 else
1158    
1159     MuonAnalyzer::LeptonOrigin MuonAnalyzer::getMother(const reco::Candidate* genParticle){
1160    
1161     int pdg_id=genParticle->pdgId();
1162     // cout << "particle id : " << pdg_id << endl;
1163    
1164     while (genParticle = genParticle->mother() ) { //izvrsava se sve dok ne dodje do nule tj. dok ne trazi majku od protona
1165     // cout << "Going up: ";
1166     // cout << "Mother id : " << genParticle->pdgId() << endl;
1167     if (abs(genParticle->pdgId())!=abs(pdg_id)) {
1168     // cout<< " good mother " <<endl;
1169     if (abs(genParticle->pdgId())==23){ // mother is Z
1170     return zboson;
1171     }
1172     if (abs(genParticle->pdgId())==24) { // mother is W
1173     return wboson;
1174     }
1175     // if (abs(genParticle->pdgId())==23 || abs(genParticle->pdgId())==24) { // mother is W or Z
1176     // WZ_matched=1;
1177     // mother=3;
1178     // }
1179     if ((((abs(genParticle->pdgId()))/100)%10 ==5)||(((abs(genParticle->pdgId()))/1000)%10 ==5 )) { // mother is b
1180     return bdecay;
1181     }
1182     if ((((abs(genParticle->pdgId()))/100)%10 ==4)||(((abs(genParticle->pdgId()))/1000)%10 ==4 )) { // mother is c
1183     return cdecay;
1184     }
1185 senka 1.2 if (abs(genParticle->pdgId()==6)) { // mother is t
1186     return tdecay;
1187     }
1188 senka 1.1 return other;
1189     }
1190     }
1191    
1192     return other;
1193     }
1194    
1195 senka 1.2 const reco::Candidate * MuonAnalyzer::getMother_particle(const reco::Candidate* genParticle){
1196    
1197     int pdg_id=genParticle->pdgId();
1198     // cout << "particle id : " << pdg_id << endl;
1199    
1200     while (genParticle = genParticle->mother() ) { //izvrsava se sve dok ne dodje do nule tj. dok ne trazi majku od protona
1201     // cout << "Going up: ";
1202     // cout << "Mother id : " << genParticle->pdgId() << endl;
1203     if (abs(genParticle->pdgId())!=abs(pdg_id)) {
1204     // cout<< " good mother " <<endl;
1205     return genParticle->mother();
1206     }
1207     }
1208    
1209     // return ;
1210     }
1211    
1212    
1213    
1214 senka 1.1 ///////////////////
1215     int MuonAnalyzer::getMother_Id(const reco::Candidate* genParticle){
1216    
1217     int mother_id=0;
1218    
1219     int pdg_id=genParticle->pdgId();
1220     // cout << "particle id : " << pdg_id << endl;
1221    
1222     while (genParticle = genParticle->mother() ) { //izvrsava se sve dok ne dodje do nule tj. dok ne trazi majku od protona
1223     if (abs(genParticle->pdgId())!=abs(pdg_id)) {
1224     mother_id=genParticle->pdgId();
1225     return mother_id;
1226     }
1227     }
1228    
1229     return mother_id;
1230     }
1231    
1232     /////////////////
1233    
1234    
1235     const reco::Candidate * MuonAnalyzer::MatchedGenParticle(const reco::Candidate * daughter){
1236     ::OverlapChecker overlap;
1237     matched_genParticle=0;
1238    
1239     for (map< int,
1240     std::vector< matchedParticles > > ::iterator iter = leptonMatching.begin();
1241     iter!=leptonMatching.end(); iter++) // entire map
1242     {
1243     for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
1244     const reco::Candidate * nesto=iter->second[j].pair.second;
1245     if (overlap(*nesto,*daughter)){
1246     matched_genParticle=iter->second[j].pair.first;
1247     }
1248     }
1249     }
1250     return matched_genParticle;
1251     }
1252    
1253     /////////////////////////////////////////
1254     // for muon fake rate
1255    
1256     const reco::Candidate * MuonAnalyzer::MatchedMuonWithJet(const reco::CaloJet * daughter){
1257     ::OverlapChecker overlap;
1258     matched_Muon=0;
1259    
1260     for (map< int,
1261     std::vector< matchedJets_and_Muons > > ::iterator iter = Muon_Jet_Matching.begin();
1262     iter!=Muon_Jet_Matching.end(); iter++) // entire map
1263     {
1264     // cout <<" ide preko cijele mape za jet" <<endl;
1265     for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
1266     const reco::CaloJet * nesto=iter->second[j].pair.first;
1267     if (overlap(*nesto,*daughter)){
1268     matched_Muon=iter->second[j].pair.second; /// tu je bila greska
1269     }
1270     }
1271     }
1272     return matched_Muon;
1273     }
1274    
1275    
1276     /////////////////////////////
1277     // Jet and B hadron matching
1278    
1279     const reco::Candidate * MuonAnalyzer::MatchedBhadronWithJet(const reco::CaloJet * daughter){
1280     ::OverlapChecker overlap;
1281     matched_B=0;
1282    
1283     for (map< int,
1284     std::vector< matchedJets_and_B > > ::iterator iter = B_Jet_Matching.begin();
1285     iter!=B_Jet_Matching.end(); iter++) // entire map
1286     {
1287     // cout <<" ide preko cijele mape za jet" <<endl;
1288     for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
1289     const reco::CaloJet * nesto=iter->second[j].pair.second;
1290     if (overlap(*nesto,*daughter)){
1291     matched_B=iter->second[j].pair.first; /// tu je bila greska
1292     }
1293     }
1294     }
1295     return matched_B;
1296     }
1297    
1298     /////////////////////////////
1299     // Jet and C hadron matching
1300    
1301     const reco::Candidate * MuonAnalyzer::MatchedChadronWithJet(const reco::CaloJet * daughter){
1302     ::OverlapChecker overlap;
1303     matched_C=0;
1304    
1305     for (map< int,
1306     std::vector< matchedJets_and_C > > ::iterator iter = C_Jet_Matching.begin();
1307     iter!=C_Jet_Matching.end(); iter++) // entire map
1308     {
1309     // cout <<" ide preko cijele mape za jet" <<endl;
1310     for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
1311     const reco::CaloJet * nesto=iter->second[j].pair.second;
1312     if (overlap(*nesto,*daughter)){
1313     matched_C=iter->second[j].pair.first; /// tu je bila greska
1314     }
1315     }
1316     }
1317     return matched_C;
1318     }
1319    
1320    
1321     float MuonAnalyzer::dR(const reco::Candidate * daughter){
1322    
1323     ::OverlapChecker overlap;
1324    
1325     float delta_r = -10;
1326     for (map< int,
1327     std::vector< matchedParticles > > ::iterator iter = leptonMatching.begin();
1328     iter!=leptonMatching.end(); iter++) // entire map
1329     {
1330     for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
1331     const reco::Candidate * nesto=iter->second[j].pair.second;
1332     if (overlap(*nesto,*daughter)){
1333     delta_r=iter->second[j].delta_R;
1334     }
1335     }
1336     }
1337     return delta_r;
1338     }
1339    
1340     float MuonAnalyzer::dR_Muon_Jet(const reco::CaloJet * daughter){
1341    
1342     ::OverlapChecker overlap;
1343    
1344     float delta_r = -10;
1345     for (map< int,
1346     std::vector<matchedJets_and_Muons > > ::iterator iter = Muon_Jet_Matching.begin();
1347     iter!=Muon_Jet_Matching.end(); iter++) // entire map
1348     {
1349     for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
1350     const reco::CaloJet * nesto=iter->second[j].pair.first;
1351     if (overlap(*nesto,*daughter)){
1352     delta_r=iter->second[j].delta_R;
1353     }
1354     }
1355     }
1356     return delta_r;
1357     }
1358    
1359     float MuonAnalyzer::dR_B_Jet(const reco::CaloJet * daughter){
1360    
1361     ::OverlapChecker overlap;
1362    
1363     float delta_r = -10;
1364     for (map< int,
1365     std::vector<matchedJets_and_B > > ::iterator iter = B_Jet_Matching.begin();
1366     iter!=B_Jet_Matching.end(); iter++) // entire map
1367     {
1368     for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
1369     const reco::CaloJet * nesto=iter->second[j].pair.second;
1370     if (overlap(*nesto,*daughter)){
1371     delta_r=iter->second[j].delta_R;
1372     }
1373     }
1374     }
1375     return delta_r;
1376     }
1377    
1378     float MuonAnalyzer::dR_C_Jet(const reco::CaloJet * daughter){
1379    
1380     ::OverlapChecker overlap;
1381    
1382     float delta_r = -10;
1383     for (map< int,
1384     std::vector<matchedJets_and_C > > ::iterator iter = C_Jet_Matching.begin();
1385     iter!=C_Jet_Matching.end(); iter++) // entire map
1386     {
1387     for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
1388     const reco::CaloJet * nesto=iter->second[j].pair.second;
1389     if (overlap(*nesto,*daughter)){
1390     delta_r=iter->second[j].delta_R;
1391     }
1392     }
1393     }
1394     return delta_r;
1395     }
1396    
1397    
1398    
1399     void MuonAnalyzer::collectMCsummary(Handle <reco::CandidateCollection> mccands) {
1400    
1401     using namespace reco;
1402    
1403     //collections
1404    
1405     MCleptZ1_pdgid = -1;
1406     MCleptZ2_pdgid = -1;
1407     MCleptW_pdgid = -1;
1408    
1409     MCleptZ1_pt = -1;
1410     MCleptZ2_pt = -1;
1411     MCleptW_pt = -1;
1412    
1413     MCleptZ1_eta = -1;
1414     MCleptZ2_eta = -1;
1415     MCleptW_eta = -1;
1416    
1417     MCleptZ1_phi = -1;
1418     MCleptZ2_phi = -1;
1419     MCleptW_phi = -1;
1420    
1421     vector<reco::GenParticleCandidate*> Tau;
1422     vector<reco::GenParticleCandidate*> StableMuons;
1423     vector<reco::GenParticleCandidate*> StableElectrons;
1424    
1425    
1426    
1427     for(CandidateCollection::const_iterator p = mccands->begin();
1428     p != mccands->end(); p++) {
1429    
1430     // reco::Candidate* p_tmp = &(*p);
1431     reco::GenParticleCandidate* ptr = (reco::GenParticleCandidate*) &(*p);
1432    
1433     if ( (ptr)->status() == 1
1434     /*&& (ptr)->momentum().perp() > 2.*/ ) { // stable particles with pt>2 only
1435     if ( abs((ptr)->pdgId()) == 11 ) {
1436    
1437     StableElectrons.push_back((ptr));
1438     //cout << "electron MCT\n";
1439     }
1440     else if ( abs((ptr)->pdgId()) == 13 ) {
1441    
1442     StableMuons.push_back((ptr)) ;
1443     //cout << "muon MCT\n";
1444     }
1445     }
1446     else if ((ptr)->status() == 2) {
1447     if ( abs((ptr)->pdgId()) == 15 ) {//tau
1448     Tau.push_back((ptr));
1449     }
1450     }
1451     }
1452    
1453     // std::cout << "# Electrons : " << StableElectrons.size()
1454     // << "# Muons : " << StableMuons.size() << std::endl
1455     // << "# Tau : " << Tau.size() << std::endl;
1456    
1457    
1458     vector<reco::GenParticleCandidate*> * StableLeptons[3] = {&StableElectrons, &StableMuons, &Tau};
1459    
1460     vector<reco::GenParticleCandidate*> TauChildLeptons;
1461    
1462    
1463     //erase electrons and muons from tau from GenparticleCandidate Vector and put to Tau Vector
1464     for (int i=1; i>=0; i--) {
1465     for ( vector<reco::GenParticleCandidate*>::iterator lepton = StableLeptons[i]->begin();
1466     lepton != StableLeptons[i]->end();lepton++) {
1467    
1468     reco::GenParticleCandidate* mcParticleRef = *lepton;
1469    
1470     int parentItr=0;
1471     while ( mcParticleRef->mother()!=0 && parentItr<2) {
1472    
1473     parentItr++;
1474    
1475     mcParticleRef = (reco::GenParticleCandidate*) mcParticleRef->mother();
1476     if (mcParticleRef==0) break;
1477    
1478     //if tau
1479     if (abs((mcParticleRef)->pdgId())==15 ) {
1480     //put into collection
1481     TauChildLeptons.push_back(*lepton);
1482     StableLeptons[i]->erase(lepton);
1483     lepton--;
1484     break;
1485     }
1486     }
1487     }
1488     }
1489    
1490    
1491     bool firstZlept = true;
1492     int MC_tauDecayTypeZ1 = 0;
1493     int MC_tauDecayTypeZ2 = 0;
1494     int MC_tauDecayTypeW = 0;
1495    
1496    
1497     for (int i=2; i>=0; i--) {
1498     while (StableLeptons[i]->size() > 0) {
1499     float maxPt = 0;
1500     vector<reco::GenParticleCandidate*>::iterator index = StableLeptons[i]->end();
1501    
1502     for ( vector<reco::GenParticleCandidate*>::iterator lepton = StableLeptons[i]->begin();
1503     lepton != StableLeptons[i]->end(); lepton++ ) {
1504    
1505     if ((*lepton)->pt() > maxPt) {
1506     maxPt = (*lepton)->pt();
1507     index = lepton;
1508     }
1509     }
1510     bool Zchild = false;
1511     bool Wchild = false;
1512     bool isTau = false;
1513    
1514    
1515     reco::GenParticleCandidate* mcParticleRef = *index;
1516    
1517     if (abs((*index)->pdgId()) == 15) isTau = true;
1518    
1519     //find W/Z mother
1520     int parentItr=0;
1521     while ( mcParticleRef->mother()!=0 && parentItr<2) {
1522    
1523     if (parentItr>=2) break;
1524     parentItr++;
1525     mcParticleRef = (reco::GenParticleCandidate*) mcParticleRef->mother();
1526     if (mcParticleRef==0) break;
1527    
1528     if (abs((mcParticleRef)->pdgId())==23 ) {
1529     Zchild=true;
1530     break;
1531     }
1532     if (abs((mcParticleRef)->pdgId())==24 ) {
1533     Wchild=true;
1534     break;
1535     }
1536     }
1537    
1538    
1539     if (maxPt >0) {
1540     int * MCtauDecayTypePtr = 0;
1541    
1542     if (Wchild) {
1543     MCleptW_pdgid=(*index)->pdgId();
1544     MCleptW_pt=(float)(*index)->pt();
1545     MCleptW_eta=(float)(*index)->eta();
1546     MCleptW_phi=(float)(*index)->phi();
1547     if (isTau) {
1548     MCtauDecayTypePtr = &MC_tauDecayTypeW;
1549     MC_tauDecayTypeW =3;//default to hadronic decay
1550     }
1551    
1552     }
1553     if (Zchild) {
1554     if (firstZlept) {
1555     firstZlept=false;
1556     MCleptZ1_pdgid=(*index)->pdgId();
1557     MCleptZ1_pt=(float)(*index)->pt();
1558     MCleptZ1_eta=(float)(*index)->eta();
1559     MCleptZ1_phi=(float)(*index)->phi();
1560     if (isTau) {
1561     MCtauDecayTypePtr = &MC_tauDecayTypeZ1;
1562     MC_tauDecayTypeZ1 =3;
1563     }
1564    
1565     }
1566     else {
1567     MCleptZ2_pdgid=(*index)->pdgId();
1568     MCleptZ2_pt=(float)(*index)->pt();
1569     MCleptZ2_eta=(float)(*index)->eta();
1570     MCleptZ2_phi=(float)(*index)->phi();
1571     if (isTau) {
1572     MCtauDecayTypePtr = &MC_tauDecayTypeZ2;
1573     MC_tauDecayTypeZ2 =3;
1574     }
1575    
1576     }
1577     }
1578     //check type of tau decay
1579     if (MCtauDecayTypePtr) {
1580     for( vector<reco::GenParticleCandidate*>::iterator p = TauChildLeptons.begin();p != TauChildLeptons.end(); p++) {
1581     int pitr = 0;
1582     reco::GenParticleCandidate* mcParticleRef = *p;
1583     while (mcParticleRef->mother() && pitr<2) {
1584     pitr++;
1585     mcParticleRef =(reco::GenParticleCandidate*) mcParticleRef->mother();
1586     if (mcParticleRef==0) break;
1587     if (mcParticleRef == *index) {
1588     if (abs((*p)->pdgId()) == 11) *MCtauDecayTypePtr = 1;
1589     if (abs((*p)->pdgId()) == 13) *MCtauDecayTypePtr = 2;
1590     }
1591     }
1592     }
1593     }
1594     }
1595     StableLeptons[i]->erase(index);
1596     }
1597     }
1598     }
1599     //define this as a plug-in
1600     // DEFINE_FWK_MODULE(WZAnalyzer);
1601