ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/Systematics.C

Comparing UserCode/cbrown/AnalysisFramework/Plotting/Modules/Systematics.C (file contents):
Revision 1.1 by buchmann, Fri Jul 15 10:40:43 2011 UTC vs.
Revision 1.5 by buchmann, Wed Jul 20 08:51:33 2011 UTC

#	Line 15 | Line 15
15		#ifndef Verbosity
16		#define Verbosity 0
17		#endif
18	–	#ifndef HUSH
19	–	#define HUSH 1
20	–	#endif
18
19		#include <TFile.h>
20		#include <TTree.h>
#	Line 37 | Line 34 | Float_t jzbMax =  243;
34		Float_t jzbSel =  100;
35		int iplot=0;
36		int verbose=0;
37	+	string geqleq;
38	+	string mcjzbexpression;
39	+	bool automatized=false;//if we're running this fully automatized we don't want each function to flood the screen
40	+
41	+	TString geq_or_leq() {
42	+	if(geqleq=="geq") return TString(">=");
43	+	if(geqleq=="leq") return TString("<=");
44	+	return TString("GEQ_OR_LEQ_ERROR");
45	+	}
46
47		//______________________________________________________________________________
48		Double_t Interpolate(Double_t x, TH1 *histo)
#	Line 81 | Line 87 | TH1F* plotEff(TTree* events, TCut kbase,
87		nBins,jzbMin,jzbMax);
88		Float_t step = (jzbMax-jzbMin)/static_cast<Float_t>(nBins);
89
90	<	events->Draw("jzb[1]","genJZBSel>-400"&&kbase,"goff");
90	>	events->Draw(mcjzbexpression.c_str(),"genJZBSel>-400"&&kbase,"goff");
91		Float_t maxEff = events->GetSelectedRows();
92	<	if(verbose>0) std::cout << hname << " (" << informalname <<") " << maxEff <<  std::endl;
92	>	if(verbose>0) dout << hname << " (" << informalname <<") " << maxEff <<  std::endl;
93
94	<	if(verbose>0) std::cout <<  "JZB max = " << jzbMax << std::endl;
94	>	if(verbose>0) dout <<  "JZB max = " << jzbMax << std::endl;
95		// Loop over steps to get efficiency curve
96		char cut[256];
97		for ( Int_t iBin = 0; iBin<nBins; ++iBin ) {
98		sprintf(cut,"genJZBSel>%3f",jzbMin+iBin*step);
99	<	events->Draw("jzb[1]",TCut(cut)&&kbase,"goff");
99	>	events->Draw(mcjzbexpression.c_str(),TCut(cut)&&kbase,"goff");
100		Float_t eff = static_cast<Float_t>(events->GetSelectedRows())/maxEff;
101	<	//     std::cout << "COUCOU " << __LINE__ << std::endl;
101	>	//     dout << "COUCOU " << __LINE__ << std::endl;
102		hJzbEff->SetBinContent(iBin+1,eff);
103		hJzbEff->SetBinError(iBin+1,TMath::Sqrt(eff*(1-eff)/maxEff));
104		}
#	Line 109 | Line 115 | float pileup(TTree *events, string infor
115		jzbMax = myJzbMax;
116
117		// Acceptance cuts
118	<	TCut kbase("abs(genMllSel-91.2)<20&&pfJetGoodNum>2&&genZPtSel>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
118	>	TCut kbase("abs(genMllSel-91.2)<20&&genNjets>2&&genZPtSel>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
119
120		TH1F* hLM4 = plotEff(events,kbase,informalname);
121		hLM4->SetMinimum(0.);
#	Line 123 | Line 129 | float pileup(TTree *events, string infor
129		// Pimped-up function
130		TF1* funcUp = (TF1*)func->Clone();
131		funcUp->SetParameter( 0., func->GetParameter(0)/1.1); // 10% systematic error (up in sigma => 0.1 in erfc)
132	<	std::cout << "  PU: " << funcUp->Eval(jzbSel) << " " <<  func->Eval(jzbSel)
132	>	if(!automatized) dout << "  PU: " << funcUp->Eval(jzbSel) << " " <<  func->Eval(jzbSel)
133		<< "(" << (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100. << "%)" << std::endl;
134
135		return (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100.;
#	Line 132 | Line 138 | float pileup(TTree *events, string infor
138
139		//____________________________________________________________________________________
140		// Total selection efficiency (MC)
141	<	void MCefficiency(TTree *events,float &res, float &reserr) {
141	>	void MCefficiency(TTree *events,float &result, float &resulterr,string mcjzb) {
142
143		char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
144		// All acceptance cuts at gen. level
145	<	TCut kbase("abs(genMllSel-91.2)<20&&pfJetGoodNum>2&&genZPt>0&&genJZB>"+TString(jzbSelStr)+"&&genId1==-genId2");
145	>	TCut kbase("abs(genMllSel-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
146		// Corresponding reco. cuts
147	<	TCut ksel("abs(mll-91.2)<20&&id1==id2&&jzb[1]+0.9>"+TString(jzbSelStr));
147	>	TCut ksel("abs(mll-91.2)<20&&id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
148
149	<	events->Draw("jzb[1]",kbase&&ksel,"goff");
149	>	events->Draw(mcjzbexpression.c_str(),kbase&&ksel,"goff");
150		Float_t sel = events->GetSelectedRows();
151	<	events->Draw("jzb[1]",kbase,"goff");
151	>	events->Draw(mcjzbexpression.c_str(),kbase,"goff");
152		Float_t tot = events->GetSelectedRows();
153
154	<	res=sel/tot;
155	<	reserr=TMath::Sqrt(sel/tot*(1-sel/tot)/tot);
156	<	std::cout << "  MC efficiency: " << res << "+-" << reserr << std::endl;
154	>	result=sel/tot;
155	>	resulterr=TMath::Sqrt(sel/tot*(1-sel/tot)/tot);
156	>	if(!automatized) dout << "  MC efficiency: " << result << "+-" << resulterr << std::endl;
157		}
158
159		float JZBefficiency(TTree *events, string informalname) {
160	<	TCut kbase("abs(genMllSel-91.2)<20&&pfJetGoodNum>2&&genZPt>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
160	>	TCut kbase("abs(genMllSel-91.2)<20&&genNjets>2&&genZPt>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
161		TH1F* hLM4 = plotEff(events,kbase,informalname);
162		Int_t bin = hLM4->FindBin(jzbSel); // To get the error
163	<	std::cout << "  Efficiency at JZB==" << jzbSel  << std::endl;
164	<	std::cout << "    " << Interpolate(jzbSel,hLM4) << "+-" << hLM4->GetBinError(bin)  << std::endl;
163	>	if(!automatized) dout << "  Efficiency at JZB==" << jzbSel  << std::endl;
164	>	if(!automatized) dout << "    " << Interpolate(jzbSel,hLM4) << "+-" << hLM4->GetBinError(bin)  << std::endl;
165		return -1;
166		}
167
#	Line 187 | Line 193 | void JZBjetScale(TTree *events, float &j
193		Float_t eff  = Interpolate(jzbSel,hist);
194		Float_t effp = Interpolate(jzbSel,histp);
195		Float_t effm = Interpolate(jzbSel,histm);
196	<	std::cout << "  Efficiency at JZB==" << jzbSel  << std::endl;
197	<	std::cout << "    JESup: " << effp << " (" << (effp-eff)/eff*100. << "%)" << std::endl;
198	<	std::cout << "    central:  " << eff << std::endl;
199	<	std::cout << "    JESdown: " << effm << " (" << (effm-eff)/eff*100. << "%)" << std::endl;
196	>	if(!automatized) dout << "  Efficiency at JZB==" << jzbSel  << std::endl;
197	>	if(!automatized) dout << "    JESup: " << effp << " (" << (effp-eff)/eff*100. << "%)" << std::endl;
198	>	if(!automatized) dout << "    central:  " << eff << std::endl;
199	>	if(!automatized) dout << "    JESdown: " << effm << " (" << (effm-eff)/eff*100. << "%)" << std::endl;
200		jesup=(effp-eff)/eff*100.;
201		jesdown=(effm-eff)/eff*100.;
202		}
#	Line 199 | Line 205 | void JZBjetScale(TTree *events, float &j
205		// Effect of energy scale on JZB efficiency
206		void doJZBscale(TTree *events, float &down, float &up, float &syst, float systematic, string informalname) {
207
208	<	TCut kbase("abs(genMllSel-91.2)<20&&genZPt>0&&pfJetGoodNum>2");
208	>	TCut kbase("abs(genMllSel-91.2)<20&&genZPt>0&&genNjets>2");
209		TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
210
211		nBins =    50;
#	Line 212 | Line 218 | void doJZBscale(TTree *events, float &do
218		Float_t eff  = Interpolate(jzbSel,hist);
219		Float_t effp = Interpolate(jzbSel*(1.+systematic),hist);
220		Float_t effm = Interpolate(jzbSel*(1.-systematic),hist);
221	<	std::cout << "  efficiency at JZB==" << jzbSel*(1.+systematic)  << "(-"<<syst*100<<"%) : " << effp << " (" << ((effp-eff)/eff)*100. << "%)"  << std::endl;
222	<	std::cout << "  efficiency at JZB==" << jzbSel  << ": " << eff << std::endl;
223	<	std::cout << "  efficiency at JZB==" << jzbSel*(1.-systematic)  << "(-"<<syst*100<<"%) : " << effm << " (" << ((effm-eff)/eff)*100. << "%)"  << std::endl;
221	>	if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.+systematic)  << "(-"<<syst*100<<"%) : " << effp << " (" << ((effp-eff)/eff)*100. << "%)"  << std::endl;
222	>	if(!automatized) dout << "  efficiency at JZB==" << jzbSel  << ": " << eff << std::endl;
223	>	if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.-systematic)  << "(-"<<syst*100<<"%) : " << effm << " (" << ((effm-eff)/eff)*100. << "%)"  << std::endl;
224		up=((effp-eff)/eff)*100;
225		down=((effm-eff)/eff)*100;
226		}
#	Line 224 | Line 230 | void doJZBscale(TTree *events, float &do
230		void JZBresponse(TTree *events, bool isMET = kFALSE, Float_t myJzbMax = 200., Int_t nPeriods = 9 ) {
231
232		jzbMin = 20;
233	<	TCut kbase("abs(genMllSel-91.2)<20&&genZPtSel>0&&pfJetGoodNum>2");
233	>	TCut kbase("abs(genMllSel-91.2)<20&&genZPtSel>0&&genNjets>2");
234		TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
235
236		TProfile* hJzbResp = new TProfile("hJzbResp","JZB response  ; JZB true (GeV/c); JZB reco. / JZB true",
237		nPeriods, jzbMin, myJzbMax, "" );
238
239	<	if (!isMET) events->Project("hJzbResp","jzb[1]/genJZBSel:genJZBSel",kbase&&ksel);
239	>	if (!isMET) events->Project("hJzbResp","("+TString(mcjzbexpression)+")/genJZBSel:genJZBSel",kbase&&ksel);
240		else events->Project("hJzbResp","met[4]/genMET:genMET",kbase&&ksel);
241
242		hJzbResp->SetMaximum(1.2);
243		hJzbResp->SetMinimum(0.2);
244		hJzbResp->Fit("pol0","Q");
245		TF1 *fittedfunction = hJzbResp->GetFunction("pol0");
246	<	cout << "  Response: " << fittedfunction->GetParameter(0) << " +/- " << fittedfunction->GetParError(0) << endl;
246	>	if(!automatized) dout << "  Response: " << fittedfunction->GetParameter(0) << " +/- " << fittedfunction->GetParError(0) << endl;
247	>	delete hJzbResp;
248		}
249
250
251	<	void do_systematics_for_one_file(TTree *events,string informalname, vector<vector<float> > &uncertainties) {
251	>	void do_systematics_for_one_file(TTree *events,string informalname, vector<vector<float> > &results,string mcjzb,string datajzb) {
252
253		float JetEnergyScaleUncert=0.1;
254		float JZBScaleUncert=0.1;
255	+	mcjzbexpression=mcjzb;
256
257		float triggereff=4;//percent!
258	<	cout << "Trigger efficiency not implemented in this script  yet, still using external one" << endl;
258	>	dout << "Trigger efficiency not implemented in this script  yet, still using external one" << endl;
259		float leptonseleff=2;//percent!
260	<	cout << "Lepton selection efficiency not implemented in this script  yet, still using external one" << endl;
260	>	dout << "Lepton selection efficiency not implemented in this script  yet, still using external one" << endl;
261
262		float mceff,mcefferr;
263	<	cout << "MC efficiencies:" << endl;
264	<	MCefficiency(events,mceff,mcefferr);
263	>	if(!automatized) dout << "MC efficiencies:" << endl;
264	>	MCefficiency(events,mceff,mcefferr,mcjzb);
265		JZBefficiency(events,informalname);
266
267	<	std::cout << "Jet energy scale: " << std::endl;
267	>	if(!automatized) dout << "Jet energy scale: " << std::endl;
268		float jesup,jesdown;
269		JZBjetScale(events,jesdown,jesup,informalname,JetEnergyScaleUncert);
270
271	<	std::cout << "JZB scale: " << std::endl;
271	>	if(!automatized) dout << "JZB scale: " << std::endl;
272		float scaleup,scaledown,scalesyst;
273		doJZBscale(events,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname);
274
275	<	std::cout << "JZB response: " << std::endl;
275	>	if(!automatized) dout << "JZB response: " << std::endl;
276		JZBresponse(events);
277
278	<	std::cout << "Pileup: " << std::endl;
278	>	if(!automatized) dout << "Pileup: " << std::endl;
279		float resolution=pileup(events,informalname);
280
281	<	cout << "_______________________________________________" << endl;
282	<	cout << "                 SUMMARY FOR " << informalname << " with JZB>" << jzbSel << endl;
283	<	cout << "Trigger efficiency: " << triggereff << endl;
284	<	cout << "Lepton Sel Eff: " << leptonseleff << endl;
285	<	cout << "For JZB>" << jzbSel << endl;
286	<	cout << "Jet energy scale: " << jesup << " " << jesdown << " --> suggesting: " << Round(0.5*(fabs(jesup)+fabs(jesdown)),1) << endl;
287	<	cout << "JZB Scale Uncert: " << scaledown << " " << scaleup << " --> suggesting: " << Round(0.5*(fabs(scaledown)+fabs(scaleup)),1) << endl;
288	<	cout << "Resolution : " << resolution << endl;
281	>	dout << "_______________________________________________" << endl;
282	>	dout << "                 SUMMARY FOR " << informalname << " with JZB>" << jzbSel << endl;
283	>	dout << "MC efficiency: " << mceff << "+/-" << mcefferr << endl;
284	>	dout << "Trigger efficiency: " << triggereff << endl;
285	>	dout << "Lepton Sel Eff: " << leptonseleff << endl;
286	>	dout << "For JZB>" << jzbSel << endl;
287	>	dout << "Jet energy scale: " << jesup << " " << jesdown << " --> suggesting: " << Round(0.5*(fabs(jesup)+fabs(jesdown)),1) << endl;
288	>	dout << "JZB Scale Uncert: " << scaledown << " " << scaleup << " --> suggesting: " << Round(0.5*(fabs(scaledown)+fabs(scaleup)),1) << endl;
289	>	dout << "Resolution : " << resolution << endl;
290	>
291
292	<	vector<float> uncert;
293	<	uncert.push_back(jzbSel);
294	<	uncert.push_back(triggereff);
295	<	uncert.push_back(leptonseleff);
296	<	uncert.push_back(0.5*(fabs(jesup)+fabs(jesdown)));
297	<	uncert.push_back(0.5*(fabs(scaledown)+fabs(scaleup)));
298	<	uncert.push_back(resolution);
292	>	float toterr=0;
293	>	toterr+=(triggereff/100)*(triggereff/100);
294	>	toterr+=(leptonseleff/100)*(leptonseleff/100);
295	>	if(fabs(jesup)>fabs(jesdown)) toterr+=(jesup/100)*(jesup/100); else toterr+=(jesdown/100)*(jesdown/100);
296	>	if(fabs(scaleup)>fabs(scaledown)) toterr+=(scaleup/100)*(scaleup/100); else toterr+=(scaledown/100)*(scaledown/100);
297	>	toterr+=(resolution/100)*(resolution/100);
298	>	toterr=TMath::Sqrt(toterr);
299	>	dout << "FINAL RESULT : " << mceff << " +/- "<< mcefferr << " (stat) +/- " << 100*toterr << " (syst)" << endl;
300	>	dout << "     we thus use the sqrt of the sum of the squares which is : " << 100*TMath::Sqrt(mcefferr*mcefferr+(toterr*toterr)) << endl;
301	>	vector<float> res;
302	>	res.push_back(jzbSel);
303	>	res.push_back(mceff);
304	>	res.push_back(mcefferr);
305	>	res.push_back(toterr);
306	>	res.push_back(TMath::Sqrt((mcefferr)*(mcefferr)+(toterr*toterr)));
307
308	<	uncertainties.push_back(uncert);
308	>	results.push_back(res);
309		}
310
311		vector<vector<float> > compute_systematics(string mcjzb, string datajzb, samplecollection &signalsamples, vector<float> bins) {
312	<	vector< vector<float> > uncertainties;
312	>	automatized=true;
313	>	vector< vector<float> > systematics;
314		for (int isignal=0; isignal<signalsamples.collection.size();isignal++) {
315	<	cout << "Looking at signal " << (signalsamples.collection)[isignal].filename << endl;
315	>	dout << "Looking at signal " << (signalsamples.collection)[isignal].filename << endl;
316		for(int ibin=0;ibin<bins.size();ibin++) {
317		jzbSel=bins[ibin];
318	<	do_systematics_for_one_file((signalsamples.collection)[isignal].events,(signalsamples.collection)[isignal].samplename,uncertainties);
318	>	geqleq="geq";
319	>	do_systematics_for_one_file((signalsamples.collection)[isignal].events,(signalsamples.collection)[isignal].samplename,systematics,mcjzb,datajzb);
320		}//end of bin loop
321		}//end of signal loop
322	<	return uncertainties;
322	>	return systematics;
323		}

Diff Legend

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