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.24 by khahn, Mon May 21 22:29:53 2012 UTC vs.
Revision 1.55 by anlevin, Wed Sep 12 12:13:27 2012 UTC

#	Line 33 | Line 33 | using namespace std;
33		#include "Electron.h"
34		#include "Muon.h"
35		#include "Vertex.h"
36	+	#include "PFMet.h"
37		#include "PFCandidate.h"
38		#include "PFCandidateCol.h"
39		#include "PileupInfoCol.h"
#	Line 40 | Line 41 | using namespace std;
41		#include "MCParticle.h"
42		#include "TriggerMask.h"
43		#include "TriggerTable.h"
44	+	#include "TriggerObject.h"
45	+	#include "TriggerObjectRel.h"
46		#include "Names.h"
47		#include "BaseMod.h"
48	+	#include "TriggerObjectsTable.h"
49
50		//
51		// our headers
#	Line 50 | Line 54 | using namespace std;
54		#include "MuonSelection.h"
55		#include "ElectronSelection.h"
56		#include "IsolationSelection.h"
57	<
54	<	//#include "Selection.h"
55	<	//#include "ReferenceSelection.h"
56	<	#include "ReferenceSelection2.h"
57	>	#include "ReferenceSelection.h"
58
59		#include "TriggerUtils.h"
60		#include "PassHLT.h"
61		#include "Angles.h"
62		#include "KinematicsStruct.h"
63		#include "InfoStruct.h"
64	<	//#include "GenInfoStruct.h"
64	>	#include "GenInfoStruct.h"
65		#include "WeightStruct.h"
65	–	//#include "GlobalDataAndFuncs.h"
66	–	#include "RunLumiRangeMap.h"
67	–
66		#include "AngleTuple.h"
67		#include "FOTuple.h"
68
69	+	#include "RunLumiRangeMap.h"
70		#include "SampleWeight.h"
71		#include "EfficiencyWeightsInterface.h"
72	+	#include "SelectionFuncs.h"
73
74		#include "SimpleLepton.h"
75
#	Line 82 | Line 82 | using namespace mithep;
82		//
83		// globals
84		//----------------------------------------------------------------------------
85	<	TH1D * hpu;
85	>	TH1D *hpu_2011, *hpu_2012;
86		RunLumiRangeMap rlrm;
87		vector<SimpleLepton> failingLeptons ; // for fake estimation
88		vector<SimpleLepton> passingLeptons;      // for fake estimation
#	Line 91 | Line 91 | vector<vector<unsigned> > zcutvec;
91		vector<vector<unsigned> > zzcutvec;
92		map<unsigned,float>       evtrhoMap;
93		vector<string> cutstrs;
94	<	bool passes_HLT_MC;
94	>	bool passes_HLT_MC,passes_HLT;
95		vector<bool>   PFnoPUflag;;
96
97		//
98		// prototypes
99		//----------------------------------------------------------------------------
100	<	bool setPV(ControlFlags,const mithep::Array<mithep::Vertex>*, mithep::Vertex& );
101	<	void initPUWeights();
102	<	double getPUWeight(unsigned npu);
103	<	void setEffiencyWeights(EventData&, WeightStruct& );
104	<	void initRunLumiRangeMap();
105	<	void initEvtRhoMap(map<unsigned,float> &);
106	<	unsigned makePFnoPUArray(mithep::Array<PFCandidate> * fPFCandidates,
107	<	vector<bool> &pfNoPileUpFlag,
108	<	mithep::Array<Vertex>      * vtxArr );
100	>	void initRunLumiRangeMap(ControlFlags &ctrl);
101		//----------------------------------------------------------------------------
102
103
#	Line 117 | Line 109 | int main(int argc, char** argv)
109		{
110		cutstrs.push_back(string("skim0"));              //0
111		cutstrs.push_back(string("skim1"));              //1
112	+	cutstrs.push_back(string("skim2"));              //1
113		cutstrs.push_back(string("trigger"));            //2
114		// -------------------------------------------------
115		cutstrs.push_back(string("Z candidate"));        //3
#	Line 128 | Line 121 | int main(int argc, char** argv)
121		cutstrs.push_back(string("ZZ 20/10"));           //8
122		cutstrs.push_back(string("resonance"));          //9
123		cutstrs.push_back(string("m4l>70"));             //10
124	<	cutstrs.push_back(string("m4l>100"));            //11
124	>	cutstrs.push_back(string("m4l>100, mZ2 > 12"));                  //11
125
126
127
#	Line 173 | Line 166 | int main(int argc, char** argv)
166		//--------------------------------------------------------------------------------------------------------------
167		TChain * chain = new TChain("Events");
168		TChain * hltchain = new TChain("HLT");
169	+	TChain * runchain = new TChain(Names::gkRunTreeName);
170
171		string fname;
172	+	unsigned total_unskimmed=0;
173		if( !ctrl.inputfiles.empty() ) {
174		ifstream f(ctrl.inputfiles.c_str());
175		while (f >> fname) {
#	Line 182 | Line 177 | int main(int argc, char** argv)
177		cout << "adding inputfile : " << fname.c_str() << endl;
178		entrymap[string(fname.c_str())] = unskimmedEntries(fname.c_str());
179		cout << "unskimmed entries: " << entrymap[string(fname.c_str())] << endl;
180	+	total_unskimmed += entrymap[string(fname.c_str())];
181		chain->AddFile(fname.c_str());
182		hltchain->AddFile(fname.c_str());
183	+	runchain->AddFile(fname.c_str());
184		}
185		} else {
186		cout << "adding inputfile : " << ctrl.inputfile.c_str() << endl;
187	<	unsigned tmpent = unskimmedEntries(ctrl.inputfile.c_str());
188	<	cout << "tmpent: " << tmpent << endl;
189	<	entrymap[string(ctrl.inputfile.c_str())] = unskimmedEntries(ctrl.inputfile.c_str());
190	<	cout << "unskimmed entries: " << entrymap[string(ctrl.inputfile.c_str())] << endl;
187	>	unsigned unsk_ents = unskimmedEntries(ctrl.inputfile.c_str());
188	>	entrymap[string(ctrl.inputfile.c_str())] = unsk_ents;
189	>	cout << "unskimmed entries: " << unsk_ents << endl;
190	>	total_unskimmed += unsk_ents;
191		chain->AddFile(ctrl.inputfile.c_str());
192		hltchain->AddFile(ctrl.inputfile.c_str());
193	+	runchain->AddFile(ctrl.inputfile.c_str());
194		}
195	+	// // write the total number of unskimmed events that went into making this output file to a text file
196	+	writeEntries(ctrl,total_unskimmed);
197
198		const char * ofname;
199		if( strcmp( ctrl.outputfile.c_str(), "") ) {
#	Line 207 | Line 207 | int main(int argc, char** argv)
207		SimpleTable xstab(xspath.c_str());
208
209
210	+
211		//
212		// Setup
213		//--------------------------------------------------------------------------------------------------------------
#	Line 214 | Line 215 | int main(int argc, char** argv)
215		KinematicsStruct kinematics;
216		InfoStruct evtinfo;
217		WeightStruct weights;
218	<	//  GenInfoStruct geninfo;
218	>	GenInfoStruct geninfo;
219
220		AngleTuple nt( (const char*)ofname, (const char*)"zznt");
221		nt.makeAngleBranch(angles);
222		nt.makeKinematicsBranch(kinematics);
223		nt.makeInfoBranch(evtinfo);
224
225	+
226		FOTuple foTree( nt.getFile(), (const char*)"FOtree");
227	+	foTree.makeInfoBranch(evtinfo);
228		foTree.makeFailedLeptonBranch(failingLeptons);
229	<	foTree.makeZLeptonBranch(passingLeptons);
229	>	foTree.makePassedLeptonBranch(passingLeptons);
230
231	<	TH1F * h_w_hpt;
231	>	TH1F * h_wf11_hpt;
232		if(ctrl.mc) {
233		nt.makeWeightBranch(weights);
234	<	//    nt.makeGenInfoBranch(geninfo);
235	<	initEfficiencyWeights();
234	>	if(ctrl.fillGen ) nt.makeGenInfoBranch(geninfo);
235	>	//initEfficiencyWeights();
236		initPUWeights();
237	<	/*
237	>
238		// Higgs only, pt reweighting
239		if( ctrl.mH <= 140 ) {
240		char wptname[256];
#	Line 239 | Line 242 | int main(int argc, char** argv)
242		TFile * f_hpt = new TFile(wptname);
243		f_hpt->Print();
244		sprintf(wptname, "weight_hqt_fehipro_fit_%i", ctrl.mH);
245	<	h_w_hpt  = (TH1F*)(f_hpt->FindObjectAny(wptname));
243	<	h_w_hpt->Print();
245	>	h_wf11_hpt  = (TH1F*)(f_hpt->FindObjectAny(wptname));
246		}
247	<	*/
247	>
248		} else {
249	<	initRunLumiRangeMap();
249	>	initRunLumiRangeMap(ctrl);
250		}
251
252		//  initMuonIDMVA();
253	<	initMuonIsoMVA();
252	<	initElectronIDMVA();
253	<	initElectronIsoMVA();
253	>	initElectronIDMVAV1();
254		initTrigger();
255	+	TrigInfo ti;
256	+	initAnalysisTriggers(ti);
257
256	–	// hack
257	–	initEvtRhoMap(evtrhoMap);
258	–
259	–
258		//
259		// Setup tree I/O
260		//--------------------------------------------------------------------------------------------------------------
#	Line 264 | Line 262 | int main(int argc, char** argv)
262		TTree *eventTree=0;
263		string currentFile("");
264
265	<	mithep::EventHeader *info    = new mithep::EventHeader();
266	<	//  mithep::TGenInfo    *ginfo  = new mithep::TGenInfo();
265	>	UInt_t fNMaxTriggers = 1024;
266	>	mithep::EventHeader                         *info          = new mithep::EventHeader();
267	>	mithep::Array<mithep::PFMet>                *metArr        = new mithep::Array<mithep::PFMet>();
268		mithep::Array<mithep::Electron>             *electronArr   = new mithep::Array<mithep::Electron>();
269		mithep::Array<mithep::Muon>                 *muonArr       = new mithep::Array<mithep::Muon>();
270		mithep::Array<mithep::Vertex>               *vtxArr        = new mithep::Array<mithep::Vertex>();
#	Line 274 | Line 273 | int main(int argc, char** argv)
273		mithep::Array<mithep::PileupEnergyDensity>  *puDArr        = new mithep::Array<mithep::PileupEnergyDensity>();
274		mithep::Array<mithep::Track>                *trkArr        = new mithep::Array<mithep::Track>();
275		mithep::Array<mithep::MCParticle>           *mcArr         = new mithep::Array<mithep::MCParticle>();
276	+	mithep::MCEventInfo                         *mcEvtInfo     = new mithep::MCEventInfo();
277		mithep::TriggerMask                         *trigMask      = new mithep::TriggerMask();
278	<	mithep::TriggerTable                        *hltTable      = new mithep::TriggerTable();
278	>	mithep::TriggerTable                        *hltTable      = new mithep::TriggerTable(fNMaxTriggers);
279		vector<string>                              *hltTableStrings  = new vector<string>();
280	+	vector<string>                              *hltLabelStrings  = new vector<string>();
281	+	mithep::Array<mithep::TriggerObject>        *hltObjArr     = new mithep::Array<mithep::TriggerObject>();
282	+	mithep::Array<mithep::TriggerObjectRel>     *hltRelsArr    = new mithep::Array<mithep::TriggerObjectRel>();
283	+	std::vector<std::string>                    *hltTab        = new vector<string>();
284	+	std::vector<std::string>                    *hltLab        = new vector<string>();
285	+	TriggerObjectsTable                         *fTrigObjs     = new TriggerObjectsTable(hltTable,fNMaxTriggers);
286	+
287	+	RunInfo                                     *runInfo       = new RunInfo();
288
289		TString fElectronName(Names::gkElectronBrn);
290		TString fMuonName(Names::gkMuonBrn);
291		TString fInfoName(Names::gkEvtHeaderBrn);
292	+	TString fPFMetName("PFMet");
293		TString fPrimVtxName(Names::gkPVBrn);
294		TString fPFCandidateName(Names::gkPFCandidatesBrn);
295		TString fPileupInfoName(Names::gkPileupInfoBrn);
296		TString fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn);
297		TString fTrackName(Names::gkTrackBrn);
298		TString fMCParticleName(Names::gkMCPartBrn);
299	+	TString fMCEvtInfoName(Names::gkMCEvtInfoBrn);
300		TString fTriggerMaskName(Names::gkHltBitBrn);
301		TString fTriggerTableName(Names::gkHltTableBrn);
302	+	TString fTriggerLabelName(Names::gkHltLabelBrn);
303	+	TString fTriggerObjectName(Names::gkHltObjBrn);
304	+	TString fTriggerObjectRelsName(Form("HLTObjectsRelation"));
305	+
306	+	//   TString fHLTLabName(Names::gkHltLabelBrn);
307	+	//   TString fHLTTabName(Names::gkHltTableBrn);
308	+	//   TString fHLTTabNamePub(Form("%sFwk",fHLTTabName.Data()));
309	+	//   TString fHLTLabNamePub(Form("%sFwk",fHLTLabName.Data()));
310	+	//   TString fObjsNamePub(Form("%sFwk",fTriggerObjectName.Data()));
311	+
312	+
313
314		chain->SetBranchAddress(fInfoName,        &info);
315	+	chain->SetBranchAddress(fPFMetName,       &metArr);
316		chain->SetBranchAddress(fElectronName,    &electronArr);
317		chain->SetBranchAddress(fMuonName,        &muonArr);
318		chain->SetBranchAddress(fPrimVtxName,     &vtxArr);
#	Line 298 | Line 320 | int main(int argc, char** argv)
320		if( ctrl.mc ) {
321		chain->SetBranchAddress(fPileupInfoName,  &puArr);
322		chain->SetBranchAddress(fMCParticleName,  &mcArr);
323	+	chain->SetBranchAddress(fMCEvtInfoName,  &mcEvtInfo);
324		}
325		chain->SetBranchAddress(fPileupEnergyDensityName, &puDArr);
326		chain->SetBranchAddress(fTrackName, &trkArr);
327		chain->SetBranchAddress(fTriggerMaskName, &trigMask);
328	+	chain->SetBranchAddress(fTriggerObjectName,  &hltObjArr);
329	+	chain->SetBranchAddress(fTriggerObjectRelsName,  &hltRelsArr);
330	+	//  chain->SetBranchAddress(fHLTTabNamePub,  &hltTab);
331	+	//  chain->SetBranchAddress(fHLTLabNamePub,  &hltLab);
332	+
333		cout << "hlttable: " << fTriggerTableName << endl;
334		hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
335	+	hltchain->SetBranchAddress(fTriggerLabelName, &hltLabelStrings);
336	+
337	+	runchain->SetBranchAddress(Names::gkRunInfoBrn, &runInfo);
338
339		mithep::Vertex              vtx;          // best primary vertex in the event
340
#	Line 318 | Line 349 | int main(int argc, char** argv)
349
350		// only 1 HLT table / file ???
351		hltchain->GetEntry(0);
321	–	cerr << "here... " << endl;
352
353	<	int imax = chain->GetEntries();
354	<	cout << "nEntries: " << imax << endl;
353	>	int imax;
354	>	//-1 means process all of the events
355	>	if(ctrl.n_events_to_process != -2 && ctrl.n_events_to_process != -1){
356	>	if(ctrl.n_events_to_process > chain->GetEntries()){
357	>	cerr << "error: n_events_to_process is larger than chain->GetEntries()" << endl;
358	>	return 1;
359	>	}
360	>	imax = ctrl.n_events_to_process;
361	>	}
362	>	else
363	>	imax = chain->GetEntries();
364
365	+	cout << "number of events to process: " << imax << endl;
366
367		//
368		// Loop !!!!!!!!!
369		//--------------------------------------------------------------------------------------------------------------
370	+	int lastHltEntry=-2;
371		for(UInt_t ientry=0; ientry<imax; ientry++)
372		{
373		chain->GetEntry(ientry);
374	+	if(!(ientry%1000)) cerr << "entry: " << ientry << endl;
375
376	<	/*
377	<	bool found_tau=false;
378	<	if( ctrl.mc) {  // && ctrl.debug
379	<	cout << "nparticles: " << mcArr->GetEntries() << endl;
380	<	cout << "-----------------------------------" << endl;
381	<	for( int i=0; i<mcArr->GetEntries(); i++ ) {
382	<	const mithep::MCParticle *mc = mcArr->At(i);
383	<	cout << "i: " << i << "\t"
384	<	<< "id: " << mc->AbsPdgId()
385	<	<< endl;
386	<	if( mc->Is(MCParticle::kTau) ) {
387	<	found_tau=true;
346	<	}
347	<	}
376	>	if( ctrl.debug ) {
377	>	cout << "-----------------------------------------------------------------" << endl;
378	>	cout << "-----------------------------------------------------------------" << endl;
379	>	cout << "Run: " << info->RunNum()
380	>	<< "\tEvt: " << info->EvtNum()
381	>	<< "\tLumi: " << info->LumiSec()
382	>	<< endl;
383	>	cerr << "Run: " << info->RunNum()
384	>	<< "\tEvt: " << info->EvtNum()
385	>	<< "\tLumi: " << info->LumiSec()
386	>	<< endl;
387	>	cout << "-----------------------------------------------------------------" << endl;
388		}
349	–	cout << "-----------------------------------" << endl;
350	–	cout << "found tau: " << found_tau << endl;
351	–	*/
352	–
353	–	// events that fail bc ele overlaps w/ reco mu
354	–	/*
355	–	if( info->EvtNum() != 67315 &&
356	–	info->EvtNum() != 288693 ) continue;
357	–	*/
358	–
359	–	// events where they apparently don't veto electrons ...
360	–	/*
361	–	if( info->EvtNum() != 141713 &&
362	–	info->EvtNum() != 178019 &&
363	–	info->EvtNum() != 60823 &&
364	–	info->EvtNum() != 192795 ) continue;
365	–	*/
366	–
367	–	if(!(ientry%1000)) cerr << "entry: " << ientry << endl;
368	–	setPV( ctrl, vtxArr, vtx);
369	–	//if (!(setPV( ctrl, vtxArr, vtx)) ) continue;
389
390	<	float rho = evtrhoMap[info->EvtNum()];
390	>	runchain->GetEvent(info->RunEntry());
391	>	hltchain->GetEntry(runInfo->HltEntry());
392	>	hltTable->Clear();
393	>	fillTriggerNames(hltTable, hltTableStrings );
394	>	if(ctrl.debug && (lastHltEntry != runInfo->HltEntry())) hltTable->Print();
395	>	lastHltEntry = runInfo->HltEntry();
396	>	fillTriggerBits( hltTable, trigMask, triggerBits );
397	>	setHLTObjectRelations( hltObjArr, hltRelsArr, hltTableStrings, hltLabelStrings, fTrigObjs );
398	>	// printTriggerObjs( hltTable, fTrigObjs);
399	>
400	>	string fname = string(chain->GetFile()->GetEndpointUrl()->GetFile());
401	>	if(ctrl.debug) cout << "era is " << ctrl.era << endl;
402	>	if( ctrl.era == 0  ) {
403	>	setEra( fname, ctrl );
404	>	if(ctrl.debug) cout << "era is now " << ctrl.era << endl;
405	>	}
406
407		//
408		// pfNoPU
#	Line 379 | Line 413 | int main(int argc, char** argv)
413		assert(pfnopu_size == pfArr->GetEntries());
414
415		//
416	+	// trigger
417	+	//
418	+	if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) {
419	+	currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile());
420	+	}
421	+	if( ctrl.debug ) cout << "file is : " << currentFile  << endl;
422	+
423	+	//
424		// data/MC
425		//
426		if(ctrl.mc) {
427		//
428	+	// gen info
429	+	//
430	+	if( ctrl.fillGen )
431	+	fillGenInfo( mcArr, mcEvtInfo, geninfo, ESampleType::kHZZ, ctrl);
432	+
433	+	//
434		// xsec & PU weights
435		//
436		weights.w = getWeight(xstab,entrymap,chain)/ctrl.totalMC;
437		cout << "xsec weight: " << weights.w << endl;
438	<	int npu = -1;
439	<	for(int i=0;i<puArr->GetEntries();i++) {
440	<	if(puArr->At(i)->GetBunchCrossing() == 0)
441	<	npu = puArr->At(i)->GetPU_NumInteractions();
442	<	}
443	<	assert(npu>=0);
396	<	evtinfo.npu;
397	<	weights.npuw = getPUWeight(evtinfo.npu);
398	<	//      cout << "weight: " << weights.w << endl;
438	>
439	>	unsigned npu = getNPU(puArr, 0);
440	>	weights.npu  = npu;
441	>	weights.npuw = getPUWeight(ctrl.era,npu);
442	>	cerr << "weights.npuw: " <<  weights.npuw << endl;
443	>
444		//
445	<	// trigger
445	>	// pt reweighting for Higgs < 140, F11
446		//
447	<	if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) {
448	<	currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile());
449	<	hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
450	<	hltchain->GetEntry(0);
451	<	hltTable->Clear();
452	<	fillTriggerBits(hltTable, hltTableStrings );
447	>	if( ctrl.fillGen ) {
448	>	if( (fname.find("f11-h115zz4l") != string::npos) ||
449	>	(fname.find("f11-h120zz4l") != string::npos) ||
450	>	(fname.find("f11-h130zz4l") != string::npos)  ) {
451	>	weights.w *= h_wf11_hpt->GetBinContent(h_wf11_hpt->FindBin(geninfo.pt_zz));
452	>	}
453		}
454	<	if( ctrl.debug ) cout << "file is : " << currentFile  << endl;
455	<	fillTriggerBits( hltTable, trigMask, triggerBits );
454	>
455	>	// trigger
456		passes_HLT_MC = true; // passed to apply as extern global, just for sync
457	<	if( !passHLTMC(triggerBits, info->RunNum(), info->EvtNum() ) ) {
457	>	if( !passHLTMC(triggerBits, info->RunNum(), info->EvtNum(), k2012_MC ) ) {
458		passes_HLT_MC = false;
459		}
460	+
461		} else {
462		//
463		// JSON
464		//
465	<	/*
466	<	RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());
467	<	if( !(rlrm.HasRunLumi(rl)) )  {
468	<	if( ctrl.debug ) cout << "\tfails JSON" << endl;
469	<	continue;
465	>	if(!(ctrl.noJSON) ) {
466	>	RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());
467	>	if( !(rlrm.HasRunLumi(rl)) )  {
468	>	if( ctrl.debug ) cout << "fails JSON" << endl;
469	>	continue;
470	>	}
471		}
472	<	*/
426	<	//
472	>
473		// trigger
474	<	//
475	<	if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) {
476	<	currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile());
477	<	hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
478	<	hltchain->GetEntry(0);
433	<	hltTable->Clear();
434	<	fillTriggerBits(hltTable, hltTableStrings );
435	<	}
436	<	if( ctrl.debug ) cout << "file is : " << currentFile  << endl;
437	<	/*
438	<	fillTriggerBits( hltTable, trigMask, triggerBits );
439	<	if( !passHLT(triggerBits, info->RunNum(), info->EvtNum() ) ) {
440	<	if( ctrl.debug ) cout << "\tfails trigger ... " << endl;
441	<	continue;
442	<	}
443	<	*/
474	>	passes_HLT = true;
475	>	//if( !passHLT(ctrl, triggerBits,ti,info->RunNum()) ) {
476	>	//  passes_HLT = false;
477	>	//}
478	>
479		}
480
481		//
#	Line 452 | Line 487 | int main(int argc, char** argv)
487
488
489		// reference
490	<	apply_HZZ4L_reference2_selection(ctrl, info,
491	<	vtx,
492	<	pfArr,
493	<	#ifdef HACKED_RHOS
494	<	rho,
495	<	#else
496	<	puDArr,
497	<	#endif
498	<	electronArr,
499	<	&electronReferencePreSelection,
500	<	&electronReferenceIDMVASelection,
501	<	&electronReferenceIsoSelection,
502	<	muonArr,
468	<	&muonReferencePreSelection,
469	<	&muonIDPFSelection,
470	<	&muonReferenceIsoSelection);
471	<
472	<	/*
473	<	// evolved
474	<	apply_HZZ4L_reference2_selection(ctrl, info,
475	<	vtx,
476	<	pfArr,
477	<	//rho,
478	<	puDArr,
479	<	electronArr,
480	<	&electronReferencePreSelection,
481	<	//&electronPreSelection,
482	<	&electronReferenceIDMVASelection,
483	<	&electronIsoMVASelection,
484	<	muonArr,
485	<	&muonReferencePreSelection,
486	<	//&muonPreSelection,
487	<	&muonIDPFSelection,
488	<	&muonIsoMVASelection);
489	<	*/
490	<
490	>	apply_HZZ4L_reference_selection(ctrl, info,
491	>	vtxArr,
492	>	pfArr,
493	>	puDArr,
494	>	electronArr,
495	>	&electronReferencePreSelection,
496	>	&electronReferenceIDMVASelectionV1,
497	>	&electronReferenceIsoSelection,
498	>	muonArr,
499	>	&muonReferencePreSelection,
500	>	&muonIDPFSelection,
501	>	&muonReferenceIsoSelection);
502	>
503		if( ctrl.debug ) cout << endl;
504
505
#	Line 501 | Line 513 | int main(int argc, char** argv)
513		if( abs(ret4l.Z2leptons[0].type) == 13  ) Z2type = 13;
514		}
515
504	–	#ifdef SYNC
505	–	cout  << "bits: " << ret4l.status.selectionBits  << "\t"
506	–	<< "Z1: " << Z1type << "\t"
507	–	<< "Z2: " << Z2type << endl;
508	–	cout << endl;
509	–	#endif
510	–
516		if( ret4l.status.pass() ) {
517		fillAngles(ret4l,angles);
518		fillKinematics(ret4l,kinematics);
519	<	fillEventInfo(info,evtinfo);
519	>	fillMassErrors(ret4l,muonArr,electronArr,kinematics);
520	>	TLorentzVector pfmet; pfmet.SetPxPyPzE(metArr->At(0)->Mex(),metArr->At(0)->Mey(),0,0);
521	>	fillEventInfo( info, pfmet, evtinfo, ctrl.mc ? getNPU(puArr) : 0);
522		if( ctrl.mc) {
523	<	// fillGenInfo(ginfo,geninfo);
524	<	setEffiencyWeights(ret4l, weights);
525	<	if(ctrl.debug)
526	<	cout << "w: " << weights.w << "\t"
527	<	<< "won: " << weights.won << "\t"
528	<	<< "woff: " << weights.woff << endl;
523	>	//setEffiencyWeights(ctrl.era, ret4l, triggerBits, hltTable, hltObjArr, fTrigObjs, weights);
524	>	if(ctrl.debug)
525	>	cout << "w: " << weights.w << "\t"
526	>	<< "npuw: " << weights.npuw << "\t"
527	>	<< "won: " << weights.won << "\t"
528	>	<< "woff: " << weights.woff << endl;
529		}
530
531	<	/*
525	<	// only for Higgs < 140
526	<	geninfo.weight *= h_w_hpt->GetBinContent(h_w_hpt->FindBin(geninfo.pt_zz));
527	<	angles.bdt = reader->EvaluateMVA("BDTG");
528	<	*/
531	>
532		nt.Fill();
533
534		cerr  << "PASS:: "
535		<< "\trun: " << evtinfo.run
536		<< "\tevt: " << evtinfo.evt
537		<< "\tlumi: " << evtinfo.lumi
538	<	<< "\tcostheta1: " << angles.costheta1
539	<	<< "\tcostheta2: " << angles.costheta2
540	<	<< "\tcostheta*: " << angles.costhetastar
541	<	<< "\tbdt: " << angles.bdt
538	>	<< "\tchannel: " << kinematics.channel
539	>	<< "\tm4l: " << kinematics.m4l
540	>	<< "\tmZ1: " << kinematics.mZ1
541	>	<< "\tmZ2: " << kinematics.mZ2
542		<< endl;
543		pass++;
544		//      if( pass > 3 ) break;
545
546	<	} else if( ret4l.status.selectionBits.test(PASS_GOODZ1) ) { // save for fakes ...
547	<	//      if(ctrl.debug) {
546	>	} else if( ret4l.status.selectionBits.test(PASS_GOODZ1) &&
547	>	ret4l.status.selectionBits.to_ulong() < (1UL<<PASS_4L)) { // save for fakes ...
548	>	if(ctrl.debug) {
549		cout << "failingLeptons : " << failingLeptons.size() << endl;
550		for( int f=0; f<failingLeptons.size(); f++ ) {
551		SimpleLepton  l = failingLeptons[f];
#	Line 550 | Line 554 | int main(int argc, char** argv)
554		<< "pT: " << l.vec.Pt()  << "\t"
555		<< "eta: " << l.vec.Eta()
556		<< endl;
557	<	//    }
557	>	}
558		}
559	+	TLorentzVector pfmet; pfmet.SetPxPyPzE(metArr->At(0)->Mex(),metArr->At(0)->Mey(),0,0);
560	+	fillEventInfo( info, pfmet, evtinfo, ctrl.mc ? getNPU(puArr) : 0);
561		foTree.Fill();
562	+	} else {
563	+	if(ctrl.debug)
564	+	cout << "failing with some other code : " << ret4l.status.selectionBits << endl;
565		}
566	<	}
566	>	if(ctrl.debug)
567	>	cout << "passes HLT: " << passes_HLT << endl;
568	>	}
569
570
571		foTree.getFile()->cd();
572		foTree.getTree()->Write();
573		nt.WriteClose();
574
575	<	cout << endl;
576	<	cout << endl;
575	>	if(ctrl.debug)
576	>	cout << endl << endl;
577		for( int i=0; i<cutvec.size(); i++ ) {
578		cout << "cut: " << i << "\t"
579		<< "pass: " << cutvec[i];
#	Line 583 | Line 594 | int main(int argc, char** argv)
594
595		}
596
586	–	//----------------------------------------------------------------------------
587	–	//
588	–	// Get primary vertices
589	–	// Assumes primary vertices are ordered by sum-pT^2 (as should be in CMSSW)
590	–	// NOTE: if no PV is found from fitting tracks, the beamspot is used
591	–	//
592	–	//----------------------------------------------------------------------------
593	–	bool setPV(ControlFlags ctrl,
594	–	const mithep::Array<mithep::Vertex> * vtxArr,
595	–	mithep::Vertex & vtx)
596	–	//----------------------------------------------------------------------------
597	–	{
598	–	const Vertex *bestPV = 0;
599	–	float best_sumpt=-1;
600	–
601	–	// good PV requirements
602	–	const UInt_t   fMinNTracksFit = 0;
603	–	const Double_t fMinNdof       = 4;
604	–	const Double_t fMaxAbsZ       = 24;
605	–	const Double_t fMaxRho        = 2;
606	–
607	–	for(int i=0; i<vtxArr->GetEntries(); ++i) {
608	–	const Vertex *pv = (Vertex*)(vtxArr->At(i));
609	–	if( ctrl.debug ) cout << "vertex :: " << i << "\tntrks: " << pv->NTracks() << endl;
610	–
611	–	// Select best PV for corrected d0; if no PV passing cuts, the first PV in the collection will be used
612	–	if(!pv->IsValid())                                continue;
613	–	if(pv->NTracksFit()       < fMinNTracksFit)       continue;
614	–	if(pv->Ndof()                 < fMinNdof)         continue;
615	–	if(fabs(pv->Z()) > fMaxAbsZ)                      continue;
616	–	if(pv->Position().Rho()   > fMaxRho)              continue;
617	–
618	–	// take the first ...
619	–	bestPV = pv;
620	–	break;
621	–
622	–	// this never reached ...
623	–	float tmp_sumpt=0;
624	–	for( int t=0; t<pv->NTracks(); t++ )
625	–	tmp_sumpt += pv->Trk(t)->Pt();
626	–
627	–	if( tmp_sumpt > best_sumpt  ) {
628	–	bestPV = pv;
629	–	best_sumpt = tmp_sumpt;
630	–	if( ctrl.debug) cout << "new PV set, pt : " << best_sumpt << endl;
631	–	}
632	–	}
633	–	// sync
634	–	if(!bestPV) bestPV = vtxArr->At(0);
635	–	//if(!bestPV) return false;
636	–	vtx.SetPosition(bestPV->X(),bestPV->Y(),bestPV->Z());
637	–	vtx.SetErrors(bestPV->XErr(),bestPV->YErr(),bestPV->ZErr());
638	–	return true;
639	–	};
640	–
641	–
642	–
643	–	//----------------------------------------------------------------------------
644	–	void setEffiencyWeights(EventData & evtdat, WeightStruct & weights )
645	–	//----------------------------------------------------------------------------
646	–	{
647	–	vector<SimpleLepton> lepvec = evtdat.Z1leptons;
648	–	lepvec.insert(lepvec.end(), evtdat.Z2leptons.begin(), evtdat.Z2leptons.end());
649	–	double w_offline=-1, werr_offline=0;
650	–	double w_online=-1, werr_online=0;
651	–
652	–	vector< pair <double,double> > wlegs; // pair here is eff & err
653	–	vector< pair <float,float> > mvec;  // pair here is eta & pt
654	–
655	–	//      vector< pair <float,float> > evec;  // pair here is eta & pt
656	–	// now deal with medium vs loose
657	–	vector< pair< bool, pair <float,float> > > evec;  // pair here is eta & pt
658	–
659	–	for( int k=0; k<lepvec.size(); k++ ) {
660	–	if( abs(lepvec[k].type) == 13 ) {
661	–	mvec.push_back( std::pair<float,float> (fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt()) );
662	–	wlegs.push_back( muonPerLegOfflineEfficiencyWeight( fabs(lepvec[k].vec.Eta()),
663	–	lepvec[k].vec.Pt() ) );
664	–	} else {
665	–
666	–	// now deal with medium vs loose
667	–	//              evec.push_back( std::pair<float,float> (fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt()) );
668	–
669	–	std::pair<float,float> tmppair(fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt());
670	–	evec.push_back( std::pair<bool, std::pair<float,float> > (lepvec[k].isTight, tmppair) );
671	–
672	–	//              wlegs.push_back( elePerLegOfflineEfficiencyWeight(  fabs(lepvec[k].vec.Eta()),
673	–	//                                                                 lepvec[k].vec.Pt() ) );
674	–
675	–	wlegs.push_back( elePerLegOfflineEfficiencyWeight(  fabs(lepvec[k].vec.Eta()),
676	–	lepvec[k].vec.Pt(),
677	–	lepvec[k].isTight ) );
678	–
679	–	}
680	–	}
681	–
682	–	pair<double,double> offpair = getOfflineEfficiencyWeight( wlegs );
683	–	weights.woff    = offpair.first;
684	–	weights.werroff = offpair.second;
685	–
686	–	pair<double,double> onpair  = getOnlineEfficiencyWeight( mvec, evec );
687	–	weights.won    = onpair.first;
688	–	weights.werron = onpair.second;
689	–	}
690	–
597
598		//----------------------------------------------------------------------------
599	<	void initRunLumiRangeMap()
599	>	void initRunLumiRangeMap(ControlFlags &ctrl)
600		//----------------------------------------------------------------------------
601		{
602		/*
#	Line 702 | Line 608 | void initRunLumiRangeMap()
608		// r11b
609		rlrm.AddJSONFile(std::string("./data/Cert_160404-180252_7TeV_PromptReco_Collisions11_JSON.txt"));
610		*/
705	–	};
706	–
707	–
708	–	//----------------------------------------------------------------------------
709	–	void initPUWeights()
710	–	//----------------------------------------------------------------------------
711	–	{
712	–	TFile * puf = new TFile("data/PileupReweighting.Summer11DYmm_To_Full2011.root");
713	–	hpu = (TH1D*)(puf->Get("puWeights"));
714	–	}
715	–
716	–	//----------------------------------------------------------------------------
717	–	double getPUWeight(unsigned npu)
718	–	//----------------------------------------------------------------------------
719	–	{
720	–	return hpu->GetBinContent(hpu->FindBin(npu));
721	–	}
722	–
723	–	//----------------------------------------------------------------------------
724	–	void initEvtRhoMap( map<unsigned, float> & evtrhoMap )
725	–	//----------------------------------------------------------------------------
726	–	{
727	–	unsigned evt;
728	–	float rho;
729	–
730	–	cout << "initialzing EvtRhoMap ";
731	–	//FILE * f = fopen("evtrho.dat", "r");
732	–	//  FILE * f = fopen("allrho.cali.uniq.dat", "r");
733	–	FILE * f = fopen("allrhoAA.cali.uniq.dat", "r");
734	–	int lcount=0;
735	–	while( fscanf( f, "%u:%f", &evt, &rho ) ) {
736	–	if( feof(f) ) break;
737	–	evtrhoMap[evt] = rho;
738	–	lcount++;
739	–	if( !(lcount%1000) ) cout << "."; flush(cout);
740	–	}
741	–	cout << " done" << endl;
742	–	}
743	–
744	–	//----------------------------------------------------------------------------
745	–	unsigned makePFnoPUArray(mithep::Array<PFCandidate> * fPFCandidates,
746	–	vector<bool> &pfNoPileUpflag,
747	–	mithep::Array<Vertex>      * vtxArr )
748	–	//----------------------------------------------------------------------------
749	–	{
750	–	for(UInt_t i = 0; i < fPFCandidates->GetEntries(); i++) {
751	–	const PFCandidate *pf = fPFCandidates->At(i);
752	–	assert(pf);
753	–
754	–	if(pf->PFType() == PFCandidate::eHadron) {
755	–	if(pf->HasTrackerTrk() &&
756	–	vtxArr->At(0)->HasTrack(pf->TrackerTrk()) &&
757	–	vtxArr->At(0)->TrackWeight(pf->TrackerTrk()) > 0) {
758	–
759	–	pfNoPileUpflag.push_back(1);
760	–
761	–	} else {
762	–
763	–	Bool_t vertexFound = kFALSE;
764	–	const Vertex *closestVtx = 0;
765	–	Double_t dzmin = 10000;
766	–
767	–	// loop over vertices
768	–	for(UInt_t j = 0; j < vtxArr->GetEntries(); j++) {
769	–	const Vertex *vtx = vtxArr->At(j);
770	–	assert(vtx);
771	–
772	–	if(pf->HasTrackerTrk() &&
773	–	vtx->HasTrack(pf->TrackerTrk()) &&
774	–	vtx->TrackWeight(pf->TrackerTrk()) > 0) {
775	–	vertexFound = kTRUE;
776	–	closestVtx = vtx;
777	–	break;
778	–	}
779	–	Double_t dz = fabs(pf->SourceVertex().Z() - vtx->Z());
780	–	if(dz < dzmin) {
781	–	closestVtx = vtx;
782	–	dzmin = dz;
783	–	}
784	–	}
785	–
786	–	//      if(fCheckClosestZVertex) {
787	–	if(1) {
788	–	// Fallback: if track is not associated with any vertex,
789	–	// associate it with the vertex closest in z
790	–	if(vertexFound || closestVtx != vtxArr->At(0)) {
791	–	//      pfPileUp->Add(pf);
792	–	pfNoPileUpflag.push_back(0);
793	–	} else {
794	–	pfNoPileUpflag.push_back(1);
795	–	}
796	–	} else {
797	–	if(vertexFound && closestVtx != vtxArr->At(0)) {
798	–	//      pfPileUp->Add(pf);
799	–	pfNoPileUpflag.push_back(0);
800	–	} else {
801	–	//      PFCandidate * pfNoPileUp->AddNew(); // Ridiculous but that's how it is
802	–	pfNoPileUpflag.push_back(1);
803	–	}
804	–	}
805	–	} //hadron & trk stuff
806	–	} else { // hadron
807	–	//      PFCandidate * ptr = pfNoPileUp->AddNew();
808	–	pfNoPileUpflag.push_back(1);
809	–	}
810	–	} // Loop over PF candidates
811	–
812	–	return pfNoPileUpflag.size();
813	–	}
611
612	+	// 2012 only for now ...
613	+	std::string jsonFile;
614	+	if(ctrl.jsonFile=="")
615	+	jsonFile = "./data/Cert_190456-194479_8TeV_PromptReco_Collisions12_JSON.txt";
616	+	else
617	+	jsonFile = ctrl.jsonFile;
618	+	rlrm.AddJSONFile(jsonFile);
619	+	};

Diff Legend

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