ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Selection/src/applySelection.cc
Revision: 1.61
Committed: Wed Oct 31 14:28:14 2012 UTC (12 years, 6 months ago) by anlevin
Content type: text/plain
Branch: MAIN
Changes since 1.60: +6 -5 lines
Log Message:
added trigger object matching for 2011

File Contents

# Content
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 <TRegexp.h>
16 #include <TTree.h>
17 #include <TChain.h>
18 #include <TBranch.h>
19 #include <TClonesArray.h>
20
21
22 //
23 // TMVA
24 //
25 #include "TMVA/Reader.h"
26 #include "TMVA/Tools.h"
27 #include "TMVA/Config.h"
28 #include "TMVA/MethodBDT.h"
29
30 #include "FactorizedJetCorrector.h"
31 #include "JetCorrectorParameters.h"
32
33 //
34 // ntuple format headers
35 //
36 #include "EventHeader.h"
37 #include "Electron.h"
38 #include "Muon.h"
39 #include "Vertex.h"
40 #include "PFMet.h"
41 #include "PFJetCol.h"
42 #include "PFCandidate.h"
43 #include "PFCandidateCol.h"
44 #include "PileupInfoCol.h"
45 #include "PileupEnergyDensity.h"
46 #include "MCParticle.h"
47 #include "TriggerMask.h"
48 #include "TriggerTable.h"
49 #include "TriggerObject.h"
50 #include "TriggerObjectRel.h"
51 #include "Names.h"
52 #include "BaseMod.h"
53 #include "TriggerObjectsTable.h"
54 #include "JetIDMVA.h"
55
56 //
57 // our headers
58 //
59 #include "ParseArgs.h"
60 #include "MuonSelection.h"
61 #include "ElectronSelection.h"
62 #include "ElectronMomentumCorrection.h"
63 #include "JetSelection.h"
64 #include "IsolationSelection.h"
65 #include "ReferenceSelection.h"
66 #include "MuonMomentumCorrection.h"
67
68 #include "TriggerUtils.h"
69 #include "PassHLT.h"
70 #include "Angles.h"
71 #include "KinematicsStruct.h"
72 #include "InfoStruct.h"
73 #include "GenInfoStruct.h"
74 #include "WeightStruct.h"
75 #include "AngleTuple.h"
76 #include "FOTuple.h"
77
78 #include "RunLumiRangeMap.h"
79 #include "SampleWeight.h"
80 #include "EfficiencyWeightsInterface.h"
81 #include "SelectionFuncs.h"
82
83 #include "SimpleLepton.h"
84
85 #ifndef CMSSW_BASE
86 #define CMSSW_BASE "../"
87 #endif
88
89 using namespace mithep;
90
91 //
92 // globals
93 //----------------------------------------------------------------------------
94 TH1D *hpu_2011, *hpu_2012;
95 RunLumiRangeMap rlrm;
96 vector<SimpleLepton> failingLeptons ; // for fake estimation
97 vector<SimpleLepton> passingLeptons; // for fake estimation
98 vector<unsigned> cutvec;
99 vector<vector<unsigned> > zcutvec;
100 vector<vector<unsigned> > zzcutvec;
101 map<unsigned,float> evtrhoMap;
102 vector<string> cutstrs;
103 bool passes_HLT;
104 vector<bool> PFnoPUflag;;
105 map<TString, map<TString,int>* > counts;
106 ElectronMomentumCorrection electron_momentum_correction;
107 MuCorr *muCorr;
108 //----------------------------------------------------------------------------
109 void initRunLumiRangeMap(ControlFlags &ctrl)
110 //----------------------------------------------------------------------------
111 {
112 if(ctrl.jsonFile==""){
113 cout << "WARNING: running with empty jsonFile, so setting noJSON to true" << endl;
114 ctrl.noJSON = true;
115 }
116 rlrm.AddJSONFile(string(ctrl.jsonFile));
117 }
118
119 //
120 // prototypes
121 //----------------------------------------------------------------------------
122
123 //
124 // MAIN
125 //----------------------------------------------------------------------------
126 int main(int argc, char** argv)
127 //----------------------------------------------------------------------------
128 {
129 vector<TString> cuts;
130 cuts.push_back("start"); // NOTE: is incremented in ReferenceSelection, so if you turn on JSONS be careful
131 cuts.push_back("trigger");
132 cuts.push_back("sfOsHiPt");
133 cuts.push_back("4<mZ2<120");
134 cuts.push_back("pt>20,10");
135 cuts.push_back("mll>4");
136 cuts.push_back("m4l>70");
137 cuts.push_back("mZ2>12");
138 cuts.push_back("m4l>100");
139
140 for(unsigned icut=0; icut<cuts.size(); icut++)
141 counts[cuts[icut]] = new map<TString,int>;
142
143 map<TString,map<TString,int>* >::iterator cit;
144 for(cit=counts.begin(); cit!=counts.end(); cit++) {
145 map<TString,int> *cts = (*cit).second;
146 (*cts)["all"] = int(0);
147 (*cts)["4e"] = int(0);
148 (*cts)["4m"] = int(0);
149 (*cts)["2e2m"] = int(0);
150 }
151
152 string cmsswpath(CMSSW_BASE);
153 cmsswpath.append("/src");
154 bitset<1024> triggerBits;
155 vector<vector<string> > inputFiles;
156 inputFiles.push_back(vector<string>());
157 map<string,unsigned> entrymap; // number of unskimmed entries in each file
158
159 //
160 // args
161 //--------------------------------------------------------------------------------------------------------------
162 ControlFlags ctrl;
163 parse_args( argc, argv, ctrl );
164 if( ctrl.inputfiles.empty() &&ctrl.inputfile.empty() )
165 {
166 cerr << "usage: applySelection.exe <flags> " << endl;
167 cerr << "\tmandoatory flags : " << endl;
168 cerr << "\t--inputfiles | file containing a list of ntuples to run over" << endl;
169 cerr << "\t--inputfile | a file to run over" << endl;
170 cerr << "\t--outputfile | file to store selected evet" << endl;
171 return 1;
172 }
173 ctrl.dump();
174 cout << "\n\nNOTE: now setting era with arguments... double check this doesnt mess something else up\n\n" << endl;
175 bool doJets = false;//TString(ctrl.mcfmfname).Contains("doJets");
176 assert(ctrl.doFSR);
177
178 //
179 // File I/O
180 //--------------------------------------------------------------------------------------------------------------
181 TChain * chain = new TChain("Events");
182 TChain * hltchain = new TChain("HLT");
183 TChain * runchain = new TChain(Names::gkRunTreeName);
184
185 string fname;
186 unsigned total_unskimmed=0;
187 if( !ctrl.inputfiles.empty() ) {
188 ifstream f(ctrl.inputfiles.c_str());
189 while (f >> fname) {
190 if( !(strncmp( fname.c_str(), "#", 1 ) ) ) continue; // skip commented lines
191 if(fname.find("/store/") != string::npos) fname = "root://eoscms//"+fname;
192 cout << "adding inputfile : " << fname.c_str() << endl;
193 entrymap[string(fname.c_str())] = unskimmedEntries(fname.c_str());
194 cout << "unskimmed entries: " << entrymap[string(fname.c_str())] << endl;
195 total_unskimmed += entrymap[string(fname.c_str())];
196 chain->AddFile(fname.c_str());
197 hltchain->AddFile(fname.c_str());
198 runchain->AddFile(fname.c_str());
199 }
200 } else {
201 if(ctrl.inputfile.find("/store/") != string::npos) ctrl.inputfile = "root://eoscms//"+ctrl.inputfile;
202 cout << "adding inputfile : " << ctrl.inputfile.c_str() << endl;
203 unsigned unsk_ents = unskimmedEntries(ctrl.inputfile.c_str());
204 entrymap[string(ctrl.inputfile.c_str())] = unsk_ents;
205 cout << "unskimmed entries: " << unsk_ents << endl;
206 total_unskimmed += unsk_ents;
207 chain->AddFile(ctrl.inputfile.c_str());
208 hltchain->AddFile(ctrl.inputfile.c_str());
209 runchain->AddFile(ctrl.inputfile.c_str());
210 }
211 // write the total number of unskimmed events that went into making this output file to a text file
212 writeEntries(ctrl,total_unskimmed);
213 ctrl.totalMC = total_unskimmed;
214
215 const char * ofname;
216 if( strcmp( ctrl.outputfile.c_str(), "") ) {
217 ofname = ctrl.outputfile.c_str();
218 } else {
219 ofname = "tmp.root";
220 }
221 // table of cross section values
222 string xspath = (cmsswpath+string("/MitPhysics/data/xs.dat"));
223 cout << "xspath: " << xspath.c_str() << endl;
224 SimpleTable xstab(xspath.c_str());
225
226 //
227 // Setup
228 //--------------------------------------------------------------------------------------------------------------
229 Angles angles;
230 KinematicsStruct kine;
231 InfoStruct evtinfo;
232 WeightStruct weights;
233 GenInfoStruct geninfo;
234 JetInfoStruct ji;
235 vector<SimpleLepton> simpleLeptonJets;
236
237 AngleTuple nt( (const char*)ofname, (const char*)"zznt");
238 nt.makeAngleBranch(angles);
239 nt.makeKinematicsBranch(kine);
240 nt.makeInfoBranch(evtinfo);
241 if(doJets)
242 nt.makeJetInfoBranch(ji);
243
244 int genSampleType(0);
245 if(ctrl.mc) {
246 nt.makeWeightBranch(weights);
247 if(ctrl.fillGen) {
248 genSampleType = getGenSampleType(ofname);
249 if(genSampleType == -1) {
250 cout << "\ngenSampleType returned as -1, setting fillGen to false" << endl;
251 ctrl.fillGen = false;
252 } else
253 nt.makeGenInfoBranch(geninfo);
254 }
255 initEfficiencyWeights();
256 // initPUWeights();
257 } else {
258 if(!(ctrl.noJSON) )
259 initRunLumiRangeMap(ctrl);
260 }
261
262 initElectronIDMVAV1();
263 initTrigger();
264 TrigInfo ti;
265 initAnalysisTriggers(ti);
266 if(ctrl.do_escale) {
267 electron_momentum_correction.initialize_regression("../EGamma/EGammaAnalysisTools/data/eleEnergyRegWeights_V1.root");
268 if(ctrl.mc) {
269 if(ctrl.era == 2011) ctrl.smearing_scaling_dataset = "Fall11";
270 else ctrl.smearing_scaling_dataset = "Summer12_DR53X_HCP2012";
271 } else {
272 if(ctrl.era == 2011) ctrl.smearing_scaling_dataset = "Jan16ReReco";
273 else ctrl.smearing_scaling_dataset = "2012Jul13ReReco";
274 }
275 }
276 if(ctrl.correct_muon_momentum) {
277 muCorr = new MuCorr;
278 muCorr->corr2011 = new rochcor;
279 muCorr->corr2012 = new rochcor2012;
280 }
281
282 vector<TString> fJetCorrParsv;
283 vector<JetCorrectorParameters> correctionParameters;
284 FactorizedJetCorrector *fJetCorrector;
285 JetIDMVA *fJetIDMVA;
286 if(doJets) {
287 initJets(ctrl, cmsswpath, fJetCorrParsv, correctionParameters, fJetCorrector, fJetIDMVA);
288 }
289
290 //
291 // Setup tree I/O
292 //--------------------------------------------------------------------------------------------------------------
293 TFile *inputFile=0;
294 TTree *eventTree=0;
295 string currentFile("");
296
297 UInt_t fNMaxTriggers = TRIGGER_BIG_NUMBER;
298 EventHeader *info = new EventHeader();
299 Array<PFMet> *metArr = new Array<PFMet>();
300 Array<Electron> *electronArr = new Array<Electron>();
301 Array<Muon> *muonArr = new Array<Muon>();
302 Array<Vertex> *vtxArr = new Array<Vertex>();
303 Array<PFCandidate> *pfArr = new Array<PFCandidate>();
304 Array<PFJet> *jetArr = new Array<PFJet>();
305 Array<PileupInfo> *puArr = new Array<PileupInfo>();
306 Array<PileupEnergyDensity> *puDArr = new Array<PileupEnergyDensity>();
307 Array<Track> *trkArr = new Array<Track>();
308 Array<MCParticle> *mcArr = new Array<MCParticle>();
309 MCEventInfo *mcEvtInfo = new MCEventInfo();
310 TriggerMask *trigMask = new TriggerMask();
311 TriggerTable *hltTable = new TriggerTable(fNMaxTriggers);
312 vector<string> *hltTableStrings = new vector<string>();
313 vector<string> *hltLabelStrings = new vector<string>();
314 Array<TriggerObject> *hltObjArr = new Array<TriggerObject>();
315 Array<TriggerObjectRel> *hltRelsArr = new Array<TriggerObjectRel>();
316 TriggerObjectsTable *fTrigObjs = new TriggerObjectsTable(hltTable,fNMaxTriggers);
317 RunInfo *runInfo = new RunInfo();
318
319 chain->SetBranchAddress(Names::gkEvtHeaderBrn, &info);
320 chain->SetBranchAddress("PFMet", &metArr);
321 chain->SetBranchAddress(Names::gkElectronBrn, &electronArr);
322 chain->SetBranchAddress(Names::gkMuonBrn, &muonArr);
323 chain->SetBranchAddress(Names::gkPVBrn, &vtxArr);
324 chain->SetBranchAddress(Names::gkPFCandidatesBrn, &pfArr);
325 chain->SetBranchAddress(Names::gkPFJetBrn, &jetArr);
326 if( ctrl.mc ) {
327 chain->SetBranchAddress(Names::gkPileupInfoBrn, &puArr);
328 chain->SetBranchAddress(Names::gkMCPartBrn, &mcArr);
329 chain->SetBranchAddress(Names::gkMCEvtInfoBrn, &mcEvtInfo);
330 }
331 chain->SetBranchAddress(Names::gkPileupEnergyDensityBrn, &puDArr);
332 chain->SetBranchAddress(Names::gkTrackBrn, &trkArr);
333 chain->SetBranchAddress(Names::gkHltBitBrn, &trigMask);
334 chain->SetBranchAddress(Names::gkHltObjBrn, &hltObjArr);
335 chain->SetBranchAddress("HLTObjectsRelation", &hltRelsArr);
336 hltchain->SetBranchAddress(Names::gkHltTableBrn, &hltTableStrings);
337 hltchain->SetBranchAddress(Names::gkHltLabelBrn, &hltLabelStrings);
338
339 runchain->SetBranchAddress(Names::gkRunInfoBrn, &runInfo);
340
341 Vertex vtx; // best primary vertex in the event
342
343 hltchain->GetEntry(0);
344
345 int imax = (ctrl.n_events_to_process < 0) ? chain->GetEntries() : ctrl.n_events_to_process;
346 cerr << "processing: " << imax << " ( / " << chain->GetEntries() << " )" << endl;
347
348 if(ctrl.mc && ctrl.debugOffMc) // necessary for batch mode -- we need debug output for, data but size is too large to store MC output
349 ctrl.debug = false;
350
351 int lastHltEntry=-2;
352 for(UInt_t ientry=0; ientry<imax; ientry++)
353 {
354 chain->GetEntry(ientry);
355 if(!(ientry%1000)) cerr << "entry: " << ientry << " ( / " << imax << " )" << endl;
356
357 if( ctrl.debug ) {
358 cout << "-----------------------------------------------------------------" << endl;
359 cout << "-----------------------------------------------------------------" << endl;
360 cout << "Run: " << info->RunNum()
361 << "\tEvt: " << info->EvtNum()
362 << "\tLumi: " << info->LumiSec()
363 << endl;
364 cerr << "Run: " << info->RunNum()
365 << "\tEvt: " << info->EvtNum()
366 << "\tLumi: " << info->LumiSec()
367 << endl;
368 cout << "-----------------------------------------------------------------" << endl;
369 }
370
371 runchain->GetEvent(info->RunEntry());
372 hltchain->GetEntry(runInfo->HltEntry());
373 hltTable->Clear();
374 fillTriggerNames(hltTable, hltTableStrings );
375 if(ctrl.debug && (lastHltEntry != runInfo->HltEntry())) hltTable->Print();
376 lastHltEntry = runInfo->HltEntry();
377 fillTriggerBits( hltTable, trigMask, triggerBits );
378 setHLTObjectRelations( hltObjArr, hltRelsArr, hltTableStrings, hltLabelStrings, fTrigObjs );
379 // printTriggerObjs( hltTable, fTrigObjs);
380
381 string fname = string(chain->GetFile()->GetEndpointUrl()->GetFile());
382 if(ctrl.debug) cout << "era is " << ctrl.era << endl;
383 if( ctrl.era == 0 ) {
384 setEra( fname, ctrl );
385 if(ctrl.debug) cout << "era is now " << ctrl.era << endl;
386 }
387
388 // pfNoPU
389 PFnoPUflag.clear();
390 int pfnopu_size = makePFnoPUArray( pfArr, PFnoPUflag, vtxArr );
391 assert(pfnopu_size == pfArr->GetEntries());
392
393 // trigger
394 if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) {
395 currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile());
396 }
397 if( ctrl.debug ) cout << "file is : " << currentFile << endl;
398 passes_HLT = passHLT(ctrl, triggerBits, ti, info->RunNum(), hltTable, fTrigObjs); // NOTE: run value is not used if it's MC
399
400
401 // data/MC
402 if(ctrl.mc) {
403 if(ctrl.fillGen) fillGenInfo( mcArr, mcEvtInfo, geninfo, genSampleType, ctrl);
404 // NOTE: xs weight and npuw are reset in merge.cc
405 weights.w = getWeight(xstab,entrymap,chain)/ctrl.totalMC;
406 if(ctrl.debug) cout ;
407 weights.npu = getNPU(puArr, 0);
408 weights.npuw = getPUWeight(ctrl.era, fname, weights.npu);
409 } else {
410 // JSON
411 if(!(ctrl.noJSON) ) {
412 RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());
413 if( !(rlrm.HasRunLumi(rl)) ) {
414 if( ctrl.debug ) cout << "fails JSON" << endl;
415 continue;
416 }
417 }
418 }
419
420 //
421 // lepton/kinematic selection ...
422 //
423 failingLeptons.clear();
424 passingLeptons.clear();
425 EventData ret4l = apply_HZZ4L_reference_selection(ctrl, info,
426 vtxArr,
427 pfArr,
428 puDArr,
429 electronArr,
430 &electronReferencePreSelection,
431 &electronReferenceIDMVASelectionV1, // not used anymore! (Blame the damn regression implementation, not me!)
432 &electronReferenceIsoSelection,
433 muonArr,
434 &muonReferencePreSelection,
435 &muonIDPFSelection,
436 &muonReferenceIsoSelection);
437
438 if( ctrl.debug ) cout << endl;
439
440 if( ret4l.status.pass() ) {
441 fillAngles(ret4l,angles);
442
443 fillKinematics(ret4l,kine);
444 // fillMassErrors(ret4l,muonArr,electronArr,kine);
445
446 if(doJets) {
447 fillJets(ctrl, jetArr, simpleLeptonJets, fJetCorrector, puDArr, fJetIDMVA, vtxArr, ret4l);
448 fillJetInfo(simpleLeptonJets,ji,ctrl);
449 }
450
451 TLorentzVector pfmet; pfmet.SetPxPyPzE(metArr->At(0)->Mex(),metArr->At(0)->Mey(),0,0);
452 fillEventInfo( info, pfmet, evtinfo, ctrl.mc ? getNPU(puArr) : 0, vtxArr->GetEntries());
453 if( ctrl.mc) {
454 setEfficiencyWeights(ctrl, ctrl.era, ret4l, triggerBits, hltTable, hltObjArr, fTrigObjs, weights);
455 if(ctrl.debug)
456 cout << "w: " << weights.w << "\t" << "npuw: " << weights.npuw << "\t" << "won: " << weights.won << "\t" << "woff: " << weights.woff << endl;
457 }
458
459 nt.Fill();
460
461 cerr << "PASS:: " << "\trun: " << evtinfo.run << "\tevt: " << evtinfo.evt << "\tlumi: " << evtinfo.lumi << "\tchannel: "
462 << kine.channel << "\tm4l: " << kine.m4l << "\tmZ1: " << kine.mZ1 << "\tmZ2: " << kine.mZ2 << endl;
463 cout << "PASS:: " << "\trun: " << evtinfo.run << "\tevt: " << evtinfo.evt << "\tlumi: " << evtinfo.lumi << "\tchannel: "
464 << kine.channel << "\tm4l: " << kine.m4l << "\tmZ1: " << kine.mZ1 << "\tmZ2: " << kine.mZ2 << endl;
465
466 } else {
467 if(ctrl.debug)
468 cout << "failing with code : " << ret4l.status.selectionBits << endl;
469 }
470 }
471
472 nt.WriteClose();
473
474 for(unsigned icut=0; icut<cuts.size(); icut++) {
475 map<TString,int> thisCount(*(counts[cuts[icut]]));
476 cout << setw(20) << cuts[icut] << "\t" << setw(9) << thisCount["all"] << " ("
477 << setw(9) << thisCount["4e"] << setw(9) << thisCount["4m"] << setw(9) << thisCount["2e2m"] << ")" << endl;
478 }
479 }