ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/NonMCBackground/src/select_fakes.cc
Revision: 1.11
Committed: Fri Jul 13 07:44:52 2012 UTC (12 years, 10 months ago) by dkralph
Content type: text/plain
Branch: MAIN
Changes since 1.10: +6 -0 lines
Log Message:
*** empty log message ***

File Contents

# User Rev Content
1 dkralph 1.2 //
2     // System headers
3     //
4     #include <vector> // STL vector class
5     #include <iostream> // standard I/O
6     #include <iomanip> // functions to format standard I/O
7     #include <bitset>
8     #include <map>
9     using namespace std;
10    
11     //
12     // root headers
13     //
14     #include <TFile.h>
15     #include <TTree.h>
16     #include <TChain.h>
17     #include <TBranch.h>
18     #include <TClonesArray.h>
19    
20     //
21     // TMVA
22     //
23     #include "TMVA/Reader.h"
24     #include "TMVA/Tools.h"
25     #include "TMVA/Config.h"
26     #include "TMVA/MethodBDT.h"
27    
28     //
29     // ntuple format headers
30     //
31     #include "EventHeader.h"
32     #include "PFMet.h"
33     #include "Electron.h"
34     #include "Muon.h"
35     #include "Vertex.h"
36     #include "PFCandidate.h"
37     #include "PFCandidateCol.h"
38     #include "PileupInfoCol.h"
39     #include "PileupEnergyDensity.h"
40     #include "MCParticle.h"
41     #include "MCEventInfo.h"
42     #include "TriggerMask.h"
43     #include "TriggerTable.h"
44     #include "Names.h"
45     #include "BaseMod.h"
46    
47     //
48     // our headers
49     //
50     #include "ParseArgs.h"
51     #include "MuonSelection.h"
52     #include "ElectronSelection.h"
53     #include "IsolationSelection.h"
54 dkralph 1.3 #include "SelectionFuncs.h"
55 dkralph 1.9 #include "CommonDefs.h"
56     #include "Angles.h"
57     #include "AngleTuple.h"
58     #include "FOTuple.h"
59 dkralph 1.2 #include "TriggerUtils.h"
60     #include "PassHLT.h"
61     #include "RunLumiRangeMap.h"
62 dkralph 1.3 #include "fake_defs.h"
63 dkralph 1.5 #include "Various.h"
64 dkralph 1.2 #include "CommonDefs.h"
65    
66     #ifndef CMSSW_BASE
67     #define CMSSW_BASE "../"
68     #endif
69    
70     using namespace mithep;
71     //
72     // globals
73     //----------------------------------------------------------------------------
74 khahn 1.6 TH1D *hpu_2011, *hpu_2012;
75 dkralph 1.2 mithep::RunLumiRangeMap rlrm;
76     vector<SimpleLepton> failingLeptons ; // for fake estimation
77     vector<SimpleLepton> passingLeptons; // for fake estimation
78     vector<unsigned> cutvec;
79     vector<vector<unsigned> > zcutvec;
80     vector<vector<unsigned> > zzcutvec;
81     map<unsigned,float> evtrhoMap;
82     vector<string> cutstrs;
83     bool passes_HLT_MC;
84     vector<bool> PFnoPUflag;;
85    
86     //
87     // prototypes
88     //----------------------------------------------------------------------------
89 dkralph 1.9 SelectionStatus analysisSelection(ControlFlags &ctrl, const mithep::Muon *mu, const mithep::Vertex *vtx, mithep::Array<mithep::PFCandidate> *pfArr,
90     mithep::Array<mithep::PileupEnergyDensity> *puDArr);
91     SelectionStatus analysisSelection(ControlFlags &ctrl, const mithep::Electron *ele, const mithep::Vertex *vtx, mithep::Array<mithep::PFCandidate> *pfArr,
92     mithep::Array<mithep::PileupEnergyDensity> *puDArr);
93 dkralph 1.7 void initRunLumiRangeMap(ControlFlags &ctrl);
94 dkralph 1.9 float getMass(const mithep::Muon *zmu1, const mithep::Muon *zmu2, const mithep::Electron *zele1, const mithep::Electron *zele2);
95 dkralph 1.2 mithep::MuonTools::EMuonEffectiveAreaTarget eraMu;
96     mithep::ElectronTools::EElectronEffectiveAreaTarget eraEle;
97     vector<const mithep::PFCandidate*> photonsToVeto;
98 dkralph 1.9 vector<SimpleLepton> findFakes(ControlFlags &ctrl,
99     mithep::Array<mithep::Muon> *muonArr,
100     mithep::Array<mithep::Electron> *electronArr,
101     mithep::Array<mithep::PFJet> *jetArr,
102     mithep::Array<mithep::PFCandidate> *pfArr,
103     mithep::Array<mithep::PileupInfo> *puArr,
104     mithep::Array<mithep::PileupEnergyDensity> *puDArr,
105     mithep::Array<mithep::Vertex> *vtxArr,
106     const mithep::Vertex *vtx,
107     const mithep::Muon *zmu1,
108     const mithep::Muon *zmu2,
109     const mithep::Electron *zele1,
110     const mithep::Electron *zele2);
111     SimpleLepton fillSimpleLepton(const Muon *mu, const Electron *ele, SelectionStatus *ctrl);
112 dkralph 1.1 //=== MAIN =================================================================================================
113     int main(int argc, char** argv)
114     {
115 dkralph 1.2
116     string cmsswpath(CMSSW_BASE);
117     cmsswpath.append("/src");
118     std::bitset<1024> triggerBits;
119 dkralph 1.1 vector<vector<string> > inputFiles;
120     inputFiles.push_back(vector<string>());
121 dkralph 1.2 map<string,unsigned> entrymap; // number of unskimmed entries in each file
122 dkralph 1.1
123     //
124     // args
125     //--------------------------------------------------------------------------------------------------------------
126     ControlFlags ctrl;
127     parse_args( argc, argv, ctrl );
128     ctrl.dump();
129    
130 dkralph 1.9 // // hack: store low and high mass values in mcfmfname
131     // float massLo,massHi;
132     // TString massWindow(ctrl.mcfmfname);
133     // massWindow.ReplaceAll("..."," ");
134     // string line(massWindow);
135     // stringstream ss(line);
136     // ss >> massLo >> massHi;
137     // assert(massLo>=0 && massLo <= 8000);
138     // assert(massHi>=0 && massHi <= 8000);
139    
140 dkralph 1.1 //
141     // File I/O
142     //--------------------------------------------------------------------------------------------------------------
143     TChain * chain = new TChain("Events");
144 dkralph 1.2 TChain * hltchain = new TChain("HLT");
145    
146 dkralph 1.1 string fname;
147 dkralph 1.2 unsigned total_unskimmed=0;
148     if( !ctrl.inputfiles.empty() ) {
149     ifstream f(ctrl.inputfiles.c_str());
150     while (f >> fname) {
151     if( !(strncmp( fname.c_str(), "#", 1 ) ) ) continue; // skip commented lines
152     cout << "adding inputfile : " << fname.c_str() << endl;
153     entrymap[string(fname.c_str())] = unskimmedEntries(fname.c_str());
154     total_unskimmed += entrymap[string(fname.c_str())];
155     chain->AddFile(fname.c_str());
156     hltchain->AddFile(fname.c_str());
157     }
158     } else {
159     cout << "adding inputfile : " << ctrl.inputfile.c_str() << endl;
160     unsigned unsk_ents = unskimmedEntries(ctrl.inputfile.c_str());
161     entrymap[string(ctrl.inputfile.c_str())] = unsk_ents;
162     total_unskimmed += unsk_ents;
163     chain->AddFile(ctrl.inputfile.c_str());
164     hltchain->AddFile(ctrl.inputfile.c_str());
165 dkralph 1.1 }
166 dkralph 1.2 // write the total number of unskimmed events that went into making this output file to a text file
167     writeEntries(ctrl,total_unskimmed);
168 dkralph 1.1
169 dkralph 1.9 UInt_t nEvents=0,nMuDenom=0,nMuNumer=0,nEleDenom=0,nEleNumer=0;
170 dkralph 1.1
171 dkralph 1.9 // move this to compute_fakes
172     // double metMax = 25;
173 dkralph 1.1
174     //
175     // Setup
176     //--------------------------------------------------------------------------------------------------------------
177     const char * ofname;
178     if( strcmp( ctrl.outputfile.c_str(), "") ) {
179     ofname = ctrl.outputfile.c_str();
180     } else {
181     ofname = "fo.root";
182     }
183    
184 dkralph 1.9 // TFile *outFile = new TFile(ofname, "RECREATE");
185    
186     Angles angles;
187     KinematicsStruct kinematics;
188     InfoStruct evtinfo;
189     // GenInfoStruct geninfo;
190    
191     AngleTuple nt( (const char*)ofname, (const char*)"zznt");
192     nt.makeAngleBranch(angles);
193     nt.makeKinematicsBranch(kinematics);
194     nt.makeInfoBranch(evtinfo);
195    
196     vector<SimpleLepton> failingLeptons,passingLeptons;
197     FOTuple foTree( nt.getFile(), (const char*)"FOtree");
198     foTree.makeFailedLeptonBranch(failingLeptons);
199     foTree.makePassedLeptonBranch(passingLeptons);
200    
201     NtuplerBase foFlatTree( nt.getFile(), (const char*)"FO");
202     const char * brname = "FO";
203     const char * varlist = VAR_LIST_FOINFO;
204     foinfo foi;
205     foFlatTree.makeBranch( brname, (void*) &foi, varlist);
206 dkralph 1.2
207     // initMuonIDMVA();
208 dkralph 1.9 // getEATargets(ctrl,eraMu,eraEle);
209 dkralph 1.2 initElectronIDMVAV1();
210     initMuonIsoMVA();
211     initElectronIDMVA();
212     initElectronIsoMVA();
213     initTrigger();
214     if(!ctrl.mc)
215 dkralph 1.7 initRunLumiRangeMap(ctrl);
216 dkralph 1.1
217     //
218     // Access samples and fill histograms
219     TFile *inputFile=0;
220     TTree *eventTree=0;
221 dkralph 1.2 string currentFile("");
222 dkralph 1.1
223 dkralph 1.2 mithep::EventHeader *info = new mithep::EventHeader();
224     mithep::Array<mithep::PFMet> *metArr = new mithep::Array<mithep::PFMet>();
225     mithep::Array<mithep::Electron> *electronArr = new mithep::Array<mithep::Electron>();
226     mithep::Array<mithep::Muon> *muonArr = new mithep::Array<mithep::Muon>();
227     mithep::Array<mithep::Vertex> *vtxArr = new mithep::Array<mithep::Vertex>();
228     mithep::Array<mithep::PFCandidate> *pfArr = new mithep::Array<mithep::PFCandidate>();
229     mithep::Array<mithep::PFJet> *jetArr = new mithep::Array<mithep::PFJet>();
230     mithep::Array<mithep::PileupInfo> *puArr = new mithep::Array<mithep::PileupInfo>();
231     mithep::Array<mithep::PileupEnergyDensity> *puDArr = new mithep::Array<mithep::PileupEnergyDensity>();
232     mithep::Array<mithep::Track> *trkArr = new mithep::Array<mithep::Track>();
233     mithep::Array<mithep::MCParticle> *mcArr = new mithep::Array<mithep::MCParticle>();
234     mithep::MCEventInfo *mcEvtInfo = new mithep::MCEventInfo();
235     mithep::TriggerMask *trigMask = new mithep::TriggerMask();
236     mithep::TriggerTable *hltTable = new mithep::TriggerTable();
237     vector<string> *hltTableStrings = new vector<string>();
238    
239     TString fElectronName(Names::gkElectronBrn);
240     TString fMuonName(Names::gkMuonBrn);
241     TString fInfoName(Names::gkEvtHeaderBrn);
242     TString fPFMetName("PFMet");
243     TString fPrimVtxName(Names::gkPVBrn);
244     TString fPFCandidateName(Names::gkPFCandidatesBrn);
245     TString fPFJetName(Names::gkPFJetBrn);
246     TString fPileupInfoName(Names::gkPileupInfoBrn);
247     TString fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn);
248     TString fTrackName(Names::gkTrackBrn);
249     TString fMCParticleName(Names::gkMCPartBrn);
250     TString fMCEvtInfoName(Names::gkMCEvtInfoBrn);
251     TString fTriggerMaskName(Names::gkHltBitBrn);
252     TString fTriggerTableName(Names::gkHltTableBrn);
253    
254     chain->SetBranchAddress(fInfoName, &info);
255     chain->SetBranchAddress(fPFMetName, &metArr);
256     chain->SetBranchAddress(fElectronName, &electronArr);
257     chain->SetBranchAddress(fMuonName, &muonArr);
258     chain->SetBranchAddress(fPrimVtxName, &vtxArr);
259     chain->SetBranchAddress(fPFCandidateName, &pfArr);
260     chain->SetBranchAddress(fPFJetName, &jetArr);
261     if( ctrl.mc ) {
262     chain->SetBranchAddress(fPileupInfoName, &puArr);
263     chain->SetBranchAddress(fMCParticleName, &mcArr);
264     chain->SetBranchAddress(fMCEvtInfoName, &mcEvtInfo);
265     }
266     chain->SetBranchAddress(fPileupEnergyDensityName, &puDArr);
267     chain->SetBranchAddress(fTrackName, &trkArr);
268     chain->SetBranchAddress(fTriggerMaskName, &trigMask);
269     cout << "hlttable: " << fTriggerTableName << endl;
270     hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
271    
272     const mithep::Vertex *vtx; // best primary vertex in the event
273 dkralph 1.1
274 dkralph 1.2 // ginfo = NULL;
275     // if( ctrl.mc ) { chain->SetBranchAddress("Gen", &ginfo);}
276    
277     // only 1 HLT table / file ???
278     hltchain->GetEntry(0);
279    
280     int imax = chain->GetEntries();
281     cout << "nEntries: " << imax << endl;
282 dkralph 1.1 for(UInt_t ientry=0; ientry<imax; ientry++) {
283     chain->GetEntry(ientry);
284 dkralph 1.2 if(!(ientry%10000)) cout << "entry: " << ientry << "\tfrac: " << setw(15) << float(ientry)/chain->GetEntries() << endl;
285    
286 dkralph 1.9 string fname = string(chain->GetFile()->GetEndpointUrl()->GetFile());
287     if( ctrl.era == 0 )
288     setEra( fname, ctrl );
289     getEATargets(ctrl,eraMu,eraEle);
290    
291 dkralph 1.2 // pfNoPU
292     PFnoPUflag.clear();
293     int pfnopu_size = makePFnoPUArray( pfArr, PFnoPUflag, vtxArr );
294     assert(pfnopu_size == pfArr->GetEntries());
295    
296 dkralph 1.7 if(!(ctrl.noJSON) ) {
297     RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());
298     if( !(rlrm.HasRunLumi(rl)) ) continue;
299     }
300 dkralph 1.2
301 dkralph 1.9 nEvents++;
302 dkralph 1.1
303 dkralph 1.10 bool goodVertex = setPV( ctrl, vtxArr, vtx );
304     if(!goodVertex) continue;
305 dkralph 1.1
306 dkralph 1.9 // move this to compute_fakes
307     // TLorentzVector pfmet; pfmet.SetPxPyPzE(metArr->At(0)->Mex(),metArr->At(0)->Mey(),0,0);
308     // if(pfmet.Pt() > metMax) continue;
309 dkralph 1.1
310 dkralph 1.2 // good lepton selection
311 dkralph 1.9 failingLeptons.clear();
312     passingLeptons.clear();
313 dkralph 1.2 vector<const mithep::Muon*> goodMuons;
314     vector<const mithep::Electron*> goodElectrons;
315 dkralph 1.9 vector<SelectionStatus> muselv,eleselv;
316 dkralph 1.2 for(int imu=0; imu<muonArr->GetEntries(); imu++) {
317     const mithep::Muon *mu = (*muonArr)[imu];
318 dkralph 1.11 SelectionStatus denomsel;
319     denomsel |= muonPreSelectionNoD0DzIP(ctrl,mu,vtx,pfArr);;
320     if( !denomsel.passPre() ) continue;
321 dkralph 1.9 SelectionStatus status = analysisSelection(ctrl,mu,vtx,pfArr,puDArr);
322     if(status.loose()) {
323     goodMuons.push_back(mu);
324     muselv.push_back(status);
325     passingLeptons.push_back(fillSimpleLepton(mu,&status));
326     } else {
327     failingLeptons.push_back(fillSimpleLepton(mu,&status));
328     }
329 dkralph 1.2 }
330 dkralph 1.1 for(Int_t i=0; i<electronArr->GetEntries(); i++) {
331 dkralph 1.2 const mithep::Electron *ele = (*electronArr)[i];
332 dkralph 1.4 bool muonMatch=false;
333     for(unsigned imu=0; imu<muonArr->GetEntries(); imu++) {
334     const mithep::Muon *mu = (*muonArr)[imu];
335     if(dr(mu,ele) < 0.05) muonMatch = true;
336     }
337 dkralph 1.9 if(muonMatch) continue;
338 dkralph 1.11 SelectionStatus denomsel;
339     denomsel |= electronPreSelectionNoD0DzIP(ctrl,ele,vtx);
340     if( !denomsel.passPre() ) continue;
341 dkralph 1.9 SelectionStatus status = analysisSelection(ctrl,ele,vtx,pfArr,puDArr);
342     if(status.loose()) {
343     goodElectrons.push_back(ele);
344     eleselv.push_back(status);
345     passingLeptons.push_back(fillSimpleLepton(ele,&status));
346     } else {
347     failingLeptons.push_back(fillSimpleLepton(ele,&status));
348     }
349 dkralph 1.1 }
350    
351 dkralph 1.2 // look for z pairs
352     int NZCandidates=0;
353 dkralph 1.9 float bestMass=-999;
354 dkralph 1.2 const mithep::Muon *zmu1=0,*zmu2=0;
355     const mithep::Electron *zele1=0,*zele2=0;
356 dkralph 1.9 SelectionStatus *status1,*status2;
357 dkralph 1.2 for(unsigned imu=0; imu<goodMuons.size(); imu++) {
358     for(unsigned jmu=imu+1; jmu<goodMuons.size(); jmu++) {
359     if(goodMuons[imu]->Charge() == goodMuons[jmu]->Charge()) continue;
360 dkralph 1.9 float mass = getMass(goodMuons[imu],goodMuons[jmu],0,0);
361     if( fabs(mass - Z_MASS) < fabs(bestMass - Z_MASS)) {
362     bestMass = mass;
363     zmu1 = goodMuons[imu]; status1 = &muselv[imu];
364     zmu2 = goodMuons[jmu]; status2 = &muselv[jmu];
365     }
366 dkralph 1.2 }
367     }
368     for(unsigned iele=0; iele<goodElectrons.size(); iele++) {
369     for(unsigned jele=iele+1; jele<goodElectrons.size(); jele++) {
370     if(goodElectrons[iele]->Charge() == goodElectrons[jele]->Charge()) continue;
371 dkralph 1.9 float mass = getMass(0,0,goodElectrons[iele],goodElectrons[jele]);
372     if( fabs(mass - Z_MASS) < fabs(bestMass - Z_MASS)) {
373     bestMass = mass;
374     zele1 = goodElectrons[iele]; status1 = &eleselv[iele];
375     zele2 = goodElectrons[jele]; status2 = &eleselv[jele];
376     zmu1 = zmu2 = 0;
377     }
378 dkralph 1.1 }
379 dkralph 1.2 }
380 dkralph 1.9 if(bestMass < 0) continue;
381     assert( (zmu1 && zmu2) || (zele1 && zele2));
382     assert( !(zmu1!=0 && zele1!=0) );
383     assert( !(zmu2!=0 && zele2!=0) );
384    
385     EventData data;
386     data.Z1leptons.push_back(fillSimpleLepton(zmu1,zele1,status1));
387     data.Z1leptons.push_back(fillSimpleLepton(zmu2,zele2,status2));
388     data.Z2leptons.push_back(fillSimpleLepton(zmu1,zele1,status1));
389     data.Z2leptons.push_back(fillSimpleLepton(zmu2,zele2,status2));
390    
391     // loop over fakes
392     vector<SimpleLepton> fakes = findFakes(ctrl,muonArr,electronArr,jetArr,pfArr,puArr,puDArr,vtxArr,vtx,zmu1,zmu2,zele1,zele2);
393 dkralph 1.10 //????????????????????????????????????????????????????????????????????????????????????????
394     // if(fakes.size() > 1) continue;
395     //????????????????????????????????????????????????????????????????????????????????????????
396 dkralph 1.9 for(unsigned ifake=0; ifake<fakes.size(); ifake++) {
397     SimpleLepton fake = fakes[ifake];
398     if(zmu1) {
399     if(dr(fake,zmu1) < 0.3) continue;
400     if(dr(fake,zmu2) < 0.3) continue;
401     }
402     if(zele1) {
403     if(dr(fake,zele1) < 0.3) continue;
404     if(dr(fake,zele2) < 0.3) continue;
405     }
406 dkralph 1.2
407 dkralph 1.9 if(abs(fake.type) == 13) {
408     nMuDenom++;
409     if(fake.isLoose) nMuNumer++;
410     } else {
411     nEleDenom++;
412     if(fake.isLoose) nEleNumer++;
413     }
414 dkralph 1.1
415 dkralph 1.9 // fill output TTree
416 dkralph 1.10 fillKinematics(data,kinematics);
417     TLorentzVector pfmet; pfmet.SetPxPyPzE(metArr->At(0)->Mex(),metArr->At(0)->Mey(),0,0);
418     fillEventInfo( info, pfmet, evtinfo, ctrl.mc ? getNPU(puArr) : 0, vtxArr->GetEntries());
419     foTree.Fill();
420     nt.Fill();
421    
422 dkralph 1.9 foi.run = info->RunNum();
423     foi.evt = info->EvtNum();
424     foi.lumi = info->LumiSec();
425     foi.type = fake.type;
426     foi.npu = ctrl.mc ? getNPU(puArr) : 0;
427     foi.npv = vtxArr->GetEntries(); // note: should maybe apply some selection here
428     foi.fopt = fake.vec.Pt();
429     foi.foeta = fake.vec.Eta();
430     foi.fophi = fake.vec.Phi();
431     foi.fopass = fake.isLoose ? 1 : 0;
432     foi.mass = getMass(zmu1,zmu2,zele1,zele2);
433     foi.l1pt = zmu1 ? zmu1->Pt() : zele1->Pt();
434     foi.l2pt = zmu2 ? zmu2->Pt() : zele2->Pt();
435     foi.met = pfmet.Pt();
436     foFlatTree.Fill();
437     }
438 dkralph 1.2 }
439    
440 dkralph 1.9 // outFile->Write();
441     // outFile->Close();
442     foTree.getFile()->cd();
443     foTree.getTree()->Write();
444     foFlatTree.getFile()->cd();
445     foFlatTree.getTree()->Write();
446     nt.WriteClose();
447 dkralph 1.2
448 dkralph 1.1 //--------------------------------------------------------------------------------------------------------------
449     // Summary print out
450     //==============================================================================================================
451     cout << endl;
452     cout << "*" << endl;
453     cout << "* SUMMARY" << endl;
454     cout << "*--------------------------------------------------" << endl;
455     cout << endl;
456    
457 dkralph 1.9 cout << " >>> No. of events processed: " << nEvents << endl;
458     cout << " >>> No. of denominator muons: " << nMuDenom << endl;
459     cout << " >>> No. of numerator muons: " << nMuNumer << endl;
460     cout << " >>> ratio: " << float(nMuNumer)/nMuDenom << endl;
461     cout << " >>> No. of denominator eles: " << nEleDenom << endl;
462     cout << " >>> No. of numerator eles: " << nEleNumer << endl;
463     cout << " >>> ratio: " << float(nEleNumer)/nEleDenom << endl;
464 dkralph 1.1 cout << endl;
465    
466 dkralph 1.2 TString outputDir(ofname);
467     outputDir = outputDir(0,outputDir.Last('/'));
468 dkralph 1.1 ofstream txtfile;
469     char txtfname[100];
470 dkralph 1.2 TString txtname(ctrl.outputfile);
471     txtname.ReplaceAll(".root",".txt");
472     cout << "writing text to: " << txtname << endl;
473     txtfile.open(txtname);
474 dkralph 1.1 txtfile << "*" << endl;
475     txtfile << "* SUMMARY" << endl;
476     txtfile << "*--------------------------------------------------" << endl;
477     txtfile << endl;
478    
479 dkralph 1.9 txtfile << " >>> No. of events processed: " << nEvents << endl;
480     txtfile << " >>> No. of denominator muons: " << nMuDenom << endl;
481     txtfile << " >>> No. of numerator muons: " << nMuNumer << endl;
482     txtfile << " >>> ratio: " << float(nMuNumer)/nMuDenom << endl;
483     txtfile << " >>> No. of denominator eles: " << nEleDenom << endl;
484     txtfile << " >>> No. of numerator eles: " << nEleNumer << endl;
485     txtfile << " >>> ratio: " << float(nEleNumer)/nEleDenom << endl;
486 dkralph 1.1 txtfile << endl;
487     txtfile.close();
488    
489     cout << " <> Output saved in " << outputDir << "/" << endl;
490     cout << endl;
491    
492     cout << "done!" << endl;
493 dkralph 1.2 }
494 dkralph 1.9 //----------------------------------------------------------------------------
495     void initRunLumiRangeMap(ControlFlags &ctrl)
496     //----------------------------------------------------------------------------
497     {
498     /*
499     rlrm.AddJSONFile(std::string("./data/Cert_136033-149442_7TeV_Apr21ReReco_Collisions10_JSON.txt"));
500     // rlrm.AddJSONFile(std::string("./data/Cert_160404-173244_7TeV_PromptReco_Collisions11_JSON_v2.txt"));
501     rlrm.AddJSONFile(std::string("./data/Cert_160404-178078_7TeV_PromptReco_Collisions11_JSON.txt"));
502     rlrm.AddJSONFile(std::string("./data/Cert_160404-163869_7TeV_May10ReReco_Collisions11_JSON_v3.txt"));
503     rlrm.AddJSONFile(std::string("./data/Cert_170249-172619_7TeV_ReReco5Aug_Collisions11_JSON.txt"));
504     // r11b
505     rlrm.AddJSONFile(std::string("./data/Cert_160404-180252_7TeV_PromptReco_Collisions11_JSON.txt"));
506     */
507    
508     // 2012 only for now ...
509     rlrm.AddJSONFile(string(ctrl.jsonFile));
510    
511     };
512     //----------------------------------------------------------------------------------------
513     float getMass(const mithep::Muon *zmu1, const mithep::Muon *zmu2, const mithep::Electron *zele1, const mithep::Electron *zele2)
514     {
515     TLorentzVector lep1,lep2;
516     if(zmu1 && zmu2) {
517     lep1.SetPtEtaPhiM(zmu1->Pt(),zmu1->Eta(),zmu1->Phi(),MUON_MASS);
518     lep2.SetPtEtaPhiM(zmu2->Pt(),zmu2->Eta(),zmu2->Phi(),MUON_MASS);
519     } else if(zele1 && zele2) {
520     lep1.SetPtEtaPhiM(zele1->Pt(),zele1->Eta(),zele1->Phi(),ELECTRON_MASS);
521     lep2.SetPtEtaPhiM(zele2->Pt(),zele2->Eta(),zele2->Phi(),ELECTRON_MASS);
522     } else assert(0);
523    
524     return (lep1 + lep2).M();
525     }
526     //----------------------------------------------------------------------------------------
527     SelectionStatus analysisSelection(ControlFlags &ctrl, const mithep::Muon *mu, const mithep::Vertex *vtx, mithep::Array<mithep::PFCandidate> *pfArr,
528     mithep::Array<mithep::PileupEnergyDensity> *puDArr)
529     {
530     SelectionStatus status;
531     status |= muonReferencePreSelection(ctrl,mu,vtx,pfArr);
532     status |= muonIDPFSelection(ctrl,mu,vtx,pfArr);
533     status |= muonReferenceIsoSelection(ctrl,mu,vtx,pfArr,puDArr,eraMu,photonsToVeto);
534     return status;
535     }
536 dkralph 1.2 //----------------------------------------------------------------------------------------
537 dkralph 1.9 SelectionStatus analysisSelection(ControlFlags &ctrl, const mithep::Electron *ele, const mithep::Vertex *vtx, mithep::Array<mithep::PFCandidate> *pfArr,
538     mithep::Array<mithep::PileupEnergyDensity> *puDArr)
539     {
540     SelectionStatus status;
541     status |= electronReferencePreSelection(ctrl,ele,vtx);
542     status |= electronReferenceIDMVASelectionV1(ctrl,ele,vtx);
543     status |= electronReferenceIsoSelection(ctrl,ele,vtx,pfArr,puDArr,eraEle,photonsToVeto);
544     return status;
545     }
546     //----------------------------------------------------------------------------------------
547     vector<SimpleLepton> findFakes(ControlFlags &ctrl,
548     mithep::Array<mithep::Muon> *muonArr,
549     mithep::Array<mithep::Electron> *electronArr,
550     mithep::Array<mithep::PFJet> *jetArr,
551     mithep::Array<mithep::PFCandidate> *pfArr,
552     mithep::Array<mithep::PileupInfo> *puArr,
553     mithep::Array<mithep::PileupEnergyDensity> *puDArr,
554     mithep::Array<mithep::Vertex> *vtxArr,
555     const mithep::Vertex *vtx,
556     const mithep::Muon *zmu1,
557     const mithep::Muon *zmu2,
558     const mithep::Electron *zele1,
559     const mithep::Electron *zele2)
560     {
561     vector<SimpleLepton> fakes;
562    
563 dkralph 1.2 for(Int_t imu=0; imu<muonArr->GetEntries(); imu++) {
564     const mithep::Muon* mu = (*muonArr)[imu];
565    
566 dkralph 1.9 if(mu==zmu1 || mu==zmu2) continue;
567 dkralph 1.2
568 dkralph 1.7 SelectionStatus denomsel;
569 dkralph 1.8 denomsel |= muonPreSelectionNoD0DzIP(ctrl,mu,vtx,pfArr);;
570     if( !denomsel.passPre() ) continue;
571 dkralph 1.2
572 dkralph 1.9 SelectionStatus status = analysisSelection(ctrl,mu,vtx,pfArr,puDArr);
573 dkralph 1.2
574 dkralph 1.9 SimpleLepton fake;
575     fake.vec.SetPtEtaPhiM(mu->Pt(),mu->Eta(),mu->Phi(),MUON_MASS);
576     fake.type = 13;
577     fake.isLoose = status.loose();
578    
579     fakes.push_back(fake);
580 dkralph 1.2 }
581     for(Int_t iele=0; iele<electronArr->GetEntries(); iele++) {
582     const mithep::Electron* ele = (*electronArr)[iele];
583    
584 dkralph 1.4 bool muonMatch=false;
585     for(unsigned imu=0; imu<muonArr->GetEntries(); imu++) {
586     const mithep::Muon *mu = (*muonArr)[imu];
587     if(dr(mu,ele) < 0.05) muonMatch = true;
588     }
589 dkralph 1.9 if(muonMatch) continue;
590 dkralph 1.4
591 dkralph 1.9 if(ele==zele1 || ele==zele2) continue;
592 dkralph 1.2
593 dkralph 1.7 SelectionStatus denomsel;
594 dkralph 1.8 denomsel |= electronPreSelectionNoD0DzIP(ctrl,ele,vtx);
595     if( !denomsel.passPre() ) continue;
596 dkralph 1.2
597 dkralph 1.9 SelectionStatus status = analysisSelection(ctrl,ele,vtx,pfArr,puDArr);
598    
599     SimpleLepton fake;
600     fake.vec.SetPtEtaPhiM(ele->Pt(),ele->Eta(),ele->Phi(),ELECTRON_MASS);
601     fake.type = 11;
602     fake.isLoose = status.loose();
603 dkralph 1.8
604 dkralph 1.9 fakes.push_back(fake);
605     }
606 dkralph 1.2
607 dkralph 1.9 return fakes;
608 dkralph 1.1 }
609 dkralph 1.9 //----------------------------------------------------------------------------------------
610     SimpleLepton fillSimpleLepton(const Muon *mu, const Electron *ele, SelectionStatus *ctrl)
611 dkralph 1.2 {
612 dkralph 1.9 if(mu) return fillSimpleLepton(mu, ctrl);
613     else if(ele) return fillSimpleLepton(ele,ctrl);
614     else assert(0);
615     return SimpleLepton();
616     }