ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Vuko/WZAnalysis/src/MuonAnalyzer.cc
Revision: 1.9
Committed: Fri Jul 4 18:52:23 2008 UTC (16 years, 9 months ago) by ymaravin
Content type: text/plain
Branch: MAIN
CVS Tags: keti_Aug25_01, keti-Aug5-2008-01, keti-Aug4-2008-01, ym-1-6-12-01, HEAD
Changes since 1.8: +6 -6 lines
Error occurred while calculating annotation data.
Log Message:
Updates for 1_6_12 release

File Contents

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