ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/NonMCBackground/src/applyEMU.cc
Revision: 1.2
Committed: Mon Jun 18 06:56:17 2012 UTC (12 years, 11 months ago) by dkralph
Content type: text/plain
Branch: MAIN
Changes since 1.1: +0 -4 lines
Log Message:
*** empty log message ***

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