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.13 by khahn, Sat May 5 11:04:01 2012 UTC vs.
Revision 1.22 by khahn, Mon May 14 17:27:34 2012 UTC

#	Line 5 | Line 5
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		//
#	Line 34 | Line 35 | using namespace std;
35		#include "Vertex.h"
36		#include "PFCandidate.h"
37		#include "PFCandidateCol.h"
38	+	#include "PileupInfoCol.h"
39	+	#include "PileupEnergyDensity.h"
40		#include "TriggerMask.h"
41		#include "TriggerTable.h"
42		#include "Names.h"
#	Line 46 | Line 49 | using namespace std;
49		#include "MuonSelection.h"
50		#include "ElectronSelection.h"
51		#include "IsolationSelection.h"
52	<	#include "Selection.h"
52	>
53	>	//#include "Selection.h"
54	>	//#include "ReferenceSelection.h"
55	>	#include "ReferenceSelection2.h"
56	>
57		#include "TriggerUtils.h"
58		#include "PassHLT.h"
59		#include "Angles.h"
#	Line 54 | Line 61 | using namespace std;
61		#include "InfoStruct.h"
62		//#include "GenInfoStruct.h"
63		#include "WeightStruct.h"
64	<	#include "GlobalDataAndFuncs.h"
64	>	//#include "GlobalDataAndFuncs.h"
65	>	#include "RunLumiRangeMap.h"
66	>
67		#include "AngleTuple.h"
68	+	#include "FOTuple.h"
69
70		#include "SampleWeight.h"
71		#include "EfficiencyWeightsInterface.h"
72
73	+	#include "SimpleLepton.h"
74	+
75		#ifndef CMSSW_BASE
76		#define CMSSW_BASE "../"
77		#endif
78
79		using namespace mithep;
80
81	+	//
82	+	// globals
83	+	//----------------------------------------------------------------------------
84	+	TH1D * hpu;
85	+	RunLumiRangeMap rlrm;
86	+	vector<SimpleLepton> failingLeptons ; // for fake estimation
87	+	vector<SimpleLepton> passingLeptons;      // for fake estimation
88	+	vector<unsigned> cutvec;
89	+	vector<vector<unsigned> > zcutvec;
90	+	vector<vector<unsigned> > zzcutvec;
91	+	map<unsigned,float>       evtrhoMap;
92	+	vector<string> cutstrs;
93	+	bool passes_HLT_MC;
94	+	vector<bool>   PFnoPUflag;;
95
96	+	//
97	+	// prototypes
98		//----------------------------------------------------------------------------
99		bool setPV(ControlFlags,const mithep::Array<mithep::Vertex>*, mithep::Vertex& );
100	+	void initPUWeights();
101	+	double getPUWeight(unsigned npu);
102	+	void setEffiencyWeights(EventData&, WeightStruct& );
103	+	void initRunLumiRangeMap();
104	+	void initEvtRhoMap(map<unsigned,float> &);
105	+	unsigned makePFnoPUArray(mithep::Array<PFCandidate> * fPFCandidates,
106	+	vector<bool> &pfNoPileUpFlag,
107	+	mithep::Array<Vertex>      * vtxArr );
108		//----------------------------------------------------------------------------
109
110	<	//=== MAIN =================================================================================================
110	>
111	>	//
112	>	// MAIN
113	>	//----------------------------------------------------------------------------
114		int main(int argc, char** argv)
115	+	//----------------------------------------------------------------------------
116		{
117	+	cutstrs.push_back(string("skim0"));              //0
118	+	cutstrs.push_back(string("skim1"));              //1
119	+	cutstrs.push_back(string("trigger"));            //2
120	+	// -------------------------------------------------
121	+	cutstrs.push_back(string("Z candidate"));        //3
122	+	cutstrs.push_back(string("good Z1"));            //4
123	+	// -------------------------------------------------
124	+	cutstrs.push_back(string("4l"));                 //5
125	+	cutstrs.push_back(string("ZZ candidate"));       //6
126	+	cutstrs.push_back(string("good Z2"));            //7
127	+	cutstrs.push_back(string("ZZ 20/10"));           //8
128	+	cutstrs.push_back(string("resonance"));          //9
129	+	cutstrs.push_back(string("m4l>70"));             //10
130	+	cutstrs.push_back(string("m4l>100"));            //11
131	+
132	+
133	+
134	+
135		string cmsswpath(CMSSW_BASE);
136		cmsswpath.append("/src");
137		std::bitset<1024> triggerBits;
138		vector<vector<string> > inputFiles;
139		inputFiles.push_back(vector<string>());
140	<	map<string,double> entrymap; // number of unskimmed entries in each file
140	>	map<string,unsigned> entrymap; // number of unskimmed entries in each file
141	>	for( int i=0; i<cutstrs.size(); i++ ) cutvec.push_back(0);
142	>	for( int i=0; i<2; i++ )  {
143	>	zcutvec.push_back(vector<unsigned>());
144	>	for( int j=0; j<cutstrs.size(); j++ ) zcutvec[i].push_back(0);
145	>	}
146	>	for( int i=0; i<3; i++ )  {
147	>	zzcutvec.push_back(vector<unsigned>());
148	>	for( int j=0; j<cutstrs.size(); j++ ) zzcutvec[i].push_back(0);
149	>	}
150	>
151
152		//
153		// args
#	Line 96 | Line 164 | int main(int argc, char** argv)
164		return 1;
165		}
166		ctrl.dump();
99	–	assert( ctrl.mH != 0 );
167
168
169
#	Line 113 | Line 180 | int main(int argc, char** argv)
180		if( !(strncmp( fname.c_str(), "#", 1 ) ) ) continue; // skip commented lines
181		cout << "adding inputfile : " << fname.c_str() << endl;
182		entrymap[string(fname.c_str())] = unskimmedEntries(fname.c_str());
183	+	cout << "unskimmed entries: " << entrymap[string(fname.c_str())] << endl;
184		chain->AddFile(fname.c_str());
185		hltchain->AddFile(fname.c_str());
186		}
187		} else {
188		cout << "adding inputfile : " << ctrl.inputfile.c_str() << endl;
189	+	unsigned tmpent = unskimmedEntries(ctrl.inputfile.c_str());
190	+	cout << "tmpent: " << tmpent << endl;
191		entrymap[string(ctrl.inputfile.c_str())] = unskimmedEntries(ctrl.inputfile.c_str());
192	+	cout << "unskimmed entries: " << entrymap[string(ctrl.inputfile.c_str())] << endl;
193		chain->AddFile(ctrl.inputfile.c_str());
194		hltchain->AddFile(ctrl.inputfile.c_str());
195		}
#	Line 149 | Line 220 | int main(int argc, char** argv)
220		nt.makeKinematicsBranch(kinematics);
221		nt.makeInfoBranch(evtinfo);
222
223	+	FOTuple foTree( nt.getFile(), (const char*)"FOtree");
224	+	foTree.makeFailedLeptonBranch(failingLeptons);
225	+	foTree.makeZLeptonBranch(passingLeptons);
226	+
227		TH1F * h_w_hpt;
228		if(ctrl.mc) {
229		nt.makeWeightBranch(weights);
230		//    nt.makeGenInfoBranch(geninfo);
231		initEfficiencyWeights();
232	+	initPUWeights();
233		/*
234		// Higgs only, pt reweighting
235		if( ctrl.mH <= 140 ) {
#	Line 170 | Line 246 | int main(int argc, char** argv)
246		initRunLumiRangeMap();
247		}
248
249	<	initMuonIDMVA();
249	>	//  initMuonIDMVA();
250		initMuonIsoMVA();
251		initElectronIDMVA();
252		initElectronIsoMVA();
253		initTrigger();
254	<
255	<
256	<	//##########################################################################
257	<	// Setup tree I/O
182	<	//##########################################################################
254	>
255	>	// hack
256	>	initEvtRhoMap(evtrhoMap);
257	>
258
259		//
260	<	// Access samples and fill histograms
260	>	// Setup tree I/O
261	>	//--------------------------------------------------------------------------------------------------------------
262		TFile *inputFile=0;
263		TTree *eventTree=0;
264		string currentFile("");
265
190	–	// Data structures to store info from TTrees
266		mithep::EventHeader *info    = new mithep::EventHeader();
267		//  mithep::TGenInfo    *ginfo  = new mithep::TGenInfo();
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>();
271		mithep::Array<mithep::PFCandidate>          *pfArr         = new mithep::Array<mithep::PFCandidate>();
272	<	mithep::Array<mithep::PileupEnergyDensity>  *puArr         = new mithep::Array<mithep::PileupEnergyDensity>();
272	>	mithep::Array<mithep::PileupInfo>           *puArr         = new mithep::Array<mithep::PileupInfo>();
273	>	mithep::Array<mithep::PileupEnergyDensity>  *puDArr        = new mithep::Array<mithep::PileupEnergyDensity>();
274		mithep::Array<mithep::Track>                *trkArr        = new mithep::Array<mithep::Track>();
275		mithep::TriggerMask                         *trigMask      = new mithep::TriggerMask();
276		mithep::TriggerTable                        *hltTable      = new mithep::TriggerTable();
#	Line 205 | Line 281 | int main(int argc, char** argv)
281		TString fInfoName(Names::gkEvtHeaderBrn);
282		TString fPrimVtxName(Names::gkPVBrn);
283		TString fPFCandidateName(Names::gkPFCandidatesBrn);
284	+	TString fPileupInfoName(Names::gkPileupInfoBrn);
285		TString fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn);
286		TString fTrackName(Names::gkTrackBrn);
287		TString fTriggerMaskName(Names::gkHltBitBrn);
#	Line 215 | Line 292 | int main(int argc, char** argv)
292		chain->SetBranchAddress(fMuonName,        &muonArr);
293		chain->SetBranchAddress(fPrimVtxName,     &vtxArr);
294		chain->SetBranchAddress(fPFCandidateName, &pfArr);
295	<	chain->SetBranchAddress(fPileupEnergyDensityName, &puArr);
295	>	if( ctrl.mc ) chain->SetBranchAddress(fPileupInfoName,  &puArr);
296	>	chain->SetBranchAddress(fPileupEnergyDensityName, &puDArr);
297		chain->SetBranchAddress(fTrackName, &trkArr);
298		chain->SetBranchAddress(fTriggerMaskName, &trigMask);
299		cout << "hlttable: " << fTriggerTableName << endl;
#	Line 240 | Line 318 | int main(int argc, char** argv)
318		cout << "nEntries: " << imax << endl;
319
320
321	<	//##########################################################################
322	<	// Loop !!!!!!!!! should alternate events here and +1 in the training ...
323	<	//##########################################################################
321	>	//
322	>	// Loop !!!!!!!!!
323	>	//--------------------------------------------------------------------------------------------------------------
324		for(UInt_t ientry=0; ientry<imax; ientry++)
325		{
326		chain->GetEntry(ientry);
249	–	if(!(ientry%1000)) cerr << "entry: " << ientry << endl;
327
328	+	// events that fail bc ele overlaps w/ reco mu
329	+	/*
330	+	if( info->EvtNum() != 67315 &&
331	+	info->EvtNum() != 288693 ) continue;
332	+	*/
333	+
334	+	// events where they apparently don't veto electrons ...
335	+	/*
336	+	if( info->EvtNum() != 141713 &&
337	+	info->EvtNum() != 178019 &&
338	+	info->EvtNum() != 60823 &&
339	+	info->EvtNum() != 192795 ) continue;
340	+	*/
341	+
342	+	if(!(ientry%1000)) cerr << "entry: " << ientry << endl;
343		setPV( ctrl, vtxArr, vtx);
344	+	//if (!(setPV( ctrl, vtxArr, vtx)) ) continue;
345	+
346	+	float rho = evtrhoMap[info->EvtNum()];
347	+
348	+	//
349	+	// pfNoPU
350	+	//
351	+	//mithep::Array<PFCandidate> pfNoPileUp;
352	+	PFnoPUflag.clear();
353	+	int pfnopu_size = makePFnoPUArray( pfArr, PFnoPUflag, vtxArr );
354	+	assert(pfnopu_size == pfArr->GetEntries());
355
356	+	//
357	+	// data/MC
358	+	//
359		if(ctrl.mc) {
360	<	weights.w = getWeight(xstab,entrymap,chain);
360	>	//
361	>	// xsec & PU weights
362	>	//
363	>	weights.w = getWeight(xstab,entrymap,chain)/ctrl.totalMC;
364	>	int npu = -1;
365	>	for(int i=0;i<puArr->GetEntries();i++) {
366	>	if(puArr->At(i)->GetBunchCrossing() == 0)
367	>	npu = puArr->At(i)->GetPU_NumInteractions();
368	>	}
369	>	assert(npu>=0);
370	>	evtinfo.npu;
371	>	weights.npuw = getPUWeight(evtinfo.npu);
372	>	//      cout << "weight: " << weights.w << endl;
373	>	//
374	>	// trigger
375	>	//
376	>	if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) {
377	>	currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile());
378	>	hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
379	>	hltchain->GetEntry(0);
380	>	hltTable->Clear();
381	>	fillTriggerBits(hltTable, hltTableStrings );
382	>	}
383	>	if( ctrl.debug ) cout << "file is : " << currentFile  << endl;
384	>	fillTriggerBits( hltTable, trigMask, triggerBits );
385	>	passes_HLT_MC = true; // passed to apply as extern global, just for sync
386	>	if( !passHLTMC(triggerBits, info->RunNum(), info->EvtNum() ) ) {
387	>	passes_HLT_MC = false;
388	>	}
389		} else {
390	+	//
391	+	// JSON
392	+	//
393	+	/*
394	+	RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec());
395	+	if( !(rlrm.HasRunLumi(rl)) )  {
396	+	if( ctrl.debug ) cout << "\tfails JSON" << endl;
397	+	continue;
398	+	}
399	+	*/
400	+	//
401	+	// trigger
402	+	//
403		if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) {
404		currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile());
405		hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings);
#	Line 261 | Line 408 | int main(int argc, char** argv)
408		fillTriggerBits(hltTable, hltTableStrings );
409		}
410		if( ctrl.debug ) cout << "file is : " << currentFile  << endl;
411	+	/*
412		fillTriggerBits( hltTable, trigMask, triggerBits );
265	–	}
266	–
267	–
268	–	//
269	–	// trigger
270	–	//
271	–	if( !ctrl.mc ) {
413		if( !passHLT(triggerBits, info->RunNum(), info->EvtNum() ) ) {
414		if( ctrl.debug ) cout << "\tfails trigger ... " << endl;
415		continue;
416		}
417	+	*/
418		}
419	<
419	>
420		//
421		// lepton/kinematic selection ...
422		//
423	+	failingLeptons.clear();
424	+	passingLeptons.clear();
425		EventData ret4l =
426	<	apply_HZZ4L_selection(ctrl, info,
427	<	vtx,
428	<	pfArr,
429	<	puArr,
430	<	electronArr,
431	<	&electronPreSelection,
432	<	//                              &electronBDTSelection,
433	<	&electronIDMVASelection,
434	<	//                              &electronIsoSelection,
435	<	&electronIsoMVASelection,
436	<	muonArr,
437	<	&muonPreSelection,
438	<	&muonIDPFSelection,
439	<	//                              &muonIDMVASelection,
440	<	//                              &passMuonSelection,
441	<	&muonIsoMVASelection);
426	>
427	>	/*
428	>	apply_HZZ4L_reference2_selection(ctrl, info,
429	>	vtx,
430	>	pfArr,
431	>	//rho,
432	>	puDArr,
433	>	electronArr,
434	>	&electronReferencePreSelection,
435	>	&electronReferenceIDMVASelection,
436	>	&electronReferenceIsoSelection,
437	>	muonArr,
438	>	&muonReferencePreSelection,
439	>	&muonIDPFSelection,
440	>	&muonReferenceIsoSelection);
441	>	*/
442	>	apply_HZZ4L_reference2_selection(ctrl, info,
443	>	vtx,
444	>	pfArr,
445	>	//rho,
446	>	puDArr,
447	>	electronArr,
448	>	&electronReferencePreSelection,
449	>	&electronReferenceIDMVASelection,
450	>	&electronIsoMVASelection,
451	>	muonArr,
452	>	&muonReferencePreSelection,
453	>	&muonIDPFSelection,
454	>	&muonIsoMVASelection);
455	>
456	>
457	>
458	>
459		if( ctrl.debug ) cout << endl;
460	+
461	+
462	+	int Z1type=0, Z2type=0;
463	+	if( ret4l.status.selectionBits.to_ulong() >= PASS_ZCANDIDATE ) {
464	+	if( abs(ret4l.Z1leptons[0].type) == 11  ) Z1type = 11;
465	+	if( abs(ret4l.Z1leptons[0].type) == 13  ) Z1type = 13;
466	+	}
467	+	if( ret4l.status.selectionBits.to_ulong() >= PASS_ZZCANDIDATE ) {
468	+	if( abs(ret4l.Z2leptons[0].type) == 11  ) Z2type = 11;
469	+	if( abs(ret4l.Z2leptons[0].type) == 13  ) Z2type = 13;
470	+	}
471	+
472	+	#ifdef SYNC
473	+	cout  << "bits: " << ret4l.status.selectionBits  << "\t"
474	+	<< "Z1: " << Z1type << "\t"
475	+	<< "Z2: " << Z2type << endl;
476	+	cout << endl;
477	+	#endif
478
479		if( ret4l.status.pass() ) {
301	–
480		fillAngles(ret4l,angles);
481		fillKinematics(ret4l,kinematics);
482		fillEventInfo(info,evtinfo);
483	<	//if(ctrl.mc) fillGenInfo(ginfo,geninfo);
483	>	if( ctrl.mc) {
484	>	// fillGenInfo(ginfo,geninfo);
485	>	setEffiencyWeights(ret4l, weights);
486	>	if(ctrl.debug)
487	>	cout << "w: " << weights.w << "\t"
488	>	<< "won: " << weights.won << "\t"
489	>	<< "woff: " << weights.woff << endl;
490	>	}
491
492		/*
493		// only for Higgs < 140
#	Line 323 | Line 508 | int main(int argc, char** argv)
508		pass++;
509		//      if( pass > 3 ) break;
510
511	<	}
511	>	} else if( ret4l.status.selectionBits.test(PASS_GOODZ1) ) { // save for fakes ...
512	>	//      if(ctrl.debug) {
513	>	cout << "failingLeptons : " << failingLeptons.size() << endl;
514	>	for( int f=0; f<failingLeptons.size(); f++ ) {
515	>	SimpleLepton  l = failingLeptons[f];
516	>	cout << "f: " << f << "\t"
517	>	<< "type: " << l.type << "\t"
518	>	<< "pT: " << l.vec.Pt()  << "\t"
519	>	<< "eta: " << l.vec.Eta()
520	>	<< endl;
521	>	//    }
522	>	}
523	>	foTree.Fill();
524	>	}
525		}
526
527	+
528	+	foTree.getFile()->cd();
529	+	foTree.getTree()->Write();
530		nt.WriteClose();
531	+
532	+	cout << endl;
533	+	cout << endl;
534	+	for( int i=0; i<cutvec.size(); i++ ) {
535	+	cout << "cut: " << i << "\t"
536	+	<< "pass: " << cutvec[i];
537	+
538	+	if( i>PASS_TRIGGER&&i<=PASS_GOODZ1 )
539	+	cout << "\t2e: " << zcutvec[0][i]
540	+	<< "\t2m: " << zcutvec[1][i];
541	+
542	+	if( i>PASS_ZCANDIDATE )
543	+	cout << "\t4e: " << zzcutvec[0][i]
544	+	<< "\t4m: " << zzcutvec[1][i]
545	+	<< "\t2e2m: " << zzcutvec[2][i];
546	+
547	+	cout << "\t" << cutstrs[i] << endl;
548	+
549	+	cout << endl;
550	+	}
551	+
552		}
553
554		//----------------------------------------------------------------------------
#	Line 355 | Line 577 | bool setPV(ControlFlags ctrl,
577		if( ctrl.debug ) cout << "vertex :: " << i << "\tntrks: " << pv->NTracks() << endl;
578
579		// Select best PV for corrected d0; if no PV passing cuts, the first PV in the collection will be used
580	<	//  if(!pv->IsValid()) continue;
580	>	if(!pv->IsValid())                                continue;
581		if(pv->NTracksFit()       < fMinNTracksFit)       continue;
582		if(pv->Ndof()                 < fMinNdof)         continue;
583		if(fabs(pv->Z()) > fMaxAbsZ)                      continue;
#	Line 376 | Line 598 | bool setPV(ControlFlags ctrl,
598		if( ctrl.debug) cout << "new PV set, pt : " << best_sumpt << endl;
599		}
600		}
601	<
601	>	// sync
602		if(!bestPV) bestPV = vtxArr->At(0);
603	+	//if(!bestPV) return false;
604		vtx.SetPosition(bestPV->X(),bestPV->Y(),bestPV->Z());
605		vtx.SetErrors(bestPV->XErr(),bestPV->YErr(),bestPV->ZErr());
606		return true;
607		};
608	+
609	+
610	+
611	+	//----------------------------------------------------------------------------
612	+	void setEffiencyWeights(EventData & evtdat, WeightStruct & weights )
613	+	//----------------------------------------------------------------------------
614	+	{
615	+	vector<SimpleLepton> lepvec = evtdat.Z1leptons;
616	+	lepvec.insert(lepvec.end(), evtdat.Z2leptons.begin(), evtdat.Z2leptons.end());
617	+	double w_offline=-1, werr_offline=0;
618	+	double w_online=-1, werr_online=0;
619	+
620	+	vector< pair <double,double> > wlegs; // pair here is eff & err
621	+	vector< pair <float,float> > mvec;  // pair here is eta & pt
622	+
623	+	//      vector< pair <float,float> > evec;  // pair here is eta & pt
624	+	// now deal with medium vs loose
625	+	vector< pair< bool, pair <float,float> > > evec;  // pair here is eta & pt
626	+
627	+	for( int k=0; k<lepvec.size(); k++ ) {
628	+	if( abs(lepvec[k].type) == 13 ) {
629	+	mvec.push_back( std::pair<float,float> (fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt()) );
630	+	wlegs.push_back( muonPerLegOfflineEfficiencyWeight( fabs(lepvec[k].vec.Eta()),
631	+	lepvec[k].vec.Pt() ) );
632	+	} else {
633	+
634	+	// now deal with medium vs loose
635	+	//              evec.push_back( std::pair<float,float> (fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt()) );
636	+
637	+	std::pair<float,float> tmppair(fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt());
638	+	evec.push_back( std::pair<bool, std::pair<float,float> > (lepvec[k].isTight, tmppair) );
639	+
640	+	//              wlegs.push_back( elePerLegOfflineEfficiencyWeight(  fabs(lepvec[k].vec.Eta()),
641	+	//                                                                 lepvec[k].vec.Pt() ) );
642	+
643	+	wlegs.push_back( elePerLegOfflineEfficiencyWeight(  fabs(lepvec[k].vec.Eta()),
644	+	lepvec[k].vec.Pt(),
645	+	lepvec[k].isTight ) );
646	+
647	+	}
648	+	}
649	+
650	+	pair<double,double> offpair = getOfflineEfficiencyWeight( wlegs );
651	+	weights.woff    = offpair.first;
652	+	weights.werroff = offpair.second;
653	+
654	+	pair<double,double> onpair  = getOnlineEfficiencyWeight( mvec, evec );
655	+	weights.won    = onpair.first;
656	+	weights.werron = onpair.second;
657	+	}
658	+
659	+
660	+	//----------------------------------------------------------------------------
661	+	void initRunLumiRangeMap()
662	+	//----------------------------------------------------------------------------
663	+	{
664	+	/*
665	+	rlrm.AddJSONFile(std::string("./data/Cert_136033-149442_7TeV_Apr21ReReco_Collisions10_JSON.txt"));
666	+	//  rlrm.AddJSONFile(std::string("./data/Cert_160404-173244_7TeV_PromptReco_Collisions11_JSON_v2.txt"));
667	+	rlrm.AddJSONFile(std::string("./data/Cert_160404-178078_7TeV_PromptReco_Collisions11_JSON.txt"));
668	+	rlrm.AddJSONFile(std::string("./data/Cert_160404-163869_7TeV_May10ReReco_Collisions11_JSON_v3.txt"));
669	+	rlrm.AddJSONFile(std::string("./data/Cert_170249-172619_7TeV_ReReco5Aug_Collisions11_JSON.txt"));
670	+	// r11b
671	+	rlrm.AddJSONFile(std::string("./data/Cert_160404-180252_7TeV_PromptReco_Collisions11_JSON.txt"));
672	+	*/
673	+	};
674	+
675	+
676	+	//----------------------------------------------------------------------------
677	+	void initPUWeights()
678	+	//----------------------------------------------------------------------------
679	+	{
680	+	TFile * puf = new TFile("data/PileupReweighting.Summer11DYmm_To_Full2011.root");
681	+	hpu = (TH1D*)(puf->Get("puWeights"));
682	+	}
683	+
684	+	//----------------------------------------------------------------------------
685	+	double getPUWeight(unsigned npu)
686	+	//----------------------------------------------------------------------------
687	+	{
688	+	return hpu->GetBinContent(hpu->FindBin(npu));
689	+	}
690	+
691	+	//----------------------------------------------------------------------------
692	+	void initEvtRhoMap( map<unsigned, float> & evtrhoMap )
693	+	//----------------------------------------------------------------------------
694	+	{
695	+	unsigned evt;
696	+	float rho;
697	+
698	+	cout << "initialzing EvtRhoMap ";
699	+	//FILE * f = fopen("evtrho.dat", "r");
700	+	//  FILE * f = fopen("allrho.cali.uniq.dat", "r");
701	+	FILE * f = fopen("allrhoAA.cali.uniq.dat", "r");
702	+	int lcount=0;
703	+	while( fscanf( f, "%u:%f", &evt, &rho ) ) {
704	+	if( feof(f) ) break;
705	+	evtrhoMap[evt] = rho;
706	+	lcount++;
707	+	if( !(lcount%1000) ) cout << "."; flush(cout);
708	+	}
709	+	cout << " done" << endl;
710	+	}
711	+
712	+	//----------------------------------------------------------------------------
713	+	unsigned makePFnoPUArray(mithep::Array<PFCandidate> * fPFCandidates,
714	+	vector<bool> &pfNoPileUpflag,
715	+	mithep::Array<Vertex>      * vtxArr )
716	+	//----------------------------------------------------------------------------
717	+	{
718	+	for(UInt_t i = 0; i < fPFCandidates->GetEntries(); i++) {
719	+	const PFCandidate *pf = fPFCandidates->At(i);
720	+	assert(pf);
721	+
722	+	if(pf->PFType() == PFCandidate::eHadron) {
723	+	if(pf->HasTrackerTrk() &&
724	+	vtxArr->At(0)->HasTrack(pf->TrackerTrk()) &&
725	+	vtxArr->At(0)->TrackWeight(pf->TrackerTrk()) > 0) {
726	+
727	+	pfNoPileUpflag.push_back(1);
728	+
729	+	} else {
730	+
731	+	Bool_t vertexFound = kFALSE;
732	+	const Vertex *closestVtx = 0;
733	+	Double_t dzmin = 10000;
734	+
735	+	// loop over vertices
736	+	for(UInt_t j = 0; j < vtxArr->GetEntries(); j++) {
737	+	const Vertex *vtx = vtxArr->At(j);
738	+	assert(vtx);
739	+
740	+	if(pf->HasTrackerTrk() &&
741	+	vtx->HasTrack(pf->TrackerTrk()) &&
742	+	vtx->TrackWeight(pf->TrackerTrk()) > 0) {
743	+	vertexFound = kTRUE;
744	+	closestVtx = vtx;
745	+	break;
746	+	}
747	+	Double_t dz = fabs(pf->SourceVertex().Z() - vtx->Z());
748	+	if(dz < dzmin) {
749	+	closestVtx = vtx;
750	+	dzmin = dz;
751	+	}
752	+	}
753	+
754	+	//      if(fCheckClosestZVertex) {
755	+	if(1) {
756	+	// Fallback: if track is not associated with any vertex,
757	+	// associate it with the vertex closest in z
758	+	if(vertexFound || closestVtx != vtxArr->At(0)) {
759	+	//      pfPileUp->Add(pf);
760	+	pfNoPileUpflag.push_back(0);
761	+	} else {
762	+	pfNoPileUpflag.push_back(1);
763	+	}
764	+	} else {
765	+	if(vertexFound && closestVtx != vtxArr->At(0)) {
766	+	//      pfPileUp->Add(pf);
767	+	pfNoPileUpflag.push_back(0);
768	+	} else {
769	+	//      PFCandidate * pfNoPileUp->AddNew(); // Ridiculous but that's how it is
770	+	pfNoPileUpflag.push_back(1);
771	+	}
772	+	}
773	+	} //hadron & trk stuff
774	+	} else { // hadron
775	+	//      PFCandidate * ptr = pfNoPileUp->AddNew();
776	+	pfNoPileUpflag.push_back(1);
777	+	}
778	+	} // Loop over PF candidates
779	+
780	+	return pfNoPileUpflag.size();
781	+	}
782	+

Diff Legend

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