ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/NonMCBackground/src/applyEMU.cc
Revision: 1.3
Committed: Wed Jun 20 19:17:19 2012 UTC (12 years, 10 months ago) by khahn
Content type: text/plain
Branch: MAIN
Changes since 1.2: +3 -21 lines
Log Message:
move PU stuff to selectionFuncs, add 2012

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