ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/ZbTools.C

Comparing UserCode/cbrown/Development/Plotting/Modules/ZbTools.C (file contents):
Revision 1.9 by buchmann, Wed Nov 14 10:46:45 2012 UTC vs.
Revision 1.17 by buchmann, Thu Nov 29 16:51:20 2012 UTC

#	Line 23 | Line 23 | using namespace PlottingSetup;
23		namespace ZbData {
24		vector<float> data_over_mc;
25
26	<	TCut ZplusBsel("pt1>20&&pt2>20&&mll>60&&mll<120");
26	>	TCut ZplusBsel("pt1>20&&pt2>20&&mll>60&&mll<120&&id1==id2");
27		TCut LeadingB("Zb3010_bTagProbCSVBP[0]>0.898");
28		TCut EtaB("abs(Zb3010_pfJetGoodEta[0])<1.3");
29	<	TCut PhiZcut("abs(pfJetDphiZ[0])>2.7");
29	>	TCut PhiZcut("abs(Zb3010_pfJetDphiZ[0])>2.7");
30
31		}
32
33	+	class ZbPointResult {
34	+	public:
35	+	ZbPointResult(string var, float ptlow, float pthigh, float alphalow, float alphahigh, Value Data, Value MC);
36	+	float PtLow;
37	+	float PtHigh;
38	+	float AlphaLow;
39	+	float AlphaHigh;
40	+	string variable;
41	+	Value Data;
42	+	Value MC;
43	+	Value Ratio;
44	+	};
45	+
46	+	ZbPointResult::ZbPointResult(string var, float ptlow, float pthigh, float alphalow, float alphahigh, Value data, Value mc) :
47	+	variable(var), PtLow(ptlow), PtHigh(pthigh), AlphaLow(alphalow), AlphaHigh(alphahigh), Data(data), MC(mc) {
48	+	Ratio = data/mc;
49	+	}
50	+
51	+	std::ostream &operator<<(std::ostream &ostr, ZbPointResult &b)
52	+	{
53	+	ostr << b.PtLow << ";" << b.PtHigh << ";" << b.AlphaLow << ";" << b.AlphaHigh << ";" << b.Data.getValue() << ";" << b.Data.getError() << ";" << b.MC.getValue() << ";" << b.MC.getError() << ";" << b.Ratio.getValue() << ";" << b.Ratio.getError();
54	+	return ostr;
55	+	}
56	+
57	+
58	+	class ZbResultContainer {
59	+	public:
60	+	string name;
61	+	vector<ZbPointResult> MPF_Rabs_location;
62	+	Value MPF_Result_FaceValue;
63	+	Value Rabs_Result_FaceValue;
64	+	Value MPF_Result_Extrapolated;
65	+	Value Rabs_Result_Extrapolated;
66	+
67	+	ZbResultContainer(string name);
68	+	ZbResultContainer(const ZbResultContainer &old);
69	+	};
70	+
71	+
72	+	ZbResultContainer::ZbResultContainer(string _name) {
73	+	name=_name;
74	+	cout<< "Result Container for \"" << name << "\" has been initialized" << endl;
75	+
76	+	}
77	+
78	+	std::ostream &operator<<(std::ostream &ostr, ZbResultContainer &b)
79	+	{
80	+	ostr << "Results for " << b.name << endl;
81	+	ostr << " MPF and Rabs locations: " << endl;
82	+	for(int i=0;i<b.MPF_Rabs_location.size();i++) ostr << "   " << b.MPF_Rabs_location[i] << endl;
83	+	ostr << " Results: " << endl;
84	+	ostr << "   Face Value: " << endl;
85	+	ostr << "     MPF " << b.MPF_Result_FaceValue << endl;
86	+	ostr << "     Rabs " << b.Rabs_Result_FaceValue << endl;
87	+	ostr << "   Extrapolated: " << endl;
88	+	ostr << "     MPF " << b.MPF_Result_Extrapolated << endl;
89	+	ostr << "             Rabs " << b.Rabs_Result_Extrapolated << endl;
90	+
91	+	return ostr;
92	+	}
93	+
94	+
95	+
96		using namespace ZbData;
97
98
#	Line 101 | Line 164 | void draw_kin_variable(string variable,
164		delete ckin;
165		}
166
104	–	Value getfrom2Dmap(TH2F *map, int ixbin, int iybin) {
105	–	Value sum(0,0);
106	–	for(int iy=1;iy<=iybin;iy++) {
107	–	sum=sum+Value(map->GetBinContent(ixbin,iy),map->GetBinError(ixbin,iy));
108	–	}
109	–	return sum;
110	–	}
111	–
167		void print_all_b_yields() {
168		cout << "Basic selection with a b jet" << endl;
169		print_yield(ZplusBsel&&"Zb3010_pfJetGoodNumBtag>0");
#	Line 120 | Line 175 | void print_all_b_yields() {
175		print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
176		cout << "10<ptZ<1000 GeV" << endl;
177		print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000");
178	+	cout << "#alpha < 0.35" << endl;
179	+	print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.35");
180		cout << "#alpha < 0.3" << endl;
181		print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.3");
182		cout << "#alpha < 0.2" << endl;
#	Line 132 | Line 189 | void print_all_b_yields() {
189		//  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.05");
190		}
191
192	<	void DrawEvilCutFlow() {
192	>	void TEST_DrawEvilCutFlow() {
193		stringstream MegaCut;
194
195		MegaCut << "   1*(" << (const char*) (ZplusBsel&&TCut("Zb3010_pfJetGoodNumBtag>0")) << ")";
#	Line 220 | Line 277 | void draw_Zb_kin_vars() {
277		draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>125&&pt<1000",40,0,2,"#alpha","Official/alpha_pt_125_to_1000",1);
278		draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",40,0,2,"#alpha","Official/alpha_pt_30_to_1000",1);
279
280	<	draw_kin_variable("pfBJetDphiZ[0]",ZplusBsel&&LeadingB&&"pt>10&&pt<1000&&pfBJetDphiZ[0]>0",30,0,3.2,"#delta#phi (Z,b lead)","DeltaPhiZBlead",0);
280	>	draw_kin_variable("Zb3010_pfBJetDphiZ[0]",ZplusBsel&&LeadingB&&"pt>10&&pt<1000&&pfBJetDphiZ[0]>0",30,0,3.2,"#delta#phi (Z,b lead)","DeltaPhiZBlead",0);
281
282		/*
283		draw_kin_variable("pt1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,200,"p_{T}^{l1}","Lep/pt1",1);
#	Line 254 | Line 311 | void draw_mpf_vars() {
311		draw_kin_variable("Zb3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.05",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz__mpf_alpha_smaller_0p05",0,"#alpha<0.05, 10 GeV < p_{T}^{Z} < 1000 GeV");
312		}
313
314	<	Value get_Zb_data_over_mc(string variable, TCut cut) {
314	>	/*double GetMedian(TH1F *histo) {
315	>	int numBins = histo->GetXaxis()->GetNbins();
316	>	Double_t x[numBins];
317	>	Double_t y[numBins];
318	>	for (int i = 1; i <= numBins; i++) {
319	>	x[i] = histo->GetBinCenter(i);
320	>	y[i] = histo->GetBinContent(i);
321	>	}
322	>	return TMath::Median(numBins, &x, &y);
323	>	}*/
324	>
325	>	Value get_Zb_data_over_mc(string ContainerName, ZbResultContainer &res, string variable, string variable2, TCut cut, string saveas, float pt1, float pt2, float a1, float a2) {
326		//  write_warning(__FUNCTION__,"Debugging this function, therefore always returning 3.1415+/-1.010101"); return Value(3.1415,1.010101);
327	<	TH1F *hdata  = allsamples.Draw("hdata",variable,1000,0,30, "#alpha [GeV]", "events", cut,data,luminosity);
328	<	TH1F *hmc    = allsamples.Draw("hmc" , variable,1000,0,30, "#alpha [GeV]", "events", cut,mc,  luminosity);
327	>	TCanvas *cn = new TCanvas("cn","cn");
328	>	string varname="MPF";
329	>	if(Contains(variable,"pfJetGood")) varname="R_{abs}";
330	>	TH1F *hdata  = allsamples.Draw("hdata",variable,1200,0,30, varname, "events",  cut,data,luminosity);
331	>	TH1F *hmc    = allsamples.Draw("hmc" , variable,1200,0,30, varname, "events", cut, mc,  luminosity);
332	>	hmc->SetLineColor(kBlue);
333
334	<	hdata->Fit("gaus","QE");
263	<	hmc->Fit("gaus","QE");
264	<	TF1 *fdata = (TF1*) hdata->GetFunction("gaus");
265	<	TF1 *fmc   = (TF1*) hmc->GetFunction("gaus");
266	<
267	<	/*  float a=hdata->GetMean();
334	>	float a=hdata->GetMean();
335		float b=hmc->GetMean();
336		float da=hdata->GetMeanError();
337	<	float db=hmc->GetMeanError();*/
271	<	float a = fdata->GetParameter(1);
272	<	float da = fdata->GetParError(1);
273	<	float b = fmc->GetParameter(1);
274	<	float db = fmc->GetParError(1);
337	>	float db=hmc->GetMeanError();
338		float factor = a / b;
339		float error = TMath::Sqrt(  (1/(b*b)) * (da*da) + ((a*a)/(b*b))*db*db);
340	+
341	+	TF1 *gaus_data = new TF1("gaus_data","gaus",0.7,1.3);
342	+	gaus_data->SetParameter(0,hdata->GetMaximum());
343	+	gaus_data->SetParameter(2,hdata->GetMean());
344	+	gaus_data->SetParameter(1,hdata->GetRMS());
345	+	hdata->Fit(gaus_data);
346	+	TF1 *gaus_mc = new TF1("gaus_mc","gaus",0.7,1.3);
347	+	gaus_mc->SetParameter(0,hmc->GetMaximum());
348	+	gaus_mc->SetParameter(1,hmc->GetMean());
349	+	gaus_mc->SetParameter(2,hmc->GetRMS());
350	+	hmc->Fit(gaus_mc);
351	+
352	+
353	+	cn->cd();
354	+	hdata->GetYaxis()->SetTitle("events / 0.1");
355	+	hdata->SetMarkerColor(kRed);
356	+	/*  hdata->Rebin(4);
357	+	hmc->Rebin(4);
358	+
359	+	gaus_mc->SetParameter(0,gaus_mc->GetParameter(0)*4);//rebinning
360	+	gaus_data->SetParameter(0,gaus_data->GetParameter(0)*4);//rebinning
361	+	*/
362	+	gaus_mc->SetLineColor(kBlue);
363	+	gaus_data->SetLineColor(kRed);
364	+
365	+	hdata->GetXaxis()->SetRangeUser(0,2);
366	+	hdata->GetYaxis()->SetTitleOffset(1.3);
367	+	hdata->Draw("e1");
368	+	hmc->Draw("histo,same");
369	+	hdata->Draw("e1,same");
370	+	gaus_data->Draw("same");
371	+	gaus_mc->Draw("same");
372	+	stringstream summary;
373	+	summary << "#splitline{#bf{" << varname << "}}{#splitline{" << pt1 << " < p_{T}^{Z} < " << pt2 << ", " << a1 << " #leq #alpha < " << a2 << "}{#splitline{data: " << std::setprecision(4) << a << " #pm " << std::setprecision(4) << da << "}{";
374	+	summary << "#splitline{mc:   " << std::setprecision(4) << b << " #pm " << std::setprecision(4) << db << "}{#splitline{ratio: " << factor;
375	+	summary << " #pm " << error << "}{#splitline{}{#splitline{data fit: " << std::setprecision(4) << gaus_data->GetParameter(1) << "+/-";
376	+	summary << std::setprecision(4) << gaus_data->GetParError(1) << "}{#splitline{";
377	+	summary << "mc fit:" << std::setprecision(4) << gaus_mc->GetParameter(1) << "#pm" << std::setprecision(4) << gaus_mc->GetParError(1) << "}{ratio: ";
378	+	summary << gaus_data->GetParameter(1)/gaus_mc->GetParameter(1) << "#pm" << gaus_data->GetParameter(1)/gaus_mc->GetParameter(1) * sqrt((gaus_data->GetParError(1) * gaus_data->GetParError(1))/(gaus_data->GetParameter(1)*gaus_data->GetParameter(1)) + (gaus_mc->GetParError(1) * gaus_mc->GetParError(1))/(gaus_mc->GetParameter(1)*gaus_mc->GetParameter(1)) ) << "}}}}}}}}";
379	+
380	+	TText *infobox = write_title(summary.str());
381	+	infobox->SetX(0.75);
382	+	infobox->SetTextSize(0.03);
383	+	infobox->SetY(0.75);
384	+	infobox->Draw();
385	+
386	+
387	+	DrawPrelim();
388	+	CompleteSave(cn,ContainerName+"/ResponsePlots/"+saveas);
389	+
390	+	res.MPF_Rabs_location.push_back(ZbPointResult(variable,pt1,pt2,a1,a2,Value(a,da),Value(b,db)));
391	+	cout << "Data;" << a << ";"<<da<<";MC"<<b<<";"<<db<<endl;
392	+	delete cn;
393		delete hdata;
394		delete hmc;
395
280	–	//    cout << a << "+/-" << da << "        " << b << "=/-" << db << endl;
281	–
396		return Value(factor,error);
397		}
398
399	+	int nFakes=1;
400	+	ZbResultContainer Fake_new_data_mc_agreement_2d(string AlphaVariable, string ContainerName, TCut BTagCut, TCut BTagWgt) {
401	+	write_warning(__FUNCTION__,"You are using the fake way to extract the correction ... ");
402	+	nFakes++;
403	+	ZbResultContainer results(ContainerName);
404	+
405
406	+	results.MPF_Result_FaceValue=Value(1+nFakes*0.02,0.03);
407	+	results.Rabs_Result_FaceValue=Value(2,0.02);
408	+	results.MPF_Result_Extrapolated=Value(3,0.03);
409	+	results.Rabs_Result_Extrapolated=Value(4,0.04);
410	+	cout << "The fake result for " << ContainerName << " for the MPF face value is " << results.MPF_Result_FaceValue << endl;
411
412	<	void new_data_mc_agreement_2d() {
412	>	return results;
413	>	}
414	>
415	>
416	>	ZbResultContainer new_data_mc_agreement_2d(string AlphaVariable, string ContainerName, TCut BTagCut, TCut BTagWgt) {
417	>	ZbResultContainer results(ContainerName);
418	>
419	>	cutWeight=TCut("(weight*(weight<1000)*(is_data+(!is_data)))")*BTagWgt;
420	>	LeadingB=BTagCut;
421	>
422		gStyle->SetOptFit(0);
423		TCut basecut(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
424	<	const int nalphacuts=5;
425	<	float alphacuts[nalphacuts] = {0.1,0.15,0.2,0.25,0.3};
292	<	int PointThree=-1;
293	<	for(int i=0;i<nalphacuts;i++) {
294	<	if((alphacuts[i]-0.3)<0.01) PointThree=i-1;
295	<	}
296	<	if(PointThree<0) {
297	<	write_warning(__FUNCTION__,"alpha=0.3 MUST BE ONE OF THE ALPHA VALUES!");
298	<	assert(PointThree<0);
299	<	}
424	>	const int nalphacuts=6;
425	>	float alphacuts[nalphacuts] = {0.1,0.15,0.2,0.25,0.3,0.35};
426
427		const int nptcuts=5;
428		float ptcuts[nptcuts]={30,50,75,125,1000};
429
430	+
431	+
432		float MPF_Results[nalphacuts][nptcuts];
433		float MPF_Errors[nalphacuts][nptcuts];
434		float RABS_Results[nalphacuts][nptcuts];
#	Line 310 | Line 438 | void new_data_mc_agreement_2d() {
438		for(int ia=0;ia<nalphacuts;ia++) {
439		for(int ipt=0;ipt<nptcuts-1;ipt++) {
440		stringstream specialcut;
441	<	specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (Zb3010_alpha<" << alphacuts[ia] << "))";
441	>	float alow,ahigh;
442	>	if(ia>0) {
443	>	specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (" << AlphaVariable << "<" << alphacuts[ia] << " && " << AlphaVariable << ">" << alphacuts[ia-1] << "))";
444	>	alow=alphacuts[ia-1];
445	>	ahigh=alphacuts[ia];
446	>	} else {
447	>	specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (" << AlphaVariable << "<" << alphacuts[ia] << "))";
448	>	alow=0;
449	>	ahigh=alphacuts[ia];
450	>	}
451		cout << specialcut.str() << endl;
452	<	Value MPF_data_over_mc = get_Zb_data_over_mc("mpf",TCut(basecut && specialcut.str().c_str()));
453	<	Value RABS_data_over_mc = get_Zb_data_over_mc("Zb3010_pfJetGoodPt[0]/pt",TCut(basecut && specialcut.str().c_str()));
452	>	Value MPF_data_over_mc = get_Zb_data_over_mc(ContainerName,results,"mpf","",TCut(basecut && specialcut.str().c_str()),"MPF___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__alpha_"+any2string(alphacuts[ia]), ptcuts[ipt], ptcuts[ipt+1],alow,ahigh);
453	>	Value RABS_data_over_mc = get_Zb_data_over_mc(ContainerName,results,"Zb3010_pfJetGoodPt[0]/pt","",TCut(basecut && specialcut.str().c_str()),"Rabs___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__alpha_"+any2string(alphacuts[ia]), ptcuts[ipt], ptcuts[ipt+1],alow,ahigh);
454
455		MPF_Results[ia][ipt]=MPF_data_over_mc.getValue();
456		MPF_Errors[ia][ipt]=MPF_data_over_mc.getError();
#	Line 329 | Line 466 | void new_data_mc_agreement_2d() {
466		float RABS_ExtraPolatedResults[nptcuts];
467
468		TCanvas *can = new TCanvas("can");
469	+	can->cd(1)->SetLeftMargin(0.15);
470
471		TGraphErrors *MPF_FinalGraph = new TGraphErrors();
472		TGraphErrors *RABS_FinalGraph = new TGraphErrors();
473	<	TGraphErrors *FaceValueAtPoint3 = new TGraphErrors();
473	>	TGraphErrors *MPF_FaceValueAtPoint3 = new TGraphErrors();
474	>	TGraphErrors *RABS_FaceValueAtPoint3 = new TGraphErrors();
475
476		for(int ipt=0;ipt<nptcuts-1;ipt++) {
477
#	Line 341 | Line 480 | void new_data_mc_agreement_2d() {
480		cout << "Going to create histo for " << filename.str() << endl;
481		TGraphErrors *mpf_gr =new TGraphErrors();
482		TGraphErrors *rabs_gr =new TGraphErrors();
483	+
484	+	int rabs_counter=0;
485	+	int mpf_counter=0;
486	+
487		for(int ia=0;ia<nalphacuts;ia++) {
488		cout << "Setting a point for an alpha cut at " << alphacuts[ia] << " with value " << MPF_Results[ia][ipt] << " (MPF) and " << RABS_Results[ia][ipt] << " (RABS)" << endl;
489	<	mpf_gr->SetPoint(ia,alphacuts[ia],MPF_Results[ia][ipt]);     // avoiding empty bins (nan's)
490	<	mpf_gr->SetPointError(ia,0,MPF_Errors[ia][ipt]);             // avoiding empty bins (nan's)
491	<	rabs_gr->SetPoint(ia,alphacuts[ia],RABS_Results[ia][ipt]); // avoiding empty bins (nan's)
492	<	rabs_gr->SetPointError(ia,0,RABS_Errors[ia][ipt]);         // avoiding empty bins (nan's)
493	<	}
494	<	cout << "About to add a point in pt (" << ptcuts[ipt]+ptcuts[ipt+1] << " ) to face value : " << MPF_Results[PointThree][ipt] << endl;
495	<	FaceValueAtPoint3->SetPoint(ipt,(ptcuts[ipt]+ptcuts[ipt+1]),MPF_Results[PointThree][ipt]);
496	<	FaceValueAtPoint3->SetPointError(ipt,0,MPF_Errors[PointThree][ipt]);
489	>	if(MPF_Results[ia][ipt]==MPF_Results[ia][ipt] && MPF_Errors[ia][ipt]==MPF_Errors[ia][ipt]) { // remove nan's
490	>	mpf_gr->SetPoint(mpf_counter,alphacuts[ia],MPF_Results[ia][ipt]);
491	>	mpf_gr->SetPointError(mpf_counter,0,MPF_Errors[ia][ipt]);
492	>	mpf_counter++;
493	>	} else {
494	>	dout << "Detected BS for MPF method (bin ignored) : " << endl;
495	>	dout << "       alpha: " << alphacuts[ia] << " , pt range: " << ptcuts[ipt] << " < pt < " << ptcuts[ipt+1] << endl;
496	>	cout << "       value (MPF) : " << MPF_Results[ia][ipt] << " +/- " << MPF_Errors[ia][ipt] << endl;
497	>	}
498	>	if(RABS_Results[ia][ipt]==RABS_Results[ia][ipt] && RABS_Errors[ia][ipt]==RABS_Errors[ia][ipt]) {
499	>	rabs_gr->SetPoint(rabs_counter,alphacuts[ia],RABS_Results[ia][ipt]);
500	>	rabs_gr->SetPointError(rabs_counter,0,RABS_Errors[ia][ipt]);
501	>	rabs_counter++;
502	>	} else {
503	>	dout << "Detected BS for RABS method (bin ignored) : " << endl;
504	>	dout << "       alpha: " << alphacuts[ia] << " , pt range: " << ptcuts[ipt] << " < pt < " << ptcuts[ipt+1] << endl;
505	>	cout << "       value (RABS) : " << RABS_Results[ia][ipt] << " +/- " << RABS_Errors[ia][ipt] << endl;
506	>	}
507	>	}
508
509	+	stringstream specialcut;
510	+	specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (" << AlphaVariable << "<0.3))";
511	+	Value MPF_data_over_mc = get_Zb_data_over_mc(ContainerName,results,"mpf","",TCut(basecut && specialcut.str().c_str()),"MPF___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__INCLUSIVEalpha_0.3",ptcuts[ipt],ptcuts[ipt+1],0,0.3);
512	+	Value RABS_data_over_mc = get_Zb_data_over_mc(ContainerName,results,"Zb3010_pfJetGoodPt[0]/pt","",TCut(basecut && specialcut.str().c_str()),"Rabs___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__INCLUSIVEalpha_0.3",ptcuts[ipt],ptcuts[ipt+1],0,0.3);
513	+	cout << "About to add a point in pt (" << (ptcuts[ipt]+ptcuts[ipt+1])/2 << " ) to face value : " << MPF_data_over_mc << " (MPF) , " <<  RABS_data_over_mc << "(RABS)" << endl;
514	+	MPF_FaceValueAtPoint3->SetPoint(ipt,(ptcuts[ipt]+ptcuts[ipt+1])/2,MPF_data_over_mc.getValue());
515	+	MPF_FaceValueAtPoint3->SetPointError(ipt,0,MPF_data_over_mc.getError());
516	+	RABS_FaceValueAtPoint3->SetPoint(ipt,(ptcuts[ipt]+ptcuts[ipt+1])/2,RABS_data_over_mc.getValue());
517	+	RABS_FaceValueAtPoint3->SetPointError(ipt,0,RABS_data_over_mc.getError());
518	+
519	+	can->cd();
520		mpf_gr->SetMarkerStyle(21);
521		mpf_gr->SetMarkerColor(kBlue);
522		mpf_gr->Draw("AP");
523	<	mpf_gr->Fit("pol1","QE"); // quiet, use minos
523	>	mpf_gr->Fit("pol1","");
524		mpf_gr->GetXaxis()->SetTitle("#alpha");
525		mpf_gr->GetXaxis()->CenterTitle();
526		mpf_gr->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
527		mpf_gr->GetYaxis()->CenterTitle();
528		TF1 *mpf_pol=(TF1*)mpf_gr->GetFunction("pol1");
529	+	mpf_pol->SetLineWidth(0);
530	+	TF1 *mpf_pol_complete = new TF1("mpf_pol_complete","[0]+[1]*x");
531	+	mpf_pol_complete->SetParameters(mpf_pol->GetParameters());
532	+	mpf_pol_complete->SetLineWidth(1);
533	+	mpf_pol_complete->SetLineColor(kBlue);
534	+	float min,max;
535	+	FindMinMax(mpf_gr,min,max);
536	+	TH1F *histo = new TH1F("histo","histo",1,0,0.4);
537	+	histo->GetXaxis()->SetTitle("#alpha");
538	+	histo->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
539	+	histo->GetXaxis()->CenterTitle();
540	+	histo->GetYaxis()->CenterTitle();
541	+	histo->GetYaxis()->SetTitleOffset(1.3);
542	+	histo->SetStats(0);
543	+	histo->GetXaxis()->SetRangeUser(0,0.4);
544	+	histo->GetYaxis()->SetRangeUser(min,max);
545	+	mpf_gr->GetYaxis()->SetRangeUser(min,max);
546	+
547	+	TPolyLine *mpf_fit_uncert = GetFitUncertaintyShape(mpf_gr, "pol1", min, max,0.0,0.4);
548	+	mpf_pol->SetLineColor(TColor::GetColor("#0101DF"));
549	+	mpf_fit_uncert->SetFillColor(TColor::GetColor("#5353FC"));
550	+	histo->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
551	+	histo->Draw();
552	+	mpf_fit_uncert->Draw("F");
553	+	mpf_fit_uncert->Draw("");
554	+	mpf_gr->Draw("P");
555	+	mpf_pol_complete->Draw("same");
556	+
557		float mpf_a=mpf_pol->GetParameter(0);
558		float mpf_b=mpf_pol->GetParameter(1);
559		MPF_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),mpf_a);
#	Line 369 | Line 562 | void new_data_mc_agreement_2d() {
562		MPF_ExtraPolatedResults[ipt]=mpf_a;
563
564		stringstream mpf_info;
565	<	mpf_info <<  "#splitline{#splitline{#splitline{" << ptcuts[ipt] << " GeV < p_{T}^{Z} < " << ptcuts[ipt+1] << " GeV }{ (MPF) }}{Extrapolated value: " << std::setprecision(3) << mpf_a << "}}{#chi^{2}/NDF : " << mpf_pol->GetChisquare() << " / " << mpf_pol->GetNDF() << "}";
565	>	mpf_info <<  "#splitline{#splitline{#splitline{" << ptcuts[ipt] << " GeV < p_{T}^{Z} < " << ptcuts[ipt+1] << " GeV }{ (MPF) }}{Extrapolated value: " << std::setprecision(4) << mpf_a << "}}{#chi^{2}/NDF : " << mpf_pol->GetChisquare() << " / " << mpf_pol->GetNDF() << "}";
566
567		TText *mpf_mark = write_title(mpf_info.str());
568		mpf_mark->SetX(0.75);
#	Line 380 | Line 573 | void new_data_mc_agreement_2d() {
573		string filenamebkp=filename.str();
574		filename << "__MPF";
575		DrawPrelim();
576	<	CompleteSave(can,filename.str());
576	>	CompleteSave(can,ContainerName+"/"+filename.str());
577		cout << "MPF : " << mpf_a << " + " << mpf_b << " * alpha " << endl;
578
579		filename.str("");
580
581		rabs_gr->SetMarkerStyle(21);
582	<	rabs_gr->SetMarkerColor(kBlue);
390	<	rabs_gr->Fit("pol1","QE"); // quiet, use minos
582	>	rabs_gr->SetMarkerColor(kRed);
583		rabs_gr->Draw("AP");
584	+	rabs_gr->Fit("pol1","");
585		rabs_gr->GetXaxis()->SetTitle("#alpha");
586		rabs_gr->GetXaxis()->CenterTitle();
587		rabs_gr->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
588		rabs_gr->GetYaxis()->CenterTitle();
396	–
589		TF1 *rabs_pol=(TF1*)rabs_gr->GetFunction("pol1");
590	+	rabs_pol->SetLineWidth(0);
591	+	TF1 *rabs_pol_complete = new TF1("rabs_pol_complete","[0]+[1]*x");
592	+	rabs_pol_complete->SetParameters(rabs_pol->GetParameters());
593	+	rabs_pol_complete->SetLineColor(kRed);
594	+	rabs_pol_complete->SetLineWidth(1);
595	+	FindMinMax(rabs_gr,min,max);
596	+	rabs_gr->GetYaxis()->SetRangeUser(min,max);
597	+	histo->GetYaxis()->SetRangeUser(min,max);
598	+	TPolyLine *rabs_fit_uncert = GetFitUncertaintyShape(rabs_gr, "pol1", min, max,0.0,0.4);
599	+	rabs_pol->SetLineColor(TColor::GetColor("#FA5858"));
600	+	rabs_pol->SetFillColor(TColor::GetColor("#FA5858"));
601	+	rabs_fit_uncert->SetLineColor(TColor::GetColor("#FA5858"));
602	+	rabs_fit_uncert->SetFillColor(TColor::GetColor("#FA5858"));
603	+	histo->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
604	+	histo->Draw();
605	+	rabs_fit_uncert->Draw("F");
606	+	rabs_fit_uncert->Draw("");
607	+	rabs_gr->Draw("P");
608	+	rabs_pol_complete->Draw("same");
609	+
610	+
611		float rabs_a=rabs_pol->GetParameter(0);
612		float rabs_b=rabs_pol->GetParameter(1);
613		RABS_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),rabs_a);
#	Line 402 | Line 615 | void new_data_mc_agreement_2d() {
615
616		cout << "!+!+!+!+!+!+!+!+!+ Just added a point to the final plot : " << rabs_a << " +/- " << rabs_pol->GetParError(0) << endl;
617
405	–
618		stringstream rabs_info;
619	<	rabs_info << "#splitline{#splitline{#splitline{" << ptcuts[ipt] << " GeV < p_{T}^{Z} < " << ptcuts[ipt+1] << " GeV }{ (R_{abs}) }}{Extrapolated value: " << std::setprecision(3) << rabs_a << "}}{#chi^{2}/NDF : " << rabs_pol->GetChisquare() << " / " << rabs_pol->GetNDF() << "}";
619	>	rabs_info << "#splitline{#splitline{#splitline{" << ptcuts[ipt] << " GeV < p_{T}^{Z} < " << ptcuts[ipt+1] << " GeV }{ (R_{abs}) }}{Extrapolated value: " << std::setprecision(4) << rabs_a << "}}{#chi^{2}/NDF : " << rabs_pol->GetChisquare() << " / " << rabs_pol->GetNDF() << "}";
620		TText *rabs_mark = write_title(rabs_info.str());
621		rabs_mark->SetX(0.75);
622		rabs_mark->SetTextSize(0.03);
#	Line 412 | Line 624 | void new_data_mc_agreement_2d() {
624		rabs_mark->Draw();
625
626
415	–
627		RABS_ExtraPolatedResults[ipt]=rabs_a;
628		filename << filenamebkp << "__RABS";
629		DrawPrelim();
630	<	CompleteSave(can,filename.str());
630	>	CompleteSave(can,ContainerName+"/"+filename.str());
631		cout << "RABS : " << rabs_a << " + " << rabs_b << " * alpha " << endl;
632		}
633
634		MPF_FinalGraph->SetMarkerStyle(21);
635	<	MPF_FinalGraph->SetMarkerColor(kRed);
636	<	MPF_FinalGraph->Fit("pol0","QE"); // quiet, use minos
635	>	MPF_FinalGraph->SetMarkerColor(kBlue);
636	>	MPF_FinalGraph->Fit("pol0",""); // quiet, use minos
637
638		TF1 *mpf_pol0=(TF1*)MPF_FinalGraph->GetFunction("pol0");
639	+	TF1 *mpf_pol0_complete = new TF1("mpf_pol0_complete","[0]+[1]*x");
640	+	mpf_pol0_complete->SetLineWidth(1);
641	+	mpf_pol0_complete->SetLineColor(kBlue);
642	+	mpf_pol0->SetLineColor(kBlue);
643		MPF_FinalGraph->Draw("AP");
644		MPF_FinalGraph->GetXaxis()->SetTitle("p_{T}^{Z}");
645		MPF_FinalGraph->GetXaxis()->CenterTitle();
646		MPF_FinalGraph->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
647		MPF_FinalGraph->GetYaxis()->CenterTitle();
648		MPF_FinalGraph->Draw("AP");
649	+	mpf_pol0_complete->Draw("same");
650
651		stringstream mpf_result;
652		mpf_result << "#splitline{C_{abs}= " << std::setprecision(5) << 1/mpf_pol0->GetParameter(0) << " #pm " << std::setprecision(2) << mpf_pol0->GetParError(0)/(mpf_pol0->GetParameter(0)*mpf_pol0->GetParameter(0)) << "}{ (MPF)}";
#	Line 445 | Line 661 | void new_data_mc_agreement_2d() {
661
662		stringstream filename2;
663		filename2 << "Extrapolation/FINAL_RESULT_MPF";
664	<	CompleteSave(can,filename2.str());
664	>	CompleteSave(can,ContainerName+"/"+filename2.str());
665
666
667		RABS_FinalGraph->SetMarkerStyle(21);
#	Line 473 | Line 689 | void new_data_mc_agreement_2d() {
689
690		filename2.str("");
691		filename2 << "Extrapolation/FINAL_RESULT_RABS";
692	<	CompleteSave(can,filename2.str());
692	>	CompleteSave(can,ContainerName+"/"+filename2.str());
693
694		can->cd();
695	<	TCanvas *FaceValueCanvas = new TCanvas("FaceValueCanvas","FaceValueCanvas");
696	<	FaceValueAtPoint3->SetMarkerStyle(21);
697	<	FaceValueAtPoint3->SetMarkerStyle(21);
698	<	FaceValueAtPoint3->Fit("pol0","QE"); // quiet, use minos
699	<	FaceValueAtPoint3->Draw("AP");
700	<	FaceValueAtPoint3->GetXaxis()->SetTitle("p_{T}^{Z}");
701	<	FaceValueAtPoint3->GetXaxis()->CenterTitle();
702	<	FaceValueAtPoint3->GetYaxis()->SetTitle("data/mc for #alpha<0.3");
703	<	FaceValueAtPoint3->GetYaxis()->CenterTitle();
704	<	FaceValueAtPoint3->Draw("AP");
705	<	TF1 *faceval_pol0=(TF1*)FaceValueAtPoint3->GetFunction("pol0");
695	>	MPF_FaceValueAtPoint3->SetMarkerStyle(21);
696	>	MPF_FaceValueAtPoint3->SetMarkerColor(kBlue);
697	>	MPF_FaceValueAtPoint3->Fit("pol0","QE"); // quiet, use minos
698	>	MPF_FaceValueAtPoint3->Draw("AP");
699	>	MPF_FaceValueAtPoint3->GetXaxis()->SetTitle("p_{T}^{Z}");
700	>	MPF_FaceValueAtPoint3->GetXaxis()->CenterTitle();
701	>	MPF_FaceValueAtPoint3->GetYaxis()->SetTitle("<MPF>__{data}/<MPF>_{mc} for #alpha<0.3");
702	>	MPF_FaceValueAtPoint3->GetYaxis()->CenterTitle();
703	>	MPF_FaceValueAtPoint3->Draw("AP");
704	>	TF1 *faceval_pol0=(TF1*)MPF_FaceValueAtPoint3->GetFunction("pol0");
705	>	faceval_pol0->SetLineColor(kBlue);
706	>	faceval_pol0->SetLineWidth(1);
707	>	faceval_pol0->Draw("same");
708	>	TF1 *faceval_pol0_complete = new TF1("faceval_pol0_complete","[0]+[1]*x");
709	>	faceval_pol0_complete->SetParameters(faceval_pol0->GetParameters());
710		stringstream faceval_result;
711		faceval_result << "#splitline{C_{abs}= " << std::setprecision(5) << 1/faceval_pol0->GetParameter(0) << " #pm " << std::setprecision(2) << faceval_pol0->GetParError(0)/(faceval_pol0->GetParameter(0)*faceval_pol0->GetParameter(0)) << "}{ (MPF at #alpha<0.3)}";;
712		TText *faceval_final_mark = write_title(faceval_result.str());
#	Line 501 | Line 721 | void new_data_mc_agreement_2d() {
721
722		filename2.str("");
723		filename2 << "Extrapolation/FINAL_RESULT_MPF_FaceValp3";
724	<	CompleteSave(can,filename2.str());
724	>	CompleteSave(can,ContainerName+"/"+filename2.str());
725	>
726	>	can->cd();
727	>	RABS_FaceValueAtPoint3->SetMarkerStyle(21);
728	>	RABS_FaceValueAtPoint3->SetMarkerColor(kRed);
729	>	RABS_FaceValueAtPoint3->Fit("pol0","QE"); // quiet, use minos
730	>	RABS_FaceValueAtPoint3->Draw("AP");
731	>	RABS_FaceValueAtPoint3->GetXaxis()->SetTitle("p_{T}^{Z}");
732	>	RABS_FaceValueAtPoint3->GetXaxis()->CenterTitle();
733	>	RABS_FaceValueAtPoint3->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc} for #alpha<0.3");
734	>	RABS_FaceValueAtPoint3->GetYaxis()->CenterTitle();
735	>	RABS_FaceValueAtPoint3->Draw("AP");
736	>	TF1 *faceval_pol0RABS=(TF1*)RABS_FaceValueAtPoint3->GetFunction("pol0");
737	>	TF1 *faceval_pol0RABS_complete = new TF1("faceval_pol0_complete","[0]+[1]*x");
738	>	faceval_pol0RABS_complete->SetParameters(faceval_pol0RABS->GetParameters());
739	>	faceval_pol0RABS->SetLineWidth(1);
740	>	faceval_pol0RABS->SetLineColor(kRed);
741	>	faceval_pol0RABS_complete->Draw("same");
742	>	stringstream RABSfaceval_result;
743	>	RABSfaceval_result << "#splitline{C_{abs}= " << std::setprecision(5) << 1/faceval_pol0RABS->GetParameter(0) << " #pm " << std::setprecision(2) << faceval_pol0RABS->GetParError(0)/(faceval_pol0RABS->GetParameter(0)*faceval_pol0RABS->GetParameter(0)) << "}{ (R_{abs} at #alpha<0.3)}";;
744	>	TText *RABSfaceval_final_mark = write_title(RABSfaceval_result.str());
745	>	RABSfaceval_final_mark->SetX(0.75);
746	>	RABSfaceval_final_mark->SetY(0.75);
747	>	RABSfaceval_final_mark->SetTextSize(0.03);
748	>	RABSfaceval_final_mark->Draw();
749	>	DrawPrelim();
750	>
751	>	float RABSFaceValResult=1/faceval_pol0RABS->GetParameter(0);
752	>	float RABSFaceValResultError=faceval_pol0RABS->GetParError(0)/(faceval_pol0RABS->GetParameter(0)*faceval_pol0RABS->GetParameter(0));
753	>
754	>	filename2.str("");
755	>	filename2 << "Extrapolation/FINAL_RESULT_RABS_FaceValp3";
756	>	CompleteSave(can,ContainerName+"/"+filename2.str());
757
758		cout << "FINAL RESULTS: " << endl;
759		cout << "MPF : " << MPFResult  << " +/- " << MPFResultError  << endl;
760		cout << "Rabs: " << RabsResult << " +/- " << RabsResultError << endl;
761	<	cout << "Face value at 0.3 :" << FaceValResult << " +/- " << FaceValResultError << endl;
761	>	cout << "Face value at 0.3 :" << FaceValResult << " +/- " << FaceValResultError << "  (MPF) " << endl;
762	>	cout << "Face value at 0.3 :" << RABSFaceValResult << " +/- " << RABSFaceValResultError << "  (RABS) " << endl;
763
764		delete can;
765
766	<
766	>	results.MPF_Result_FaceValue=Value(FaceValResult,FaceValResultError);
767	>	results.Rabs_Result_FaceValue=Value(RABSFaceValResult,RABSFaceValResultError);
768	>	results.MPF_Result_Extrapolated=Value(MPFResult,MPFResultError);
769	>	results.Rabs_Result_Extrapolated=Value(RabsResult,RabsResultError);
770	>	cout << results << endl;
771	>	return results;
772		}
773
774
775	<	void DoPUStudy(string identifier) {
775	>	void TEST_DoPUStudy(string identifier) {
776
777		float numVtxcuts[7]={0,10,15,20};
778
#	Line 564 | Line 822 | void DoPUStudy(string identifier) {
822
823		}
824
825	<	void ScenarioComparison() {
826	<
825	>	void TEST_ScenarioComparison() {
826	>	//very dumb way to do this.
827		TGraphAsymmErrors *gr[5];
828		TF1 *fit[5];
829		bool dofit=false;
#	Line 598 | Line 856 | void ScenarioComparison() {
856		double dyl[4] = {0.011,0.009,0.013,0.020};
857		double dyh[4] = {0.011,0.009,0.013,0.020};
858		gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
859	<	name="30/10";
859	>	name="30/15";
860		color="#2464bc";
861		}
862		if(i==2) {
#	Line 654 | Line 912 | void ScenarioComparison() {
912
913		}
914
915	<	void ScenarioComparisonInclusive() {
916	<
915	>	void TEST_ScenarioComparisonInclusive() {
916	>	//dumbest way ever to do this. but ok we only need to do it once.
917		TGraphAsymmErrors *gr[5];
918		TF1 *fit[5];
919		bool dofit=true;
#	Line 688 | Line 946 | void ScenarioComparisonInclusive() {
946		double dyl[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
947		double dyh[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
948		gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
949	<	name="30/10";
949	>	name="30/15";
950		color="#2464bc";
951		}
952		if(i==2) {
#	Line 744 | Line 1002 | void ScenarioComparisonInclusive() {
1002
1003		}
1004
1005	<	void compare_selection(string identifier) {
1005	>	void TEST_compare_selection(string identifier) {
1006		bool recognized_scenario=false;
1007
1008		cout << "Running with identifier " << identifier << endl;
#	Line 874 | Line 1132 | void compare_selection(string identifier
1132
1133		}
1134
1135	<	void do_basic_ZB_analysis() {
1135	>	void GetNumberEventsInsideOutsideAlphaWindow(TCut specialcut, string title) {
1136	>
1137	>	//  TCut cut =  ZplusBsel&&LeadingB&&PhiZcut&&specialcut&&TCut("Zb3010_alpha<0.3");
1138	>	TCut cut = specialcut;
1139	>	TCut in  =  EtaB;
1140	>	TCut out =  TCut("abs(Zb3010_pfJetGoodEta[0])>1.3");
1141	>
1142	>	TH1F *data_IN  = allsamples.Draw("data_IN",    "mll",1,0,150, "nothing [GeV]", "events", cut&&in, data,luminosity);
1143	>	TH1F *data_OUT = allsamples.Draw("data_OUT",   "mll",1,0,150, "nothing [GeV]", "events", cut&&out,data,luminosity);
1144	>	TH1F *mc_IN    = allsamples.Draw("mc_IN",      "mll",1,0,150, "nothing [GeV]", "events", cut&&in, mc  ,luminosity);
1145	>	TH1F *mc_OUT   = allsamples.Draw("mc_OUT",     "mll",1,0,150, "nothing [GeV]", "events", cut&&out,mc  ,luminosity);
1146	>
1147	>	dout << title << " : " << endl;
1148	>	dout << "  Data: In " << data_IN->Integral() << " vs out " << data_OUT->Integral() << endl;
1149	>	dout << "  MC  : In " << mc_IN->Integral()   << " vs out " << mc_OUT->Integral()   << endl;
1150	>	dout << "       ratio in  : " << data_IN->Integral()/mc_IN->Integral() << " +/- " << (data_IN->Integral()/mc_IN->Integral()) * sqrt(1.0/mc_IN->Integral()+ 1.0/data_IN->Integral()) << endl;
1151	>	dout << "       ratio out : " << data_OUT->Integral()/mc_OUT->Integral() << " +/- " << (data_OUT->Integral()/mc_OUT->Integral()) * sqrt(1.0/mc_OUT->Integral()+ 1.0/data_OUT->Integral()) << endl;
1152	>
1153	>	delete data_IN;
1154	>	delete data_OUT;
1155	>	delete mc_IN;
1156	>	delete mc_OUT;
1157	>	}
1158	>
1159	>	void TEST_GetNumberEventsInsideOutsideAlphaWindow() {
1160	>	TCanvas *ca = new TCanvas("ca","ca");
1161	>	//   GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==0"), "ee");
1162	>	//   GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==1"), "mm");
1163	>	//   GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2"), "ee/mm");
1164	>
1165	>	cout << "BASE SELECTION" << endl;
1166	>	GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==0"), "ee");
1167	>	GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==1"), "mm");
1168	>	GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2"), "ee/mm");
1169	>
1170	>	cout << "BASE SELECTION: Z+b" << endl;
1171	>	GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
1172	>	GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
1173	>	GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2"), "ee/mm");
1174	>
1175	>	cout << "Leading B" << endl;
1176	>	GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
1177	>	GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
1178	>	GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2"), "ee/mm");
1179	>
1180	>	cout << "PhiZcut" << endl;
1181	>	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
1182	>	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
1183	>	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2"), "ee/mm");
1184	>
1185	>	cout << "alpha cut" << endl;
1186	>	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0")&&TCut("Zb3010_alpha<0.3"), "ee");
1187	>	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1")&&TCut("Zb3010_alpha<0.3"), "mm");
1188	>	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2")&&TCut("Zb3010_alpha<0.3"), "ee/mm");
1189	>
1190	>	delete ca;
1191	>
1192	>
1193	>	}
1194	>
1195	>
1196	>
1197	>
1198	>	void do_basic_ZB_analysis(bool DoCompleteAnalysis) {
1199	>
1200	>	//https://twiki.cern.ch/twiki/bin/viewauth/CMS/BtagPOG#2011_Data_and_MC
1201	>	//https://twiki.cern.ch/twiki/pub/CMS/BtagPOG/SFb-ttbar_payload.txt
1202
1203	+	cout << "The lepton requirement is " << (const char*) leptoncut << endl;
1204	+	bool doquick=true;
1205
1206		TCanvas *zbcanvas = new TCanvas("zbcanvas","zbcanvas");
1207
#	Line 893 | Line 1219 | void do_basic_ZB_analysis() {
1219		//     compare_selection("Zb3010");
1220		//     compare_selection("Zb40");
1221
1222	<
1223	<	//  print_all_b_yields();
1224	<	//  draw_mpf_vars();
1225	<	//  DrawEvilCutFlow();
1226	<
1227	<	//  draw_Zb_kin_vars();
1228	<
1229	<	new_data_mc_agreement_2d();
1222	>	if(!doquick) {
1223	>	print_all_b_yields();
1224	>	draw_mpf_vars();
1225	>	//  DrawEvilCutFlow();
1226	>
1227	>	draw_Zb_kin_vars();
1228	>
1229	>	}
1230	>
1231	>	//  GetNumberEventsInsideOutsideAlphaWindow();
1232	>
1233	>	dout <<"THIS IS WHERE IT GETS REALLY INTERESTING!!!!" << endl;
1234	>	if(DoCompleteAnalysis) {
1235	>	/* need to run the following scenarios:
1236	>	* - normal
1237	>	* - inclusive
1238	>	* - medium WP
1239	>	* - alpha up
1240	>	* - alpha down
1241	>	*/
1242	>
1243	>	ZbResultContainer inclusive = new_data_mc_agreement_2d("Zb3010_alpha","inclusive",TCut("Zb3010_bTagProbCSVBP[0]>-10000"),TCut("1.0"));
1244	>	ZbResultContainer Reference = new_data_mc_agreement_2d("Zb3010_alpha","reference",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"));
1245	>	ZbResultContainer effmistagUp = new_data_mc_agreement_2d("Zb3010_alpha","effmistagUp",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtTUp"));
1246	>	ZbResultContainer effmistagDown = new_data_mc_agreement_2d("Zb3010_alpha","effmistagDown",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtTDown"));
1247	>	ZbResultContainer medium = new_data_mc_agreement_2d("Zb3010_alpha","medium",TCut("Zb3010_bTagProbCSVBP[0]>0.679"),TCut("ZbCHS3010_BTagWgtM"));
1248	>	ZbResultContainer loose = new_data_mc_agreement_2d("Zb3010_alpha","loose",TCut("Zb3010_bTagProbCSVBP[0]>0.244"),TCut("ZbCHS3010_BTagWgtL"));
1249	>	ZbResultContainer AlphaUp = new_data_mc_agreement_2d("Zb3010_alphaUp","AlphaUp",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"));
1250	>	ZbResultContainer AlphaDown = new_data_mc_agreement_2d("Zb3010_alphaDown","AlphaDown",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"));
1251	>
1252	>	//and now let's compute some systematics!
1253	>	//1 alpha variation
1254	>	float SysMPF_Alpha_down = Reference.MPF_Result_FaceValue.getValue()-AlphaDown.MPF_Result_FaceValue.getValue();
1255	>	float SysMPF_Alpha_up = Reference.MPF_Result_FaceValue.getValue()-AlphaUp.MPF_Result_FaceValue.getValue();
1256	>
1257	>	//2 btagging efficiency/mistag correction
1258	>	float SysMPF_EFF_down = Reference.MPF_Result_FaceValue.getValue()-effmistagDown.MPF_Result_FaceValue.getValue();
1259	>	float SysMPF_EFF_up = Reference.MPF_Result_FaceValue.getValue()-effmistagUp.MPF_Result_FaceValue.getValue();
1260	>
1261	>	//3 purity
1262	>	zbcanvas->cd();
1263	>	gStyle->SetOptFit(0);
1264	>	write_info(__FUNCTION__,"Don't know the exact purity at the moment, still needs to be determined!");
1265	>	TH1F *purityhisto = new TH1F("purityhisto","purityhisto",1,0,1);
1266	>	purityhisto->GetXaxis()->SetTitle("Purity");
1267	>	purityhisto->GetYaxis()->SetTitle("C_{abs}");
1268	>	purityhisto->GetXaxis()->CenterTitle();
1269	>	purityhisto->GetYaxis()->CenterTitle();
1270	>
1271	>	TGraphErrors *gSysMPF_Purity = new TGraphErrors();
1272	>	gSysMPF_Purity->SetPoint(0,0.05,inclusive.MPF_Result_FaceValue.getValue());
1273	>	gSysMPF_Purity->SetPointError(0,0.0,inclusive.MPF_Result_FaceValue.getError());
1274	>
1275	>	gSysMPF_Purity->SetPoint(1,0.20,loose.MPF_Result_FaceValue.getValue());//x value?
1276	>	gSysMPF_Purity->SetPointError(1,0.0,loose.MPF_Result_FaceValue.getError());
1277	>	gSysMPF_Purity->SetPoint(2,0.54,medium.MPF_Result_FaceValue.getValue());//x value?
1278	>	gSysMPF_Purity->SetPointError(2,0.0,medium.MPF_Result_FaceValue.getError());
1279	>	gSysMPF_Purity->SetPoint(3,0.82,Reference.MPF_Result_FaceValue.getValue());//x value?
1280	>	gSysMPF_Purity->SetPointError(3,0.0,Reference.MPF_Result_FaceValue.getError());
1281	>
1282	>
1283	>	float min,max;
1284	>	FindMinMax(gSysMPF_Purity,min,max);
1285	>	purityhisto->SetStats(0);
1286	>	purityhisto->GetYaxis()->SetRangeUser(min,max);
1287	>
1288	>	gSysMPF_Purity->Fit("pol1");
1289	>	TPolyLine *purity_fit_uncert = GetFitUncertaintyShape(gSysMPF_Purity, "pol1", min, max,0.0,1.0);
1290	>	TF1* pufit = (TF1*)gSysMPF_Purity->GetFunction("pol1");
1291	>
1292	>	TF1* pufit_c = new TF1("pufit_c","[0]+[1]*x");
1293	>	pufit_c->SetParameters(pufit->GetParameters());
1294	>	pufit_c->SetLineColor(TColor::GetColor("#0101DF"));
1295	>	purity_fit_uncert->SetFillColor(TColor::GetColor("#5353FC"));
1296	>	pufit->SetLineColor(TColor::GetColor("#5353FC"));
1297	>
1298	>	purityhisto->Draw();
1299	>	purity_fit_uncert->Draw("f");
1300	>	gSysMPF_Purity->Draw("P*");
1301	>	pufit_c->Draw("same");
1302	>	DrawPrelim();
1303	>
1304	>	float ExtrapolatedResult=pufit->GetParameter(0)+1.0*pufit->GetParameter(1);
1305	>	float dExtrapolatedResult=sqrt(pufit->GetParError(0)*pufit->GetParError(0)+1.0*1.0*pufit->GetParError(1)*pufit->GetParError(1));
1306	>
1307	>	stringstream summary;
1308	>	summary << "#splitline{Reference: " << std::setprecision(4) << Reference.MPF_Result_FaceValue.getValue() << "+/-" << std::setprecision(4) << Reference.MPF_Result_FaceValue.getError() << "}{";
1309	>	summary << "Extrapolation: " << std::setprecision(4) << ExtrapolatedResult << " +/- " << std::setprecision(4) << dExtrapolatedResult << "}";
1310	>
1311	>	dout << "Have found an extrapolated result at 1.0 of " << ExtrapolatedResult << " ;  will use the difference between this result and the face value as ";
1312	>	dout << "purity-related systematic uncertainty" << endl;
1313	>
1314	>	TText *infobox = write_title(summary.str());
1315	>	infobox->SetX(0.75);
1316	>	infobox->SetTextSize(0.03);
1317	>	infobox->SetY(0.75);
1318	>	infobox->Draw();
1319	>
1320	>	CompleteSave(zbcanvas,"Systematics/Purity");
1321	>
1322	>
1323	>
1324	>	float SysMPF_Purity = abs(ExtrapolatedResult-Reference.MPF_Result_FaceValue.getValue());
1325	>
1326	>	float sys_alpha  = abs(SysMPF_Alpha_down)>abs(SysMPF_Alpha_up)?abs(SysMPF_Alpha_down):abs(SysMPF_Alpha_up);
1327	>	float sys_eff    = abs(SysMPF_EFF_down)>abs(SysMPF_EFF_up)?abs(SysMPF_EFF_down):abs(SysMPF_EFF_up);
1328	>	float sys_purity = abs(SysMPF_Purity);
1329	>	float sys        = sqrt(sys_alpha*sys_alpha+sys_eff*sys_eff+sys_purity*sys_purity);
1330	>
1331	>	dout << "MPF METHOD:" << endl;
1332	>	dout << "   Systematics :    down    up     (selected)" << endl;
1333	>	dout << "      Alpha variation      : " << SysMPF_Alpha_down << "   ,  " << SysMPF_Alpha_up << "   ( " << sys_alpha << ") " << endl;
1334	>	dout << "      eff/mistag variation : " << SysMPF_EFF_down << "   ,  " << SysMPF_EFF_up << "   ( " << sys_eff << ") " << endl;
1335	>	dout << "      purity               : " << SysMPF_Purity << "   ( " << sys_purity << " )" << endl;
1336	>	dout << "      total                : " << sys << endl;
1337	>	dout << endl << endl;
1338	>	dout << "   FINAL RESULTS : " << endl;
1339	>	dout << "     C_{abs}^{b}   = " << Reference.MPF_Result_FaceValue << " (stat) +/- " << sys << " (sys) " << endl;
1340	>	dout << "     C_{abs}^{inc} = " << inclusive.MPF_Result_FaceValue << " (stat) " << endl;
1341	>
1342	>	float sysfinal = Reference.MPF_Result_FaceValue.getValue()/inclusive.MPF_Result_FaceValue.getValue();
1343	>	sysfinal *= sqrt(sys*sys/(Reference.MPF_Result_FaceValue.getValue()*Reference.MPF_Result_FaceValue.getValue()));
1344	>	//yes, this could be simplified but it would look like a bug.
1345	>
1346	>	dout << "     C_{corr}      = " << Reference.MPF_Result_FaceValue/inclusive.MPF_Result_FaceValue << " (stat) +/- " << sysfinal << " (sys) " << endl;
1347	>
1348	>	LeadingB=TCut("Zb3010_bTagProbCSVBP[0]>0.898");
1349	>
1350	>	} else {
1351	>	ZbResultContainer Reference = new_data_mc_agreement_2d("Zb3010_alpha","reference",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgt"));
1352	>	}
1353
1354
1355		delete zbcanvas;

Diff Legend

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