ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/NonMCBackground/src/select_fakes.cc
Revision: 1.9
Committed: Fri Jul 6 15:01:29 2012 UTC (12 years, 10 months ago) by dkralph
Content type: text/plain
Branch: MAIN
Changes since 1.8: +271 -242 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.2 setPV( ctrl, vtxArr, vtx);
304 dkralph 1.1
305 dkralph 1.9 // move this to compute_fakes
306     // TLorentzVector pfmet; pfmet.SetPxPyPzE(metArr->At(0)->Mex(),metArr->At(0)->Mey(),0,0);
307     // if(pfmet.Pt() > metMax) continue;
308 dkralph 1.1
309 dkralph 1.2 // good lepton selection
310 dkralph 1.9 failingLeptons.clear();
311     passingLeptons.clear();
312 dkralph 1.2 vector<const mithep::Muon*> goodMuons;
313     vector<const mithep::Electron*> goodElectrons;
314 dkralph 1.9 vector<SelectionStatus> muselv,eleselv;
315 dkralph 1.2 for(int imu=0; imu<muonArr->GetEntries(); imu++) {
316     const mithep::Muon *mu = (*muonArr)[imu];
317 dkralph 1.9 SelectionStatus status = analysisSelection(ctrl,mu,vtx,pfArr,puDArr);
318     if(status.loose()) {
319     goodMuons.push_back(mu);
320     muselv.push_back(status);
321     passingLeptons.push_back(fillSimpleLepton(mu,&status));
322     } else {
323     failingLeptons.push_back(fillSimpleLepton(mu,&status));
324     }
325 dkralph 1.2 }
326 dkralph 1.1 for(Int_t i=0; i<electronArr->GetEntries(); i++) {
327 dkralph 1.2 const mithep::Electron *ele = (*electronArr)[i];
328 dkralph 1.4 bool muonMatch=false;
329     for(unsigned imu=0; imu<muonArr->GetEntries(); imu++) {
330     const mithep::Muon *mu = (*muonArr)[imu];
331     if(dr(mu,ele) < 0.05) muonMatch = true;
332     }
333 dkralph 1.9 if(muonMatch) continue;
334     SelectionStatus status = analysisSelection(ctrl,ele,vtx,pfArr,puDArr);
335     if(status.loose()) {
336     goodElectrons.push_back(ele);
337     eleselv.push_back(status);
338     passingLeptons.push_back(fillSimpleLepton(ele,&status));
339     } else {
340     failingLeptons.push_back(fillSimpleLepton(ele,&status));
341     }
342 dkralph 1.1 }
343    
344 dkralph 1.2 // look for z pairs
345     int NZCandidates=0;
346 dkralph 1.9 float bestMass=-999;
347 dkralph 1.2 const mithep::Muon *zmu1=0,*zmu2=0;
348     const mithep::Electron *zele1=0,*zele2=0;
349 dkralph 1.9 SelectionStatus *status1,*status2;
350 dkralph 1.2 for(unsigned imu=0; imu<goodMuons.size(); imu++) {
351     for(unsigned jmu=imu+1; jmu<goodMuons.size(); jmu++) {
352     if(goodMuons[imu]->Charge() == goodMuons[jmu]->Charge()) continue;
353 dkralph 1.9 float mass = getMass(goodMuons[imu],goodMuons[jmu],0,0);
354     if( fabs(mass - Z_MASS) < fabs(bestMass - Z_MASS)) {
355     bestMass = mass;
356     zmu1 = goodMuons[imu]; status1 = &muselv[imu];
357     zmu2 = goodMuons[jmu]; status2 = &muselv[jmu];
358     }
359 dkralph 1.2 }
360     }
361     for(unsigned iele=0; iele<goodElectrons.size(); iele++) {
362     for(unsigned jele=iele+1; jele<goodElectrons.size(); jele++) {
363     if(goodElectrons[iele]->Charge() == goodElectrons[jele]->Charge()) continue;
364 dkralph 1.9 float mass = getMass(0,0,goodElectrons[iele],goodElectrons[jele]);
365     if( fabs(mass - Z_MASS) < fabs(bestMass - Z_MASS)) {
366     bestMass = mass;
367     zele1 = goodElectrons[iele]; status1 = &eleselv[iele];
368     zele2 = goodElectrons[jele]; status2 = &eleselv[jele];
369     zmu1 = zmu2 = 0;
370     }
371 dkralph 1.1 }
372 dkralph 1.2 }
373 dkralph 1.9 if(bestMass < 0) continue;
374     assert( (zmu1 && zmu2) || (zele1 && zele2));
375     assert( !(zmu1!=0 && zele1!=0) );
376     assert( !(zmu2!=0 && zele2!=0) );
377    
378     EventData data;
379     data.Z1leptons.push_back(fillSimpleLepton(zmu1,zele1,status1));
380     data.Z1leptons.push_back(fillSimpleLepton(zmu2,zele2,status2));
381     data.Z2leptons.push_back(fillSimpleLepton(zmu1,zele1,status1));
382     data.Z2leptons.push_back(fillSimpleLepton(zmu2,zele2,status2));
383     fillKinematics(data,kinematics);
384     TLorentzVector pfmet; pfmet.SetPxPyPzE(metArr->At(0)->Mex(),metArr->At(0)->Mey(),0,0);
385     fillEventInfo( info, pfmet, evtinfo, ctrl.mc ? getNPU(puArr) : 0, vtxArr->GetEntries());
386     foTree.Fill();
387     nt.Fill();
388    
389     // loop over fakes
390     vector<SimpleLepton> fakes = findFakes(ctrl,muonArr,electronArr,jetArr,pfArr,puArr,puDArr,vtxArr,vtx,zmu1,zmu2,zele1,zele2);
391     if(fakes.size() > 1) continue;
392     for(unsigned ifake=0; ifake<fakes.size(); ifake++) {
393     SimpleLepton fake = fakes[ifake];
394     if(zmu1) {
395     if(dr(fake,zmu1) < 0.3) continue;
396     if(dr(fake,zmu2) < 0.3) continue;
397     }
398     if(zele1) {
399     if(dr(fake,zele1) < 0.3) continue;
400     if(dr(fake,zele2) < 0.3) continue;
401     }
402 dkralph 1.2
403 dkralph 1.9 if(abs(fake.type) == 13) {
404     nMuDenom++;
405     if(fake.isLoose) nMuNumer++;
406     } else {
407     nEleDenom++;
408     if(fake.isLoose) nEleNumer++;
409     }
410 dkralph 1.1
411 dkralph 1.9 // fill output TTree
412     foi.run = info->RunNum();
413     foi.evt = info->EvtNum();
414     foi.lumi = info->LumiSec();
415     foi.type = fake.type;
416     foi.npu = ctrl.mc ? getNPU(puArr) : 0;
417     foi.npv = vtxArr->GetEntries(); // note: should maybe apply some selection here
418     foi.fopt = fake.vec.Pt();
419     foi.foeta = fake.vec.Eta();
420     foi.fophi = fake.vec.Phi();
421     foi.fopass = fake.isLoose ? 1 : 0;
422     foi.mass = getMass(zmu1,zmu2,zele1,zele2);
423     foi.l1pt = zmu1 ? zmu1->Pt() : zele1->Pt();
424     foi.l2pt = zmu2 ? zmu2->Pt() : zele2->Pt();
425     foi.met = pfmet.Pt();
426     foFlatTree.Fill();
427     }
428 dkralph 1.2 }
429    
430 dkralph 1.9 // outFile->Write();
431     // outFile->Close();
432     foTree.getFile()->cd();
433     foTree.getTree()->Write();
434     foFlatTree.getFile()->cd();
435     foFlatTree.getTree()->Write();
436     nt.WriteClose();
437 dkralph 1.2
438 dkralph 1.1 //--------------------------------------------------------------------------------------------------------------
439     // Summary print out
440     //==============================================================================================================
441     cout << endl;
442     cout << "*" << endl;
443     cout << "* SUMMARY" << endl;
444     cout << "*--------------------------------------------------" << endl;
445     cout << endl;
446    
447 dkralph 1.9 cout << " >>> No. of events processed: " << nEvents << endl;
448     cout << " >>> No. of denominator muons: " << nMuDenom << endl;
449     cout << " >>> No. of numerator muons: " << nMuNumer << endl;
450     cout << " >>> ratio: " << float(nMuNumer)/nMuDenom << endl;
451     cout << " >>> No. of denominator eles: " << nEleDenom << endl;
452     cout << " >>> No. of numerator eles: " << nEleNumer << endl;
453     cout << " >>> ratio: " << float(nEleNumer)/nEleDenom << endl;
454 dkralph 1.1 cout << endl;
455    
456 dkralph 1.2 TString outputDir(ofname);
457     outputDir = outputDir(0,outputDir.Last('/'));
458 dkralph 1.1 ofstream txtfile;
459     char txtfname[100];
460 dkralph 1.2 TString txtname(ctrl.outputfile);
461     txtname.ReplaceAll(".root",".txt");
462     cout << "writing text to: " << txtname << endl;
463     txtfile.open(txtname);
464 dkralph 1.1 txtfile << "*" << endl;
465     txtfile << "* SUMMARY" << endl;
466     txtfile << "*--------------------------------------------------" << endl;
467     txtfile << endl;
468    
469 dkralph 1.9 txtfile << " >>> No. of events processed: " << nEvents << endl;
470     txtfile << " >>> No. of denominator muons: " << nMuDenom << endl;
471     txtfile << " >>> No. of numerator muons: " << nMuNumer << endl;
472     txtfile << " >>> ratio: " << float(nMuNumer)/nMuDenom << endl;
473     txtfile << " >>> No. of denominator eles: " << nEleDenom << endl;
474     txtfile << " >>> No. of numerator eles: " << nEleNumer << endl;
475     txtfile << " >>> ratio: " << float(nEleNumer)/nEleDenom << endl;
476 dkralph 1.1 txtfile << endl;
477     txtfile.close();
478    
479     cout << " <> Output saved in " << outputDir << "/" << endl;
480     cout << endl;
481    
482     cout << "done!" << endl;
483 dkralph 1.2 }
484 dkralph 1.9 //----------------------------------------------------------------------------
485     void initRunLumiRangeMap(ControlFlags &ctrl)
486     //----------------------------------------------------------------------------
487     {
488     /*
489     rlrm.AddJSONFile(std::string("./data/Cert_136033-149442_7TeV_Apr21ReReco_Collisions10_JSON.txt"));
490     // rlrm.AddJSONFile(std::string("./data/Cert_160404-173244_7TeV_PromptReco_Collisions11_JSON_v2.txt"));
491     rlrm.AddJSONFile(std::string("./data/Cert_160404-178078_7TeV_PromptReco_Collisions11_JSON.txt"));
492     rlrm.AddJSONFile(std::string("./data/Cert_160404-163869_7TeV_May10ReReco_Collisions11_JSON_v3.txt"));
493     rlrm.AddJSONFile(std::string("./data/Cert_170249-172619_7TeV_ReReco5Aug_Collisions11_JSON.txt"));
494     // r11b
495     rlrm.AddJSONFile(std::string("./data/Cert_160404-180252_7TeV_PromptReco_Collisions11_JSON.txt"));
496     */
497    
498     // 2012 only for now ...
499     rlrm.AddJSONFile(string(ctrl.jsonFile));
500    
501     };
502     //----------------------------------------------------------------------------------------
503     float getMass(const mithep::Muon *zmu1, const mithep::Muon *zmu2, const mithep::Electron *zele1, const mithep::Electron *zele2)
504     {
505     TLorentzVector lep1,lep2;
506     if(zmu1 && zmu2) {
507     lep1.SetPtEtaPhiM(zmu1->Pt(),zmu1->Eta(),zmu1->Phi(),MUON_MASS);
508     lep2.SetPtEtaPhiM(zmu2->Pt(),zmu2->Eta(),zmu2->Phi(),MUON_MASS);
509     } else if(zele1 && zele2) {
510     lep1.SetPtEtaPhiM(zele1->Pt(),zele1->Eta(),zele1->Phi(),ELECTRON_MASS);
511     lep2.SetPtEtaPhiM(zele2->Pt(),zele2->Eta(),zele2->Phi(),ELECTRON_MASS);
512     } else assert(0);
513    
514     return (lep1 + lep2).M();
515     }
516     //----------------------------------------------------------------------------------------
517     SelectionStatus analysisSelection(ControlFlags &ctrl, const mithep::Muon *mu, const mithep::Vertex *vtx, mithep::Array<mithep::PFCandidate> *pfArr,
518     mithep::Array<mithep::PileupEnergyDensity> *puDArr)
519     {
520     SelectionStatus status;
521     status |= muonReferencePreSelection(ctrl,mu,vtx,pfArr);
522     status |= muonIDPFSelection(ctrl,mu,vtx,pfArr);
523     status |= muonReferenceIsoSelection(ctrl,mu,vtx,pfArr,puDArr,eraMu,photonsToVeto);
524     return status;
525     }
526 dkralph 1.2 //----------------------------------------------------------------------------------------
527 dkralph 1.9 SelectionStatus analysisSelection(ControlFlags &ctrl, const mithep::Electron *ele, const mithep::Vertex *vtx, mithep::Array<mithep::PFCandidate> *pfArr,
528     mithep::Array<mithep::PileupEnergyDensity> *puDArr)
529     {
530     SelectionStatus status;
531     status |= electronReferencePreSelection(ctrl,ele,vtx);
532     status |= electronReferenceIDMVASelectionV1(ctrl,ele,vtx);
533     status |= electronReferenceIsoSelection(ctrl,ele,vtx,pfArr,puDArr,eraEle,photonsToVeto);
534     return status;
535     }
536     //----------------------------------------------------------------------------------------
537     vector<SimpleLepton> findFakes(ControlFlags &ctrl,
538     mithep::Array<mithep::Muon> *muonArr,
539     mithep::Array<mithep::Electron> *electronArr,
540     mithep::Array<mithep::PFJet> *jetArr,
541     mithep::Array<mithep::PFCandidate> *pfArr,
542     mithep::Array<mithep::PileupInfo> *puArr,
543     mithep::Array<mithep::PileupEnergyDensity> *puDArr,
544     mithep::Array<mithep::Vertex> *vtxArr,
545     const mithep::Vertex *vtx,
546     const mithep::Muon *zmu1,
547     const mithep::Muon *zmu2,
548     const mithep::Electron *zele1,
549     const mithep::Electron *zele2)
550     {
551     vector<SimpleLepton> fakes;
552    
553 dkralph 1.2 for(Int_t imu=0; imu<muonArr->GetEntries(); imu++) {
554     const mithep::Muon* mu = (*muonArr)[imu];
555    
556 dkralph 1.9 if(mu==zmu1 || mu==zmu2) continue;
557 dkralph 1.2
558 dkralph 1.7 SelectionStatus denomsel;
559 dkralph 1.8 denomsel |= muonPreSelectionNoD0DzIP(ctrl,mu,vtx,pfArr);;
560     if( !denomsel.passPre() ) continue;
561 dkralph 1.2
562 dkralph 1.9 SelectionStatus status = analysisSelection(ctrl,mu,vtx,pfArr,puDArr);
563 dkralph 1.2
564 dkralph 1.9 SimpleLepton fake;
565     fake.vec.SetPtEtaPhiM(mu->Pt(),mu->Eta(),mu->Phi(),MUON_MASS);
566     fake.type = 13;
567     fake.isLoose = status.loose();
568    
569     fakes.push_back(fake);
570 dkralph 1.2 }
571     for(Int_t iele=0; iele<electronArr->GetEntries(); iele++) {
572     const mithep::Electron* ele = (*electronArr)[iele];
573    
574 dkralph 1.4 bool muonMatch=false;
575     for(unsigned imu=0; imu<muonArr->GetEntries(); imu++) {
576     const mithep::Muon *mu = (*muonArr)[imu];
577     if(dr(mu,ele) < 0.05) muonMatch = true;
578     }
579 dkralph 1.9 if(muonMatch) continue;
580 dkralph 1.4
581 dkralph 1.9 if(ele==zele1 || ele==zele2) continue;
582 dkralph 1.2
583 dkralph 1.7 SelectionStatus denomsel;
584 dkralph 1.8 denomsel |= electronPreSelectionNoD0DzIP(ctrl,ele,vtx);
585     if( !denomsel.passPre() ) continue;
586 dkralph 1.2
587 dkralph 1.9 SelectionStatus status = analysisSelection(ctrl,ele,vtx,pfArr,puDArr);
588    
589     SimpleLepton fake;
590     fake.vec.SetPtEtaPhiM(ele->Pt(),ele->Eta(),ele->Phi(),ELECTRON_MASS);
591     fake.type = 11;
592     fake.isLoose = status.loose();
593 dkralph 1.8
594 dkralph 1.9 fakes.push_back(fake);
595     }
596 dkralph 1.2
597 dkralph 1.9 return fakes;
598 dkralph 1.1 }
599 dkralph 1.9 //----------------------------------------------------------------------------------------
600     SimpleLepton fillSimpleLepton(const Muon *mu, const Electron *ele, SelectionStatus *ctrl)
601 dkralph 1.2 {
602 dkralph 1.9 if(mu) return fillSimpleLepton(mu, ctrl);
603     else if(ele) return fillSimpleLepton(ele,ctrl);
604     else assert(0);
605     return SimpleLepton();
606     }