ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/MitHzz4l/Selection/src/applySelection.cc

Comparing UserCode/MitHzz4l/Selection/src/applySelection.cc (file contents):
Revision 1.61 by anlevin, Wed Oct 31 14:28:14 2012 UTC vs.
Revision 1.62 by dkralph, Mon Dec 17 17:12:21 2012 UTC

#	Line 73 | Line 73 | using namespace std;
73		#include "GenInfoStruct.h"
74		#include "WeightStruct.h"
75		#include "AngleTuple.h"
76	+	#include "JetInfoStruct.h"
77		#include "FOTuple.h"
78
79		#include "RunLumiRangeMap.h"
#	Line 99 | Line 100 | vector<unsigned> cutvec;
100		vector<vector<unsigned> > zcutvec;
101		vector<vector<unsigned> > zzcutvec;
102		map<unsigned,float>       evtrhoMap;
103	+	TrigInfo ti;
104		vector<string> cutstrs;
105		bool passes_HLT;
106		vector<bool>   PFnoPUflag;;
107		map<TString, map<TString,int>* > counts;
108		ElectronMomentumCorrection electron_momentum_correction;
109		MuCorr *muCorr;
110	+	void addDummyZ2Lepotons(EventData &ret4l);
111	+	void setNtupleVals(ControlFlags &ctrl, EventData &ret4l, InfoStruct &evtinfo, Angles &angles, KinematicsStruct &kine,
112	+	JetInfoStruct &ji, WeightStruct &weights, EventHeader *info, Array<PFJet> *jetArr, Array<PFMet> *metArr,
113	+	vector<SimpleLepton> &simpleLeptonJets, FactorizedJetCorrector *fJetCorrector, Array<PileupInfo> *puArr,
114	+	Array<PileupEnergyDensity> *puDArr, JetIDMVA *fJetIDMVA, Array<Vertex> *vtxArr, bitset<1024> &triggerBits,
115	+	TriggerTable *hltTable, Array<TriggerObject> *hltObjArr, TriggerObjectsTable *fTrigObjs, AngleTuple *nt);
116		//----------------------------------------------------------------------------
117		void initRunLumiRangeMap(ControlFlags &ctrl)
118		//----------------------------------------------------------------------------
#	Line 161 | Line 169 | int main(int argc, char** argv)
169		//--------------------------------------------------------------------------------------------------------------
170		ControlFlags ctrl;
171		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	–	}
172		ctrl.dump();
173		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");
174		assert(ctrl.doFSR);
175
176		//
#	Line 191 | Line 189 | int main(int argc, char** argv)
189		if(fname.find("/store/") != string::npos) fname = "root://eoscms//"+fname;
190		cout << "adding inputfile : " << fname.c_str() << endl;
191		entrymap[string(fname.c_str())] = unskimmedEntries(fname.c_str());
192	<	cout << "unskimmed entries: " << entrymap[string(fname.c_str())] << endl;
192	>	cerr << "unskimmed entries: " << entrymap[string(fname.c_str())] << endl;
193		total_unskimmed += entrymap[string(fname.c_str())];
194		chain->AddFile(fname.c_str());
195		hltchain->AddFile(fname.c_str());
#	Line 202 | Line 200 | int main(int argc, char** argv)
200		cout << "adding inputfile : " << ctrl.inputfile.c_str() << endl;
201		unsigned unsk_ents = unskimmedEntries(ctrl.inputfile.c_str());
202		entrymap[string(ctrl.inputfile.c_str())] = unsk_ents;
203	<	cout << "unskimmed entries: " << unsk_ents << endl;
203	>	cerr << "unskimmed entries: " << unsk_ents << endl;
204		total_unskimmed += unsk_ents;
205		chain->AddFile(ctrl.inputfile.c_str());
206		hltchain->AddFile(ctrl.inputfile.c_str());
#	Line 222 | Line 220 | int main(int argc, char** argv)
220		string xspath = (cmsswpath+string("/MitPhysics/data/xs.dat"));
221		cout << "xspath: " << xspath.c_str() << endl;
222		SimpleTable xstab(xspath.c_str());
223	+	if(ctrl.mc) {
224	+	getWeight(xstab,entrymap,chain); // test *now* whether we've got the xs in here, so it fails before it creates the output
225	+	}
226
227		//
228		// Setup
#	Line 238 | Line 239 | int main(int argc, char** argv)
239		nt.makeAngleBranch(angles);
240		nt.makeKinematicsBranch(kine);
241		nt.makeInfoBranch(evtinfo);
242	<	if(doJets)
243	<	nt.makeJetInfoBranch(ji);
242	>	nt.makeJetInfoBranch(ji);
243	>	nt.makeJetAngleBranches("DataStruct/data/jetAngleVars.txt");
244	>
245	>	FOTuple *foTree=0;
246	>	if(ctrl.fakeScheme != "none") {
247	>	foTree = new FOTuple( nt.getFile(), (const char*)"FOtree");
248	>	foTree->makeInfoBranch(evtinfo);
249	>	foTree->makeFailedLeptonBranch(failingLeptons);
250	>	foTree->makePassedLeptonBranch(passingLeptons);
251	>	}
252
253		int genSampleType(0);
254		if(ctrl.mc) {
#	Line 261 | Line 270 | int main(int argc, char** argv)
270
271		initElectronIDMVAV1();
272		initTrigger();
264	–	TrigInfo ti;
273		initAnalysisTriggers(ti);
274	+	ti.dump();
275	+
276		if(ctrl.do_escale) {
277	<	electron_momentum_correction.initialize_regression("../EGamma/EGammaAnalysisTools/data/eleEnergyRegWeights_V1.root");
277	>	string regPath("EGamma/EGammaAnalysisTools/data/eleEnergyRegWeights_V1.root");
278	>	// if(TString(getenv("HOSTNAME")).Contains(TRegexp("t3[bd][te][cs][hk]")))
279	>	regPath = cmsswpath + "/" + regPath;
280	>	electron_momentum_correction.initialize_regression(regPath);
281		if(ctrl.mc) {
282		if(ctrl.era == 2011) ctrl.smearing_scaling_dataset = "Fall11";
283		else                 ctrl.smearing_scaling_dataset = "Summer12_DR53X_HCP2012";
#	Line 283 | Line 296 | int main(int argc, char** argv)
296		vector<JetCorrectorParameters> correctionParameters;
297		FactorizedJetCorrector *fJetCorrector;
298		JetIDMVA *fJetIDMVA;
299	<	if(doJets) {
287	<	initJets(ctrl, cmsswpath, fJetCorrParsv, correctionParameters, fJetCorrector, fJetIDMVA);
288	<	}
299	>	initJets(ctrl, cmsswpath, fJetCorrParsv, correctionParameters, fJetCorrector, fJetIDMVA);
300
301		//
302		// Setup tree I/O
#	Line 349 | Line 360 | int main(int argc, char** argv)
360		ctrl.debug = false;
361
362		int lastHltEntry=-2;
363	<	for(UInt_t ientry=0; ientry<imax; ientry++)
364	<	{
365	<	chain->GetEntry(ientry);
366	<	if(!(ientry%1000)) cerr << "entry: " << ientry << " ( / " << imax << " )" << endl;
367	<
368	<	if( ctrl.debug ) {
369	<	cout << "-----------------------------------------------------------------" << endl;
370	<	cout << "-----------------------------------------------------------------" << endl;
371	<	cout << "Run: " << info->RunNum()
372	<	<< "\tEvt: " << info->EvtNum()
373	<	<< "\tLumi: " << info->LumiSec()
374	<	<< endl;
375	<	cerr << "Run: " << info->RunNum()
376	<	<< "\tEvt: " << info->EvtNum()
377	<	<< "\tLumi: " << info->LumiSec()
378	<	<< endl;
379	<	cout << "-----------------------------------------------------------------" << endl;
380	<	}
363	>	for(UInt_t ientry=0; ientry<imax; ientry++) {
364	>	chain->GetEntry(ientry);
365	>	if(!(ientry%1000)) cerr << "entry: " << ientry << " ( / " << imax << " )" << endl;
366	>
367	>	increment(counts,"start");
368	>	if( ctrl.debug ) {
369	>	cout << "-----------------------------------------------------------------" << endl;
370	>	cout << "-----------------------------------------------------------------" << endl;
371	>	cout << "Run: " << info->RunNum()
372	>	<< "\tEvt: " << info->EvtNum()
373	>	<< "\tLumi: " << info->LumiSec()
374	>	<< endl;
375	>	cerr << "Run: " << info->RunNum()
376	>	<< "\tEvt: " << info->EvtNum()
377	>	<< "\tLumi: " << info->LumiSec()
378	>	<< endl;
379	>	cout << "-----------------------------------------------------------------" << endl;
380	>	}
381
382	<	runchain->GetEvent(info->RunEntry());
383	<	hltchain->GetEntry(runInfo->HltEntry());
384	<	hltTable->Clear();
385	<	fillTriggerNames(hltTable, hltTableStrings );
386	<	if(ctrl.debug && (lastHltEntry != runInfo->HltEntry())) hltTable->Print();
387	<	lastHltEntry = runInfo->HltEntry();
388	<	fillTriggerBits( hltTable, trigMask, triggerBits );
389	<	setHLTObjectRelations( hltObjArr, hltRelsArr, hltTableStrings, hltLabelStrings, fTrigObjs );
390	<	// printTriggerObjs( hltTable, fTrigObjs);
382	>	runchain->GetEvent(info->RunEntry());
383	>	hltchain->GetEntry(runInfo->HltEntry());
384	>	hltTable->Clear();
385	>	fillTriggerNames(hltTable, hltTableStrings );
386	>	// if(ctrl.debug && (lastHltEntry != runInfo->HltEntry())) hltTable->Print();
387	>	lastHltEntry = runInfo->HltEntry();
388	>	fillTriggerBits( hltTable, trigMask, triggerBits );
389	>	setHLTObjectRelations( hltObjArr, hltRelsArr, hltTableStrings, hltLabelStrings, fTrigObjs );
390	>	// printTriggerObjs( hltTable, fTrigObjs);
391
392	<	string fname = string(chain->GetFile()->GetEndpointUrl()->GetFile());
393	<	if(ctrl.debug) cout << "era is " << ctrl.era << endl;
394	<	if( ctrl.era == 0  ) {
395	<	setEra( fname, ctrl );
396	<	if(ctrl.debug) cout << "era is now " << ctrl.era << endl;
397	<	}
392	>	string fname = string(chain->GetFile()->GetEndpointUrl()->GetFile());
393	>	if(ctrl.debug) cout << "era is " << ctrl.era << endl;
394	>	if( ctrl.era == 0  ) {
395	>	setEra( fname, ctrl );
396	>	if(ctrl.debug) cout << "era is now " << ctrl.era << endl;
397	>	}
398
399	<	// pfNoPU
400	<	PFnoPUflag.clear();
401	<	int pfnopu_size = makePFnoPUArray( pfArr, PFnoPUflag, vtxArr );
402	<	assert(pfnopu_size == pfArr->GetEntries());
403	<
404	<	// trigger
405	<	if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) {
406	<	currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile());
407	<	}
408	<	if( ctrl.debug ) cout << "file is : " << currentFile  << endl;
409	<	passes_HLT = passHLT(ctrl, triggerBits, ti, info->RunNum(), hltTable, fTrigObjs); // NOTE: run value is not used if it's MC
399	>	// pfNoPU
400	>	PFnoPUflag.clear();
401	>	int pfnopu_size = makePFnoPUArray( pfArr, PFnoPUflag, vtxArr );
402	>	assert(pfnopu_size == pfArr->GetEntries());
403	>
404	>	// trigger
405	>	if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) {
406	>	currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile());
407	>	}
408	>	if( ctrl.debug ) cout << "file is : " << currentFile  << endl;
409	>	passes_HLT = passHLT(ctrl, triggerBits, ti, info->RunNum(), hltTable, fTrigObjs); // NOTE: run value is not used if it's MC
410
411	<
412	<	// data/MC
413	<	if(ctrl.mc) {
414	<	if(ctrl.fillGen) fillGenInfo( mcArr, mcEvtInfo, geninfo, genSampleType, ctrl);
415	<	// NOTE: xs weight and npuw are reset in merge.cc
416	<	weights.w = getWeight(xstab,entrymap,chain)/ctrl.totalMC;
417	<	if(ctrl.debug) cout ;
418	<	weights.npu = getNPU(puArr, 0);
419	<	weights.npuw = getPUWeight(ctrl.era, fname, weights.npu);
420	<	} else {
421	<	// JSON
422	<	if(!(ctrl.noJSON) ) {
423	<	RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());
424	<	if( !(rlrm.HasRunLumi(rl)) )  {
425	<	if( ctrl.debug ) cout << "fails JSON" << endl;
415	<	continue;
416	<	}
411	>	// data/MC
412	>	if(ctrl.mc) {
413	>	if(ctrl.fillGen) fillGenInfo( mcArr, mcEvtInfo, geninfo, genSampleType, ctrl);
414	>	// NOTE: xs weight and npuw are reset in merge.cc
415	>	weights.w = getWeight(xstab,entrymap,chain)/ctrl.totalMC;
416	>	if(ctrl.debug) cout ;
417	>	weights.npu = getNPU(puArr, 0);
418	>	weights.npuw = getPUWeight(ctrl.era, fname, weights.npu);
419	>	} else {
420	>	// JSON
421	>	if(!(ctrl.noJSON) ) {
422	>	RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());
423	>	if( !(rlrm.HasRunLumi(rl)) )  {
424	>	if( ctrl.debug ) cout << "fails JSON" << endl;
425	>	continue;
426		}
427		}
428	+	}
429
430	<	//
431	<	// lepton/kinematic selection ...
432	<	//
433	<	failingLeptons.clear();
434	<	passingLeptons.clear();
435	<	EventData ret4l = apply_HZZ4L_reference_selection(ctrl, info,
436	<	vtxArr,
437	<	pfArr,
438	<	puDArr,
439	<	electronArr,
440	<	&electronReferencePreSelection,
441	<	&electronReferenceIDMVASelectionV1, // not used anymore! (Blame the damn regression implementation, not me!)
442	<	&electronReferenceIsoSelection,
443	<	muonArr,
444	<	&muonReferencePreSelection,
445	<	&muonIDPFSelection,
446	<	&muonReferenceIsoSelection);
430	>	//
431	>	// lepton/kinematic selection ...
432	>	//
433	>	failingLeptons.clear();
434	>	passingLeptons.clear();
435	>	EventData ret4l = apply_HZZ4L_reference_selection(ctrl, info,
436	>	vtxArr,
437	>	pfArr,
438	>	puDArr,
439	>	electronArr,
440	>	&electronReferencePreSelection,
441	>	&electronReferenceIDMVASelectionV1, // not used anymore! (Blame the damn regression implementation, not me!)
442	>	&electronReferenceIsoSelection,
443	>	muonArr,
444	>	&muonReferencePreSelection,
445	>	&muonIDPFSelection,
446	>	&muonReferenceIsoSelection);
447
448	<	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);
448	>	if( ctrl.debug ) cout << endl;
449
450	<	if(doJets) {
451	<	fillJets(ctrl, jetArr, simpleLeptonJets, fJetCorrector, puDArr, fJetIDMVA, vtxArr, ret4l);
452	<	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	<
450	>	if(ctrl.fakeScheme == "none") { // fill angletuple if event passes full selection
451	>	if( ret4l.status.pass() ) {
452	>	setNtupleVals(ctrl, ret4l, evtinfo, angles, kine, ji, weights, info, jetArr, metArr, simpleLeptonJets, fJetCorrector, puArr, puDArr, fJetIDMVA, vtxArr, triggerBits, hltTable, hltObjArr, fTrigObjs, &nt);
453		nt.Fill();
454	<
455	<	cerr  << "PASS:: " << "\trun: " << evtinfo.run << "\tevt: " << evtinfo.evt << "\tlumi: " << evtinfo.lumi << "\tchannel: "
456	<	<< kine.channel << "\tm4l: " << kine.m4l << "\tmZ1: " << kine.mZ1 << "\tmZ2: " << kine.mZ2 << endl;
457	<	cout  << "PASS:: " << "\trun: " << evtinfo.run << "\tevt: " << evtinfo.evt << "\tlumi: " << evtinfo.lumi << "\tchannel: "
458	<	<< kine.channel << "\tm4l: " << kine.m4l << "\tmZ1: " << kine.mZ1 << "\tmZ2: " << kine.mZ2 << endl;
459	<
460	<	} else {
461	<	if(ctrl.debug)
462	<	cout << "failing with code : " << ret4l.status.selectionBits << endl;
454	>	} else {
455	>	if(ctrl.debug) cout << "failing with code : " << ret4l.status.selectionBits << endl;
456	>	}
457	>	} else { // fill FOtuple and angletuple if event has a good Z and one or two extra leptons
458	>	assert(ctrl.fakeScheme=="ZPlusF" || ctrl.fakeScheme=="ZPlusFF");
459	>	if(ret4l.status.selectionBits.test(PASS_GOODZ1) ) {
460	>	unsigned nleps = passingLeptons.size() + failingLeptons.size();
461	>	if( (nleps>=3 && ctrl.fakeScheme=="ZPlusF" ) || (nleps>=4 && ctrl.fakeScheme=="ZPlusFF" ) ) {
462	>	addDummyZ2Lepotons(ret4l);
463	>	setNtupleVals(ctrl, ret4l, evtinfo, angles, kine, ji, weights, info, jetArr, metArr, simpleLeptonJets, fJetCorrector, puArr, puDArr, fJetIDMVA, vtxArr, triggerBits, hltTable, hltObjArr, fTrigObjs, &nt);
464	>	nt.Fill();
465	>	foTree->Fill();
466	>	}
467	>	} else {
468	>	if(ctrl.debug) cout << "failing PASS_GOODZ1" << endl;
469		}
470		}
471	+	}
472
473	+	if(ctrl.fakeScheme != "none") {
474	+	foTree->getFile()->cd();
475	+	foTree->getTree()->Write();
476	+	}
477		nt.WriteClose();
478
479		for(unsigned icut=0; icut<cuts.size(); icut++) {
480		map<TString,int> thisCount(*(counts[cuts[icut]]));
481	<	cout << setw(20) << cuts[icut] << "\t" << setw(9) << thisCount["all"] << "  ("
481	>	cerr << setw(20) << cuts[icut] << "\t" << setw(9) << thisCount["all"] << "  ("
482		<< setw(9) << thisCount["4e"] << setw(9) << thisCount["4m"] << setw(9) << thisCount["2e2m"] << ")" << endl;
483		}
484	<	}
484	>	}
485	>	//----------------------------------------------------------------------------------------
486	>	void addDummyZ2Lepotons(EventData &ret4l)
487	>	// for fake ntuples, we don't know which are the Z2 leptons yet, but need to fill something so we don't seg fault
488	>	{
489	>	SimpleLepton dum1,dum2;
490	>	dum1.vec.SetPtEtaPhiM(1,2,3,4);
491	>	dum2.vec.SetPtEtaPhiM(5,6,7,8);
492	>	ret4l.Z2leptons.push_back(dum1);
493	>	ret4l.Z2leptons.push_back(dum2);
494	>	}
495	>	//----------------------------------------------------------------------------------------
496	>	void setNtupleVals(ControlFlags &ctrl, EventData &ret4l, InfoStruct &evtinfo, Angles &angles, KinematicsStruct &kine,
497	>	JetInfoStruct &ji, WeightStruct &weights, EventHeader *info, Array<PFJet> *jetArr, Array<PFMet> *metArr,
498	>	vector<SimpleLepton> &simpleLeptonJets, FactorizedJetCorrector *fJetCorrector, Array<PileupInfo> *puArr,
499	>	Array<PileupEnergyDensity> *puDArr, JetIDMVA *fJetIDMVA, Array<Vertex> *vtxArr, bitset<1024> &triggerBits,
500	>	TriggerTable *hltTable, Array<TriggerObject> *hltObjArr, TriggerObjectsTable *fTrigObjs, AngleTuple *nt)
501	>	{
502	>	fillKinematics(ret4l,kine);
503	>	// fillMassErrors(ret4l,muonArr,electronArr,kine);
504	>	TLorentzVector pfmet; pfmet.SetPxPyPzE(metArr->At(0)->Mex(),metArr->At(0)->Mey(),0,0);
505	>	fillEventInfo( info, pfmet, evtinfo, ctrl.mc ? getNPU(puArr) : 0, vtxArr->GetEntries());
506	>	fillJets(ctrl, jetArr, simpleLeptonJets, fJetCorrector, puDArr, fJetIDMVA, vtxArr, ret4l);
507	>	fillJetInfo(simpleLeptonJets, ji, ctrl);
508	>	fillJetAngleBranches(ji, ret4l, nt, ctrl);
509	>	evtinfo.status = 0;
510	>	evtinfo.status = setHltBits(ti, triggerBits);
511	>
512	>	if(ctrl.fakeScheme == "none") {
513	>	fillAngles(ret4l,angles);
514	>
515	>	if( ctrl.mc) {
516	>	setEfficiencyWeights(ctrl, ctrl.era, ret4l, triggerBits, hltTable, hltObjArr, fTrigObjs, weights);
517	>	if(ctrl.debug)
518	>	cout << "w: " << weights.w << "\t" << "npuw: " << weights.npuw << "\t" << "won: " << weights.won << "\t" << "woff: " << weights.woff << endl;
519	>	}
520	>
521	>	if(ctrl.debug) {
522	>	stringstream ss;
523	>	ss << "PASS:: " << "\trun: " << evtinfo.run << "\tevt: " << evtinfo.evt << "\tlumi: " << evtinfo.lumi << "\tchannel: "
524	>	<< kine.channel << "\tm4l: " << kine.m4l << "\tmZ1: " << kine.mZ1 << "\tmZ2: " << kine.mZ2 << endl;
525	>	cerr << ss.str();
526	>	cout << ss.str();
527	>	}
528	>	} else {
529	>	// fillAngles(ret4l,angles); // at this point we don't know which four leptons will constitute the ZZ system, so doesn't make sense to fill the angles
530	>	}
531	>	}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines