ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

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

Comparing UserCode/MitHzz4l/Selection/src/applySelectionFillAnglesRunBDT.cc (file contents):
Revision 1.1 by khahn, Mon Nov 7 05:39:59 2011 UTC vs.
Revision 1.7 by khahn, Mon Apr 30 21:42:16 2012 UTC

#	Line 2 | Line 2
2		// System headers
3		//
4		#include <vector>                   // STL vector class
5	–	//#include <vector>                   // STL pair
5		#include <iostream>                 // standard I/O
6		#include <iomanip>                  // functions to format standard I/O
8	–	#include <fstream>                  // functions for file I/O
9	–	#include <string>                   // C++ string class
10	–	#include <sstream>                  // class for parsing strings
11	–	#include <assert.h>
12	–	#include <stdlib.h>
13	–	#include <getopt.h>
7		using namespace std;
8
9		//
10	<	// ROOT headers
10	>	// root headers
11		//
12	<	#include <TROOT.h>                  // access to gROOT, entry point to ROOT system
13	<	#include <TSystem.h>                // interface to OS
14	<	#include <TFile.h>                  // file handle class
15	<	#include <TNtuple.h>
16	<	#include <TTree.h>                  // class to access ntuples
17	<	#include <TChain.h>                 //
25	<	#include <TBranch.h>                // class to access branches in TTree
26	<	#include <TClonesArray.h>           // ROOT array class
27	<	#include <TCanvas.h>                // class for drawing
28	<	#include <TH1F.h>                   // 1D histograms
29	<	#include <TBenchmark.h>             // class to track macro running statistics
30	<	#include <TLorentzVector.h>         // 4-vector class
31	<	#include <TVector3.h>               // 3D vector class
12	>	#include <TFile.h>
13	>	#include <TTree.h>
14	>	#include <TChain.h>
15	>	#include <TBranch.h>
16	>	#include <TClonesArray.h>
17	>
18
19		//
20		// TMVA
#	Line 41 | Line 27 | using namespace std;
27		//
28		// ntuple format headers
29		//
30	<	#include "HiggsAnaDefs.hh"
31	<	#include "TEventInfo.hh"
32	<	#include "TElectron.hh"
33	<	#include "TMuon.hh"
34	<	#include "TJet.hh"
35	<	#include "SiMVAElectronSelection.h"
36	<	#include "RunLumiRangeMap.h"
37	<	#include "EfficiencyWeightsInterface.h"
30	>	#include "EventHeader.h"
31	>	#include "Electron.h"
32	>	#include "Muon.h"
33	>	#include "Vertex.h"
34	>	#include "PFCandidate.h"
35	>	#include "PFCandidateCol.h"
36	>	#include "Names.h"
37	>
38
39		//
40	<	// utility headers
40	>	// our headers
41		//
42		#include "ParseArgs.h"
43	<	#include "SampleWeight.h"
43	>	#include "MuonSelection.h"
44	>	#include "ElectronSelection.h"
45	>	#include "IsolationSelection.h"
46		#include "Selection.h"
47	<	#include "HZZCiCElectronSelection.h"
60	<	#include "HZZLikelihoodElectronSelection.h"
61	<	#include "HZZBDTElectronSelection.h"
62	<	#include "SampleWeight.h"
47	>
48		#include "Angles.h"
49	+	#include "KinematicsStruct.h"
50	+	#include "InfoStruct.h"
51	+	//#include "GenInfoStruct.h"
52	+	#include "WeightStruct.h"
53	+	#include "GlobalDataAndFuncs.h"
54	+	#include "AngleTuple.h"
55
56	<	//#define BDTFAIL_THEIR_EVENTS
57	<	//#define THEIR_EVENTS
56	>	#include "SampleWeight.h"
57	>	#include "EfficiencyWeightsInterface.h"
58
59	<	Double_t etaFromY(Double_t pt, Double_t y, Double_t phi, Double_t m)
69	<	{
70	<	Double_t a  = (1+exp(2*y))/(exp(2*y)-1); // intermediate term
71	<	if(a*a<1) { cout << "a too small" << endl; assert(0); }
72	<	Double_t E  = sqrt( a*a*(pt*pt+m*m)/(a*a-1) );
73	<	Double_t pz = E*E - pt*pt - m*m;
74	<	if(pz<0) { cout << "imag. pz" << endl; assert(0); }
75	<	pz = sqrt(pz);
76	<	if(y<0) pz *= -1;
77	<	TLorentzVector v;
78	<	v.SetPxPyPzE(pt*cos(phi),pt*sin(phi),pz,E);
79	<	Double_t th = v.Theta();
80	<	return -log(tan(th/2));
81	<	};
82	<
83	<	int getChannel(mithep::TGenInfo * ginfo) {
84	<
85	<	int gchannel=-1;
86	<
87	<	if( abs(ginfo->id_1_a) == EGenType::kElectron && abs(ginfo->id_1_b) == EGenType::kElectron )  gchannel=0;
88	<	else if( abs(ginfo->id_1_a) == EGenType::kMuon && abs(ginfo->id_1_b) == EGenType::kMuon ) gchannel=1;
89	<	else if( (abs(ginfo->id_1_a) == EGenType::kElectron && abs(ginfo->id_1_b) == EGenType::kMuon) ||
90	<	(abs(ginfo->id_1_a) == EGenType::kMuon && abs(ginfo->id_1_b) == EGenType::kElectron) ) gchannel=2;
91	<
92	<	return gchannel;
93	<	}
59	>	using namespace mithep;
60
95	–
61		//=== MAIN =================================================================================================
62		int main(int argc, char** argv)
63		{
#	Line 134 | Line 99 | int main(int argc, char** argv)
99		}
100
101		// table of cross section values
102	<	SimpleTable xstab("./data/xs.dat");
102	>	//  SimpleTable xstab("./data/xs.dat");
103	>
104	>
105	>
106
107		//
108		// Setup
#	Line 145 | Line 113 | int main(int argc, char** argv)
113		} else {
114		ofname = "tmp.root";
115		}
116	<	TFile *file = new TFile(ofname, "RECREATE");
149	<	TH1F * h_evt = new TH1F("hevt", "hevt", 2, 0, 2 );
150	<	TH1F * h_evtfail = new TH1F("hevtfail", "hevtfail", 1024, 0, 1024 );
151	<
152	<	unsigned run, evt, lumi, channel, gchannel;
153	<	float mZ1, mZ2, m4l, pt4l, fchannel;
154	<	unsigned tZ1, tZ2;
116	>
117		Angles angles;
118	<	double bdt;
119	<	double w, woff, werroff, won, werron;
120	<	unsigned npu; double npuw;
121	<
122	<	TTree * angletuple = new TTree("angletuple", "angletuple" );
123	<	angletuple->Branch("run", &run);
124	<	angletuple->Branch("evt", &evt);
125	<	angletuple->Branch("lumi", &lumi );
126	<	angletuple->Branch("channel", &channel );
127	<	angletuple->Branch("costheta1", &(angles.costheta1) );
128	<	angletuple->Branch("costheta2", &(angles.costheta2) );
129	<	angletuple->Branch("costhetastar", &(angles.costhetastar) );
130	<	angletuple->Branch("Phi", &(angles.Phi) );
131	<	angletuple->Branch("Phi1", &(angles.Phi1) );
132	<	angletuple->Branch("mZ2", &(angles.mz2) );
133	<	angletuple->Branch("mZ1", &(angles.mz1) );
134	<	angletuple->Branch("m4l", &(angles.m4l) );
135	<	angletuple->Branch("bdt", &bdt );
136	<	angletuple->Branch("w",     &w );
137	<	angletuple->Branch("won",   &won );
138	<	angletuple->Branch("woff",  &woff );
139	<	angletuple->Branch("npuw",  &npuw );
140	<
141	<	TTree * passtuple = new TTree("passtuple", "passtuple" );
142	<	if( ctrl.mc )   {
143	<	passtuple->Branch("gchannel", &gchannel);
144	<	passtuple->Branch("woff", &woff);
145	<	passtuple->Branch("werroff", &werroff);
146	<	passtuple->Branch("won", &won);
185	<	passtuple->Branch("werron", &werron);
186	<	passtuple->Branch("npu",  &npu);
187	<	passtuple->Branch("npuw", &npuw);
118	>	KinematicsStruct kinematics;
119	>	InfoStruct evtinfo;
120	>	WeightStruct weights;
121	>	//  GenInfoStruct geninfo;
122	>
123	>	AngleTuple nt( (const char*)ofname, (const char*)"zznt");
124	>	nt.makeAngleBranch(angles);
125	>	nt.makeKinematicsBranch(kinematics);
126	>	nt.makeInfoBranch(evtinfo);
127	>
128	>	TH1F * h_w_hpt;
129	>	if(ctrl.mc) {
130	>	nt.makeWeightBranch(weights);
131	>	//    nt.makeGenInfoBranch(geninfo);
132	>	initEfficiencyWeights();
133	>	/*
134	>	// Higgs only, pt reweighting
135	>	if( ctrl.mH <= 140 ) {
136	>	char wptname[256];
137	>	sprintf( wptname, "/data/blue/pharris/Flat/ntupler/root/weight_ptH_%i.root", ctrl.mH );
138	>	TFile * f_hpt = new TFile(wptname);
139	>	f_hpt->Print();
140	>	sprintf(wptname, "weight_hqt_fehipro_fit_%i", ctrl.mH);
141	>	h_w_hpt  = (TH1F*)(f_hpt->FindObjectAny(wptname));
142	>	h_w_hpt->Print();
143	>	}
144	>	*/
145	>	} else {
146	>	initRunLumiRangeMap();
147		}
189	–	passtuple->Branch("channel", &channel);
190	–	passtuple->Branch("run", &run);
191	–	passtuple->Branch("evt", &evt);
192	–	passtuple->Branch("lumi", &lumi);
193	–	passtuple->Branch("mZ1", &mZ1);
194	–	passtuple->Branch("mZ2", &mZ2);
195	–	passtuple->Branch("tZ1", &tZ1);
196	–	passtuple->Branch("tZ2", &tZ2);
197	–	passtuple->Branch("m4l", &m4l);
198	–	passtuple->Branch("pt4l", &pt4l);
199	–	passtuple->Branch("w", &w);
200	–
201	–
202	–	initCiCSelection();
203	–	if(ctrl.eleSele=="lik") initLikSelection();
204	–	if(ctrl.eleSele=="bdt") initBDTSelection();
205	–	// if(ctrl.eleSele=="si" ) initSiMVAElectronSelection();
206	–	initRunLumiRangeMap();
207	–	initEfficiencyWeights();
148
149	+	//  initBDTSelection();
150	+	initMuonIDMVA();
151	+	initMuonIsoMVA();
152	+	initElectronIDMVA();
153	+	initElectronIsoMVA();
154
155		//
156		// tmva
157		//
158	+	int channel;
159		TMVA::Reader * reader = new TMVA::Reader( "V" );
160	<	float f_costheta1, f_costheta2, f_costhetastar, f_Phi, f_Phi1;
161	<	float f_mz2, f_m4l;
162	<	int i_channel;
163	<
164	<	reader->AddVariable( "costheta1",     &f_costheta1);
165	<	reader->AddVariable( "costheta2",     &f_costheta2);
166	<	reader->AddVariable( "costhetastar",  &f_costhetastar);
167	<	reader->AddVariable( "Phi",           &f_Phi);
168	<	reader->AddVariable( "Phi1",          &f_Phi1);
169	<	reader->AddSpectator( "m4l",          &f_m4l);
224	<	reader->AddVariable( "mZ2",           &f_mz2);
225	<	reader->AddSpectator( "channel",      &i_channel);
226	<
227	<	char wfilename[512]; sprintf(wfilename, "./weights/hzz4lTMVA_mH%d_BDTG.weights.xml", ctrl.mH);
228	<	reader->BookMVA("BDTG", wfilename);
229	<	//   char wfilename[512]; sprintf(wfilename, "./weights/hzz4lTMVA_mH%d_Likelihood.weights.xml", ctrl.mH);
230	<	//   reader->BookMVA("Likelihood", wfilename);
160	>	reader->AddVariable( "costheta1",     &angles.costheta1);
161	>	reader->AddVariable( "costheta2",     &angles.costheta2);
162	>	reader->AddVariable( "costhetastar",  &angles.costhetastar);
163	>	reader->AddVariable( "Phi",           &angles.Phi);
164	>	reader->AddVariable( "Phi1",          &angles.Phi1);
165	>	reader->AddSpectator( "m4l",          &kinematics.m4l);
166	>	reader->AddVariable( "mZ2",           &kinematics.mZ2);
167	>	reader->AddSpectator( "channel",          &channel);
168	>	//   reader->AddVariable( "ZZpt",           &kinematics.ZZpt);
169	>	//   reader->AddVariable( "ZZeta",           &kinematics.ZZeta);
170
171	<	//
172	<	// Loop !!!!!!!!! should alternate events here and +1 in the training ...
173	<	//--------------------------------------------------------------------------------------------------------------
174	<	if(ctrl.fakeScheme!="none") {cout << "wrong executable!" << endl; assert(0);}
171	>
172	>	char wfilename[512]; sprintf(wfilename, "./weights/hzz4lTMVA_mH%d_BDTG.weights.xml", ctrl.mH);
173	>	reader->BookMVA("BDTG", wfilename);
174	>
175	>
176	>	//##########################################################################
177	>	// Setup tree I/O
178	>	//##########################################################################
179
180		//
181		// Access samples and fill histograms
182		TFile *inputFile=0;
183		TTree *eventTree=0;
184	<
184	>
185		// Data structures to store info from TTrees
186	<	mithep::TEventInfo *info    = new mithep::TEventInfo();
187	<	mithep::TGenInfo    *ginfo  = new mithep::TGenInfo();
188	<	TClonesArray *electronArr   = new TClonesArray("mithep::TElectron");
189	<	TClonesArray *muonArr       = new TClonesArray("mithep::TMuon");
186	>	mithep::EventHeader *info    = new mithep::EventHeader();
187	>	//  mithep::TGenInfo    *ginfo  = new mithep::TGenInfo();
188	>	mithep::Array<mithep::Electron>             *electronArr   = new mithep::Array<mithep::Electron>();
189	>	mithep::Array<mithep::Muon>                 *muonArr       = new mithep::Array<mithep::Muon>();
190	>	mithep::Array<mithep::Vertex>               *vtxArr        = new mithep::Array<mithep::Vertex>();
191	>	mithep::Array<mithep::PFCandidate>          *pfArr         = new mithep::Array<mithep::PFCandidate>();
192	>	mithep::Array<mithep::PileupEnergyDensity>  *puArr         = new mithep::Array<mithep::PileupEnergyDensity>();
193	>	mithep::Array<mithep::Track>                *trkArr        = new mithep::Array<mithep::Track>();
194	>
195	>	TString fElectronName(Names::gkElectronBrn);
196	>	TString fMuonName(Names::gkMuonBrn);
197	>	TString fInfoName(Names::gkEvtHeaderBrn);
198	>	TString fPrimVtxName(Names::gkPVBrn);
199	>	TString fPFCandidateName(Names::gkPFCandidatesBrn);
200	>	TString fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn);
201	>	TString fTrackName(Names::gkTrackBrn);
202
203	+	chain->SetBranchAddress(fInfoName,        &info);
204	+	chain->SetBranchAddress(fElectronName,    &electronArr);
205	+	chain->SetBranchAddress(fMuonName,        &muonArr);
206	+	chain->SetBranchAddress(fPrimVtxName,     &vtxArr);
207	+
208	+	chain->SetBranchAddress(fPFCandidateName, &pfArr);
209	+	chain->SetBranchAddress(fPileupEnergyDensityName, &puArr);
210	+	chain->SetBranchAddress(fTrackName, &trkArr);
211	+
212	+	Vertex              vtx;          // best primary vertex in the event
213
214	<	chain->SetBranchAddress("Info",       &info);
215	<	chain->SetBranchAddress("Electron", &electronArr);
251	<	chain->SetBranchAddress("Muon", &muonArr);
252	<	if( ctrl.mc ) { chain->SetBranchAddress("Gen",        &ginfo);}
253	<	//ginfo = NULL;
214	>	//  ginfo = NULL;
215	>	// if( ctrl.mc ) { chain->SetBranchAddress("Gen",        &ginfo);}
216
217		int count=0, pass=0;
218	+	float passcorr=0., passcorr_errup=0., passcorr_errdown=0.;
219		float denom[3]={0.,0.,0.};
220		float numer[3]={0.,0.,0.};
221		float numercorr[3]={0.,0.,0.};
#	Line 260 | Line 223 | int main(int argc, char** argv)
223		UInt_t imax = chain->GetEntries();
224		cout << "nEntries: " << imax << endl;
225
263	–	for(UInt_t ientry=0; ientry<imax; ientry++) {
264	–	chain->GetEntry(ientry);
265	–	if(!(ientry%100000)) cerr << "entry: " << ientry << endl;
266	–
267	–	// get event weight for this file
268	–	double eventweight=1;
269	–	if(ctrl.mc) eventweight = getWeight(xstab,entrymap,chain);
270	–	w = eventweight;
271	–
272	–	//
273	–	// corresponds to 2124.8
274	–	//
275	–	if( !ctrl.mc ) {
276	–	if(  !( (info->runNum >= 160404 && info->runNum <= 163869) ||
277	–	(info->runNum >= 165088 && info->runNum <= 167913) ||
278	–	(info->runNum >= 170249 && info->runNum <= 172619) ||
279	–	(info->runNum >= 172620 && info->runNum <= 173692) ) ) continue;
280	–	std::pair<unsigned,unsigned> tmppair(info->runNum,info->evtNum);
281	–	if( find( runevtvec.begin(), runevtvec.end(), tmppair ) != runevtvec.end() )
282	–	continue;
283	–	}
284	–
285	–	LabVectors labvectors;
286	–	unsigned evtfail = fails_HZZ4L_selection(ctrl, info, electronArr, muonArr, eventweight, passtuple, &labvectors, NULL);
226
227	<	if( !evtfail ) {
228	<
229	<	if( !ctrl.mc ) runevtvec.push_back(std::pair<unsigned,unsigned> (info->runNum,info->evtNum) );
227	>	//##########################################################################
228	>	// Loop !!!!!!!!! should alternate events here and +1 in the training ...
229	>	//##########################################################################
230	>	for(UInt_t ientry=0; ientry<imax; ientry++)
231	>	{
232	>	chain->GetEntry(ientry);
233	>	if(!(ientry%1000)) cerr << "entry: " << ientry << endl;
234	>
235	>	//      if( !(info->EvtNum() == 324923 ) && !(info->EvtNum() == 325156)) continue;
236	>
237	>	//
238	>	// Get primary vertices
239	>	// Assumes primary vertices are ordered by sum-pT^2 (as should be in CMSSW)
240	>	// NOTE: if no PV is found from fitting tracks, the beamspot is used
241	>	//
242	>	const Vertex *bestPV = 0;
243	>	float best_sumpt=-1;
244	>	// good PV requirements
245	>	const UInt_t   fMinNTracksFit = 0;
246	>	const Double_t fMinNdof       = 4;
247	>	const Double_t fMaxAbsZ       = 24;
248	>	const Double_t fMaxRho        = 2;
249	>	//       const UInt_t   fMinNTracksFit = 0;
250	>	//       const Double_t fMinNdof       = 5;
251	>	//       const Double_t fMaxAbsZ       = 15;
252	>	//       const Double_t fMaxRho        = 2;
253
254	<	unsigned thechannel, therun, theevt, thelumi;
255	<	double thew, thewoff, thewon, thenpuw;
294	<
295	<	passtuple->ResetBranchAddresses();
296	<	passtuple->SetBranchAddress("channel", &thechannel);
297	<	passtuple->SetBranchAddress("run", &therun);
298	<	passtuple->SetBranchAddress("evt", &theevt);
299	<	passtuple->SetBranchAddress("lumi", &thelumi);
300	<	if( ctrl.mc ) {
301	<	passtuple->SetBranchAddress("w",       &thew); // must put this here ...
302	<	passtuple->SetBranchAddress("woff",    &thewoff); // must put this here ...
303	<	passtuple->SetBranchAddress("won",     &thewon); // must put this here ...
304	<	passtuple->SetBranchAddress("npuw",    &thenpuw); // must put this here ...
305	<	}
306	<	int ret = passtuple->GetEntry(passtuple->GetEntries()-1);
307	<	getAngles( labvectors, angles );
308	<
309	<	// for tmva
310	<	f_costheta1 = angles.costheta1;
311	<	f_costheta2 = angles.costheta2;
312	<	f_costhetastar = angles.costhetastar;
313	<	f_Phi = angles.Phi;
314	<	f_Phi1 = angles.Phi1;
315	<	f_m4l = angles.m4l;
316	<	f_mz2 = angles.mz2;
317	<
318	<	bdt = reader->EvaluateMVA("BDTG");
319	<
320	<	// for angletuple
321	<	channel = thechannel;
322	<	run = therun;
323	<	evt = theevt;
324	<	lumi = thelumi;
325	<	w    = thew;
326	<	won  = thewon;
327	<	woff = thewoff;
328	<	npuw = thenpuw;
329	<
330	<	angletuple->Fill();
331	<	}
332	<
254	>	for(UInt_t i=0; i<vtxArr->GetEntries(); ++i) {
255	>	const Vertex *pv = (Vertex*)(vtxArr->At(i));
256
257	<
258	<	#ifdef THEIR_EVENTS
259	<	if(evtfail!=0) cout << "Failed!" << endl;
260	<	#endif
261	<
262	<
263	<	if( ctrl.mc && ginfo ) {
264	<	// BR is implicit here ...
265	<	int gchannel = getChannel(ginfo);
266	<	if( gchannel >= 0 ) {
267	<	//        denom[gchannel] += eventweight;
268	<	//        if( !evtfail ) numer[gchannel] += eventweight;
269	<	denom[gchannel] += 1;
270	<	if( !evtfail ) numer[gchannel] += 1;
257	>	if( ctrl.debug ) cout << "vertex :: " << i << "\tntrks: " << pv->NTracks() << endl;
258	>
259	>	// Select best PV for corrected d0; if no PV passing cuts, the first PV in the collection will be used
260	>	//      if(!pv->IsValid()) continue;
261	>	if(pv->NTracksFit()       < fMinNTracksFit)   continue;
262	>	if(pv->Ndof()             < fMinNdof)         continue;
263	>	if(fabs(pv->Z()) > fMaxAbsZ)                  continue;
264	>	if(pv->Position().Rho()   > fMaxRho)          continue;
265	>
266	>	// tmp, take the first ...
267	>	bestPV = pv;
268	>	break;
269	>
270	>	float tmp_sumpt=0;
271	>	for( int t=0; t<pv->NTracks(); t++ )
272	>	tmp_sumpt += pv->Trk(t)->Pt();
273
274	<	// get eff weight for current evt
275	<	if( !evtfail ) {
276	<	numercorr[gchannel] += woff*won*npuw;
274	>	if( tmp_sumpt > best_sumpt  ) {
275	>	bestPV = pv;
276	>	best_sumpt = tmp_sumpt;
277	>	if( ctrl.debug) cout << "new PV set, pt : " << best_sumpt << endl;
278		}
279	+
280		}
354	–	}
355	–
356	–	count++;
357	–	if( !evtfail ) pass++;
281
282	<	} //end loop over data
283	<
284	<	if( ctrl.mc ) {
285	<	cout << "--------------" << endl;
363	<	cout << "denom: " << count << endl;
364	<	cout << "pass: " << pass << endl;
365	<	cout << "axe: " << (float)pass/count << endl;
366	<	cout << "--------------" << endl;
367	<	}
368	<
369	<	delete info;
370	<	delete electronArr;
371	<	delete muonArr;
282	>	if(!bestPV) continue;
283	>	vtx.SetPosition(bestPV->X(),bestPV->Y(),bestPV->Z());
284	>	vtx.SetErrors(bestPV->XErr(),bestPV->YErr(),bestPV->ZErr());
285	>
286
287
288	+	// get event weight for this process
289	+	if(ctrl.mc) {
290	+	//weights.w = getWeight(xstab,entrymap,chain);
291	+	} else {
292	+	// skip duplicates
293	+	std::pair<unsigned,unsigned> tmppair(info->RunNum(),info->EvtNum());
294	+	if( find( runevtvec.begin(), runevtvec.end(), tmppair ) != runevtvec.end() )
295	+	continue;
296	+	}
297	+
298	+
299	+	//      ctrl.debug = false;
300	+	if( info->RunNum() == 171578 && info->EvtNum() == 651944556  )  ctrl.debug = true;
301	+	if( ctrl.debug ) cout << "file is : " << chain->GetFile()->GetEndpointUrl()->GetFile()  << endl;
302	+
303	+	EventData ret4l =
304	+	apply_HZZ4L_selection(ctrl, info,
305	+	vtx,
306	+	pfArr,
307	+	puArr,
308	+	electronArr,
309	+	&electronPreSelection,
310	+	//                              &electronBDTSelection,
311	+	&electronIDMVASelection,
312	+	//                              &electronIsoSelection,
313	+	&electronIsoMVASelection,
314	+	muonArr,
315	+	&muonPreSelection,
316	+	&muonIDMVASelection,
317	+	//                              &passMuonSelection,
318	+	&muonIsoMVASelection);
319	+	if( ctrl.debug ) cout << endl;
320	+
321	+	if( ret4l.status.pass() ) {
322	+
323	+	runevtvec.push_back(pair<unsigned,unsigned> (info->RunNum(),info->EvtNum()) );
324	+
325	+	fillAngles(ret4l,angles);
326	+	fillKinematics(ret4l,kinematics);
327	+	fillEventInfo(info,evtinfo);
328	+	//if(ctrl.mc) fillGenInfo(ginfo,geninfo);
329	+
330	+	/*
331	+	// only for Higgs < 140
332	+	geninfo.weight *= h_w_hpt->GetBinContent(h_w_hpt->FindBin(geninfo.pt_zz));
333	+	angles.bdt = reader->EvaluateMVA("BDTG");
334	+	*/
335	+	nt.Fill();
336	+
337	+	cerr  << "PASS:: "
338	+	<< "\trun: " << evtinfo.run
339	+	<< "\tevt: " << evtinfo.evt
340	+	<< "\tlumi: " << evtinfo.lumi
341	+	<< "\tcostheta1: " << angles.costheta1
342	+	<< "\tcostheta2: " << angles.costheta2
343	+	<< "\tcostheta*: " << angles.costhetastar
344	+	<< "\tbdt: " << angles.bdt
345	+	<< endl;
346	+	pass++;
347	+	//      if( pass > 3 ) break;
348
349	<	//--------------------------------------------------------------------------------------------------------------
350	<	// Save Histograms;
377	<	//==============================================================================================================
378	<	// const char * ofname;
379	<	// if( strcmp( ctrl.outputfile.c_str(), "") ) {
380	<	//   ofname = ctrl.outputfile.c_str();
381	<	// } else {
382	<	//   ofname = "tmp.root";
383	<	// }
384	<
385	<	// TFile *file = new TFile(ofname, "RECREATE");
386	<	gROOT->cd(file->GetPath());
387	<	angletuple->Write();
388	<	passtuple->Write();
389	<	file->Close();
390	<	delete file;
391	<
392	<	// map<TString,TMVA::Reader*>::iterator it;
393	<	// for(it=readers.begin(); it!=readers.end(); it++) delete (*it).second;
394	<
395	<
396	<	float acc_4e        = numer[0]/denom[0];
397	<	float acc_4mu       = numer[1]/denom[1];
398	<	float acc_2e2mu     = numer[2]/denom[2];
399	<	float acc_err_4e    = (1./denom[0])*sqrt(numer[0]*(1-acc_4e));
400	<	float acc_err_4mu   = (1./denom[1])*sqrt(numer[1]*(1-acc_4mu));
401	<	float acc_err_2e2mu = (1./denom[2])*sqrt(numer[2]*(1-acc_2e2mu));
402	<
403	<	float acc_c_4e        = numercorr[0]/denom[0];
404	<	float acc_c_4mu       = numercorr[1]/denom[1];
405	<	float acc_c_2e2mu     = numercorr[2]/denom[2];
406	<	float acc_c_err_4e    = (1./denom[0])*sqrt(numer[0]*(1-acc_4e));
407	<	float acc_c_err_4mu   = (1./denom[1])*sqrt(numer[1]*(1-acc_4mu));
408	<	float acc_c_err_2e2mu = (1./denom[2])*sqrt(numer[2]*(1-acc_2e2mu));
409	<
410	<	if( ctrl.mc ) {
411	<	cerr << "---------------------" << endl;
412	<	//   cerr << "acc_4e: " << fixed << setw(8) << setprecision(7) << numer[0] << endl; // numer already weighted ...
413	<	//   cerr << "acc_4mu: " << fixed << setw(8) << setprecision(7) << numer[1] <<  endl;
414	<	//   cerr << "acc_2e2mu: " << fixed << setw(8) << setprecision(7) << numer[2] << endl;
415	<	cerr << "acc_4e: " << fixed << setw(8) << setprecision(7) << acc_4e << " +/- " << acc_err_4e << endl;
416	<	cerr << "acc_4mu: " << fixed << setw(8) << setprecision(7) << acc_4mu << " +/- " << acc_err_4mu << endl;
417	<	cerr << "acc_2e2mu: " << fixed << setw(8) << setprecision(7) << acc_2e2mu << " +/- " << acc_err_2e2mu << endl;
418	<	cerr << "-------corrected------" << endl;
419	<	cerr << "acc_c_4e: " << fixed << setw(8) << setprecision(7) << acc_c_4e << " +/- " << acc_c_err_4e << endl;
420	<	cerr << "acc_c_4mu: " << fixed << setw(8) << setprecision(7) << acc_c_4mu << " +/- " << acc_c_err_4mu << endl;
421	<	cerr << "acc_c_2e2mu: " << fixed << setw(8) << setprecision(7) << acc_c_2e2mu << " +/- " << acc_c_err_2e2mu << endl;
422	<	cerr << "---------------------" << endl;
423	<	}
349	>	}
350	>	}
351
352	<	cerr << "done!" << endl;
352	>	nt.WriteClose();
353		}
354
355

Diff Legend

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