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
Log Message:
muon eff. analysis

File Contents

# User Rev Content
1 vuko 1.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