ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Vuko/WZAnalysis/src/MuonAnalyzer.cc
Revision: 1.7
Committed: Fri Apr 4 10:22:44 2008 UTC (17 years ago) by senka
Content type: text/plain
Branch: MAIN
CVS Tags: keti_June6_02, keti_June6_01, SB_23Apr08, keti_Apr22_01, vuko-10apr08, srecko_10_Apr_v2, srecko_10_Apr, keti_Apr10_01, keti_Apr9_02, vuko-9apr08-2, vuko-9apr08
Changes since 1.6: +40 -40 lines
Log Message:
check if track pointer exist

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