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.4 by buchmann, Fri Aug 10 11:47:53 2012 UTC vs.
Revision 1.21 by buchmann, Thu Jan 17 21:58:22 2013 UTC

#	Line 21 | Line 21 | using namespace std;
21		using namespace PlottingSetup;
22
23		namespace ZbData {
24	<	vector<float> data_over_mc;
25	<
26	<	TCut ZplusBsel("pt1>20&&pt2>20&&mll>60&&mll<120");
27	<	TCut LeadingB("Zb30_bTagProbCSVBP[0]>0.679");
28	<	TCut EtaB("abs(Zb30_pfJetGoodEta[0])<1.3");
29	<	TCut PhiZcut("Zb30_pfJetDphiZ[0]>2.7");
24	>	TCut ZplusBsel("pt1>20&&pt2>20&&mll>60&&mll<120&&id1==id2");
25	>	TCut LeadingB("ZbCHS3010_bTagProbCSVBP[0]>0.898");
26	>	TCut EtaB("abs(ZbCHS3010_pfJetGoodEta[0])<1.3");
27	>	TCut PhiZcut("abs(ZbCHS3010_pfJetDphiZ[0])>2.7");
28
29		}
30
31	+	TCut STDWEIGHT;
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	<	void print_yield(TCut cut) {
98	>
99	>	//*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*  DELETE EVERYTHING BETWEEN THIS LINE /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*//
100	>	void SpecialCutFlow(string identifier);
101	>
102	>	//*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*            AND THIS LINE             /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*//
103	>
104	>
105	>	class CutYields{
106	>	public:
107	>	float data;
108	>	float Zb;
109	>	float Zc;
110	>	float Zl;
111	>	float tt;
112	>	float SiTo;
113	>	float Dibosons;
114	>	float WJets;
115	>	float Other;
116	>
117	>	CutYields();
118	>	};
119	>
120	>	CutYields::CutYields() {
121	>	this->data=0;
122	>	this->Zb=0;
123	>	this->Zc=0;
124	>	this->Zl=0;
125	>	this->tt=0;
126	>	this->SiTo=0;
127	>	this->Dibosons=0;
128	>	this->WJets=0;
129	>	this->Other=0;
130	>	}
131	>
132	>	std::ostream &operator<<(std::ostream &ostr, CutYields y)
133	>	{
134	>	return ostr << y.data << ";" << y.Zb << ";" << y.Zc << ";" << y.Zl << ";" << y.tt << ";" << y.SiTo << ";" << y.Dibosons << ";" << y.WJets << ";" << endl;
135	>	}
136	>
137	>
138	>	CutYields print_yield(TCut cut) {
139	>
140	>	CutYields Yield;
141		TH1F *data  = allsamples.Draw("data",    "mll",1,0,10000000, "m_{ll} [GeV]", "events", cut,data,luminosity);
142	<	cout << "  Yields for " << (const char*) cut << endl;
143	<	cout << "     data: " << data->Integral() << endl;
142	>	dout << "  Yields for " << (const char*) cut << endl;
143	>	//  dout << "     data: " << data->Integral() << endl;
144	>
145	>	Yield.data=data->Integral();
146		THStack mcm = allsamples.DrawStack("mc", "mll",1,0,10000000, "m_{ll} [GeV]", "events", cut,mc,luminosity);
147		int nhists = mcm.GetHists()->GetEntries();
41	–	cout << "Going to loop over " << nhists << " histograms " << endl;
148		TFile *test = new TFile("test.root","RECREATE");
149		mcm.Write();
150		test->Close();
151		for (int i = 0; i < nhists; i++) {
152		TH1* hist = static_cast<TH1*>(mcm.GetHists()->At(i));
153	<	cout << "     " << hist->GetName() << ": " << hist->Integral() << endl;
153	>	string name=hist->GetName();
154	>	if(Contains(name,"t_bar_t")) Yield.tt+=hist->Integral();
155	>	if(Contains(name,"W_Jets")) Yield.WJets+=hist->Integral();
156	>	if(Contains(name,"Single_top")) Yield.SiTo+=hist->Integral();
157	>	if(Contains(name,"Dibosons")) Yield.Dibosons+=hist->Integral();
158	>	if(Contains(name,"Z_l_")) Yield.Zl+=hist->Integral();
159	>	if(Contains(name,"Z_c_")) Yield.Zc+=hist->Integral();
160	>	if(Contains(name,"Z_b_")) Yield.Zb+=hist->Integral();
161		}
162	+
163	+	cout << Yield << endl;
164		cout << endl << endl;
165	+
166		delete data;
167	+	CleanLegends();
168	+	DeleteStack(mcm);
169	+
170	+	return Yield;
171		}
172
173		void draw_kin_variable(string variable, TCut cut, int nbins, float min, float max, string xlabel, string saveas, bool dologscale, string speciallabel="") {
54	–
174		TCanvas *ckin = new TCanvas("ckin","kin variable canvas");
175	<	ckin->SetLogy(dologscale);
175	>	if(dologscale) ckin->SetLogy(1);
176		TH1F *datahisto = allsamples.Draw("datahisto",variable,nbins,min,max, xlabel, "events",cut,data,luminosity);
177		datahisto->SetMarkerSize(DataMarkerSize);
178		THStack mcstack = allsamples.DrawStack("mcstack",variable,nbins,min,max, xlabel, "events",cut,mc,luminosity);
179		if(dologscale) datahisto->SetMaximum(3*datahisto->GetMaximum());
180		else datahisto->SetMaximum(1.5*datahisto->GetMaximum());
181	<	TText *speciall;
182	<	if(speciallabel!="") {
183	<	speciall = new TText(0.2,0.8,speciallabel.c_str());
184	<	speciall->Draw();
181	>
182	>	float SumMC=0.0;
183	>	TIter nextH(mcstack.GetHists());
184	>	TH1F* h;
185	>	while ( h = (TH1F*)nextH() ) {
186	>	SumMC += h->Integral();
187		}
188	+	float scale = datahisto->Integral()/SumMC;
189	+
190	+	// Re-scale stack to data integral
191	+	nextH = TIter(mcstack.GetHists());
192	+
193	+	while ( h = (TH1F*)nextH() ) {
194	+	h->Scale(scale);
195	+	}
196	+
197	+	TText *speciall;
198		datahisto->Draw("e1");
199		ckin->Update();
200	<	mcstack.Draw("same");
200	>	mcstack.Draw("histo,same");
201		datahisto->Draw("same,e1");
202		TLegend *kinleg = allsamples.allbglegend();
203	+
204		kinleg->Draw();
205
206		TPad *kinpad = new TPad("kinpad","kinpad",0,0,1,1);
207		kinpad->cd();
208	<	kinpad->SetLogy(dologscale);
208	>	if(dologscale) kinpad->SetLogy(1);
209		datahisto->Draw("e1");
210	<	mcstack.Draw("same");
210	>	mcstack.Draw("histo,same");
211		datahisto->Draw("same,e1");
212		datahisto->Draw("same,axis");
213		kinleg->Draw();
214		DrawPrelim();
215		saveas="kin/"+saveas;
216	<	save_with_ratio(datahisto,mcstack,kinpad->cd(),saveas);
85	<
86	<
216	>	kinpad->cd();
217
218	+	if(speciallabel!="") {
219	+	speciall = write_title(speciallabel);
220	+	speciall->SetX(0.18);
221	+	speciall->SetTextAlign(11);
222	+	speciall->SetTextSize(0.03);
223	+	speciall->SetY(0.85);
224	+	speciall->Draw();
225	+	}
226
227	<	ZbData::data_over_mc.push_back(datahisto->Integral()/CollapseStack(mcstack)->Integral());
227	>	save_with_ratio(datahisto,mcstack,kinpad->cd(),saveas);
228
229		datahisto->Delete();
230	+	delete kinleg;
231		delete ckin;
232	<	}
233	<
95	<	Value getfrom2Dmap(TH2F *map, int ixbin, int iybin) {
96	<	Value sum(0,0);
97	<	for(int iy=1;iy<=iybin;iy++) {
98	<	sum=sum+Value(map->GetBinContent(ixbin,iy),map->GetBinError(ixbin,iy));
99	<	}
100	<	return sum;
101	<	}
102	<
103	<
104	<	void data_mc_agreement_2d() {
105	<
106	<	vector<int> emptyvector;
107	<	float binX[5]={10,50,100,200,1000};
108	<	float binY[5]={0.0,0.1,0.15,0.2,0.3};
109	<	vector<float>binningx;
110	<	vector<float>binningy;
111	<	binningx.assign(binX,binX+5);
112	<	binningy.assign(binY,binY+5);
113	<
114	<	TCut cut(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
115	<
116	<	TCanvas *agrcan = new TCanvas("agrcan","agrcan");
117	<
118	<	cout << "Drawing 2d for data " << endl;
119	<	TH2F *alphazptplane_data = allsamples.Draw("alphazptplane_data","alpha:pt",binningx,binningy,"Z p_{T}","#alpha",cut,data,luminosity,emptyvector,false);
120	<	cout << "Drawing 2d for MC " << endl;
121	<	TH2F *alphazptplane_mc   = allsamples.Draw("alphazptplane_mc",  "alpha:pt",binningx,binningy,"Z p_{T}","#alpha",cut,mc,  luminosity,emptyvector,false);
122	<
123	<
124	<	TFile *f = new TFile("test.root","RECREATE");
125	<	alphazptplane_data->Write();
126	<	f->Close();
127	<
128	<	for(int ix=1;ix<binningx.size();ix++) {
129	<	cout << "pt in " << binningx[ix-1] << " , " << binningx[ix] << endl;
130	<	TGraphErrors gr;
131	<	for(int iy=1;iy<binningy.size();iy++) {
132	<	Value mc   = Value((getfrom2Dmap(alphazptplane_mc,ix,iy)).getValue(),0.0);// not taking MC statistics into account
133	<	Value data = getfrom2Dmap(alphazptplane_data,ix,iy);
134	<	cout << " For an alpha cut at " << binningy[iy] << " we get : mc (" << mc << ") , data (" << data << ") --> data/mc = " << data/mc << endl;
135	<	gr.SetPoint(iy-1,binningy[iy],(data/mc).getValue());
136	<	gr.SetPointError(iy-1,0,(data/mc).getError());
137	<	}
138	<
139	<	gr.Draw("AP*");
140	<	gr.GetYaxis()->SetTitle("data / mc");
141	<	gr.GetYaxis()->CenterTitle();
142	<	gr.GetXaxis()->SetTitle("#alpha");
143	<	gr.GetXaxis()->CenterTitle();
144	<	gr.GetXaxis()->SetRangeUser(0,0.4);
145	<	gr.GetYaxis()->SetRangeUser(0.5,1.5);
146	<	gr.Draw("AP*");
147	<	DrawPrelim();
148	<	stringstream ptbin;
149	<	ptbin << binningx[ix-1] << " < p_{T}^{Z} < " << binningx[ix] << " GeV";
150	<	TText *ptt = write_title(ptbin.str());
151	<	ptt->SetX(0.7);
152	<	ptt->SetY(0.7);
153	<	ptt->SetTextSize(0.03);
154	<	ptt->Draw();
155	<	TLine *l = new TLine(0.0,1.0,0.4,1.0);
156	<	l->SetLineStyle(2);
157	<	l->SetLineColor(kBlue);
158	<	l->Draw();
159	<
160	<	stringstream saveas;
161	<	saveas << "Zplusb_data_over_mc___" << binningx[ix-1] << "to" << binningx[ix];
162	<	DrawPrelim();
163	<	CompleteSave(agrcan,saveas.str());
164	<	//    agrcan->SaveAs("tester.png");
165	<	}
166	<
167	<	delete agrcan;
168	<
232	>	CleanLegends();
233	>	DeleteStack(mcstack);
234		}
235
236		void print_all_b_yields() {
237	+	vector<CutYields> yields;
238		cout << "Basic selection with a b jet" << endl;
239	<	print_yield(ZplusBsel&&"Zb30_pfJetGoodNumBtag>0");
239	>	yields.push_back(print_yield(ZplusBsel&&"ZbCHS3010_pfJetGoodNumBtag>0"));
240		cout << "Leading jet is a b-jet " << endl;
241	<	print_yield(ZplusBsel&&LeadingB);
241	>	yields.push_back(print_yield(ZplusBsel&&LeadingB));
242		cout << "|#eta_{b}|<1.3 " << endl;
243	<	print_yield(ZplusBsel&&LeadingB&&EtaB);
243	>	yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB));
244		cout << "|#delta#phi(b<Z)|>2.7" << endl;
245	<	print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
246	<	cout << "10<ptZ<1000 GeV" << endl;
247	<	print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000");
245	>	yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut));
246	>	cout << "30<ptZ<1000 GeV" << endl;
247	>	yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"));
248	>	cout << "#alpha < 0.35" << endl;
249	>	yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.35"));
250		cout << "#alpha < 0.3" << endl;
251	<	print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.3");
252	<	cout << "#alpha < 0.2" << endl;
253	<	print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.2");
251	>	yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.3"));
252	>	/*  cout << "#alpha < 0.2" << endl;
253	>	yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.2"));
254		cout << "#alpha < 0.15" << endl;
255	<	print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.15");
255	>	yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.15"));
256		cout << "#alpha < 0.1" << endl;
257	<	print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.1");
258	<	cout << "#alpha < 0.05" << endl;
259	<	print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.05");
257	>	yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.1"));*/
258	>
259	>	for(int iy=0;iy<yields.size();iy++) {
260	>	cout << yields[iy] << endl;
261	>	}
262		}
263
264	+	void DrawEvilCutFlow() {
265	+	stringstream MegaCut;
266	+
267	+	MegaCut << "   1*(" << (const char*) (ZplusBsel&&TCut("ZbCHS3010_pfJetGoodNumBtag>0")) << ")";
268	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB) << ")";
269	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB) << ")";
270	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut) << ")";
271	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000")) << ")";
272	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.3")) << ")";
273	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.2")) << ")";
274	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.15")) << ")";
275	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.1")) << ")";
276	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.05")) << ")";
277	+	MegaCut << " ";
278	+
279	+	TCanvas *can = new TCanvas("can","can");
280	+	can->SetLogy(1);
281	+	TCut basecut("mll>6||mll<10");
282	+
283	+
284	+
285	+	TLegend *leg = allsamples.allbglegend();
286	+
287	+
288	+	TH1F *data  = allsamples.Draw("data",    MegaCut.str(),11,-0.5,10.5, "", "events", basecut,data,luminosity);
289	+	THStack mcm = allsamples.DrawStack("mc", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,mc,luminosity);
290	+
291	+	float runningsum=0;
292	+	for(int i=data->GetNbinsX();i>0;i--) {
293	+	runningsum+=data->GetBinContent(i);
294	+	data->SetBinContent(i,runningsum);
295	+	}
296	+
297	+	float minimum=data->GetMaximum();
298	+
299	+	TH1F* h;
300	+	TIter nextOF(mcm.GetHists());
301	+	while ( h = (TH1F*)nextOF() ) {
302	+	float runningsum=0;
303	+	minimum=data->GetMaximum();//done deliberately at each step so get the minimum of the last (!) contribution
304	+	for(int i=h->GetNbinsX();i>0;i--) {
305	+	runningsum+=h->GetBinContent(i);
306	+	h->SetBinContent(i,runningsum);
307	+	if(runningsum<minimum&&runningsum>0) minimum=runningsum;
308	+	}
309	+	}
310	+
311	+
312	+
313	+	data->SetMinimum(0.05*minimum);
314	+	can->cd(1)->SetBottomMargin(0.25);
315	+	data->GetXaxis()->SetBinLabel(1,"Pre-selection");
316	+	data->GetXaxis()->SetBinLabel(2,"Contains b jet");
317	+	data->GetXaxis()->SetBinLabel(3,"Leading jet is b-jet");
318	+	data->GetXaxis()->SetBinLabel(4,"|#eta_{b}|<1.3 ");
319	+	data->GetXaxis()->SetBinLabel(5,"|#delta#phi(b,Z)|>2.7");
320	+	data->GetXaxis()->SetBinLabel(6,"30<p_{T}^{Z}<1000 GeV");
321	+	data->GetXaxis()->SetBinLabel(7,"#alpha < 0.3");
322	+	data->GetXaxis()->SetBinLabel(8,"#alpha < 0.2");
323	+	data->GetXaxis()->SetBinLabel(9,"#alpha < 0.15");
324	+	data->GetXaxis()->SetBinLabel(10,"#alpha < 0.1");
325	+	data->GetXaxis()->SetBinLabel(11,"#alpha < 0.05");
326	+	data->GetXaxis()->LabelsOption("v");
327	+
328	+	data->Draw("e1");
329	+	mcm.Draw("histo,same");
330	+	data->Draw("same,e1,axis");
331	+	data->Draw("same,e1");
332	+	//  data->Draw("same,TEXT");
333	+
334	+	leg->Draw();
335	+
336	+	DrawPrelim();
337	+
338	+	CompleteSave(can,"CutFlow");
339	+	CleanLegends();
340	+	DeleteStack(mcm);
341	+
342	+	}
343	+
344		void draw_Zb_kin_vars() {
345	<	draw_kin_variable("Zb30_pfBJetDphiZ",ZplusBsel&&LeadingB,100,0,3.2,"#delta#phi (Z,b lead)","DeltaPhiZBlead",1);
346	<	draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",50,0,200,"Z p_{T}","Zpt",1);
347	<	draw_kin_variable("Zb30_pfJetGoodPt[1]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",50,0,200,"Sub-Leading Jet Pt","Jet2Pt",1);
348	<	draw_kin_variable("Zb30_pfJetGoodPtlBtag[0]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",50,0,200,"Leading Jet Pt (B)","LeadingJetPt",1);
349	<	draw_kin_variable("Zb30_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,2,"#alpha","alpha",1);
345	>
346	>	draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000"),25,0,200,"Z p_{T}","Official/Zpt",1);
347	>	draw_kin_variable("ZbCHS3010_pfJetGoodPt[1]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000"),25,0,200,"Sub-Leading Jet Pt","Official/Jet2Pt",1);
348	>	draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000"),25,0,200,"Leading Jet Pt (B)","Official/LeadingJetPt",1);
349	>	draw_kin_variable("ZbCHS3010_pfJetGoodPt[1]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000"),20,0,2,"p_{T}^{J2} / p_{T}^{Z}","Official/SubLeadingJetPt_Over_ZPt",1);
350	>	draw_kin_variable("ZbCHS3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<50"),   40,0,2,"#alpha","Official/alpha_pt_30_to_50",1);
351	>	draw_kin_variable("ZbCHS3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>50&&pt<75"),  40,0,2,"#alpha","Official/alpha_50_to_75",1);
352	>	draw_kin_variable("ZbCHS3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>75&&pt<125"), 40,0,2,"#alpha","Official/alpha_pt_75_to_125",1);
353	>	draw_kin_variable("ZbCHS3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>125&&pt<1000"),40,0,2,"#alpha","Official/alpha_pt_125_to_1000",1);
354	>	draw_kin_variable("ZbCHS3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000"),40,0,2,"#alpha","Official/alpha_pt_30_to_1000",1);
355	>
356	>	draw_kin_variable("ZbCHS3010_pfBJetDphiZ[0]",ZplusBsel&&LeadingB&&TCut("pt>30&&pt<1000&&pfBJetDphiZ[0]>0"),30,0,3.2,"#delta#phi (Z,b lead)","DeltaPhiZBlead",0);
357		}
358
359		void draw_mpf_vars() {
360	<	ZbData::data_over_mc.clear();
361	<	draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.3",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p3",0,"#alpha<0.3, 10 GeV < p_{T}^{Z} < 1000 GeV");
362	<	draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.2",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p2",0,"#alpha<0.2, 10 GeV < p_{T}^{Z} < 1000 GeV");
363	<	draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.15",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p15",0,"#alpha<0.15, 10 GeV < p_{T}^{Z} < 1000 GeV");
364	<	draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.1",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p1",0,"#alpha<0.1, 10 GeV < p_{T}^{Z} < 1000 GeV");
365	<	draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.05",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p05",0,"#alpha<0.05, 10 GeV < p_{T}^{Z} < 1000 GeV");
366	<
367	<	draw_kin_variable("Zb30_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.3",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p3",0,"#alpha<0.3, 10 GeV < p_{T}^{Z} < 1000 GeV");
368	<	draw_kin_variable("Zb30_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.2",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p2",0,"#alpha<0.2, 10 GeV < p_{T}^{Z} < 1000 GeV");
369	<	draw_kin_variable("Zb30_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.15",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p15",0,"#alpha<0.15, 10 GeV < p_{T}^{Z} < 1000 GeV");
370	<	draw_kin_variable("Zb30_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_alpha<0.1",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p1",0,"#alpha<0.1, 10 GeV < p_{T}^{Z} < 1000 GeV");
214	<	draw_kin_variable("Zb30_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb30_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");
360	>	draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.3"), 20,0.0,2.0,"Z+b MPF","mpf_alpha_smaller_0p3",0,"#alpha<0.3, 30 GeV < p_{T}^{Z} < 1000 GeV");
361	>	draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.2"), 20,0.0,2.0,"Z+b MPF","mpf_alpha_smaller_0p2",0,"#alpha<0.2, 30 GeV < p_{T}^{Z} < 1000 GeV");
362	>	draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.15"),20,0.0,2.0,"Z+b MPF","mpf_alpha_smaller_0p15",0,"#alpha<0.15, 30 GeV < p_{T}^{Z} < 1000 GeV");
363	>	draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.1"), 20,0.0,2.0,"Z+b MPF","mpf_alpha_smaller_0p1",0,"#alpha<0.1, 30 GeV < p_{T}^{Z} < 1000 GeV");
364	>	draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.05"),20,0.0,2.0,"Z+b MPF","mpf_alpha_smaller_0p05",0,"#alpha<0.05, 30 GeV < p_{T}^{Z} < 1000 GeV");
365	>
366	>	draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.3"), 20,0.0,2.0,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p3",    0,"#alpha<0.3, 30 GeV < p_{T}^{Z} < 1000 GeV");
367	>	draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.2"), 20,0.0,2.0,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p2",    0,"#alpha<0.2, 30 GeV < p_{T}^{Z} < 1000 GeV");
368	>	draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.15"),20,0.0,2.0,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p15",   0,"#alpha<0.15, 30 GeV < p_{T}^{Z} < 1000 GeV");
369	>	draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.1"), 20,0.0,2.0,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p1",    0,"#alpha<0.1, 30 GeV < p_{T}^{Z} < 1000 GeV");
370	>	draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.05"),20,0.0,2.0,"p_{T} b-jet / p_{T} Z","ptb_over_ptz__mpf_alpha_smaller_0p05",0,"#alpha<0.05, 30 GeV < p_{T}^{Z} < 1000 GeV");
371		}
372
373	<	Value get_Zb_data_over_mc(string variable, TCut cut) {
373	>	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) {
374		//  write_warning(__FUNCTION__,"Debugging this function, therefore always returning 3.1415+/-1.010101"); return Value(3.1415,1.010101);
375	<	TH1F *hdata  = allsamples.Draw("hdata",variable,1000,0,30, "#alpha [GeV]", "events", cut,data,luminosity);
376	<	TH1F *hmc    = allsamples.Draw("hmc" , variable,1000,0,30, "#alpha [GeV]", "events", cut,mc,  luminosity);
375	>	TCanvas *cn = new TCanvas("cn","cn");
376	>	string varname="MPF";
377	>
378	>	bool DoFit=false;
379	>	bool UseFit=false;
380	>
381	>
382	>	if(Contains(variable,"pfJetGood")) varname="R_{abs}";
383	>	//  TH1F *hdata  = allsamples.Draw("hdata",variable,1200,0,30, varname, "events",  cut,data,luminosity);
384	>	//  TH1F *hmc    = allsamples.Draw("hmc" , variable,1200,0,30, varname, "events", cut, mc,  luminosity);
385	>	TH1F *hdata  = allsamples.Draw("hdata",variable,1200,0.5,1.5, varname, "events",  cut,data,luminosity); // making it more precise
386	>	TH1F *hmc    = allsamples.Draw("hmc" , variable,1200,0.5,1.5, varname, "events", cut, mc,  luminosity); // making it more precise
387	>	hmc->SetLineColor(kBlue);
388
389		float a=hdata->GetMean();
390		float b=hmc->GetMean();
#	Line 225 | Line 392 | Value get_Zb_data_over_mc(string variabl
392		float db=hmc->GetMeanError();
393		float factor = a / b;
394		float error = TMath::Sqrt(  (1/(b*b)) * (da*da) + ((a*a)/(b*b))*db*db);
395	+
396	+	TF1 *gaus_data = new TF1("gaus_data","gaus",0.7,1.3);
397	+	TF1 *gaus_mc = new TF1("gaus_mc","gaus",0.7,1.3);
398	+	if(DoFit) {
399	+	gaus_data->SetParameter(0,hdata->GetMaximum());
400	+	gaus_data->SetParameter(2,hdata->GetMean());
401	+	gaus_data->SetParameter(1,hdata->GetRMS());
402	+	hdata->Fit(gaus_data);
403	+
404	+	gaus_mc->SetParameter(0,hmc->GetMaximum());
405	+	gaus_mc->SetParameter(1,hmc->GetMean());
406	+	gaus_mc->SetParameter(2,hmc->GetRMS());
407	+	hmc->Fit(gaus_mc);
408	+	}
409	+
410	+	cn->cd();
411	+	hdata->GetYaxis()->SetTitle("events / 0.1");
412	+	hdata->SetMarkerColor(kRed);
413	+	if(DoFit) {
414	+	gaus_mc->SetLineColor(kBlue);
415	+	gaus_data->SetLineColor(kRed);
416	+	}
417	+
418	+	hdata->GetXaxis()->SetRangeUser(0,2);
419	+	hdata->GetYaxis()->SetTitleOffset(1.3);
420	+	hdata->Draw("e1");
421	+	hmc->Draw("histo,same");
422	+	hdata->Draw("e1,same");
423	+	if(DoFit) {
424	+	gaus_data->Draw("same");
425	+	gaus_mc->Draw("same");
426	+	}
427	+
428	+	stringstream summary;
429	+	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 << "}{";
430	+	summary << "#splitline{mc:   " << std::setprecision(4) << b << " #pm " << std::setprecision(4) << db << "}{#splitline{ratio: " << factor;
431	+	summary << " #pm " << error;
432	+	if(DoFit) {
433	+	summary << "}{#splitline{}{#splitline{data fit: " << std::setprecision(4) << gaus_data->GetParameter(1) << "+/-";
434	+	summary << std::setprecision(4) << gaus_data->GetParError(1) << "}{#splitline{";
435	+	summary << "mc fit:" << std::setprecision(4) << gaus_mc->GetParameter(1) << "#pm" << std::setprecision(4) << gaus_mc->GetParError(1) << "}{ratio: ";
436	+	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)) ) << "}}}}";
437	+	summary << "}}}}";
438	+	} else {
439	+	summary << "}{}}}}}";
440	+	}
441	+
442	+	if(UseFit) {
443	+	write_error(__FUNCTION__,"Requested using fit result for final result - not implemented yet!");
444	+	assert(0);
445	+	}
446	+
447	+
448	+	TText *infobox = write_title(summary.str());
449	+	infobox->SetX(0.75);
450	+	infobox->SetTextSize(0.03);
451	+	infobox->SetY(0.75);
452	+	infobox->Draw();
453	+
454	+
455	+	DrawPrelim();
456	+	CompleteSave(cn,ContainerName+"/ResponsePlots/"+saveas);
457	+
458	+	res.MPF_Rabs_location.push_back(ZbPointResult(variable,pt1,pt2,a1,a2,Value(a,da),Value(b,db)));
459	+	cout << "Data;" << a << ";"<<da<<";MC;"<<b<<";"<<db<<endl;
460	+	delete cn;
461		delete hdata;
462		delete hmc;
463	+
464		return Value(factor,error);
465		}
466
467	+	int nFakes=1;
468	+
469	+	float PurityInclusive=0.05;
470	+	float PurityLoose=0.20;
471	+	float PurityMedium=0.54;
472	+	float PurityTight=0.82;
473	+
474	+
475	+	void DeterminePurity(TCut basecut, int nptcuts, float ptcuts[], string AlphaVariable, string ContainerName, float FaceValueToBeTakenAt) {
476	+
477	+	bool ComputePurity=false;
478	+
479	+	if(ContainerName=="inclusive") ComputePurity=true;
480	+	if(ContainerName=="loose") ComputePurity=true;
481	+	if(ContainerName=="medium") ComputePurity=true;
482	+	if(ContainerName=="reference") ComputePurity=true;
483	+
484	+	if(!ComputePurity) return;
485	+
486	+	stringstream specialcut;
487	+	specialcut << "((pt>30&& pt< 1000) && (" << AlphaVariable << "<" << FaceValueToBeTakenAt << "))";
488	+
489	+	THStack mcm    = allsamples.DrawStack("mcm" , "mpf",1200,0,30, "MPF", "events", TCut(basecut && specialcut.str().c_str()), mc,  luminosity);
490	+	int nhists = mcm.GetHists()->GetEntries();
491	+
492	+	float total=0;
493	+	float Zb=0;
494
495	+	for (int i = 0; i < nhists; i++) {
496	+	TH1* hist = static_cast<TH1*>(mcm.GetHists()->At(i));
497	+	string name=hist->GetName();
498	+	if(Contains(name,"t_bar_t")) total+=hist->Integral();
499	+	if(Contains(name,"W_Jets")) total+=hist->Integral();
500	+	if(Contains(name,"Single_top")) total+=hist->Integral();
501	+	if(Contains(name,"Dibosons")) total+=hist->Integral();
502	+	if(Contains(name,"Z_l_")) total+=hist->Integral();
503	+	if(Contains(name,"Z_c_")) total+=hist->Integral();
504	+	if(Contains(name,"Z_b_")) {
505	+	total+=hist->Integral();
506	+	Zb+=hist->Integral();
507	+	}
508	+	}
509	+
510	+	float purity=Zb/total;
511	+	dout << "Determined a purity of " << 100*purity << " % for the container " << ContainerName << endl;
512	+
513	+	if(ContainerName=="inclusive") PurityInclusive=purity;
514	+	if(ContainerName=="loose") PurityLoose=purity;
515	+	if(ContainerName=="medium") PurityMedium=purity;
516	+	if(ContainerName=="reference") PurityTight=purity;
517	+
518	+	DeleteStack(mcm);
519	+	}
520	+
521
522	<	void new_data_mc_agreement_2d() {
522	>	ZbResultContainer new_data_mc_agreement_2d(bool gofast, string AlphaVariable, string ContainerName, TCut BTagCut, TCut BTagWgt, float FaceValueToBeTakenAt) {
523	>
524	>
525	>	if(gofast) write_info(__FUNCTION__,"Fast mode activated for "+ContainerName);
526	>	ZbResultContainer results(ContainerName);
527	>
528	>	cutWeight=STDWEIGHT*BTagWgt;
529	>	LeadingB=BTagCut;
530	>
531	>
532		gStyle->SetOptFit(0);
533		TCut basecut(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
534	+	const int nalphacuts=6;
535	+	float alphacuts[nalphacuts] = {0.1,0.15,0.2,0.25,0.3,0.35};
536	+
537
239	–	const int nalphacuts=5;
240	–	float alphacuts[nalphacuts] = {0.1,0.15,0.2,0.25,0.3};
538		const int nptcuts=5;
539	<	float ptcuts[nptcuts]={10,50,100,200,1000};
539	>	float ptcuts[nptcuts]={30,50,75,125,1000};
540	>
541	>	DeterminePurity(basecut, nptcuts,ptcuts, AlphaVariable, ContainerName, FaceValueToBeTakenAt);
542
543		float MPF_Results[nalphacuts][nptcuts];
544		float MPF_Errors[nalphacuts][nptcuts];
#	Line 247 | Line 546 | void new_data_mc_agreement_2d() {
546		float RABS_Errors[nalphacuts][nptcuts];
547
548
549	<	for(int ia=0;ia<nalphacuts;ia++) {
549	>	for(int ia=0;ia<nalphacuts && !gofast;ia++) {
550		for(int ipt=0;ipt<nptcuts-1;ipt++) {
551		stringstream specialcut;
552	<	specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (alpha<" << alphacuts[ia] << "))";
553	<	Value MPF_data_over_mc = get_Zb_data_over_mc("mpf",TCut(basecut && specialcut.str().c_str()));
554	<	Value RABS_data_over_mc = get_Zb_data_over_mc("Zb30_pfJetGoodPt[0]/pt",TCut(basecut && specialcut.str().c_str()));
552	>	float alow,ahigh;
553	>	if(ia>0) {
554	>	specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (" << AlphaVariable << "<" << alphacuts[ia] << " && " << AlphaVariable << ">" << alphacuts[ia-1] << "))";
555	>	alow=alphacuts[ia-1];
556	>	ahigh=alphacuts[ia];
557	>	} else {
558	>	specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (" << AlphaVariable << "<" << alphacuts[ia] << "))";
559	>	alow=0;
560	>	ahigh=alphacuts[ia];
561	>	}
562	>	cout << specialcut.str() << endl;
563	>	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);
564	>	Value RABS_data_over_mc = get_Zb_data_over_mc(ContainerName,results,"ZbCHS3010_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);
565
566		MPF_Results[ia][ipt]=MPF_data_over_mc.getValue();
567		MPF_Errors[ia][ipt]=MPF_data_over_mc.getError();
#	Line 260 | Line 569 | void new_data_mc_agreement_2d() {
569		RABS_Results[ia][ipt]=RABS_data_over_mc.getValue();
570		RABS_Errors[ia][ipt]=RABS_data_over_mc.getError();
571
572	<	cout << "alpha cut at " << alphacuts[ia] << " and pt range " << ptcuts[ipt] << " , " << ptcuts[ipt+1] << " \t : \t" << specialcut.str() << " \t \t --> " << MPF_data_over_mc << " (MPF) " << RABS_data_over_mc << " (RABS)" << endl;
572	>	cout << "alpha cut at " << alphacuts[ia+1] << " and pt range " << ptcuts[ipt] << " , " << ptcuts[ipt+1] << " \t : \t" << specialcut.str() << " \t \t --> " << MPF_data_over_mc << " (MPF) " << RABS_data_over_mc << " (RABS)" << endl;
573		}
574		}
575
#	Line 268 | Line 577 | void new_data_mc_agreement_2d() {
577		float RABS_ExtraPolatedResults[nptcuts];
578
579		TCanvas *can = new TCanvas("can");
580	+	can->cd(1)->SetLeftMargin(0.15);
581
582		TGraphErrors *MPF_FinalGraph = new TGraphErrors();
583		TGraphErrors *RABS_FinalGraph = new TGraphErrors();
584	+	TGraphErrors *MPF_FaceValueAtPoint3 = new TGraphErrors();
585	+	TGraphErrors *RABS_FaceValueAtPoint3 = new TGraphErrors();
586
587		for(int ipt=0;ipt<nptcuts-1;ipt++) {
588
#	Line 279 | Line 591 | void new_data_mc_agreement_2d() {
591		cout << "Going to create histo for " << filename.str() << endl;
592		TGraphErrors *mpf_gr =new TGraphErrors();
593		TGraphErrors *rabs_gr =new TGraphErrors();
282	–	for(int ia=0;ia<nalphacuts;ia++) {
283	–	mpf_gr->SetPoint(ia,alphacuts[ia],MPF_Results[ia][ipt]);
284	–	mpf_gr->SetPointError(ia,0,MPF_Errors[ia][ipt]);
285	–	rabs_gr->SetPoint(ia,alphacuts[ia],RABS_Results[ia][ipt]);
286	–	rabs_gr->SetPointError(ia,0,RABS_Errors[ia][ipt]);
287	–	}
288	–
289	–	mpf_gr->SetMarkerStyle(21);
290	–	mpf_gr->SetMarkerColor(kBlue);
291	–	mpf_gr->Draw("AP");
292	–	mpf_gr->Fit("pol1","QE"); // quiet, use minos
293	–	mpf_gr->GetXaxis()->SetTitle("#alpha");
294	–	mpf_gr->GetXaxis()->CenterTitle();
295	–	mpf_gr->GetYaxis()->SetTitle("data/mc");
296	–	mpf_gr->GetYaxis()->CenterTitle();
297	–	TF1 *mpf_pol=(TF1*)mpf_gr->GetFunction("pol1");
298	–	float mpf_a=mpf_pol->GetParameter(0);
299	–	float mpf_b=mpf_pol->GetParameter(1);
300	–	MPF_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),mpf_a);
301	–	MPF_FinalGraph->SetPointError(ipt,0,mpf_pol->GetParError(0));
302	–
303	–	MPF_ExtraPolatedResults[ipt]=mpf_a;
304	–
305	–	stringstream mpf_info;
306	–	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() << "}";
307	–
308	–	TText *mpf_mark = write_title(mpf_info.str());
309	–	mpf_mark->SetX(0.75);
310	–	mpf_mark->SetTextSize(0.03);
311	–	mpf_mark->SetY(0.75);
312	–	mpf_mark->Draw();
594
595	<	string filenamebkp=filename.str();
596	<	filename << "__MPF";
316	<	DrawPrelim();
317	<	CompleteSave(can,filename.str());
318	<	cout << "MPF : " << mpf_a << " + " << mpf_b << " * alpha " << endl;
595	>	int rabs_counter=0;
596	>	int mpf_counter=0;
597
598	<	filename.str("");
598	>	for(int ia=0;ia<nalphacuts && !gofast;ia++) {
599	>	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;
600	>	if(MPF_Results[ia][ipt]==MPF_Results[ia][ipt] && MPF_Errors[ia][ipt]==MPF_Errors[ia][ipt]) { // remove nan's
601	>	mpf_gr->SetPoint(mpf_counter,alphacuts[ia],MPF_Results[ia][ipt]);
602	>	mpf_gr->SetPointError(mpf_counter,0,MPF_Errors[ia][ipt]);
603	>	mpf_counter++;
604	>	} else {
605	>	dout << "Detected BS for MPF method (bin ignored) : " << endl;
606	>	dout << "       alpha: " << alphacuts[ia] << " , pt range: " << ptcuts[ipt] << " < pt < " << ptcuts[ipt+1] << endl;
607	>	cout << "       value (MPF) : " << MPF_Results[ia][ipt] << " +/- " << MPF_Errors[ia][ipt] << endl;
608	>	}
609	>	if(RABS_Results[ia][ipt]==RABS_Results[ia][ipt] && RABS_Errors[ia][ipt]==RABS_Errors[ia][ipt]) {
610	>	rabs_gr->SetPoint(rabs_counter,alphacuts[ia],RABS_Results[ia][ipt]);
611	>	rabs_gr->SetPointError(rabs_counter,0,RABS_Errors[ia][ipt]);
612	>	rabs_counter++;
613	>	} else {
614	>	dout << "Detected BS for RABS method (bin ignored) : " << endl;
615	>	dout << "       alpha: " << alphacuts[ia] << " , pt range: " << ptcuts[ipt] << " < pt < " << ptcuts[ipt+1] << endl;
616	>	cout << "       value (RABS) : " << RABS_Results[ia][ipt] << " +/- " << RABS_Errors[ia][ipt] << endl;
617	>	}
618	>	}
619	>
620	>	stringstream specialcut;
621	>	specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (" << AlphaVariable << "<" << FaceValueToBeTakenAt << "))";
622	>	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_"+any2string(FaceValueToBeTakenAt),ptcuts[ipt],ptcuts[ipt+1],0,FaceValueToBeTakenAt);
623	>	Value RABS_data_over_mc = get_Zb_data_over_mc(ContainerName,results,"ZbCHS3010_pfJetGoodPt[0]/pt","",TCut(basecut && specialcut.str().c_str()),"Rabs___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__INCLUSIVEalpha_"+any2string(FaceValueToBeTakenAt),ptcuts[ipt],ptcuts[ipt+1],0,FaceValueToBeTakenAt);
624	>	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;
625	>	MPF_FaceValueAtPoint3->SetPoint(ipt,(ptcuts[ipt]+ptcuts[ipt+1])/2,MPF_data_over_mc.getValue());
626	>	MPF_FaceValueAtPoint3->SetPointError(ipt,0,MPF_data_over_mc.getError());
627	>	RABS_FaceValueAtPoint3->SetPoint(ipt,(ptcuts[ipt]+ptcuts[ipt+1])/2,RABS_data_over_mc.getValue());
628	>	RABS_FaceValueAtPoint3->SetPointError(ipt,0,RABS_data_over_mc.getError());
629
322	–	rabs_gr->SetMarkerStyle(21);
323	–	rabs_gr->SetMarkerColor(kBlue);
324	–	rabs_gr->Fit("pol1","QE"); // quiet, use minos
325	–	rabs_gr->Draw("AP");
326	–	rabs_gr->GetXaxis()->SetTitle("#alpha");
327	–	rabs_gr->GetXaxis()->CenterTitle();
328	–	rabs_gr->GetYaxis()->SetTitle("data/mc");
329	–	rabs_gr->GetYaxis()->CenterTitle();
630
331	–	TF1 *rabs_pol=(TF1*)rabs_gr->GetFunction("pol1");
332	–	float rabs_a=rabs_pol->GetParameter(0);
333	–	float rabs_b=rabs_pol->GetParameter(1);
334	–	RABS_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),rabs_a);
335	–	RABS_FinalGraph->SetPointError(ipt,0,rabs_pol->GetParError(0));
631
632	<	cout << "!+!+!+!+!+!+!+!+!+ Just added a point to the final plot : " << rabs_a << " +/- " << rabs_pol->GetParError(0) << endl;
632	>	if(!gofast) {
633	>	can->cd();
634	>	mpf_gr->SetMarkerStyle(21);
635	>	mpf_gr->SetMarkerColor(kBlue);
636	>	mpf_gr->Draw("AP");
637	>	mpf_gr->Fit("pol1","");
638	>	mpf_gr->GetXaxis()->SetTitle("#alpha");
639	>	mpf_gr->GetXaxis()->CenterTitle();
640	>	mpf_gr->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
641	>	mpf_gr->GetYaxis()->CenterTitle();
642	>	TF1 *mpf_pol=(TF1*)mpf_gr->GetFunction("pol1");
643	>	mpf_pol->SetLineWidth(0);
644	>	TF1 *mpf_pol_complete = new TF1("mpf_pol_complete","[0]+[1]*x");
645	>	mpf_pol_complete->SetParameters(mpf_pol->GetParameters());
646	>	mpf_pol_complete->SetLineWidth(1);
647	>	mpf_pol_complete->SetLineColor(kBlue);
648	>
649	>	float min,max;
650	>	FindMinMax(mpf_gr,min,max);
651	>	TH1F *histo = new TH1F("histo","histo",1,0,0.4);
652	>	histo->GetXaxis()->SetTitle("#alpha");
653	>	histo->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
654	>	histo->GetXaxis()->CenterTitle();
655	>	histo->GetYaxis()->CenterTitle();
656	>	histo->GetYaxis()->SetTitleOffset(1.3);
657	>	histo->SetStats(0);
658	>	histo->GetXaxis()->SetRangeUser(0,0.4);
659	>	histo->GetYaxis()->SetRangeUser(min,max);
660	>	mpf_gr->GetYaxis()->SetRangeUser(min,max);
661	>
662	>	TPolyLine *mpf_fit_uncert = GetFitUncertaintyShape(mpf_gr, "pol1", min, max,0.0,0.4);
663	>	mpf_pol->SetLineColor(TColor::GetColor("#0101DF"));
664	>	mpf_fit_uncert->SetFillColor(TColor::GetColor("#5353FC"));
665	>	histo->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
666	>	histo->Draw();
667	>	mpf_fit_uncert->Draw("F");
668	>	mpf_fit_uncert->Draw("");
669	>	mpf_gr->Draw("P");
670	>	mpf_pol_complete->Draw("same");
671	>
672	>	float mpf_a=mpf_pol->GetParameter(0);
673	>	float mpf_b=mpf_pol->GetParameter(1);
674	>	MPF_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),mpf_a);
675	>	MPF_FinalGraph->SetPointError(ipt,0,mpf_pol->GetParError(0));
676	>
677	>	MPF_ExtraPolatedResults[ipt]=mpf_a;
678	>
679	>	stringstream mpf_info;
680	>	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() << "}";
681	>
682	>	TText *mpf_mark = write_title(mpf_info.str());
683	>	mpf_mark->SetX(0.75);
684	>	mpf_mark->SetTextSize(0.03);
685	>	mpf_mark->SetY(0.75);
686	>	mpf_mark->Draw();
687	>
688	>	string filenamebkp=filename.str();
689	>	filename << "__MPF";
690	>	DrawPrelim();
691	>	CompleteSave(can,ContainerName+"/"+filename.str());
692	>	cout << "MPF : " << mpf_a << " + " << mpf_b << " * alpha " << endl;
693	>
694	>	filename.str("");
695	>
696	>	rabs_gr->SetMarkerStyle(21);
697	>	rabs_gr->SetMarkerColor(kRed);
698	>	rabs_gr->Draw("AP");
699	>	rabs_gr->Fit("pol1","");
700	>	rabs_gr->GetXaxis()->SetTitle("#alpha");
701	>	rabs_gr->GetXaxis()->CenterTitle();
702	>	rabs_gr->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
703	>	rabs_gr->GetYaxis()->CenterTitle();
704	>	TF1 *rabs_pol=(TF1*)rabs_gr->GetFunction("pol1");
705	>	rabs_pol->SetLineWidth(0);
706	>	TF1 *rabs_pol_complete = new TF1("rabs_pol_complete","[0]+[1]*x");
707	>	rabs_pol_complete->SetParameters(rabs_pol->GetParameters());
708	>	rabs_pol_complete->SetLineColor(kRed);
709	>	rabs_pol_complete->SetLineWidth(1);
710	>	FindMinMax(rabs_gr,min,max);
711	>	rabs_gr->GetYaxis()->SetRangeUser(min,max);
712	>	histo->GetYaxis()->SetRangeUser(min,max);
713	>	TPolyLine *rabs_fit_uncert = GetFitUncertaintyShape(rabs_gr, "pol1", min, max,0.0,0.4);
714	>	rabs_pol->SetLineColor(TColor::GetColor("#FA5858"));
715	>	rabs_pol->SetFillColor(TColor::GetColor("#FA5858"));
716	>	rabs_fit_uncert->SetLineColor(TColor::GetColor("#FA5858"));
717	>	rabs_fit_uncert->SetFillColor(TColor::GetColor("#FA5858"));
718	>	histo->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
719	>	histo->Draw();
720	>	rabs_fit_uncert->Draw("F");
721	>	rabs_fit_uncert->Draw("");
722	>	rabs_gr->Draw("P");
723	>	rabs_pol_complete->Draw("same");
724	>
725	>	float rabs_a=rabs_pol->GetParameter(0);
726	>	float rabs_b=rabs_pol->GetParameter(1);
727	>	RABS_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),rabs_a);
728	>	RABS_FinalGraph->SetPointError(ipt,0,rabs_pol->GetParError(0));
729	>
730	>	cout << "!+!+!+!+!+!+!+!+!+ Just added a point to the final plot : " << rabs_a << " +/- " << rabs_pol->GetParError(0) << endl;
731	>
732	>	stringstream rabs_info;
733	>	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() << "}";
734	>	TText *rabs_mark = write_title(rabs_info.str());
735	>	rabs_mark->SetX(0.75);
736	>	rabs_mark->SetTextSize(0.03);
737	>	rabs_mark->SetY(0.75);
738	>	rabs_mark->Draw();
739	>
740	>	RABS_ExtraPolatedResults[ipt]=rabs_a;
741	>	filename << filenamebkp << "__RABS";
742	>	DrawPrelim();
743	>	CompleteSave(can,ContainerName+"/"+filename.str());
744	>	cout << "RABS : " << rabs_a << " + " << rabs_b << " * alpha " << endl;
745	>	}
746	>	}
747	>
748	>	float MPFResult;
749	>	float MPFResultError;
750	>
751	>	float RabsResult;
752	>	float RabsResultError;
753	>
754	>	if(!gofast) {
755	>	MPF_FinalGraph->SetMarkerStyle(21);
756	>	MPF_FinalGraph->SetMarkerColor(kBlue);
757	>	MPF_FinalGraph->Fit("pol0",""); // quiet, use minos
758
759	+	TF1 *mpf_pol0=(TF1*)MPF_FinalGraph->GetFunction("pol0");
760	+	TF1 *mpf_pol0_complete = new TF1("mpf_pol0_complete","[0]+[1]*x");
761	+	mpf_pol0_complete->SetLineWidth(1);
762	+	mpf_pol0_complete->SetLineColor(kBlue);
763	+	mpf_pol0->SetLineColor(kBlue);
764	+	MPF_FinalGraph->Draw("AP");
765	+	MPF_FinalGraph->GetXaxis()->SetTitle("p_{T}^{Z}");
766	+	MPF_FinalGraph->GetXaxis()->CenterTitle();
767	+	MPF_FinalGraph->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
768	+	MPF_FinalGraph->GetYaxis()->CenterTitle();
769	+	MPF_FinalGraph->Draw("AP");
770	+	mpf_pol0_complete->Draw("same");
771
772	<	stringstream rabs_info;
773	<	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() << "}";
774	<	TText *rabs_mark = write_title(rabs_info.str());
775	<	rabs_mark->SetX(0.75);
776	<	rabs_mark->SetTextSize(0.03);
777	<	rabs_mark->SetY(0.75);
778	<	rabs_mark->Draw();
772	>	stringstream mpf_result;
773	>	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)}";
774	>	TText *rmpf_final_mark = write_title(mpf_result.str());
775	>	rmpf_final_mark->SetX(0.75);
776	>	rmpf_final_mark->SetY(0.75);
777	>	rmpf_final_mark->SetTextSize(0.03);
778	>	rmpf_final_mark->Draw();
779	>	DrawPrelim();
780	>	MPFResult=1/mpf_pol0->GetParameter(0);
781	>	MPFResultError=mpf_pol0->GetParError(0)/(mpf_pol0->GetParameter(0)*mpf_pol0->GetParameter(0));
782
783	+	stringstream filename2;
784	+	filename2 << "Extrapolation/FINAL_RESULT_MPF";
785	+	CompleteSave(can,ContainerName+"/"+filename2.str());
786
787	<
788	<	RABS_ExtraPolatedResults[ipt]=rabs_a;
789	<	filename << filenamebkp << "__RABS";
787	>	RABS_FinalGraph->SetMarkerStyle(21);
788	>	RABS_FinalGraph->SetMarkerStyle(21);
789	>	RABS_FinalGraph->SetMarkerColor(kRed);
790	>	RABS_FinalGraph->Fit("pol0","QE"); // quiet, use minos
791	>	RABS_FinalGraph->Draw("AP");
792	>	RABS_FinalGraph->GetXaxis()->SetTitle("p_{T}^{Z}");
793	>	RABS_FinalGraph->GetXaxis()->CenterTitle();
794	>	RABS_FinalGraph->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
795	>	RABS_FinalGraph->GetYaxis()->CenterTitle();
796	>	RABS_FinalGraph->Draw("AP");
797	>	TF1 *rabs_pol0=(TF1*)RABS_FinalGraph->GetFunction("pol0");
798	>	stringstream rabs_result;
799	>	rabs_result << "#splitline{C_{abs}= " << std::setprecision(5) << 1/rabs_pol0->GetParameter(0) << " #pm " << std::setprecision(2) << rabs_pol0->GetParError(0)/(rabs_pol0->GetParameter(0)*rabs_pol0->GetParameter(0)) << "}{ (R_{abs})}";;
800	>	TText *rabs_final_mark = write_title(rabs_result.str());
801	>	rabs_final_mark->SetX(0.75);
802	>	rabs_final_mark->SetY(0.75);
803	>	rabs_final_mark->SetTextSize(0.03);
804	>	rabs_final_mark->Draw();
805		DrawPrelim();
806	<	CompleteSave(can,filename.str());
807	<	cout << "RABS : " << rabs_a << " + " << rabs_b << " * alpha " << endl;
806	>
807	>	RabsResult=1/rabs_pol0->GetParameter(0);
808	>	RabsResultError=rabs_pol0->GetParError(0)/(rabs_pol0->GetParameter(0)*rabs_pol0->GetParameter(0));
809	>
810	>	filename2.str("");
811	>	filename2 << "Extrapolation/FINAL_RESULT_RABS";
812	>	CompleteSave(can,ContainerName+"/"+filename2.str());
813		}
814	<
815	<	MPF_FinalGraph->SetMarkerStyle(21);
816	<	MPF_FinalGraph->SetMarkerColor(kRed);
817	<	MPF_FinalGraph->Fit("pol0","QE"); // quiet, use minos
818	<
819	<	TF1 *mpf_pol0=(TF1*)MPF_FinalGraph->GetFunction("pol0");
820	<	MPF_FinalGraph->Draw("AP");
821	<	MPF_FinalGraph->GetXaxis()->SetTitle("p_{T}^{Z}");
822	<	MPF_FinalGraph->GetXaxis()->CenterTitle();
823	<	MPF_FinalGraph->GetYaxis()->SetTitle("data/mc");
824	<	MPF_FinalGraph->GetYaxis()->CenterTitle();
825	<	MPF_FinalGraph->Draw("AP");
826	<
827	<	stringstream mpf_result;
828	<	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)}";
829	<	TText *rmpf_final_mark = write_title(mpf_result.str());
830	<	rmpf_final_mark->SetX(0.75);
831	<	rmpf_final_mark->SetY(0.75);
832	<	rmpf_final_mark->SetTextSize(0.03);
833	<	rmpf_final_mark->Draw();
814	>
815	>	can->cd();
816	>	MPF_FaceValueAtPoint3->SetMarkerStyle(21);
817	>	MPF_FaceValueAtPoint3->SetMarkerColor(kBlue);
818	>	MPF_FaceValueAtPoint3->Fit("pol0","QE"); // quiet, use minos
819	>	MPF_FaceValueAtPoint3->Draw("AP");
820	>	MPF_FaceValueAtPoint3->GetXaxis()->SetTitle("p_{T}^{Z}");
821	>	MPF_FaceValueAtPoint3->GetXaxis()->CenterTitle();
822	>	MPF_FaceValueAtPoint3->GetYaxis()->SetTitle(((string)"<MPF>__{data}/<MPF>_{mc} for #alpha<"+any2string(FaceValueToBeTakenAt)).c_str());
823	>	MPF_FaceValueAtPoint3->GetYaxis()->CenterTitle();
824	>	MPF_FaceValueAtPoint3->Draw("AP");
825	>	TF1 *faceval_pol0=(TF1*)MPF_FaceValueAtPoint3->GetFunction("pol0");
826	>	faceval_pol0->SetLineColor(kBlue);
827	>	faceval_pol0->SetLineWidth(1);
828	>	faceval_pol0->Draw("same");
829	>	TF1 *faceval_pol0_complete = new TF1("faceval_pol0_complete","[0]+[1]*x");
830	>	faceval_pol0_complete->SetParameters(faceval_pol0->GetParameters());
831	>	stringstream faceval_result;
832	>	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<" << FaceValueToBeTakenAt << ")}";;
833	>	TText *faceval_final_mark = write_title(faceval_result.str());
834	>	faceval_final_mark->SetX(0.75);
835	>	faceval_final_mark->SetY(0.75);
836	>	faceval_final_mark->SetTextSize(0.03);
837	>	faceval_final_mark->Draw();
838		DrawPrelim();
839
840	<	stringstream filename2;
841	<	filename2 << "Extrapolation/FINAL_RESULT_MPF";
380	<	CompleteSave(can,filename2.str());
840	>	float FaceValResult=1/faceval_pol0->GetParameter(0);
841	>	float FaceValResultError=faceval_pol0->GetParError(0)/(faceval_pol0->GetParameter(0)*faceval_pol0->GetParameter(0));
842
843	+	stringstream filename2;
844	+	filename2.str("");
845	+	filename2 << "Extrapolation/FINAL_RESULT_MPF_FaceValp3";
846	+	CompleteSave(can,ContainerName+"/"+filename2.str());
847
848	<	RABS_FinalGraph->SetMarkerStyle(21);
849	<	RABS_FinalGraph->SetMarkerStyle(21);
850	<	RABS_FinalGraph->SetMarkerColor(kRed);
851	<	RABS_FinalGraph->Fit("pol0","QE"); // quiet, use minos
852	<	RABS_FinalGraph->Draw("AP");
853	<	RABS_FinalGraph->GetXaxis()->SetTitle("p_{T}^{Z}");
854	<	RABS_FinalGraph->GetXaxis()->CenterTitle();
855	<	RABS_FinalGraph->GetYaxis()->SetTitle("data/mc");
856	<	RABS_FinalGraph->GetYaxis()->CenterTitle();
857	<	RABS_FinalGraph->Draw("AP");
858	<	TF1 *rabs_pol0=(TF1*)RABS_FinalGraph->GetFunction("pol0");
859	<	stringstream rabs_result;
860	<	rabs_result << "#splitline{C_{abs}= " << std::setprecision(5) << 1/rabs_pol0->GetParameter(0) << " #pm " << std::setprecision(2) << rabs_pol0->GetParError(0)/(rabs_pol0->GetParameter(0)*rabs_pol0->GetParameter(0)) << "}{ (R_{abs})}";;
861	<	TText *rabs_final_mark = write_title(rabs_result.str());
862	<	rabs_final_mark->SetX(0.75);
863	<	rabs_final_mark->SetY(0.75);
864	<	rabs_final_mark->SetTextSize(0.03);
865	<	rabs_final_mark->Draw();
848	>	can->cd();
849	>	RABS_FaceValueAtPoint3->SetMarkerStyle(21);
850	>	RABS_FaceValueAtPoint3->SetMarkerColor(kRed);
851	>	RABS_FaceValueAtPoint3->Fit("pol0","QE"); // quiet, use minos
852	>	RABS_FaceValueAtPoint3->Draw("AP");
853	>	RABS_FaceValueAtPoint3->GetXaxis()->SetTitle("p_{T}^{Z}");
854	>	RABS_FaceValueAtPoint3->GetXaxis()->CenterTitle();
855	>	RABS_FaceValueAtPoint3->GetYaxis()->SetTitle(((string)"<R_{abs}>_{data}/<R_{abs}>_{mc} for #alpha<"+any2string(FaceValueToBeTakenAt)).c_str());
856	>	RABS_FaceValueAtPoint3->GetYaxis()->CenterTitle();
857	>	RABS_FaceValueAtPoint3->Draw("AP");
858	>	TF1 *faceval_pol0RABS=(TF1*)RABS_FaceValueAtPoint3->GetFunction("pol0");
859	>	TF1 *faceval_pol0RABS_complete = new TF1("faceval_pol0_complete","[0]+[1]*x");
860	>	faceval_pol0RABS_complete->SetParameters(faceval_pol0RABS->GetParameters());
861	>	faceval_pol0RABS->SetLineWidth(1);
862	>	faceval_pol0RABS->SetLineColor(kRed);
863	>	faceval_pol0RABS_complete->Draw("same");
864	>	stringstream RABSfaceval_result;
865	>	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<" << FaceValueToBeTakenAt << ")}";;
866	>	TText *RABSfaceval_final_mark = write_title(RABSfaceval_result.str());
867	>	RABSfaceval_final_mark->SetX(0.75);
868	>	RABSfaceval_final_mark->SetY(0.75);
869	>	RABSfaceval_final_mark->SetTextSize(0.03);
870	>	RABSfaceval_final_mark->Draw();
871		DrawPrelim();
872
873	+	float RABSFaceValResult=1/faceval_pol0RABS->GetParameter(0);
874	+	float RABSFaceValResultError=faceval_pol0RABS->GetParError(0)/(faceval_pol0RABS->GetParameter(0)*faceval_pol0RABS->GetParameter(0));
875	+
876		filename2.str("");
877	<	filename2 << "Extrapolation/FINAL_RESULT_RABS";
878	<	CompleteSave(can,filename2.str());
877	>	filename2 << "Extrapolation/FINAL_RESULT_RABS_FaceValp3";
878	>	CompleteSave(can,ContainerName+"/"+filename2.str());
879	>
880	>	cout << "FINAL RESULTS: " << endl;
881	>	if(!gofast) cout << "MPF : " << MPFResult  << " +/- " << MPFResultError  << endl;
882	>	if(!gofast) cout << "Rabs: " << RabsResult << " +/- " << RabsResultError << endl;
883	>	cout << "Face value at " << FaceValueToBeTakenAt << " :" << FaceValResult << " +/- " << FaceValResultError << "  (MPF) " << endl;
884	>	cout << "Face value at " << FaceValueToBeTakenAt << " :" << RABSFaceValResult << " +/- " << RABSFaceValResultError << "  (RABS) " << endl;
885	>
886	>	delete can;
887	>
888	>	results.MPF_Result_FaceValue=Value(FaceValResult,FaceValResultError);
889	>	results.Rabs_Result_FaceValue=Value(RABSFaceValResult,RABSFaceValResultError);
890	>	if(!gofast) results.MPF_Result_Extrapolated=Value(MPFResult,MPFResultError);
891	>	if(!gofast) results.Rabs_Result_Extrapolated=Value(RabsResult,RabsResultError);
892	>	cout << results << endl;
893	>	return results;
894	>	}
895	>
896	>
897	>	void DoPUStudy(string identifier) {
898	>
899	>	float numVtxcuts[7]={0,10,15,20};
900
901	+	TCanvas *can = new TCanvas("can","can");
902	+	TH1F *malpha[8];
903	+	TH1F *dalpha[8];
904	+	cout << identifier << ";";
905	+	cout << endl;
906	+	for(int i=1;i<5;i++) {
907	+	stringstream sSpecialCut;
908	+	if(i<4) {
909	+	sSpecialCut << "numVtx>" << numVtxcuts[i-1] << " && numVtx<" << numVtxcuts[i];
910	+	cout << numVtxcuts[i-1] << "  < nVtx < " << numVtxcuts[i] << ";";
911	+	} else {
912	+	sSpecialCut << "numVtx>" << numVtxcuts[i-1];
913	+	cout << numVtxcuts[i-1] << ";";
914	+	}
915	+
916	+	sSpecialCut << " && " << identifier << "_alpha<0.3";
917	+	TCut SpecialCut(sSpecialCut.str().c_str());
918	+
919	+
920	+	malpha[i] = allsamples.Draw("malpha"+any2string(i),"mpf",1000,0,2,"MPF","events",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&SpecialCut,mc,luminosity);
921	+	dalpha[i] = allsamples.Draw("dalpha"+any2string(i),"mpf",1000,0,2,"MPF","events",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&SpecialCut,data,luminosity);
922	+
923	+
924	+	float a=dalpha[i]->GetMean();
925	+	float b=malpha[i]->GetMean();
926	+	float da=dalpha[i]->GetMeanError();
927	+	float db=malpha[i]->GetMeanError();
928	+	float factor = a / b;
929	+	float error = TMath::Sqrt(  (1/(b*b)) * (da*da) + ((a*a)/(b*b))*db*db);
930	+
931	+	cout << malpha[i]->GetMean() << ";" << malpha[i]->GetMeanError() << ";" << dalpha[i]->GetMean() << ";" << dalpha[i]->GetMeanError() << ";" << malpha[i]->Integral()<<";"<<dalpha[i]->Integral() << ";"<<factor << ";" << error <<  endl;
932	+	delete malpha[i];
933	+	delete dalpha[i];
934	+	}
935
936		delete can;
937
938	+	}
939	+
940	+	void ScenarioComparison() {
941	+	//very dumb way to do this.
942	+
943	+	TGraphAsymmErrors *gr[5];
944	+	TF1 *fit[5];
945	+	bool dofit=true;
946	+
947	+	TCanvas *can = new TCanvas("can","can",500,500);
948	+
949	+	TLegend *leg = new TLegend(0.2,0.2,0.5,0.45);
950	+	leg->SetBorderSize(0);
951	+	leg->SetFillColor(kWhite);
952	+
953	+	for(int i=0;i<6;i++) {
954	+	gr[i] = new TGraphAsymmErrors();
955	+	double x[4]   = {5,12.5,17.5,25};
956	+	double dxl[4] = {5,2.5,2.5,5};
957	+	double dxh[4] = {5,2.5,2.5,10};
958	+
959	+	string name="";
960	+	string color="";
961	+
962	+	if(i==0) {
963	+	double y[4]   = {1.01935,1.01558,0.996438,1.01901};
964	+	double dyl[4] = {0.0132973,0.0131661,0.017629,0.0319767};
965	+	double dyh[4] = {0.0132973,0.0131661,0.017629,0.0319767};
966	+	name="30/30";
967	+	gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
968	+	color="#a01d99";
969	+	}
970	+	if(i==1) {
971	+	double y[4]   = {1.00577,1.00706,0.991029,0.954098};
972	+	double dyl[4] = {0.0135575,0.0136548,0.0170778,0.0263919};
973	+	double dyh[4] = {0.0135575,0.0136548,0.0170778,0.0263919};
974	+	gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
975	+	name="30/15";
976	+	color="#2464bc";
977	+	}
978	+	if(i==2) {
979	+	double y[4]   = {1.023,1.00648,1.02217,1.03181};
980	+	double dyl[4] = {0.0157372,0.0187778,0.0261512,0.0422295};
981	+	double dyh[4] = {0.0157372,0.0187778,0.0261512,0.0422295};
982	+	gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
983	+	name="15/15";
984	+	color="#f49230";
985	+	}
986	+	if(i==3) {
987	+	double y[4]   = {1.02029,1.01338,0.999769,1.00357};
988	+	double dyl[4] = {0.0129509,0.0120127,0.0160537,0.0268972};
989	+	double dyh[4] = {0.0129509,0.0120127,0.0160537,0.0268972};
990	+	gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
991	+	name="CHS 30/30";
992	+	color="#c72f3c";
993	+	}
994	+	if(i==4) {
995	+	double y[4]   = {1.00274,1.01056,0.993214,0.997963};
996	+	double dyl[4] = {0.0127179,0.0140344,0.0162814,0.0278914};
997	+	double dyh[4] = {0.0127179,0.0140344,0.0162814,0.0278914};
998	+	gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
999	+	name="CHS 30/15";
1000	+	color="#7cb433";
1001	+	}
1002	+	if(i==5) {
1003	+	double y[4]   = {1.04103,1.03726,1.04649,1.05383};
1004	+	double dyl[4] = {0.0156465,0.0168194,0.0231982,0.0376114};
1005	+	double dyh[4] = {0.0156465,0.0168194,0.0231982,0.0376114};
1006	+	gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
1007	+	name="CHS 15/15";
1008	+	color="#DF013A";
1009	+	}
1010	+
1011	+
1012	+	gr[i]->SetTitle();
1013	+	gr[i]->GetXaxis()->SetTitle("N(Vtx)");
1014	+	gr[i]->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
1015	+	gr[i]->GetXaxis()->CenterTitle();
1016	+	gr[i]->GetYaxis()->CenterTitle();
1017	+	gr[i]->SetName(name.c_str());
1018	+	gr[i]->SetMarkerSize(0.00001);
1019	+	gr[i]->GetYaxis()->SetRangeUser(0.90,1.10);
1020	+	gr[i]->SetLineColor(TColor::GetColor(color.c_str()));
1021	+	gr[i]->SetMarkerColor(TColor::GetColor(color.c_str()));
1022	+	if(dofit) {
1023	+	gr[i]->Fit("pol1");
1024	+	fit[i] = (TF1*)gr[i]->GetFunction("pol1");
1025	+	fit[i]->SetLineColor(TColor::GetColor(color.c_str()));
1026	+	}
1027	+	leg->AddEntry(gr[i],name.c_str(),"l");
1028	+	if(i==0) gr[i]->Draw("AP*");
1029	+	else gr[i]->Draw("P");
1030	+	}
1031	+	leg->Draw();
1032	+	DrawPrelim();
1033	+
1034	+	CompleteSave(can,"ScenarioComparison");
1035	+	delete can;
1036
1037		}
1038
1039	+	void ScenarioComparisonInclusive() {
1040	+	//dumbest way ever to do this. but ok we only need to do it once.
1041	+	TGraphAsymmErrors *gr[5];
1042	+	TF1 *fit[5];
1043	+	bool dofit=true;
1044	+
1045	+	TCanvas *can = new TCanvas("can","can",500,500);
1046	+
1047	+	TLegend *leg = new TLegend(0.2,0.2,0.65,0.45);
1048	+	leg->SetBorderSize(0);
1049	+	leg->SetFillColor(kWhite);
1050	+
1051	+	for(int i=0;i<6;i++) {
1052	+	gr[i] = new TGraphAsymmErrors();
1053	+	double x[4]   = {5,12.5,17.5,25};
1054	+	double dxl[4] = {5,2.5,2.5,5};
1055	+	double dxh[4] = {5,2.5,2.5,5};
1056	+
1057	+	string name="";
1058	+	string color="";
1059	+
1060	+	if(i==0) {
1061	+	double y[4]   = {1.01093,1.01379,1.0147,1.02766};
1062	+	double dyl[4] = {0.00243596,0.00242534,0.00341941,0.00600663};
1063	+	double dyh[4] = {0.00243596,0.00242534,0.00341941,0.00600663};
1064	+	name="30/30";
1065	+	gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
1066	+	color="#a01d99";
1067	+	}
1068	+	if(i==1) {
1069	+	double y[4]   = {1.0068,1.00779,0.999692,0.994484};
1070	+	double dyl[4] = {0.00249233,0.00257169,0.00371399,0.00631138};
1071	+	double dyh[4] = {0.00249233,0.00257169,0.00371399,0.00631138};
1072	+	gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
1073	+	name="30/15";
1074	+	color="#2464bc";
1075	+	}
1076	+	if(i==2) {
1077	+	double y[4]   = {1.01208,1.01991,1.02185,1.03479};
1078	+	double dyl[4] = {0.00366745,0.00400001,0.00584054,0.0104661};
1079	+	double dyh[4] = {0.00366745,0.00400001,0.00584054,0.0104661};
1080	+	gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
1081	+	name="15/15";
1082	+	color="#f49230";
1083	+	}
1084	+	if(i==3) {
1085	+	double y[4]   = {1.01177,1.01568,1.01616,1.03211};
1086	+	double dyl[4] = {0.00217317,0.00203014,0.00276838,0.00465643};
1087	+	double dyh[4] = {0.00217317,0.00203014,0.00276838,0.00465643};
1088	+	gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
1089	+	name="CHS 30/30";
1090	+	color="#c72f3c";
1091	+	}
1092	+	if(i==4) {
1093	+	double y[4]   = {1.00751,1.00964,1.00804,1.01073};
1094	+	double dyl[4] = {0.00224851,0.00215563,0.0029592,0.00499805};
1095	+	double dyh[4] = {0.00224851,0.00215563,0.0029592,0.00499805};
1096	+	gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
1097	+	name="CHS 30/15";
1098	+	color="#7cb433";
1099	+	}
1100	+	if(i==5) {
1101	+	double y[4]   = {1.01557,1.0318,1.03535,1.06331};
1102	+	double dyl[4] = {0.00300783,0.00302486,0.00421064,0.00739416};
1103	+	double dyh[4] = {0.00300783,0.00302486,0.00421064,0.00739416};
1104	+	gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
1105	+	name="CHS 15/15";
1106	+	color="#DF013A";
1107	+	}
1108	+
1109	+	gStyle->SetOptFit(0);
1110	+	gr[i]->SetTitle();
1111	+	gr[i]->GetXaxis()->SetTitle("N(Vtx)");
1112	+	gr[i]->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
1113	+	gr[i]->GetXaxis()->CenterTitle();
1114	+	gr[i]->GetYaxis()->CenterTitle();
1115	+	gr[i]->SetName(name.c_str());
1116	+	gr[i]->SetMarkerSize(0.00001);
1117	+	gr[i]->GetYaxis()->SetRangeUser(0.90,1.10);
1118	+	gr[i]->SetLineColor(TColor::GetColor(color.c_str()));
1119	+	gr[i]->SetMarkerColor(TColor::GetColor(color.c_str()));
1120	+	if(dofit) {
1121	+	gr[i]->Fit("pol1");
1122	+	fit[i] = (TF1*)gr[i]->GetFunction("pol1");
1123	+	fit[i]->SetLineColor(TColor::GetColor(color.c_str()));
1124	+	}
1125	+
1126	+	stringstream nameplus;
1127	+	nameplus << name << " (slope: " << std::setprecision(3) << fit[i]->GetParameter(1) << ")";
1128	+
1129	+	leg->AddEntry(gr[i],nameplus.str().c_str(),"l");
1130	+	if(i==0) gr[i]->Draw("AP*");
1131	+	else gr[i]->Draw("P");
1132	+	}
1133	+	leg->Draw();
1134	+	DrawPrelim();
1135	+
1136	+	CompleteSave(can,"ScenarioComparisonInclusive");
1137	+	delete can;
1138	+
1139	+	}
1140
1141	+	void compare_selection(string identifier) {
1142	+	bool recognized_scenario=false;
1143	+
1144	+	cout << "Running with identifier " << identifier << endl;
1145	+	if(identifier=="Zb30") {
1146	+	LeadingB=TCut ("Zb30_bTagProbCSVBP[0]>0.898");
1147	+	EtaB=TCut("abs(Zb30_pfJetGoodEta[0])<1.3");
1148	+	PhiZcut=TCut("abs(Zb30_pfJetDphiZ[0])>2.7");
1149	+	recognized_scenario=true;
1150	+	}
1151	+	if(identifier=="Zb3010") {
1152	+	LeadingB=TCut ("Zb3010_bTagProbCSVBP[0]>0.898");
1153	+	EtaB=TCut("abs(Zb3010_pfJetGoodEta[0])<1.3");
1154	+	PhiZcut=TCut("abs(Zb3010_pfJetDphiZ[0])>2.7");
1155	+	recognized_scenario=true;
1156	+	}
1157	+
1158	+	if(identifier=="ZbCHS3010") {
1159	+	LeadingB=TCut ("ZbCHS1510_bTagProbCSVBP[0]>0.898");
1160	+	EtaB=TCut("abs(ZbCHS1510_pfJetGoodEta[0])<1.3");
1161	+	PhiZcut=TCut("abs(ZbCHS1510_pfJetDphiZ[0])>2.7");
1162	+	recognized_scenario=true;
1163	+	}
1164	+
1165	+	if(identifier=="Zb40") {
1166	+	LeadingB=TCut ("Zb40_bTagProbCSVBP[0]>0.898");
1167	+	EtaB=TCut("abs(Zb40_pfJetGoodEta[0])<1.3");
1168	+	PhiZcut=TCut("abs(Zb40_pfJetDphiZ[0])>2.7");
1169	+	recognized_scenario=true;
1170	+	}
1171	+
1172	+	if(identifier=="ZbCHS301030") {
1173	+	LeadingB=TCut ("ZbCHS301030_bTagProbCSVBP[0]>0.898");
1174	+	EtaB=TCut("abs(ZbCHS301030_pfJetGoodEta[0])<1.3");
1175	+	PhiZcut=TCut("abs(ZbCHS301030_pfJetDphiZ[0])>2.7");
1176	+	recognized_scenario=true;
1177	+	}
1178	+
1179	+	if(identifier=="Zb1530") {
1180	+	LeadingB=TCut ("Zb1530_bTagProbCSVBP[0]>0.898");
1181	+	EtaB=TCut("abs(Zb1530_pfJetGoodEta[0])<1.3");
1182	+	PhiZcut=TCut("abs(Zb1530_pfJetDphiZ[0])>2.7");
1183	+	recognized_scenario=true;
1184	+	}
1185	+
1186	+	if(identifier=="ZbCHS301010") {
1187	+	LeadingB=TCut ("ZbCHS301010_bTagProbCSVBP[0]>0.898");
1188	+	EtaB=TCut("abs(ZbCHS301010_pfJetGoodEta[0])<1.3");
1189	+	PhiZcut=TCut("abs(ZbCHS301010_pfJetDphiZ[0])>2.7");
1190	+	recognized_scenario=true;
1191	+	}
1192	+
1193	+	if(identifier=="Zb1510") {
1194	+	LeadingB=TCut ("Zb1510_bTagProbCSVBP[0]>0.898");
1195	+	EtaB=TCut("abs(Zb1510_pfJetGoodEta[0])<1.3");
1196	+	PhiZcut=TCut("abs(Zb1510_pfJetDphiZ[0])>2.7");
1197	+	recognized_scenario=true;
1198	+	}
1199	+
1200	+	if(identifier=="ZbCHS301010") {
1201	+	LeadingB=TCut ("ZbCHS301010_bTagProbCSVBP[0]>0.898");
1202	+	EtaB=TCut("abs(ZbCHS301010_pfJetGoodEta[0])<1.3");
1203	+	PhiZcut=TCut("abs(ZbCHS301010_pfJetDphiZ[0])>2.7");
1204	+	recognized_scenario=true;
1205	+	}
1206	+
1207	+	if(identifier=="ZbMikko") {
1208	+	LeadingB=TCut ("ZbMikko_bTagProbCSVBP[0]>0.898 && ZbMikko__IsClean");
1209	+	EtaB=TCut("abs(ZbMikko_pfJetGoodEta[0])<1.3 && ZbMikko__IsClean");
1210	+	PhiZcut=TCut("abs(ZbMikko_pfJetDphiZ[0])>2.7 && ZbMikko__IsClean");
1211	+	recognized_scenario=true;
1212	+	}
1213	+
1214	+	if(identifier=="ZbCHS3010_SecEta3") {
1215	+	LeadingB=TCut ("ZbCHS3010_SecEta3_bTagProbCSVBP[0]>0.898");
1216	+	EtaB=TCut("abs(ZbCHS3010_SecEta3_pfJetGoodEta[0])<1.3");
1217	+	PhiZcut=TCut("abs(ZbCHS3010_SecEta3_pfJetDphiZ[0])>2.7");
1218	+	recognized_scenario=true;
1219	+	}
1220	+
1221	+	if(identifier=="ZbCHS3010_SecEta5") {
1222	+	LeadingB=TCut ("ZbCHS3010_SecEta5_bTagProbCSVBP[0]>0.898");
1223	+	EtaB=TCut("abs(ZbCHS3010_SecEta5_pfJetGoodEta[0])<1.3");
1224	+	PhiZcut=TCut("abs(ZbCHS3010_SecEta5_pfJetDphiZ[0])>2.7");
1225	+	recognized_scenario=true;
1226	+	}
1227	+
1228	+	if(identifier=="ZbCHS3010_p5Clean") {
1229	+	LeadingB=TCut ("ZbCHS3010_p5Clean_bTagProbCSVBP[0]>0.898 && ZbCHS3010_p5Clean_IsClean");
1230	+	EtaB=TCut("abs(ZbCHS3010_p5Clean_pfJetGoodEta[0])<1.3 && ZbCHS3010_p5Clean_IsClean");
1231	+	PhiZcut=TCut("abs(ZbCHS3010_p5Clean_pfJetDphiZ[0])>2.7 && ZbCHS3010_p5Clean_IsClean");
1232	+	recognized_scenario=true;
1233	+	}
1234	+
1235	+	if(identifier=="ZbCHS3010") {
1236	+	LeadingB=TCut ("ZbCHS3010_bTagProbCSVBP[0]>0.898");
1237	+	EtaB=TCut("abs(ZbCHS3010_pfJetGoodEta[0])<1.3");
1238	+	PhiZcut=TCut("abs(ZbCHS3010_pfJetDphiZ[0])>2.7");
1239	+	recognized_scenario=true;
1240	+	}
1241	+	if(identifier=="ZbCHS3010") {
1242	+	LeadingB=TCut ("ZbCHS3010_bTagProbCSVBP[0]>0.898");
1243	+	EtaB=TCut("abs(ZbCHS3010_pfJetGoodEta[0])<1.3");
1244	+	PhiZcut=TCut("abs(ZbCHS3010_pfJetDphiZ[0])>2.7");
1245	+	recognized_scenario=true;
1246	+	}
1247	+	if(identifier=="ZbCHS30") {
1248	+	LeadingB=TCut ("ZbCHS30_bTagProbCSVBP[0]>0.898");
1249	+	EtaB=TCut("abs(ZbCHS30_pfJetGoodEta[0])<1.3");
1250	+	PhiZcut=TCut("abs(ZbCHS30_pfJetDphiZ[0])>2.7");
1251	+	recognized_scenario=true;
1252	+	}
1253	+
1254	+	assert(recognized_scenario);
1255	+
1256	+	cout << "Cuts: " << endl;
1257	+	cout << "   " << (const char*) LeadingB << endl;
1258	+	cout << "   " << (const char*) EtaB << endl;
1259	+	cout << "   " << (const char*) PhiZcut << endl;
1260	+
1261	+	cout << endl << endl;
1262	+
1263	+	//  ScenarioComparison();
1264	+	//  ScenarioComparisonInclusive();
1265	+
1266	+	//    DoPUStudy(identifier);
1267	+	//    SpecialCutFlow(identifier);
1268	+
1269	+	//    draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",17000,30,200,"Z p_{T}","Official/"+identifier+"/Zpt__"+identifier,1);
1270	+	//    draw_kin_variable(identifier+"_pfJetGoodPt[0]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",40,0,200,"Leading Jet Pt (B)","Official/"+identifier+"/LeadingJetPt__"+identifier,1);
1271	+	//    draw_kin_variable(identifier+"_pfJetGoodPt[1]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",40,0,200,"Sub-Leading Jet Pt [GeV]","Official/"+identifier+"/Jet2Pt__"+identifier,1);
1272	+	//    draw_kin_variable(identifier+"_pfJetGoodPt[1]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",20,0,2,"p_{T}^{J2} / p_{T}^{Z}","Official/"+identifier+"/SubLeadingJetPt_Over_ZPt__"+identifier,1);
1273	+	draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<50"),20,0,2,"#alpha","Official/"+identifier+"/alpha_pt_30_to_50__"+identifier,1);
1274	+	draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>50&&pt<75"),20,0,2,"#alpha","Official/"+identifier+"/alpha_50_to_75__"+identifier,1);
1275	+	draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>75&&pt<125"),20,0,2,"#alpha","Official/"+identifier+"/alpha_pt_75_to_125__"+identifier,1);
1276	+	draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>125&&pt<1000"),20,0,2,"#alpha","Official/"+identifier+"/alpha_pt_125_to_1000__"+identifier,1);
1277	+	draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut,20,0,2,"#alpha","Official/"+identifier+"/alpha_full__"+identifier,1);
1278	+	draw_kin_variable(identifier+"_pfBJetDphiZ[0]",ZplusBsel&&LeadingB&&"pt>30&&pt<1000&&pfBJetDphiZ[0]>0",30,0,3.2,"#delta#phi (Z,b lead)","Official/"+identifier+"/DeltaPhiZBlead__"+identifier,0);
1279	+
1280	+	draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&pt<1000&&"+identifier+"_alpha<0.3").c_str()),20,0,2,"Z+b MPF","Official/"+identifier+"/mpf_alpha_smaller_0p3___"+identifier,0,"#alpha<0.3, 30 GeV < p_{T}^{Z} < 1000 GeV");
1281	+	draw_kin_variable(identifier+"_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&pt<1000&&"+identifier+"_alpha<0.3").c_str()),20,0,2,"p_{T} b-jet / p_{T} Z","Official/"+identifier+"/ptb_over_ptz___alpha_smaller_0p3______"+identifier,0,"#alpha<0.3, 30 GeV < p_{T}^{Z} < 1000 GeV");
1282	+	draw_kin_variable(identifier+"_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&pt<50&&"+identifier+"_alpha<0.3").c_str()),20,0,2,"p_{T} b-jet / p_{T} Z","Official/"+identifier+"/ptb_over_ptz___alpha_smaller_0p3___pt3050___"+identifier,0,"#alpha<0.3, 30 GeV < p_{T}^{Z} < 1000 GeV}{30<p_{T}^{Z}<50 GeV}");
1283	+
1284	+	}
1285	+
1286	+	void GetNumberEventsInsideOutsideAlphaWindow(TCut specialcut, string title) {
1287	+
1288	+	//  TCut cut =  ZplusBsel&&LeadingB&&PhiZcut&&specialcut&&TCut("ZbCHS3010_alpha<0.3");
1289	+	TCut cut = specialcut;
1290	+	TCut in  =  EtaB;
1291	+	TCut out =  TCut("abs(ZbCHS3010_pfJetGoodEta[0])>1.3");
1292	+
1293	+	TH1F *data_IN  = allsamples.Draw("data_IN",    "mll",1,0,150, "nothing [GeV]", "events", cut&&in, data,luminosity);
1294	+	TH1F *data_OUT = allsamples.Draw("data_OUT",   "mll",1,0,150, "nothing [GeV]", "events", cut&&out,data,luminosity);
1295	+	TH1F *mc_IN    = allsamples.Draw("mc_IN",      "mll",1,0,150, "nothing [GeV]", "events", cut&&in, mc  ,luminosity);
1296	+	TH1F *mc_OUT   = allsamples.Draw("mc_OUT",     "mll",1,0,150, "nothing [GeV]", "events", cut&&out,mc  ,luminosity);
1297	+
1298	+	dout << title << " : " << endl;
1299	+	dout << "  Data: In " << data_IN->Integral() << " vs out " << data_OUT->Integral() << endl;
1300	+	dout << "  MC  : In " << mc_IN->Integral()   << " vs out " << mc_OUT->Integral()   << endl;
1301	+	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;
1302	+	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;
1303	+
1304	+	delete data_IN;
1305	+	delete data_OUT;
1306	+	delete mc_IN;
1307	+	delete mc_OUT;
1308	+	}
1309	+
1310	+	void TEST_GetNumberEventsInsideOutsideAlphaWindow() {
1311	+	TCanvas *ca = new TCanvas("ca","ca");
1312	+	//   GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==0"), "ee");
1313	+	//   GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==1"), "mm");
1314	+	//   GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2"), "ee/mm");
1315	+
1316	+	cout << "BASE SELECTION" << endl;
1317	+	GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==0"), "ee");
1318	+	GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==1"), "mm");
1319	+	GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2"), "ee/mm");
1320
1321	<	void do_basic_ZB_analysis() {
1322	<	TCanvas *zbcanvas = new TCanvas("zbcanvas","zbcanvas");
1323	<	write_warning(__FUNCTION__,"TESTING NEW DATA MC AGREEMENT");
1324	<
1325	<	/*
1326	<	//  print_all_b_yields();
1327	<	//  draw_mpf_vars();
1328	<	draw_Zb_kin_vars();
1329	<
1330	<	data_mc_agreement_2d();
1331	<	*/
1332	<	new_data_mc_agreement_2d();
1321	>	cout << "BASE SELECTION: Z+b" << endl;
1322	>	GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
1323	>	GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
1324	>	GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2"), "ee/mm");
1325	>
1326	>	cout << "Leading B" << endl;
1327	>	GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
1328	>	GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
1329	>	GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2"), "ee/mm");
1330	>
1331	>	cout << "PhiZcut" << endl;
1332	>	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
1333	>	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
1334	>	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2"), "ee/mm");
1335	>
1336	>	cout << "alpha cut" << endl;
1337	>	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0")&&TCut("ZbCHS3010_alpha<0.3"), "ee");
1338	>	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1")&&TCut("ZbCHS3010_alpha<0.3"), "mm");
1339	>	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2")&&TCut("ZbCHS3010_alpha<0.3"), "ee/mm");
1340	>
1341	>	delete ca;
1342	>
1343	>
1344	>	}
1345	>
1346	>
1347	>
1348	>
1349	>	void do_basic_ZB_analysis(bool DoCompleteAnalysis) {
1350	>	STDWEIGHT=TCut(cutWeight);
1351	>	cutWeight=TCut(STDWEIGHT*TCut("ZbCHS3010_BTagWgtT"));
1352	>	cout << "The lepton requirement is " << (const char*) leptoncut << endl;
1353	>	bool doquick=false;
1354	>
1355	>	TCanvas *zbcanvas = new TCanvas("zbcanvas","zbcanvas");
1356	>
1357	>	//   ScenarioComparison();
1358	>	//   ScenarioComparisonInclusive();
1359	>
1360	>	//   compare_selection("ZbCHS3010_p5Clean");
1361	>	//   compare_selection("ZbCHS3010");
1362	>	//   compare_selection("Zb30");
1363	>	//   compare_selection("Zb3010");
1364	>	//   compare_selection("Zb1510");
1365	>	//   compare_selection("ZbCHS30");
1366	>	//   compare_selection("ZbCHS3010");
1367	>	//   compare_selection("ZbCHS3010");
1368	>	//   compare_selection("ZbCHS301010");
1369	>	//   compare_selection("ZbCHS3010_SecEta3");
1370	>	//   compare_selection("ZbCHS3010_SecEta5");
1371	>	//   compare_selection("Zb1510");
1372	>	//   compare_selection("ZbCHS301010");
1373	>	//   compare_selection("ZbMikko");
1374	>	//   compare_selection("ZbCHS3010");
1375	>	//   compare_selection("Zb40");
1376	>
1377	>	if(!doquick) {
1378	>	//draw_kin_variable("mll","",28,60.0,200.0,"m_{ll}","mll",1,"");//nice and quick test to see if the normalization is ok (i.e. all data is processed etc.)
1379	>	//print_all_b_yields();
1380	>	draw_mpf_vars();
1381	>	//DrawEvilCutFlow();
1382	>
1383	>	draw_Zb_kin_vars();
1384	>
1385	>	}
1386	>
1387	>	//  GetNumberEventsInsideOutsideAlphaWindow();
1388	>
1389	>	if(DoCompleteAnalysis) {
1390	>	/* need to run the following scenarios:
1391	>	* - normal
1392	>	* - inclusive
1393	>	* - medium WP
1394	>	* - alpha up
1395	>	* - alpha down
1396	>	*/
1397	>
1398	>	bool fast=true;bool slow=false;//don't change these
1399	>
1400	>	ZbResultContainer inclusive = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","inclusive",TCut("mll>0"),TCut("1.0"),0.3);
1401	>	ZbResultContainer Reference = new_data_mc_agreement_2d(slow,"ZbCHS3010_alpha","reference",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"),0.3);
1402	>
1403	>	ZbResultContainer NeutrinoQ = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","NeutrinoQ",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtTUp*(1.0*softMuon)"),0.3);
1404	>	ZbResultContainer ANeutrinoQ = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","ANeutrinoQ",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtTUp*(1.0*(1.0-softMuon))"),0.3);
1405	>
1406	>	ZbResultContainer effmistagUp = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","effmistagUp",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtTUp"),0.3);
1407	>	ZbResultContainer effmistagDown = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","effmistagDown",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtTDown"),0.3);
1408	>	ZbResultContainer medium = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","medium",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.679"),TCut("ZbCHS3010_BTagWgtM"),0.3);
1409	>	ZbResultContainer loose = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","loose",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.244"),TCut("ZbCHS3010_BTagWgtL"),0.3);
1410	>	ZbResultContainer AlphaUp = new_data_mc_agreement_2d(fast,"ZbCHS3010_alphaUp","AlphaUp",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"),0.3);
1411	>	ZbResultContainer AlphaDown = new_data_mc_agreement_2d(fast,"ZbCHS3010_alphaDown","AlphaDown",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"),0.3);
1412	>
1413	>	ZbResultContainer L5inclusive = new_data_mc_agreement_2d(fast,"ZbCHS3010_alphaL5","L5inclusive",TCut("mll>0"),TCut("1.0"),0.3);
1414	>	ZbResultContainer L5Reference = new_data_mc_agreement_2d(fast,"ZbCHS3010_alphaL5","L5reference",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"),0.3);
1415	>
1416	>	ZbResultContainer AlphaThresholdVariationDown10 =    new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationDown10",   TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),   TCut("ZbCHS3010_BTagWgtT"),0.3-0.1*0.3);
1417	>	ZbResultContainer AlphaThresholdVariationUp10 =      new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationUp10",     TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),   TCut("ZbCHS3010_BTagWgtT"),0.3+0.1*0.3);
1418	>	ZbResultContainer AlphaThresholdVariationDown30 =    new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationDown30",   TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),   TCut("ZbCHS3010_BTagWgtT"),0.3-0.3*0.3);
1419	>	ZbResultContainer AlphaThresholdVariationUp30 =      new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationUp30",     TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),   TCut("ZbCHS3010_BTagWgtT"),0.3+0.3*0.3);
1420	>	ZbResultContainer AlphaThresholdVariationDown10inc = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationDown10inc",TCut("mll>0"),                              TCut("1.0"),0.3-0.1*0.3);
1421	>	ZbResultContainer AlphaThresholdVariationUp10inc =   new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationUp10inc",  TCut("mll>0"),                              TCut("1.0"),0.3+0.1*0.3);
1422	>	ZbResultContainer AlphaThresholdVariationDown30inc = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationDown30inc",TCut("mll>0"),                              TCut("1.0"),0.3-0.3*0.3);
1423	>	ZbResultContainer AlphaThresholdVariationUp30inc =   new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationUp30inc",  TCut("mll>0"),                              TCut("1.0"),0.3+0.3*0.3);
1424	>
1425	>	//and now let's compute some systematics!
1426	>	//1 alpha variation
1427	>	float SysMPF_Alpha_down = Reference.MPF_Result_FaceValue.getValue()-AlphaDown.MPF_Result_FaceValue.getValue();
1428	>	float SysMPF_Alpha_up = Reference.MPF_Result_FaceValue.getValue()-AlphaUp.MPF_Result_FaceValue.getValue();
1429	>
1430	>	//2 btagging efficiency/mistag correction
1431	>	float SysMPF_EFF_down = Reference.MPF_Result_FaceValue.getValue()-effmistagDown.MPF_Result_FaceValue.getValue();
1432	>	float SysMPF_EFF_up = Reference.MPF_Result_FaceValue.getValue()-effmistagUp.MPF_Result_FaceValue.getValue();
1433	>
1434	>	//3 Presence of soft muons (->neutrino question)
1435	>	float SysMPF_Neutrino_up = Reference.MPF_Result_FaceValue.getValue()-NeutrinoQ.MPF_Result_FaceValue.getValue();
1436	>	float SysMPF_Neutrino_down = Reference.MPF_Result_FaceValue.getValue()-ANeutrinoQ.MPF_Result_FaceValue.getValue();
1437	>
1438	>	//4 purity
1439	>	zbcanvas->cd();
1440	>	gStyle->SetOptFit(0);
1441	>	write_info(__FUNCTION__,"Don't know the exact purity at the moment, still needs to be determined!");
1442	>	TH1F *purityhisto = new TH1F("purityhisto","purityhisto",1,0,1);
1443	>	purityhisto->GetXaxis()->SetTitle("Purity");
1444	>	purityhisto->GetYaxis()->SetTitle("C_{abs}");
1445	>	purityhisto->GetXaxis()->CenterTitle();
1446	>	purityhisto->GetYaxis()->CenterTitle();
1447	>
1448	>	TGraphErrors *gSysMPF_Purity = new TGraphErrors();
1449	>
1450	>	gSysMPF_Purity->SetPoint(0,PurityInclusive,inclusive.MPF_Result_FaceValue.getValue());
1451	>	gSysMPF_Purity->SetPointError(0,0.0,inclusive.MPF_Result_FaceValue.getError());
1452	>
1453	>	gSysMPF_Purity->SetPoint(1,PurityLoose,loose.MPF_Result_FaceValue.getValue());//x value?
1454	>	gSysMPF_Purity->SetPointError(1,0.0,loose.MPF_Result_FaceValue.getError());
1455	>	gSysMPF_Purity->SetPoint(2,PurityMedium,medium.MPF_Result_FaceValue.getValue());//x value?
1456	>	gSysMPF_Purity->SetPointError(2,0.0,medium.MPF_Result_FaceValue.getError());
1457	>	gSysMPF_Purity->SetPoint(3,PurityTight,Reference.MPF_Result_FaceValue.getValue());//x value?
1458	>	gSysMPF_Purity->SetPointError(3,0.0,Reference.MPF_Result_FaceValue.getError());
1459	>
1460	>
1461	>	float min,max;
1462	>	FindMinMax(gSysMPF_Purity,min,max);
1463	>	purityhisto->SetStats(0);
1464	>	purityhisto->GetYaxis()->SetRangeUser(min,max);
1465	>
1466	>	gSysMPF_Purity->Fit("pol1");
1467	>	TPolyLine *purity_fit_uncert = GetFitUncertaintyShape(gSysMPF_Purity, "pol1", min, max,0.0,1.0);
1468	>	TF1* pufit = (TF1*)gSysMPF_Purity->GetFunction("pol1");
1469	>
1470	>	TF1* pufit_c = new TF1("pufit_c","[0]+[1]*x");
1471	>	pufit_c->SetParameters(pufit->GetParameters());
1472	>	pufit_c->SetLineColor(TColor::GetColor("#0101DF"));
1473	>	purity_fit_uncert->SetFillColor(TColor::GetColor("#5353FC"));
1474	>	pufit->SetLineColor(TColor::GetColor("#5353FC"));
1475	>
1476	>	purityhisto->Draw();
1477	>	purity_fit_uncert->Draw("f");
1478	>	gSysMPF_Purity->Draw("P*");
1479	>	pufit_c->Draw("same");
1480	>	DrawPrelim();
1481	>
1482	>	float ExtrapolatedResult=pufit->GetParameter(0)+1.0*pufit->GetParameter(1);
1483	>	float dExtrapolatedResult=sqrt(pufit->GetParError(0)*pufit->GetParError(0)+1.0*1.0*pufit->GetParError(1)*pufit->GetParError(1));
1484	>
1485	>	stringstream summary;
1486	>	summary << "#splitline{Reference: " << std::setprecision(4) << Reference.MPF_Result_FaceValue.getValue() << "+/-" << std::setprecision(4) << Reference.MPF_Result_FaceValue.getError() << "}{";
1487	>	summary << "Extrapolation: " << std::setprecision(4) << ExtrapolatedResult << " +/- " << std::setprecision(4) << dExtrapolatedResult << "}";
1488	>
1489	>	dout << "Have found an extrapolated result at 1.0 of " << ExtrapolatedResult << " ;  will use the difference between this result and the face value as ";
1490	>	dout << "purity-related systematic uncertainty" << endl;
1491	>
1492	>	TText *infobox = write_title(summary.str());
1493	>	infobox->SetX(0.75);
1494	>	infobox->SetTextSize(0.03);
1495	>	infobox->SetY(0.75);
1496	>	infobox->Draw();
1497	>
1498	>	CompleteSave(zbcanvas,"Systematics/Purity");
1499	>
1500	>
1501	>
1502	>	float SysMPF_Purity = abs(ExtrapolatedResult-Reference.MPF_Result_FaceValue.getValue());
1503	>
1504	>	float sys_alpha    = abs(SysMPF_Alpha_down)>abs(SysMPF_Alpha_up)?abs(SysMPF_Alpha_down):abs(SysMPF_Alpha_up);
1505	>	float sys_eff      = abs(SysMPF_EFF_down)>abs(SysMPF_EFF_up)?abs(SysMPF_EFF_down):abs(SysMPF_EFF_up);
1506	>	float sys_purity   = abs(SysMPF_Purity);
1507	>	float sys_neutrino = abs(SysMPF_Neutrino_down)>abs(SysMPF_Neutrino_up)?abs(SysMPF_Neutrino_down):abs(SysMPF_Neutrino_up);
1508	>	float sys          = sqrt(sys_alpha*sys_alpha+sys_eff*sys_eff+sys_purity*sys_purity+sys_neutrino*sys_neutrino);
1509	>
1510	>	dout << "MPF METHOD:" << endl;
1511	>	dout << "   Systematics :    down    up     (selected)" << endl;
1512	>	dout << "      Oversmearing         : " << SysMPF_Alpha_down << "   ,  " << SysMPF_Alpha_up << "   ( " << sys_alpha << ") " << endl;
1513	>	dout << "      eff/mistag variation : " << SysMPF_EFF_down << "   ,  " << SysMPF_EFF_up << "   ( " << sys_eff << ") " << endl;
1514	>	dout << "      purity               : " << SysMPF_Purity << "   ( " << sys_purity << " )" << endl;
1515	>	dout << "      neutrino             : " << SysMPF_Neutrino_down << "   ,  " << SysMPF_Neutrino_up << "   ( " << sys_neutrino << ") " << endl;
1516	>	dout << "      total                : " << sys << "     (note that this is the sys on C_{abs}^{b}, not yet on C_{corr} - there will also be an additional source!)" << endl;
1517	>	dout << endl << endl;
1518	>	dout << "   FINAL RESULTS : " << endl;
1519	>	dout << "     C_{abs}^{b}   = " << Reference.MPF_Result_FaceValue << " (stat) +/- " << sys << " (sys) " << endl;
1520	>	dout << "     C_{abs}^{inc} = " << inclusive.MPF_Result_FaceValue << " (stat) " << endl;
1521	>
1522	>	float sysfinal = Reference.MPF_Result_FaceValue.getValue()/inclusive.MPF_Result_FaceValue.getValue();
1523	>	sysfinal *= sqrt(sys*sys/(Reference.MPF_Result_FaceValue.getValue()*Reference.MPF_Result_FaceValue.getValue()));
1524	>	//yes, this could be simplified but it would look like a bug.
1525	>
1526	>
1527	>
1528	>
1529	>	dout << "Note: Varying the alpha cut (->face value) we get the following results: " << endl;
1530	>	dout << "-30%: " << AlphaThresholdVariationDown30.MPF_Result_FaceValue/AlphaThresholdVariationDown30inc.MPF_Result_FaceValue << " (stat)" << endl;
1531	>	dout << "-10%: " << AlphaThresholdVariationDown10.MPF_Result_FaceValue/AlphaThresholdVariationDown10inc.MPF_Result_FaceValue << " (stat)" << endl;
1532	>	dout << "  0%: " << Reference.MPF_Result_FaceValue/inclusive.MPF_Result_FaceValue << " (stat) +/- " << endl;
1533	>	dout << "+10%: " << AlphaThresholdVariationUp10.MPF_Result_FaceValue/AlphaThresholdVariationUp10inc.MPF_Result_FaceValue << " (stat)" << endl;
1534	>	dout << "+30%: " << AlphaThresholdVariationUp30.MPF_Result_FaceValue/AlphaThresholdVariationUp30inc.MPF_Result_FaceValue << " (stat)" << endl;
1535	>
1536	>	dout << "          XCheck for Rabs, varying the alpha cut (->face value) we get the following results: " << endl;
1537	>	dout << "             -30%: " << AlphaThresholdVariationDown30.Rabs_Result_FaceValue/AlphaThresholdVariationDown30inc.Rabs_Result_FaceValue << " (stat)" << endl;
1538	>	dout << "             -10%: " << AlphaThresholdVariationDown10.Rabs_Result_FaceValue/AlphaThresholdVariationDown10inc.Rabs_Result_FaceValue << " (stat)" << endl;
1539	>	dout << "               0%: " << Reference.Rabs_Result_FaceValue/inclusive.Rabs_Result_FaceValue << " (stat) +/- " << endl;
1540	>	dout << "             +10%: " << AlphaThresholdVariationUp10.Rabs_Result_FaceValue/AlphaThresholdVariationUp10inc.Rabs_Result_FaceValue << " (stat)" << endl;
1541	>	dout << "             +30%: " << AlphaThresholdVariationUp30.Rabs_Result_FaceValue/AlphaThresholdVariationUp30inc.Rabs_Result_FaceValue << " (stat)" << endl;
1542	>
1543	>
1544	>	float sys_alphavar_up   = abs(AlphaThresholdVariationUp30.MPF_Result_FaceValue.getValue()/AlphaThresholdVariationUp30inc.MPF_Result_FaceValue.getValue() - Reference.MPF_Result_FaceValue.getValue()/inclusive.MPF_Result_FaceValue.getValue());
1545	>	float sys_alphavar_down = abs(AlphaThresholdVariationDown30.MPF_Result_FaceValue.getValue()/AlphaThresholdVariationDown30inc.MPF_Result_FaceValue.getValue() - Reference.MPF_Result_FaceValue.getValue()/inclusive.MPF_Result_FaceValue.getValue());
1546	>	float sys_alphavar = sys_alphavar_up>sys_alphavar_down?sys_alphavar_up:sys_alphavar_down;
1547	>
1548	>	dout << " -> Additional systematic (from 30%) : " << sys_alphavar << endl;
1549	>	sysfinal=sqrt(sysfinal*sysfinal+sys_alphavar*sys_alphavar);
1550	>	dout << "Final systematic is : " << sysfinal << endl;
1551	>
1552	>	dout << "   ******************************************** <  FINAL RESULT > ******************************************** " << endl;
1553	>	dout << "     C_{corr}      = " << Reference.MPF_Result_FaceValue/inclusive.MPF_Result_FaceValue << " (stat) +/- " << sysfinal << " (sys) " << endl;
1554	>	dout << "   ******************************************** < /FINAL RESULT > ******************************************** " << endl;
1555	>
1556	>	dout << "Note that if L5 corrections are applied, we get " << L5Reference.MPF_Result_FaceValue/L5inclusive.MPF_Result_FaceValue << " (stat) " << endl;
1557	>
1558	>	LeadingB=TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898");
1559	>
1560	>	} else {
1561	>	ZbResultContainer Reference = new_data_mc_agreement_2d(false,"ZbCHS3010_alpha","reference",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"),0.3);
1562	>	}
1563	>
1564
1565		delete zbcanvas;
1566		}
1567	+
1568	+
1569	+
1570	+
1571	+	//*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*  DELETE EVERYTHING BELOW THIS LINE /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*//
1572	+
1573	+
1574	+	const char* concatenate(string bla1, string bla2) {
1575	+	stringstream bla;
1576	+	bla << bla1 << bla2;
1577	+	return bla.str().c_str();
1578	+	}
1579	+
1580	+	void SpecialCutFlow(string identifier) {
1581	+	stringstream MegaCut;
1582	+
1583	+
1584	+
1585	+	MegaCut << "   1*(" << (const char*) (ZplusBsel&&TCut(concatenate(identifier,"_pfJetGoodNumBtag>0"))) << ")";
1586	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB) << ")";
1587	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB) << ")";
1588	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut) << ")";
1589	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&pt<1000").c_str())) << ")";
1590	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&pt<1000&&"+identifier+"_alpha<0.3").c_str())) << ")";
1591	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&pt<1000&&"+identifier+"_alpha<0.2").c_str())) << ")";
1592	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&pt<1000&&"+identifier+"_alpha<0.15").c_str())) << ")";
1593	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&pt<1000&&"+identifier+"_alpha<0.1").c_str())) << ")";
1594	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&pt<1000&&"+identifier+"_alpha<0.05").c_str())) << ")";
1595	+	MegaCut << " ";
1596	+
1597	+	TCanvas *can = new TCanvas("can","can");
1598	+	can->SetLogy(1);
1599	+	TCut basecut(ZplusBsel);
1600	+
1601	+
1602	+
1603	+	TLegend *leg = allsamples.allbglegend();
1604	+
1605	+
1606	+	TH1F *data  = allsamples.Draw("data",    MegaCut.str(),11,-0.5,10.5, "", "events", basecut,data,luminosity);
1607	+	THStack mcm = allsamples.DrawStack("mc", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,mc,luminosity);
1608	+
1609	+	float runningsum=0;
1610	+	for(int i=data->GetNbinsX();i>0;i--) {
1611	+	runningsum+=data->GetBinContent(i);
1612	+	data->SetBinContent(i,runningsum);
1613	+	}
1614	+
1615	+	float minimum=data->GetMaximum();
1616	+
1617	+	TH1F* h;
1618	+	TIter nextOF(mcm.GetHists());
1619	+
1620	+	float mcsum=0;
1621	+	float zb=0;
1622	+	while ( h = (TH1F*)nextOF() ) {
1623	+	float runningsum=0;
1624	+	minimum=data->GetMaximum();//done deliberately at each step so get the minimum of the last (!) contribution
1625	+	for(int i=h->GetNbinsX();i>0;i--) {
1626	+	runningsum+=h->GetBinContent(i);
1627	+	h->SetBinContent(i,runningsum);
1628	+	if(runningsum<minimum&&runningsum>0) minimum=runningsum;
1629	+	}
1630	+	if(Contains(h->GetName(),"Z_b_")) {
1631	+	zb=h->GetBinContent(h->GetNbinsX());
1632	+	}
1633	+	mcsum+=h->GetBinContent(h->GetNbinsX());
1634	+	}
1635	+
1636	+
1637	+
1638	+	data->SetMinimum(0.05*minimum);
1639	+	can->cd(1)->SetBottomMargin(0.25);
1640	+	data->GetXaxis()->SetBinLabel(1,"Pre-selection");
1641	+	data->GetXaxis()->SetBinLabel(2,"Contains b jet");
1642	+	data->GetXaxis()->SetBinLabel(3,"Leading jet is b-jet");
1643	+	data->GetXaxis()->SetBinLabel(4,"|#eta_{b}|<1.3 ");
1644	+	data->GetXaxis()->SetBinLabel(5,"|#delta#phi(b,Z)|>2.7");
1645	+	data->GetXaxis()->SetBinLabel(6,"30<p_{T}^{Z}<1000 GeV");
1646	+	data->GetXaxis()->SetBinLabel(7,"#alpha < 0.3");
1647	+	data->GetXaxis()->SetBinLabel(8,"#alpha < 0.2");
1648	+	data->GetXaxis()->SetBinLabel(9,"#alpha < 0.15");
1649	+	data->GetXaxis()->SetBinLabel(10,"#alpha < 0.1");
1650	+	data->GetXaxis()->SetBinLabel(11,"#alpha < 0.05");
1651	+	data->GetXaxis()->LabelsOption("v");
1652	+
1653	+	stringstream purityinfo;
1654	+	purityinfo << "Purity: " << std::setprecision(3) << 100*zb/mcsum << " %";
1655	+	stringstream neventsinfo;
1656	+	neventsinfo << "Nevents: " << data->GetBinContent(data->GetNbinsX());
1657	+	TH1F *crap = new TH1F("crap","",1,0,1);
1658	+	crap->SetLineColor(kWhite);
1659	+	leg->AddEntry(crap,purityinfo.str().c_str(),"l");
1660	+	leg->AddEntry(crap,neventsinfo.str().c_str(),"l");
1661	+
1662	+	data->Draw("e1");
1663	+	mcm.Draw("histo,same");
1664	+	data->Draw("same,e1,axis");
1665	+	data->Draw("same,e1");
1666	+	//  data->Draw("same,TEXT");
1667	+
1668	+	leg->Draw();
1669	+
1670	+	DrawPrelim();
1671	+
1672	+	CompleteSave(can,"CutFlow___"+identifier);
1673	+
1674	+	delete crap;
1675	+	delete data;
1676	+	CleanLegends();
1677	+	DeleteStack(mcm);
1678	+
1679	+
1680	+	}

Diff Legend

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