ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/NonMCBackground/src/applyEMU.cc
Revision: 1.4
Committed: Thu Jun 21 21:26:33 2012 UTC (12 years, 10 months ago) by dkralph
Content type: text/plain
Branch: MAIN
Changes since 1.3: +7 -14 lines
Log Message:
make json file an arg; add npu and met to event info ofiller

File Contents

# User Rev Content
1 dkralph 1.1 //
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 "TriggerMask.h"
42     #include "TriggerTable.h"
43     #include "Names.h"
44     #include "BaseMod.h"
45    
46     //
47     // our headers
48     //
49     #include "ParseArgs.h"
50     #include "MuonSelection.h"
51     #include "ElectronSelection.h"
52     #include "IsolationSelection.h"
53     #include "SelectionEMU.h"
54    
55     #include "TriggerUtils.h"
56     #include "PassHLT.h"
57     #include "Angles.h"
58     #include "KinematicsStruct.h"
59     #include "InfoStruct.h"
60     #include "GenInfoStruct.h"
61     #include "WeightStruct.h"
62     //#include "GlobalDataAndFuncs.h"
63     #include "RunLumiRangeMap.h"
64     #include "SelectionFuncs.h"
65    
66     #include "AngleTuple.h"
67     #include "FOTuple.h"
68    
69     #include "SampleWeight.h"
70     #include "EfficiencyWeightsInterface.h"
71    
72     #include "SimpleLepton.h"
73    
74     #ifndef CMSSW_BASE
75     #define CMSSW_BASE "../"
76     #endif
77    
78     using namespace mithep;
79    
80     //
81     // globals
82     //----------------------------------------------------------------------------
83 khahn 1.3 TH1D *hpu_2011, *hpu_2012;
84 dkralph 1.1 mithep::RunLumiRangeMap rlrm;
85     vector<SimpleLepton> failingLeptons ; // for fake estimation
86     vector<SimpleLepton> passingLeptons; // for fake estimation
87     vector<unsigned> cutvec;
88     vector<vector<unsigned> > zcutvec;
89     vector<vector<unsigned> > zzcutvec;
90     map<unsigned,float> evtrhoMap;
91     vector<string> cutstrs;
92     bool passes_HLT_MC;
93     vector<bool> PFnoPUflag;;
94    
95     //
96     // prototypes
97     //----------------------------------------------------------------------------
98 dkralph 1.4 void initRunLumiRangeMap(ControlFlags &ctrl);
99 dkralph 1.1
100     //
101     // MAIN
102     //----------------------------------------------------------------------------
103     int main(int argc, char** argv)
104     //----------------------------------------------------------------------------
105     {
106    
107     string cmsswpath(CMSSW_BASE);
108     cmsswpath.append("/src");
109     std::bitset<1024> triggerBits;
110     map<string,unsigned> entrymap; // number of unskimmed entries in each file
111    
112     //
113     // args
114     //--------------------------------------------------------------------------------------------------------------
115     ControlFlags ctrl;
116     parse_args( argc, argv, ctrl );
117     if( ctrl.inputfiles.empty() &&ctrl.inputfile.empty() )
118     {
119     cerr << "usage: applySelection.exe <flags> " << endl;
120     cerr << "\tmandoatory flags : " << endl;
121     cerr << "\t--inputfiles | file containing a list of ntuples to run over" << endl;
122     cerr << "\t--inputfile | a file to run over" << endl;
123     cerr << "\t--outputfile | file to store selected evet" << endl;
124     return 1;
125     }
126     ctrl.dump();
127    
128    
129    
130     //
131     // File I/O
132     //--------------------------------------------------------------------------------------------------------------
133     TChain * chain = new TChain("Events");
134     TChain * hltchain = new TChain("HLT");
135    
136     string fname;
137     unsigned total_unskimmed=0;
138     if( !ctrl.inputfiles.empty() ) {
139     ifstream f(ctrl.inputfiles.c_str());
140     while (f >> fname) {
141     if( !(strncmp( fname.c_str(), "#", 1 ) ) ) continue; // skip commented lines
142     cout << "adding inputfile : " << fname.c_str() << endl;
143     entrymap[string(fname.c_str())] = unskimmedEntries(fname.c_str());
144     total_unskimmed += entrymap[string(fname.c_str())];
145     chain->AddFile(fname.c_str());
146     hltchain->AddFile(fname.c_str());
147     }
148     } else {
149     cout << "adding inputfile : " << ctrl.inputfile.c_str() << endl;
150     unsigned unsk_ents = unskimmedEntries(ctrl.inputfile.c_str());
151     entrymap[string(ctrl.inputfile.c_str())] = unsk_ents;
152     total_unskimmed += unsk_ents;
153     chain->AddFile(ctrl.inputfile.c_str());
154     hltchain->AddFile(ctrl.inputfile.c_str());
155     }
156     // write the total number of unskimmed events that went into making this output file to a text file
157     writeEntries(ctrl,total_unskimmed);
158    
159     const char * ofname;
160     if( strcmp( ctrl.outputfile.c_str(), "") ) {
161     ofname = ctrl.outputfile.c_str();
162     } else {
163     ofname = "tmp.root";
164     }
165     // table of cross section values
166     string xspath = (cmsswpath+string("/MitPhysics/data/xs.dat"));
167     cout << "xspath: " << xspath.c_str() << endl;
168     SimpleTable xstab(xspath.c_str());
169    
170     //
171     // Setup
172     //--------------------------------------------------------------------------------------------------------------
173     Angles angles;
174     KinematicsStruct kinematics;
175     InfoStruct evtinfo;
176     WeightStruct weights;
177     GenInfoStruct geninfo;
178    
179     AngleTuple nt( (const char*)ofname, (const char*)"zznt");
180     nt.makeAngleBranch(angles);
181     nt.makeKinematicsBranch(kinematics);
182     nt.makeInfoBranch(evtinfo);
183    
184    
185     FOTuple foTree( nt.getFile(), (const char*)"FOtree");
186     foTree.makeFailedLeptonBranch(failingLeptons);
187     foTree.makePassedLeptonBranch(passingLeptons);
188    
189     TH1F * h_wf11_hpt;
190     if(ctrl.mc) {
191     nt.makeWeightBranch(weights);
192     if(ctrl.fillGen ) nt.makeGenInfoBranch(geninfo);
193     initEfficiencyWeights();
194     initPUWeights();
195    
196     // Higgs only, pt reweighting
197     if( ctrl.mH <= 140 ) {
198     char wptname[256];
199     sprintf( wptname, "/scratch/dkralph/pharris/Flat/ntupler/root/weight_ptH_%i.root", ctrl.mH );
200     TFile * f_hpt = new TFile(wptname);
201     f_hpt->Print();
202     sprintf(wptname, "weight_hqt_fehipro_fit_%i", ctrl.mH);
203     h_wf11_hpt = (TH1F*)(f_hpt->FindObjectAny(wptname));
204     }
205    
206     } else {
207 dkralph 1.4 initRunLumiRangeMap(ctrl);
208 dkralph 1.1 }
209    
210     // initMuonIDMVA();
211     initElectronIDMVAV1();
212     initMuonIsoMVA();
213     initElectronIDMVA();
214     initElectronIsoMVA();
215     initTrigger();
216    
217     // hack
218     initEvtRhoMap(evtrhoMap);
219    
220    
221     //
222     // Setup tree I/O
223     //--------------------------------------------------------------------------------------------------------------
224     string currentFile("");
225    
226     mithep::EventHeader *info = new mithep::EventHeader();
227     mithep::Array<mithep::PFMet> *metArr = new mithep::Array<mithep::PFMet>();
228     mithep::Array<mithep::Electron> *electronArr = new mithep::Array<mithep::Electron>();
229     mithep::Array<mithep::Muon> *muonArr = new mithep::Array<mithep::Muon>();
230     mithep::Array<mithep::Vertex> *vtxArr = new mithep::Array<mithep::Vertex>();
231     mithep::Array<mithep::PFCandidate> *pfArr = new mithep::Array<mithep::PFCandidate>();
232     mithep::Array<mithep::PileupInfo> *puArr = new mithep::Array<mithep::PileupInfo>();
233     mithep::Array<mithep::PileupEnergyDensity> *puDArr = new mithep::Array<mithep::PileupEnergyDensity>();
234     mithep::Array<mithep::Track> *trkArr = new mithep::Array<mithep::Track>();
235     mithep::Array<mithep::MCParticle> *mcArr = new mithep::Array<mithep::MCParticle>();
236     mithep::MCEventInfo *mcEvtInfo = new mithep::MCEventInfo();
237     mithep::TriggerMask *trigMask = new mithep::TriggerMask();
238     mithep::TriggerTable *hltTable = new mithep::TriggerTable();
239     vector<string> *hltTableStrings = new vector<string>();
240    
241     TString fElectronName(Names::gkElectronBrn);
242     TString fMuonName(Names::gkMuonBrn);
243     TString fInfoName(Names::gkEvtHeaderBrn);
244     TString fPFMetName("PFMet");
245     TString fPrimVtxName(Names::gkPVBrn);
246     TString fPFCandidateName(Names::gkPFCandidatesBrn);
247     TString fPileupInfoName(Names::gkPileupInfoBrn);
248     TString fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn);
249     TString fTrackName(Names::gkTrackBrn);
250     TString fMCParticleName(Names::gkMCPartBrn);
251     TString fMCEvtInfoName(Names::gkMCEvtInfoBrn);
252     TString fTriggerMaskName(Names::gkHltBitBrn);
253     TString fTriggerTableName(Names::gkHltTableBrn);
254    
255     chain->SetBranchAddress(fInfoName, &info);
256     chain->SetBranchAddress(fPFMetName, &metArr);
257     chain->SetBranchAddress(fElectronName, &electronArr);
258     chain->SetBranchAddress(fMuonName, &muonArr);
259     chain->SetBranchAddress(fPrimVtxName, &vtxArr);
260     chain->SetBranchAddress(fPFCandidateName, &pfArr);
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     mithep::Vertex vtx; // best primary vertex in the event
273    
274     // 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    
283    
284     //
285     // Loop !!!!!!!!!
286     //--------------------------------------------------------------------------------------------------------------
287     for(UInt_t ientry=0; ientry<imax; ientry++) {
288     chain->GetEntry(ientry);
289     if(!(ientry%1000)) cerr << "entry: " << ientry << "\t" << float(ientry)/imax << endl;
290    
291     string fname = string(chain->GetFile()->GetEndpointUrl()->GetFile());
292     setEra( fname, ctrl );
293    
294    
295     // pfNoPU
296     PFnoPUflag.clear();
297     int pfnopu_size = makePFnoPUArray( pfArr, PFnoPUflag, vtxArr );
298     assert(pfnopu_size == pfArr->GetEntries());
299    
300     // trigger
301     if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) {
302     currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile());
303     hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
304     hltchain->GetEntry(0);
305     hltTable->Clear();
306     fillTriggerNames(hltTable, hltTableStrings );
307     }
308     fillTriggerBits( hltTable, trigMask, triggerBits );
309    
310     //
311     // data/MC
312     //
313     if(ctrl.mc) {
314     if( ctrl.fillGen )
315     fillGenInfo( mcArr, mcEvtInfo, geninfo, ESampleType::kHZZ, ctrl);
316     } else {
317     if(!(ctrl.noJSON) ) {
318 dkralph 1.4 RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());
319     if( !(rlrm.HasRunLumi(rl)) ) continue;
320 dkralph 1.1 }
321     // if( !passHLTEMU(ctrl,triggerBits) ) {
322     // continue;
323     // }
324     }
325    
326     // lepton/kinematic selection ...
327     failingLeptons.clear();
328     passingLeptons.clear();
329     EventData ret4l =
330    
331     apply_HZZ4L_EMU_selection(ctrl, info,
332     vtxArr,
333     pfArr,
334     puDArr,
335     electronArr,
336     &electronReferencePreSelection,
337     &electronReferenceIDMVASelectionV1,
338     &electronReferenceIsoSelection,
339     muonArr,
340     &muonReferencePreSelection,
341     &muonIDPFSelection,
342     &muonReferenceIsoSelection);
343    
344    
345     if( ret4l.status.pass() ) {
346    
347     TLorentzVector pfmet; pfmet.SetPxPyPzE(metArr->At(0)->Mex(),metArr->At(0)->Mey(),0,0);
348 dkralph 1.4 fillEventInfo( info, pfmet, evtinfo, ctrl.mc ? getNPU(puArr) : 0);
349 dkralph 1.1 foTree.Fill();
350    
351     if( ctrl.mc)
352 khahn 1.3 setEffiencyWeights(ctrl.era, ret4l, weights);
353 dkralph 1.1
354     nt.Fill();
355     }
356    
357     }
358     foTree.getFile()->cd();
359     foTree.getTree()->Write();
360     nt.WriteClose();
361     }
362     //----------------------------------------------------------------------------
363 dkralph 1.4 void initRunLumiRangeMap(ControlFlags &ctrl)
364 dkralph 1.1 //----------------------------------------------------------------------------
365     {
366     /*
367     rlrm.AddJSONFile(std::string("./data/Cert_136033-149442_7TeV_Apr21ReReco_Collisions10_JSON.txt"));
368     // rlrm.AddJSONFile(std::string("./data/Cert_160404-173244_7TeV_PromptReco_Collisions11_JSON_v2.txt"));
369     rlrm.AddJSONFile(std::string("./data/Cert_160404-178078_7TeV_PromptReco_Collisions11_JSON.txt"));
370     rlrm.AddJSONFile(std::string("./data/Cert_160404-163869_7TeV_May10ReReco_Collisions11_JSON_v3.txt"));
371     rlrm.AddJSONFile(std::string("./data/Cert_170249-172619_7TeV_ReReco5Aug_Collisions11_JSON.txt"));
372     // r11b
373     rlrm.AddJSONFile(std::string("./data/Cert_160404-180252_7TeV_PromptReco_Collisions11_JSON.txt"));
374     */
375    
376     // 2012 only for now ...
377 dkralph 1.4 rlrm.AddJSONFile(string(ctrl.jsonFile));
378 dkralph 1.1
379     };