ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Vuko/WZAnalysis/src/MuonEfficiencyAnalyzer.cc
Revision: 1.1
Committed: Mon Jun 23 00:31:11 2008 UTC (16 years, 10 months ago) by vuko
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
Error occurred while calculating annotation data.
Log Message:
muon eff. analysis

File Contents

# Content
1 // -*- C++ -*-
2 //
3 // Package: MuonEfficiencyAnalyzer
4 // Class: MuonEfficiencyAnalyzer
5 //
6 /**\class MuonEfficiencyAnalyzer MuonEfficiencyAnalyzer.cc Vuko/WZAnalysis/src/MuonEfficiencyAnalyzer.cc
7
8 Description: <one line class summary>
9
10 Implementation:
11 <Notes on implementation>
12 */
13 //
14 // Original Author: Vuko Brigljevic
15 //
16 // $Id:
17 //
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/MuonEfficiencyAnalyzer.h"
35 #include "DataFormats/Candidate/interface/OverlapChecker.h"
36
37 #include "Vuko/WZAnalysis/interface/ElectronProperties.h"
38 #include "Vuko/WZAnalysis/interface/MuonProperties.h"
39 #include "Vuko/WZAnalysis/interface/AnglesUtil.h"
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 //#include "TVector3.h"
61
62 #include <map>
63
64 //
65 // constants, enums and typedefs
66 //
67
68 //
69 // static data member definitions
70 //
71
72 //
73 // constructors and destructor
74 //
75 MuonEfficiencyAnalyzer::MuonEfficiencyAnalyzer(const edm::ParameterSet& iConfig)
76
77 {
78 //now do what ever initialization is needed
79 // Recreate = iConfig.getUntrackedParameter<bool>("SecondRound");
80 fOutputFileName = iConfig.getUntrackedParameter<string>("HistOutFile");
81 theMCTruthCollection = iConfig.getParameter<edm::InputTag>("genParticles");
82
83 // theMCTruthCollection = iConfig.getParameter<edm::InputTag>("genParticleCollection");
84
85 theGoodMuonCollection = iConfig.getParameter<edm::InputTag>("GoodMuons");
86
87 theLooseElectronCollection = iConfig.getParameter<edm::InputTag>("LooseElectrons");
88
89 storeHLTresults=false;
90 storeHLTresults = iConfig.getParameter<bool>("storeHLTresults");
91
92 triggerResultsTag = iConfig.getParameter<edm::InputTag>("TriggerResults");
93 firstTriggerResult = true;
94
95 }
96
97
98 MuonEfficiencyAnalyzer::~MuonEfficiencyAnalyzer()
99 {
100
101 // do anything here that needs to be done at desctruction time
102 // (e.g. close files, deallocate resources etc.)
103
104 }
105
106
107 //
108 // member functions
109 //
110
111 // ------------ method called to for each event ------------
112 void
113 MuonEfficiencyAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
114 {
115 using namespace edm;
116 using namespace reco;
117 using namespace std;
118
119 ///////////////////////////////////////
120 //
121 /// EVENT INITIALISATION
122 ///
123
124 //read run number and event ID number
125 RunNumber = iEvent.id().run();
126 EventID = iEvent.id().event();
127
128 // Reset a few things
129 for (map<string,wzana::LeptonProperties* >::iterator iter = leptonProperties.begin();
130 iter!=leptonProperties.end(); iter++)
131 {
132 iter->second->ResetValues();
133 }
134
135
136 const Candidate * theZCandidate = 0;
137 const Candidate * theWlepton = 0;
138
139
140 ////////////////////////////////////////
141 // Get lists
142
143 // loose electrons
144 Handle<CandidateCollection> looseElectronCands;
145 iEvent.getByLabel( theLooseElectronCollection , looseElectronCands);
146
147
148 // good muons
149 Handle<CandidateCollection> goodMuonCands;
150 iEvent.getByLabel( theGoodMuonCollection , goodMuonCands);
151
152
153
154 // get hold of TriggerResults
155 Handle<TriggerResults> HLTR; //moja promjena: ovo prebaceno 2 reda nize!
156 if (storeHLTresults) {
157 // Handle<TriggerResults> HLTR;
158 iEvent.getByLabel(triggerResultsTag,HLTR);
159 if (firstTriggerResult) {
160 firstTriggerResult = false;
161 trigNames.init((edm::TriggerResults&)*HLTR);
162 numTriggers = trigNames.size();
163 }
164 }
165
166 if (storeHLTresults) {
167
168 //def: 1:single electron 2:single relaxed 3:double electron 4: double relaxed 5,6,7,8: muon e+mu: 9
169
170 triggerBitmask=0;
171 triggerBitmaskQCD=0;
172
173 for (size_t i=0; i<numTriggers; ++i) {
174
175 if (!HLTR->wasrun(i)) {
176 continue;
177 }
178 if (HLTR->error(i) ) {
179 continue;
180 }
181 if (HLTR->accept(i)) {
182 //if (triggerNames[i].compare("mcpath")!=0)
183 //TRaccept=true;
184 }else {
185 continue;
186 }
187 if (trigNames.triggerName(i).compare("HLT1Electron")==0) triggerBitmask |= 1;
188 if (trigNames.triggerName(i).compare("HLT1ElectronRelaxed")==0) triggerBitmask |= 2;
189 if (trigNames.triggerName(i).compare("HLT2Electron")==0) triggerBitmask |= 4;
190 if (trigNames.triggerName(i).compare("HLT2ElectronRelaxed")==0) triggerBitmask |= 8;
191
192 if (trigNames.triggerName(i).compare("HLT1MuonIso")==0) triggerBitmask |= 16;
193 if (trigNames.triggerName(i).compare("HLT1MuonNonIso")==0) triggerBitmask |= 32;
194 if (trigNames.triggerName(i).compare("HLT2MuonNonIso")==0) triggerBitmask |= 64;
195 if (trigNames.triggerName(i).compare("HLT2MuonZ")==0) triggerBitmask |= 128;
196
197 if (trigNames.triggerName(i).compare("HLTXElectronMuon")==0) triggerBitmask |= 256;
198 if (trigNames.triggerName(i).compare("HLTXElectronMuonRelaxed")==0) triggerBitmask |= 512;
199
200 if (trigNames.triggerName(i).compare("CandHLT2MuonIso")==0) triggerBitmask |= 1024;
201
202 //for QCD control sample
203 if (trigNames.triggerName(i).compare("HLT1Jet")==0) triggerBitmaskQCD |= 1;
204 if (trigNames.triggerName(i).compare("HLT2Jet")==0) triggerBitmaskQCD |= 2;
205 if (trigNames.triggerName(i).compare("HLT3Jet")==0) triggerBitmaskQCD |= 4;
206 if (trigNames.triggerName(i).compare("HLT4Jet")==0) triggerBitmaskQCD |= 8;
207 if (trigNames.triggerName(i).compare("HLT2JetAco")==0) triggerBitmaskQCD |= 16;
208 if (trigNames.triggerName(i).compare("HLTXMuonJets")==0) triggerBitmaskQCD |= 32;
209 if (trigNames.triggerName(i).compare("HLT1MET")==0) triggerBitmaskQCD |= 64;
210 if (trigNames.triggerName(i).compare("HLT1Jet1METAco")==0) triggerBitmaskQCD |= 128;
211 if (trigNames.triggerName(i).compare("HLT1Jet1MET")==0) triggerBitmaskQCD |= 256;
212 if (trigNames.triggerName(i).compare("HLT2Jet1MET")==0) triggerBitmaskQCD |= 512;
213 if (trigNames.triggerName(i).compare("HLT3Jet1MET")==0) triggerBitmaskQCD |= 1024;
214 if (trigNames.triggerName(i).compare("HLT4Jet1MET")==0) triggerBitmaskQCD |= 2048;
215 if (trigNames.triggerName(i).compare("HLT1MET1HT")==0) triggerBitmaskQCD |= 4096;
216 if (trigNames.triggerName(i).compare("HLTS2JetAco")==0) triggerBitmaskQCD |= 8192;
217
218 }
219 }
220
221 // One and only one electron in the event
222 if (looseElectronCands->size() != 1) return;
223
224
225 matching(iEvent,goodMuonCands, 13);
226 matching(iEvent,looseElectronCands, 11);
227
228
229
230 for(CandidateCollection::const_iterator electron = looseElectronCands->begin();
231 electron != looseElectronCands->end(); electron++) {
232
233 leptonProperties["electron"]->FillProperties(&(*electron),
234 iEvent, iSetup,
235 MatchedGenParticle(&(*electron)),
236 dR(&(*electron)));
237
238
239 }
240
241
242 for(CandidateCollection::const_iterator muon = goodMuonCands->begin();
243 muon != goodMuonCands->end(); muon++) {
244
245 leptonProperties["muon"]->FillProperties(&(*muon),
246 iEvent, iSetup,
247 MatchedGenParticle(&(*muon)),
248 dR(&(*muon)));
249
250 muonTree->Fill();
251
252 }
253 EventCounter++;
254 }
255
256
257 // ------------ method called once each job just before starting event loop ------------
258 void
259 MuonEfficiencyAnalyzer::beginJob(const edm::EventSetup&)
260 {
261
262 using namespace wzana;
263 theFile = new TFile( fOutputFileName.c_str(), "RECREATE" ) ;
264 // theFile = new TFile( "wz.root","UPDATE") ;
265
266 muonTree = new TTree("MuonTree","Muon information");
267
268
269 // jetTree = new TTree("JetTree","jet deltaR information");
270
271 leptonProperties["electron"] = new ElectronProperties();
272 leptonProperties["muon"] = new MuonProperties(); // SAME COMMENT AS BELOW: TEMPORARY FIX
273
274
275 leptonProperties["electron"]->CreateBranches(muonTree, "Ele_");
276 leptonProperties["muon"]->CreateBranches(muonTree, "Mu_");
277
278
279 initialiseTree();
280
281
282 //prepare HLT TriggerResults branch
283 if (storeHLTresults) {
284 muonTree->Branch("WZ_hltBitmask",&triggerBitmask,"WZ_hltBitmask/I");
285 muonTree->Branch("QCD_hltBitmask",&triggerBitmaskQCD,"QCD_hltBitmask/I");
286 }
287
288 EventCounter=0;
289 }
290
291 // ------------ method called once each job just after ending the event loop ------------
292 void
293 MuonEfficiencyAnalyzer::endJob() {
294 // cout << "In the end job with treeNumber = " << treeNumber << endl;
295 theFile->Write();
296 theFile->Close();
297
298 }
299
300
301 void MuonEfficiencyAnalyzer::initialiseTree() {
302
303
304
305 }
306
307
308
309 ////////////////////////
310 // matching
311 //
312
313 class MatchingInfo {
314 public:
315
316 MatchingInfo(float dr, int gen, int reco) : deltaR(dr), genMuon(gen),recoMuon(reco) {}
317
318 bool operator<(const MatchingInfo m) { return deltaR < m.deltaR;}
319
320 float deltaR;
321 int genMuon;
322 int recoMuon;
323
324 };
325
326
327 bool betterGenMatch(MatchingInfo m1, MatchingInfo m2)
328 { return m1.deltaR < m2.deltaR;}
329
330
331 // method that gives the vector<MatchingInfo> matching_Infos containing (dR,i,j); i=genParticle index, j=ParticleCollection index:
332
333 void MuonEfficiencyAnalyzer::matching(const edm::Event& iEvent, Handle<reco::CandidateCollection> cands, int pdgid){
334
335 using namespace edm;
336 using namespace reco;
337 using namespace std;
338
339 //-------------------------
340 Handle<CandidateCollection> mccands;
341 iEvent.getByLabel( theMCTruthCollection,mccands );
342
343 vector <GenParticleCandidate*> genparticles;
344
345 for(CandidateCollection::const_iterator p = mccands->begin();
346 p != mccands->end(); p++) {
347
348 // reco::Candidate* p_tmp = &(*p);
349 reco::GenParticleCandidate* ptr = (reco::GenParticleCandidate*) &(*p);
350
351 if ( (ptr)->status() == 1
352 /*&& (ptr)->momentum().perp() > 2.*/ ) { // stable particles with pt>2 only
353 if ( abs((ptr)->pdgId()) == pdgid ) {
354 genparticles.push_back((ptr));
355 //cout << "electron MCT\n";
356 }
357 }
358 }
359
360
361 int n1=0;
362 vector<MatchingInfo> matching_Infos;
363 for ( unsigned int i=0; i<genparticles.size(); i++ ) // po svim MC muonima
364 {
365
366 for (unsigned int j = 0; j < cands->size(); j++)// po svim reco muonima
367 {
368
369 double d_eta2=(*cands)[j].eta()-(genparticles)[i]->eta();
370 //double d_phi=fabs((*cands)[j].phi()-(genparticles)[i]->phi());
371 double d_phi= kinem::delta_phi((*cands)[j].phi(),(genparticles)[i]->phi()) ;
372 double dR_byhand= sqrt(d_eta2*d_eta2+d_phi*d_phi);
373 double dR=dR_byhand;
374
375 if (dR<0.15){
376 n1++;
377 matching_Infos.push_back(MatchingInfo(dR,i,j));
378 }
379
380 }
381 }
382
383 sort(matching_Infos.begin(),matching_Infos.end(),betterGenMatch); // sorting (dR,i,j)
384
385 if (genparticles.size()!=0 && cands->size()!=0){
386 // Now exclude combinations using same gen or rec muons
387 for (std::vector<MatchingInfo>::iterator it=matching_Infos.begin();
388 it != matching_Infos.end(); it++) {
389 for (std::vector<MatchingInfo>::iterator it2=matching_Infos.begin();
390 it2!=it; it2++) {
391 if ( (it->genMuon == it2->genMuon) || (it->recoMuon == it2->recoMuon) ) {
392 matching_Infos.erase(it);
393 it=matching_Infos.begin();
394 break;
395 }
396 }
397 }
398 }
399
400 // Now put result in vector of pairs
401
402 leptonMatching[pdgid].clear();
403
404 for (vector<MatchingInfo>::iterator match=matching_Infos.begin();
405 match !=matching_Infos.end(); match++) { // po tablici dr,i,j
406
407 pair<const GenParticleCandidate*, const reco::Candidate *> pair;
408 // vector<pair,float > matchedParticles;
409 pair.first = genparticles[match->genMuon];
410 pair.second = & (*cands)[match->recoMuon];
411
412 matchedParticles cestice(genparticles[match->genMuon] ,& (*cands)[match->recoMuon] , match->deltaR);
413
414 leptonMatching[pdgid].push_back(cestice);
415
416 }
417
418 }
419
420
421 // get mother of particle
422 // returns mother = 1 if mother is Z boson
423 // returns mother = 2 if mother is W boson
424 // returns mother = 4 if mother is b
425 // returns mother = 0 else
426
427 MuonEfficiencyAnalyzer::LeptonOrigin MuonEfficiencyAnalyzer::getMother(const reco::Candidate* genParticle){
428
429 int pdg_id=genParticle->pdgId();
430 // cout << "particle id : " << pdg_id << endl;
431
432 while (genParticle = genParticle->mother() ) { //izvrsava se sve dok ne dodje do nule tj. dok ne trazi majku od protona
433 // cout << "Going up: ";
434 // cout << "Mother id : " << genParticle->pdgId() << endl;
435 if (abs(genParticle->pdgId())!=abs(pdg_id)) {
436 // cout<< " good mother " <<endl;
437 if (abs(genParticle->pdgId())==23){ // mother is Z
438 return zboson;
439 }
440 if (abs(genParticle->pdgId())==24) { // mother is W
441 return wboson;
442 }
443 // if (abs(genParticle->pdgId())==23 || abs(genParticle->pdgId())==24) { // mother is W or Z
444 // WZ_matched=1;
445 // mother=3;
446 // }
447 if ((((abs(genParticle->pdgId()))/100)%10 ==5)||(((abs(genParticle->pdgId()))/1000)%10 ==5 )) { // mother is b
448 return bdecay;
449 }
450 if ((((abs(genParticle->pdgId()))/100)%10 ==4)||(((abs(genParticle->pdgId()))/1000)%10 ==4 )) { // mother is c
451 return cdecay;
452 }
453 return other;
454 }
455 }
456
457 return other;
458 }
459
460 const reco::Candidate * MuonEfficiencyAnalyzer::MatchedGenParticle(const reco::Candidate * daughter){
461 ::OverlapChecker overlap;
462 matched_genParticle=0;
463
464 for (map< int,
465 std::vector< matchedParticles > > ::iterator iter = leptonMatching.begin();
466 iter!=leptonMatching.end(); iter++) // entire map
467 {
468 for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
469 const reco::Candidate * nesto=iter->second[j].pair.second;
470 if (overlap(*nesto,*daughter)){
471 matched_genParticle=iter->second[j].pair.first;
472 }
473 }
474 }
475 return matched_genParticle;
476 }
477
478
479 float MuonEfficiencyAnalyzer::dR(const reco::Candidate * daughter){
480
481 ::OverlapChecker overlap;
482
483 float delta_r = -10;
484 for (map< int,
485 std::vector< matchedParticles > > ::iterator iter = leptonMatching.begin();
486 iter!=leptonMatching.end(); iter++) // entire map
487 {
488 for (unsigned int j=0; j<iter->second.size(); j++){ // entire vector<class<>>
489 const reco::Candidate * nesto=iter->second[j].pair.second;
490 if (overlap(*nesto,*daughter)){
491 delta_r=iter->second[j].delta_R;
492 }
493 }
494 }
495 return delta_r;
496 }
497
498
499
500
501 //define this as a plug-in
502 // DEFINE_FWK_MODULE(MuonEfficiencyAnalyzer);
503