ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Selection/src/applySelection.cc
Revision: 1.18
Committed: Thu May 10 03:01:17 2012 UTC (13 years ago) by khahn
Content type: text/plain
Branch: MAIN
Changes since 1.17: +13 -20 lines
Log Message:
*** empty log message ***

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 using namespace std;
9
10 //
11 // root headers
12 //
13 #include <TFile.h>
14 #include <TTree.h>
15 #include <TChain.h>
16 #include <TBranch.h>
17 #include <TClonesArray.h>
18
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 "Electron.h"
33 #include "Muon.h"
34 #include "Vertex.h"
35 #include "PFCandidate.h"
36 #include "PFCandidateCol.h"
37 #include "PileupInfoCol.h"
38 #include "PileupEnergyDensity.h"
39 #include "TriggerMask.h"
40 #include "TriggerTable.h"
41 #include "Names.h"
42 #include "BaseMod.h"
43
44 //
45 // our headers
46 //
47 #include "ParseArgs.h"
48 #include "MuonSelection.h"
49 #include "ElectronSelection.h"
50 #include "IsolationSelection.h"
51
52 //#include "Selection.h"
53 #include "ReferenceSelection.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
65 #include "AngleTuple.h"
66 #include "FOTuple.h"
67
68 #include "SampleWeight.h"
69 #include "EfficiencyWeightsInterface.h"
70
71 #include "SimpleLepton.h"
72
73 #ifndef CMSSW_BASE
74 #define CMSSW_BASE "../"
75 #endif
76
77 using namespace mithep;
78
79 //
80 // globals
81 //----------------------------------------------------------------------------
82 TH1D * hpu;
83 RunLumiRangeMap rlrm;
84 vector<SimpleLepton> failingLeptons ; // for fake estimation
85 vector<SimpleLepton> passingLeptons; // for fake estimation
86
87 //
88 // prototypes
89 //----------------------------------------------------------------------------
90 bool setPV(ControlFlags,const mithep::Array<mithep::Vertex>*, mithep::Vertex& );
91 void initPUWeights();
92 double getPUWeight(unsigned npu);
93 void setEffiencyWeights(EventData&, WeightStruct& );
94 void initRunLumiRangeMap();
95 //----------------------------------------------------------------------------
96
97
98 //
99 // MAIN
100 //----------------------------------------------------------------------------
101 int main(int argc, char** argv)
102 //----------------------------------------------------------------------------
103 {
104 string cmsswpath(CMSSW_BASE);
105 cmsswpath.append("/src");
106 std::bitset<1024> triggerBits;
107 vector<vector<string> > inputFiles;
108 inputFiles.push_back(vector<string>());
109 map<string,unsigned> entrymap; // number of unskimmed entries in each file
110
111 //
112 // args
113 //--------------------------------------------------------------------------------------------------------------
114 ControlFlags ctrl;
115 parse_args( argc, argv, ctrl );
116 if( ctrl.inputfiles.empty() &&ctrl.inputfile.empty() )
117 {
118 cerr << "usage: applySelection.exe <flags> " << endl;
119 cerr << "\tmandoatory flags : " << endl;
120 cerr << "\t--inputfiles | file containing a list of ntuples to run over" << endl;
121 cerr << "\t--inputfile | a file to run over" << endl;
122 cerr << "\t--outputfile | file to store selected evet" << endl;
123 return 1;
124 }
125 ctrl.dump();
126
127
128
129 //
130 // File I/O
131 //--------------------------------------------------------------------------------------------------------------
132 TChain * chain = new TChain("Events");
133 TChain * hltchain = new TChain("HLT");
134
135 string fname;
136 if( !ctrl.inputfiles.empty() ) {
137 ifstream f(ctrl.inputfiles.c_str());
138 while (f >> fname) {
139 if( !(strncmp( fname.c_str(), "#", 1 ) ) ) continue; // skip commented lines
140 cout << "adding inputfile : " << fname.c_str() << endl;
141 entrymap[string(fname.c_str())] = unskimmedEntries(fname.c_str());
142 cout << "unskimmed entries: " << entrymap[string(fname.c_str())] << endl;
143 chain->AddFile(fname.c_str());
144 hltchain->AddFile(fname.c_str());
145 }
146 } else {
147 cout << "adding inputfile : " << ctrl.inputfile.c_str() << endl;
148 unsigned tmpent = unskimmedEntries(ctrl.inputfile.c_str());
149 cout << "tmpent: " << tmpent << endl;
150 entrymap[string(ctrl.inputfile.c_str())] = unskimmedEntries(ctrl.inputfile.c_str());
151 cout << "unskimmed entries: " << entrymap[string(ctrl.inputfile.c_str())] << endl;
152 chain->AddFile(ctrl.inputfile.c_str());
153 hltchain->AddFile(ctrl.inputfile.c_str());
154 }
155
156 const char * ofname;
157 if( strcmp( ctrl.outputfile.c_str(), "") ) {
158 ofname = ctrl.outputfile.c_str();
159 } else {
160 ofname = "tmp.root";
161 }
162 // table of cross section values
163 string xspath = (cmsswpath+string("/MitPhysics/data/xs.dat"));
164 cout << "xspath: " << xspath.c_str() << endl;
165 SimpleTable xstab(xspath.c_str());
166
167
168 //
169 // Setup
170 //--------------------------------------------------------------------------------------------------------------
171 Angles angles;
172 KinematicsStruct kinematics;
173 InfoStruct evtinfo;
174 WeightStruct weights;
175 // GenInfoStruct geninfo;
176
177 AngleTuple nt( (const char*)ofname, (const char*)"zznt");
178 nt.makeAngleBranch(angles);
179 nt.makeKinematicsBranch(kinematics);
180 nt.makeInfoBranch(evtinfo);
181
182 FOTuple foTree( nt.getFile(), (const char*)"FOtree");
183 foTree.makeFailedLeptonBranch(failingLeptons);
184 foTree.makeZLeptonBranch(passingLeptons);
185
186 TH1F * h_w_hpt;
187 if(ctrl.mc) {
188 nt.makeWeightBranch(weights);
189 // nt.makeGenInfoBranch(geninfo);
190 initEfficiencyWeights();
191 initPUWeights();
192 /*
193 // Higgs only, pt reweighting
194 if( ctrl.mH <= 140 ) {
195 char wptname[256];
196 sprintf( wptname, "/data/blue/pharris/Flat/ntupler/root/weight_ptH_%i.root", ctrl.mH );
197 TFile * f_hpt = new TFile(wptname);
198 f_hpt->Print();
199 sprintf(wptname, "weight_hqt_fehipro_fit_%i", ctrl.mH);
200 h_w_hpt = (TH1F*)(f_hpt->FindObjectAny(wptname));
201 h_w_hpt->Print();
202 }
203 */
204 } else {
205 initRunLumiRangeMap();
206 }
207
208 // initMuonIDMVA();
209 initMuonIsoMVA();
210 initElectronIDMVA();
211 initElectronIsoMVA();
212 initTrigger();
213
214
215
216 //
217 // Setup tree I/O
218 //--------------------------------------------------------------------------------------------------------------
219 TFile *inputFile=0;
220 TTree *eventTree=0;
221 string currentFile("");
222
223 mithep::EventHeader *info = new mithep::EventHeader();
224 // mithep::TGenInfo *ginfo = new mithep::TGenInfo();
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::PileupInfo> *puArr = new mithep::Array<mithep::PileupInfo>();
230 mithep::Array<mithep::PileupEnergyDensity> *puDArr = new mithep::Array<mithep::PileupEnergyDensity>();
231 mithep::Array<mithep::Track> *trkArr = new mithep::Array<mithep::Track>();
232 mithep::TriggerMask *trigMask = new mithep::TriggerMask();
233 mithep::TriggerTable *hltTable = new mithep::TriggerTable();
234 vector<string> *hltTableStrings = new vector<string>();
235
236 TString fElectronName(Names::gkElectronBrn);
237 TString fMuonName(Names::gkMuonBrn);
238 TString fInfoName(Names::gkEvtHeaderBrn);
239 TString fPrimVtxName(Names::gkPVBrn);
240 TString fPFCandidateName(Names::gkPFCandidatesBrn);
241 TString fPileupInfoName(Names::gkPileupInfoBrn);
242 TString fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn);
243 TString fTrackName(Names::gkTrackBrn);
244 TString fTriggerMaskName(Names::gkHltBitBrn);
245 TString fTriggerTableName(Names::gkHltTableBrn);
246
247 chain->SetBranchAddress(fInfoName, &info);
248 chain->SetBranchAddress(fElectronName, &electronArr);
249 chain->SetBranchAddress(fMuonName, &muonArr);
250 chain->SetBranchAddress(fPrimVtxName, &vtxArr);
251 chain->SetBranchAddress(fPFCandidateName, &pfArr);
252 if( ctrl.mc ) chain->SetBranchAddress(fPileupInfoName, &puArr);
253 chain->SetBranchAddress(fPileupEnergyDensityName, &puDArr);
254 chain->SetBranchAddress(fTrackName, &trkArr);
255 chain->SetBranchAddress(fTriggerMaskName, &trigMask);
256 cout << "hlttable: " << fTriggerTableName << endl;
257 hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
258
259 mithep::Vertex vtx; // best primary vertex in the event
260
261 // ginfo = NULL;
262 // if( ctrl.mc ) { chain->SetBranchAddress("Gen", &ginfo);}
263
264 int count=0, pass=0;
265 float passcorr=0., passcorr_errup=0., passcorr_errdown=0.;
266 float denom[3]={0.,0.,0.};
267 float numer[3]={0.,0.,0.};
268 float numercorr[3]={0.,0.,0.};
269
270 // only 1 HLT table / file ???
271 hltchain->GetEntry(0);
272 cerr << "here... " << endl;
273
274 int imax = chain->GetEntries();
275 cout << "nEntries: " << imax << endl;
276
277
278 //
279 // Loop !!!!!!!!!
280 //--------------------------------------------------------------------------------------------------------------
281 for(UInt_t ientry=0; ientry<imax; ientry++)
282 {
283 chain->GetEntry(ientry);
284 if(!(ientry%1000)) cerr << "entry: " << ientry << endl;
285 setPV( ctrl, vtxArr, vtx);
286
287
288 //
289 // data/MC
290 //
291 if(ctrl.mc) {
292 //
293 // xsec & PU weights
294 //
295 weights.w = getWeight(xstab,entrymap,chain)/ctrl.totalMC;
296 int npu = -1;
297 for(int i=0;i<puArr->GetEntries();i++) {
298 if(puArr->At(i)->GetBunchCrossing() == 0)
299 npu = puArr->At(i)->GetPU_NumInteractions();
300 }
301 assert(npu>=0);
302 evtinfo.npu;
303 weights.npuw = getPUWeight(evtinfo.npu);
304 // cout << "weight: " << weights.w << endl;
305 } else {
306 //
307 // JSON
308 //
309 /*
310 RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());
311 if( !(rlrm.HasRunLumi(rl)) ) {
312 if( ctrl.debug ) cout << "\tfails JSON" << endl;
313 continue;
314 }
315 */
316 //
317 // trigger
318 //
319 if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) {
320 currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile());
321 hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
322 hltchain->GetEntry(0);
323 hltTable->Clear();
324 fillTriggerBits(hltTable, hltTableStrings );
325 }
326 if( ctrl.debug ) cout << "file is : " << currentFile << endl;
327 fillTriggerBits( hltTable, trigMask, triggerBits );
328 if( !passHLT(triggerBits, info->RunNum(), info->EvtNum() ) ) {
329 if( ctrl.debug ) cout << "\tfails trigger ... " << endl;
330 continue;
331 }
332 }
333
334 //
335 // lepton/kinematic selection ...
336 //
337 failingLeptons.clear();
338 passingLeptons.clear();
339 EventData ret4l =
340 apply_HZZ4L_reference_selection(ctrl, info,
341 vtx,
342 pfArr,
343 puDArr,
344 electronArr,
345 &electronReferencePreSelection,
346 &electronReferenceIDMVASelection,
347 &electronReferenceIsoSelection,
348 muonArr,
349 &muonReferencePreSelection,
350 &muonIDPFSelection,
351 &muonReferenceIsoSelection);
352 /*
353 apply_HZZ4L_reference_selection(ctrl, info,
354 vtx,
355 pfArr,
356 puArr,
357 electronArr,
358 &electronPreSelection,
359 &electronIDMVASelection,
360 &electronIsoMVASelection,
361 muonArr,
362 &muonPreSelection,
363 &muonIDPFSelection,
364 &muonIsoMVASelection);
365 */
366
367 if( ctrl.debug ) cout << endl;
368
369 // cout << "bits: " << ret4l.status.selectionBits << endl;
370
371 if( ret4l.status.pass() ) {
372 fillAngles(ret4l,angles);
373 fillKinematics(ret4l,kinematics);
374 fillEventInfo(info,evtinfo);
375 if( ctrl.mc) {
376 // fillGenInfo(ginfo,geninfo);
377 setEffiencyWeights(ret4l, weights);
378 if(ctrl.debug)
379 cout << "w: " << weights.w << "\t"
380 << "won: " << weights.won << "\t"
381 << "woff: " << weights.woff << endl;
382 }
383
384 /*
385 // only for Higgs < 140
386 geninfo.weight *= h_w_hpt->GetBinContent(h_w_hpt->FindBin(geninfo.pt_zz));
387 angles.bdt = reader->EvaluateMVA("BDTG");
388 */
389 nt.Fill();
390
391 cerr << "PASS:: "
392 << "\trun: " << evtinfo.run
393 << "\tevt: " << evtinfo.evt
394 << "\tlumi: " << evtinfo.lumi
395 << "\tcostheta1: " << angles.costheta1
396 << "\tcostheta2: " << angles.costheta2
397 << "\tcostheta*: " << angles.costhetastar
398 << "\tbdt: " << angles.bdt
399 << endl;
400 pass++;
401 // if( pass > 3 ) break;
402
403 } else if( ret4l.status.selectionBits.test(PASS_ZCANDIDATE) ) { // save for fakes ...
404 if(ctrl.debug) {
405 cout << "failingLeptons : " << failingLeptons.size() << endl;
406 for( int f=0; f<failingLeptons.size(); f++ ) {
407 SimpleLepton l = failingLeptons[f];
408 cout << "f: " << f << "\t"
409 << "type: " << l.type << "\t"
410 << "pT: " << l.vec.Pt() << "\t"
411 << "eta: " << l.vec.Eta()
412 << endl;
413 }
414 }
415 foTree.Fill();
416 }
417 }
418
419
420 foTree.getFile()->cd();
421 foTree.getTree()->Write();
422 nt.WriteClose();
423 }
424
425 //----------------------------------------------------------------------------
426 //
427 // Get primary vertices
428 // Assumes primary vertices are ordered by sum-pT^2 (as should be in CMSSW)
429 // NOTE: if no PV is found from fitting tracks, the beamspot is used
430 //
431 //----------------------------------------------------------------------------
432 bool setPV(ControlFlags ctrl,
433 const mithep::Array<mithep::Vertex> * vtxArr,
434 mithep::Vertex & vtx)
435 //----------------------------------------------------------------------------
436 {
437 const Vertex *bestPV = 0;
438 float best_sumpt=-1;
439
440 // good PV requirements
441 const UInt_t fMinNTracksFit = 0;
442 const Double_t fMinNdof = 4;
443 const Double_t fMaxAbsZ = 24;
444 const Double_t fMaxRho = 2;
445
446 for(int i=0; i<vtxArr->GetEntries(); ++i) {
447 const Vertex *pv = (Vertex*)(vtxArr->At(i));
448 if( ctrl.debug ) cout << "vertex :: " << i << "\tntrks: " << pv->NTracks() << endl;
449
450 // Select best PV for corrected d0; if no PV passing cuts, the first PV in the collection will be used
451 // if(!pv->IsValid()) continue;
452 if(pv->NTracksFit() < fMinNTracksFit) continue;
453 if(pv->Ndof() < fMinNdof) continue;
454 if(fabs(pv->Z()) > fMaxAbsZ) continue;
455 if(pv->Position().Rho() > fMaxRho) continue;
456
457 // take the first ...
458 bestPV = pv;
459 break;
460
461 // this never reached ...
462 float tmp_sumpt=0;
463 for( int t=0; t<pv->NTracks(); t++ )
464 tmp_sumpt += pv->Trk(t)->Pt();
465
466 if( tmp_sumpt > best_sumpt ) {
467 bestPV = pv;
468 best_sumpt = tmp_sumpt;
469 if( ctrl.debug) cout << "new PV set, pt : " << best_sumpt << endl;
470 }
471 }
472
473 if(!bestPV) bestPV = vtxArr->At(0);
474 vtx.SetPosition(bestPV->X(),bestPV->Y(),bestPV->Z());
475 vtx.SetErrors(bestPV->XErr(),bestPV->YErr(),bestPV->ZErr());
476 return true;
477 };
478
479
480
481 //----------------------------------------------------------------------------
482 void setEffiencyWeights(EventData & evtdat, WeightStruct & weights )
483 //----------------------------------------------------------------------------
484 {
485 vector<SimpleLepton> lepvec = evtdat.Z1leptons;
486 lepvec.insert(lepvec.end(), evtdat.Z2leptons.begin(), evtdat.Z2leptons.end());
487 double w_offline=-1, werr_offline=0;
488 double w_online=-1, werr_online=0;
489
490 vector< pair <double,double> > wlegs; // pair here is eff & err
491 vector< pair <float,float> > mvec; // pair here is eta & pt
492
493 // vector< pair <float,float> > evec; // pair here is eta & pt
494 // now deal with medium vs loose
495 vector< pair< bool, pair <float,float> > > evec; // pair here is eta & pt
496
497 for( int k=0; k<lepvec.size(); k++ ) {
498 if( abs(lepvec[k].type) == 13 ) {
499 mvec.push_back( std::pair<float,float> (fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt()) );
500 wlegs.push_back( muonPerLegOfflineEfficiencyWeight( fabs(lepvec[k].vec.Eta()),
501 lepvec[k].vec.Pt() ) );
502 } else {
503
504 // now deal with medium vs loose
505 // evec.push_back( std::pair<float,float> (fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt()) );
506
507 std::pair<float,float> tmppair(fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt());
508 evec.push_back( std::pair<bool, std::pair<float,float> > (lepvec[k].isTight, tmppair) );
509
510 // wlegs.push_back( elePerLegOfflineEfficiencyWeight( fabs(lepvec[k].vec.Eta()),
511 // lepvec[k].vec.Pt() ) );
512
513 wlegs.push_back( elePerLegOfflineEfficiencyWeight( fabs(lepvec[k].vec.Eta()),
514 lepvec[k].vec.Pt(),
515 lepvec[k].isTight ) );
516
517 }
518 }
519
520 pair<double,double> offpair = getOfflineEfficiencyWeight( wlegs );
521 weights.woff = offpair.first;
522 weights.werroff = offpair.second;
523
524 pair<double,double> onpair = getOnlineEfficiencyWeight( mvec, evec );
525 weights.won = onpair.first;
526 weights.werron = onpair.second;
527 }
528
529
530 //----------------------------------------------------------------------------
531 void initRunLumiRangeMap()
532 //----------------------------------------------------------------------------
533 {
534 /*
535 rlrm.AddJSONFile(std::string("./data/Cert_136033-149442_7TeV_Apr21ReReco_Collisions10_JSON.txt"));
536 // rlrm.AddJSONFile(std::string("./data/Cert_160404-173244_7TeV_PromptReco_Collisions11_JSON_v2.txt"));
537 rlrm.AddJSONFile(std::string("./data/Cert_160404-178078_7TeV_PromptReco_Collisions11_JSON.txt"));
538 rlrm.AddJSONFile(std::string("./data/Cert_160404-163869_7TeV_May10ReReco_Collisions11_JSON_v3.txt"));
539 rlrm.AddJSONFile(std::string("./data/Cert_170249-172619_7TeV_ReReco5Aug_Collisions11_JSON.txt"));
540 // r11b
541 rlrm.AddJSONFile(std::string("./data/Cert_160404-180252_7TeV_PromptReco_Collisions11_JSON.txt"));
542 */
543 };
544
545
546 //----------------------------------------------------------------------------
547 void initPUWeights()
548 //----------------------------------------------------------------------------
549 {
550 TFile * puf = new TFile("data/PileupReweighting.Summer11DYmm_To_Full2011.root");
551 hpu = (TH1D*)(puf->Get("puWeights"));
552 }
553
554 //----------------------------------------------------------------------------
555 double getPUWeight(unsigned npu)
556 //----------------------------------------------------------------------------
557 {
558 return hpu->GetBinContent(hpu->FindBin(npu));
559 }