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 |
|