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.1 by buchmann, Mon Jun 25 13:34:21 2012 UTC vs.
Revision 1.16 by buchmann, Wed Nov 28 17:51:14 2012 UTC

#	Line 23 | Line 23 | using namespace PlottingSetup;
23		namespace ZbData {
24		vector<float> data_over_mc;
25
26	<	TCut ZplusBsel("pt1>20&&pt2>20&&mll>60&&mll<120");
27	<	TCut LeadingB("bTagProbCSVBP[0]>0.679");
28	<	TCut EtaB("abs(pfJetGoodEta[0])<1.3");
29	<	TCut PhiZcut("pfJetDphiZ[0]>2.7");
26	>	TCut ZplusBsel("pt1>20&&pt2>20&&mll>60&&mll<120&&id1==id2");
27	>	TCut LeadingB("Zb3010_bTagProbCSVBP[0]>0.898");
28	>	TCut EtaB("abs(Zb3010_pfJetGoodEta[0])<1.3");
29	>	TCut PhiZcut("abs(Zb3010_pfJetDphiZ[0])>2.7");
30
31		}
32
33	+	class ZbPointResult {
34	+	public:
35	+	ZbPointResult(string var, float ptlow, float pthigh, float alphalow, float alphahigh, Value Data, Value MC);
36	+	float PtLow;
37	+	float PtHigh;
38	+	float AlphaLow;
39	+	float AlphaHigh;
40	+	string variable;
41	+	Value Data;
42	+	Value MC;
43	+	Value Ratio;
44	+	};
45	+
46	+	ZbPointResult::ZbPointResult(string var, float ptlow, float pthigh, float alphalow, float alphahigh, Value data, Value mc) :
47	+	variable(var), PtLow(ptlow), PtHigh(pthigh), AlphaLow(alphalow), AlphaHigh(alphahigh), Data(data), MC(mc) {
48	+	Ratio = data/mc;
49	+	}
50	+
51	+	std::ostream &operator<<(std::ostream &ostr, ZbPointResult &b)
52	+	{
53	+	ostr << b.PtLow << ";" << b.PtHigh << ";" << b.AlphaLow << ";" << b.AlphaHigh << ";" << b.Data.getValue() << ";" << b.Data.getError() << ";" << b.MC.getValue() << ";" << b.MC.getError() << ";" << b.Ratio.getValue() << ";" << b.Ratio.getError();
54	+	return ostr;
55	+	}
56	+
57	+
58	+	class ZbResultContainer {
59	+	public:
60	+	string name;
61	+	vector<ZbPointResult> MPF_Rabs_location;
62	+	Value MPF_Result_FaceValue;
63	+	Value Rabs_Result_FaceValue;
64	+	Value MPF_Result_Extrapolated;
65	+	Value Rabs_Result_Extrapolated;
66	+
67	+	ZbResultContainer(string name);
68	+	ZbResultContainer(const ZbResultContainer &old);
69	+	};
70	+
71	+
72	+	ZbResultContainer::ZbResultContainer(string _name) {
73	+	name=_name;
74	+	cout<< "Result Container for \"" << name << "\" has been initialized" << endl;
75	+
76	+	}
77	+
78	+	std::ostream &operator<<(std::ostream &ostr, ZbResultContainer &b)
79	+	{
80	+	ostr << "Results for " << b.name << endl;
81	+	ostr << " MPF and Rabs locations: " << endl;
82	+	for(int i=0;i<b.MPF_Rabs_location.size();i++) ostr << "   " << b.MPF_Rabs_location[i] << endl;
83	+	ostr << " Results: " << endl;
84	+	ostr << "   Face Value: " << endl;
85	+	ostr << "     MPF " << b.MPF_Result_FaceValue << endl;
86	+	ostr << "     Rabs " << b.Rabs_Result_FaceValue << endl;
87	+	ostr << "   Extrapolated: " << endl;
88	+	ostr << "     MPF " << b.MPF_Result_Extrapolated << endl;
89	+	ostr << "             Rabs " << b.Rabs_Result_Extrapolated << endl;
90	+
91	+	return ostr;
92	+	}
93	+
94	+
95	+
96		using namespace ZbData;
97
98	+
99	+	//*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*  DELETE EVERYTHING BETWEEN THIS LINE /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*//
100	+	void SpecialCutFlow(string identifier);
101	+
102	+	//*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*            AND THIS LINE             /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*//
103	+
104	+
105		void print_yield(TCut cut) {
106		TH1F *data  = allsamples.Draw("data",    "mll",1,0,10000000, "m_{ll} [GeV]", "events", cut,data,luminosity);
107	<	cout << "  Yields for " << (const char*) cut << endl;
108	<	cout << "     data: " << data->Integral() << endl;
109	<	THStack mcm = allsamples.DrawStack("mc", "mll",1,0,10000000, "m_{ll} [GeV}", "events", cut,mc,luminosity);
107	>	dout << "  Yields for " << (const char*) cut << endl;
108	>	dout << "     data: " << data->Integral() << endl;
109	>	THStack mcm = allsamples.DrawStack("mc", "mll",1,0,10000000, "m_{ll} [GeV]", "events", cut,mc,luminosity);
110		int nhists = mcm.GetHists()->GetEntries();
111	<	cout << "Going to loop over " << nhists << " histograms " << endl;
111	>	dout << "Going to loop over " << nhists << " histograms " << endl;
112		TFile *test = new TFile("test.root","RECREATE");
113		mcm.Write();
114		test->Close();
#	Line 46 | Line 116 | void print_yield(TCut cut) {
116		TH1* hist = static_cast<TH1*>(mcm.GetHists()->At(i));
117		cout << "     " << hist->GetName() << ": " << hist->Integral() << endl;
118		}
119	+
120		cout << endl << endl;
121	+
122		delete data;
123		}
124
125	<	void draw_kin_variable(string variable, TCut cut, int nbins, float min, float max, string xlabel, string saveas, bool dologscale) {
125	>	void draw_kin_variable(string variable, TCut cut, int nbins, float min, float max, string xlabel, string saveas, bool dologscale, string speciallabel="") {
126
127		TCanvas *ckin = new TCanvas("ckin","kin variable canvas");
128		ckin->SetLogy(dologscale);
129		TH1F *datahisto = allsamples.Draw("datahisto",variable,nbins,min,max, xlabel, "events",cut,data,luminosity);
130		datahisto->SetMarkerSize(DataMarkerSize);
131		THStack mcstack = allsamples.DrawStack("mcstack",variable,nbins,min,max, xlabel, "events",cut,mc,luminosity);
132	<	datahisto->SetMaximum(5.3*datahisto->GetMaximum());
132	>	if(dologscale) datahisto->SetMaximum(3*datahisto->GetMaximum());
133	>	else datahisto->SetMaximum(1.5*datahisto->GetMaximum());
134	>	TText *speciall;
135	>	if(speciallabel!="") {
136	>	speciall = new TText(0.2,0.8,speciallabel.c_str());
137	>	speciall->Draw();
138	>	}
139		datahisto->Draw("e1");
140		ckin->Update();
141		mcstack.Draw("same");
#	Line 84 | Line 162 | void draw_kin_variable(string variable,
162
163		datahisto->Delete();
164		delete ckin;
165	+	}
166
167	+	void print_all_b_yields() {
168	+	cout << "Basic selection with a b jet" << endl;
169	+	print_yield(ZplusBsel&&"Zb3010_pfJetGoodNumBtag>0");
170	+	cout << "Leading jet is a b-jet " << endl;
171	+	print_yield(ZplusBsel&&LeadingB);
172	+	cout << "|#eta_{b}|<1.3 " << endl;
173	+	print_yield(ZplusBsel&&LeadingB&&EtaB);
174	+	cout << "|#delta#phi(b<Z)|>2.7" << endl;
175	+	print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
176	+	cout << "10<ptZ<1000 GeV" << endl;
177	+	print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000");
178	+	cout << "#alpha < 0.35" << endl;
179	+	print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.35");
180	+	cout << "#alpha < 0.3" << endl;
181	+	print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.3");
182	+	cout << "#alpha < 0.2" << endl;
183	+	print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.2");
184	+	cout << "#alpha < 0.15" << endl;
185	+	print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.15");
186	+	cout << "#alpha < 0.1" << endl;
187	+	print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.1");
188	+	//  cout << "#alpha < 0.05" << endl;
189	+	//  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.05");
190	+	}
191	+
192	+	void TEST_DrawEvilCutFlow() {
193	+	stringstream MegaCut;
194
195	+	MegaCut << "   1*(" << (const char*) (ZplusBsel&&TCut("Zb3010_pfJetGoodNumBtag>0")) << ")";
196	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB) << ")";
197	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB) << ")";
198	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut) << ")";
199	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000")) << ")";
200	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.3")) << ")";
201	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.2")) << ")";
202	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.15")) << ")";
203	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.1")) << ")";
204	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.05")) << ")";
205	+	MegaCut << " ";
206	+
207	+	TCanvas *can = new TCanvas("can","can");
208	+	can->SetLogy(1);
209	+	TCut basecut("mll>6||mll<10");
210	+
211	+
212	+
213	+	TLegend *leg = allsamples.allbglegend();
214
90	–	write_warning(__FUNCTION__,"Not implemented yet");
91	–	write_warning(__FUNCTION__,"Need to store overall data_over_mc");
92	–	}
215
216	<	Value getfrom2Dmap(TH2F *map, int ixbin, int iybin) {
217	<	Value sum(0,0);
218	<	for(int iy=1;iy<=iybin;iy++) {
219	<	sum=sum+Value(map->GetBinContent(ixbin,iy),map->GetBinError(ixbin,iy));
216	>	TH1F *data  = allsamples.Draw("data",    MegaCut.str(),11,-0.5,10.5, "", "events", basecut,data,luminosity);
217	>	THStack mcm = allsamples.DrawStack("mc", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,mc,luminosity);
218	>
219	>	float runningsum=0;
220	>	for(int i=data->GetNbinsX();i>0;i--) {
221	>	runningsum+=data->GetBinContent(i);
222	>	data->SetBinContent(i,runningsum);
223	>	}
224	>
225	>	float minimum=data->GetMaximum();
226	>
227	>	TH1F* h;
228	>	TIter nextOF(mcm.GetHists());
229	>	while ( h = (TH1F*)nextOF() ) {
230	>	float runningsum=0;
231	>	minimum=data->GetMaximum();//done deliberately at each step so get the minimum of the last (!) contribution
232	>	for(int i=h->GetNbinsX();i>0;i--) {
233	>	runningsum+=h->GetBinContent(i);
234	>	h->SetBinContent(i,runningsum);
235	>	if(runningsum<minimum&&runningsum>0) minimum=runningsum;
236	>	}
237		}
99	–	return sum;
100	–	}
238
239
240	<	void data_mc_agreement_2d() {
240	>
241	>	data->SetMinimum(0.05*minimum);
242	>	can->cd(1)->SetBottomMargin(0.25);
243	>	data->GetXaxis()->SetBinLabel(1,"Pre-selection");
244	>	data->GetXaxis()->SetBinLabel(2,"Contains b jet");
245	>	data->GetXaxis()->SetBinLabel(3,"Leading jet is b-jet");
246	>	data->GetXaxis()->SetBinLabel(4,"|#eta_{b}|<1.3 ");
247	>	data->GetXaxis()->SetBinLabel(5,"|#delta#phi(b,Z)|>2.7");
248	>	data->GetXaxis()->SetBinLabel(6,"10<p_{T}^{Z}<1000 GeV");
249	>	data->GetXaxis()->SetBinLabel(7,"#alpha < 0.3");
250	>	data->GetXaxis()->SetBinLabel(8,"#alpha < 0.2");
251	>	data->GetXaxis()->SetBinLabel(9,"#alpha < 0.15");
252	>	data->GetXaxis()->SetBinLabel(10,"#alpha < 0.1");
253	>	data->GetXaxis()->SetBinLabel(11,"#alpha < 0.05");
254	>	data->GetXaxis()->LabelsOption("v");
255	>
256	>	data->Draw("e1");
257	>	mcm.Draw("same");
258	>	data->Draw("same,e1,axis");
259	>	data->Draw("same,e1");
260	>	//  data->Draw("same,TEXT");
261	>
262	>	leg->Draw();
263	>
264	>	DrawPrelim();
265	>
266	>	CompleteSave(can,"CutFlow");
267	>	}
268	>
269	>	void draw_Zb_kin_vars() {
270	>	draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",25,0,200,"Z p_{T}","Official/Zpt",1);
271	>	draw_kin_variable("Zb3010_pfJetGoodPt[1]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",25,0,200,"Sub-Leading Jet Pt","Official/Jet2Pt",1);
272	>	draw_kin_variable("Zb3010_pfJetGoodPt[0]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",25,0,200,"Leading Jet Pt (B)","Official/LeadingJetPt",1);
273	>	draw_kin_variable("Zb3010_pfJetGoodPt[1]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,2,"p_{T}^{J2} / p_{T}^{Z}","Official/SubLeadingJetPt_Over_ZPt",1);
274	>	draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<50",   40,0,2,"#alpha","Official/alpha_pt_30_to_50",1);
275	>	draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>50&&pt<75",  40,0,2,"#alpha","Official/alpha_50_to_75",1);
276	>	draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>75&&pt<125", 40,0,2,"#alpha","Official/alpha_pt_75_to_125",1);
277	>	draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>125&&pt<1000",40,0,2,"#alpha","Official/alpha_pt_125_to_1000",1);
278	>	draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",40,0,2,"#alpha","Official/alpha_pt_30_to_1000",1);
279	>
280	>	draw_kin_variable("Zb3010_pfBJetDphiZ[0]",ZplusBsel&&LeadingB&&"pt>10&&pt<1000&&pfBJetDphiZ[0]>0",30,0,3.2,"#delta#phi (Z,b lead)","DeltaPhiZBlead",0);
281	>
282	>	/*
283	>	draw_kin_variable("pt1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,200,"p_{T}^{l1}","Lep/pt1",1);
284	>	draw_kin_variable("pt1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id1==0",20,0,200,"p_{T}^{e1}","Lep/pt1_e",1);
285	>	draw_kin_variable("pt1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id1==1",20,0,200,"p_{T}^{#mu1}","Lep/pt1_m",1);
286	>	draw_kin_variable("pt2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,200,"p_{T}^{l2}","Lep/pt2",1);
287	>	draw_kin_variable("pt2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id2==0",20,0,200,"p_{T}^{e2}","Lep/pt2_e",1);
288	>	draw_kin_variable("pt2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id2==1",20,0,200,"p_{T}^{#mu2}","Lep/pt2_m",1);
289	>	draw_kin_variable("eta1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,-3.1,3.1,"#eta^{l1}","Lep/eta1",1);
290	>	draw_kin_variable("eta1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id1==0",20,-3.1,3.1,"#eta^{e1}","Lep/eta1_e",1);
291	>	draw_kin_variable("eta1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id1==1",20,-3.1,3.1,"#eta^{m1}","Lep/eta1_m",1);
292	>	draw_kin_variable("eta2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,-3.1,3.1,"#eta^{l1}","Lep/eta2",1);
293	>	draw_kin_variable("eta2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id2==0",20,-3.1,3.1,"#eta^{e2}","Lep/eta2_e",1);
294	>	draw_kin_variable("eta2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id2==1",20,-3.1,3.1,"#eta^{#mu2}","Lep/eta2_m",1);
295	>	draw_kin_variable("numVtx",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",30,0,30,"N(Vtx)","Lep/nVtx",1);
296	>	*/
297	>	}
298	>
299	>	void draw_mpf_vars() {
300	>	ZbData::data_over_mc.clear();
301	>	draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.3",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p3",0,"#alpha<0.3, 10 GeV < p_{T}^{Z} < 1000 GeV");
302	>	draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.2",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p2",0,"#alpha<0.2, 10 GeV < p_{T}^{Z} < 1000 GeV");
303	>	draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.15",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p15",0,"#alpha<0.15, 10 GeV < p_{T}^{Z} < 1000 GeV");
304	>	draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.1",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p1",0,"#alpha<0.1, 10 GeV < p_{T}^{Z} < 1000 GeV");
305	>	draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.05",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p05",0,"#alpha<0.05, 10 GeV < p_{T}^{Z} < 1000 GeV");
306	>
307	>	draw_kin_variable("Zb3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_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");
308	>	draw_kin_variable("Zb3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_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");
309	>	draw_kin_variable("Zb3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_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");
310	>	draw_kin_variable("Zb3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_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");
311	>	draw_kin_variable("Zb3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.05",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz__mpf_alpha_smaller_0p05",0,"#alpha<0.05, 10 GeV < p_{T}^{Z} < 1000 GeV");
312	>	}
313	>
314	>	/*double GetMedian(TH1F *histo) {
315	>	int numBins = histo->GetXaxis()->GetNbins();
316	>	Double_t x[numBins];
317	>	Double_t y[numBins];
318	>	for (int i = 1; i <= numBins; i++) {
319	>	x[i] = histo->GetBinCenter(i);
320	>	y[i] = histo->GetBinContent(i);
321	>	}
322	>	return TMath::Median(numBins, &x, &y);
323	>	}*/
324	>
325	>	Value get_Zb_data_over_mc(string ContainerName, ZbResultContainer &res, string variable, string variable2, TCut cut, string saveas, float pt1, float pt2, float a1, float a2) {
326	>	//  write_warning(__FUNCTION__,"Debugging this function, therefore always returning 3.1415+/-1.010101"); return Value(3.1415,1.010101);
327	>	TCanvas *cn = new TCanvas("cn","cn");
328	>	string varname="MPF";
329	>	if(Contains(variable,"pfJetGood")) varname="R_{abs}";
330	>	TH1F *hdata  = allsamples.Draw("hdata",variable,1200,0,30, varname, "events",  cut,data,luminosity);
331	>	TH1F *hmc    = allsamples.Draw("hmc" , variable,1200,0,30, varname, "events", cut, mc,  luminosity);
332	>	hmc->SetLineColor(kBlue);
333	>
334	>	float a=hdata->GetMean();
335	>	float b=hmc->GetMean();
336	>	float da=hdata->GetMeanError();
337	>	float db=hmc->GetMeanError();
338	>	float factor = a / b;
339	>	float error = TMath::Sqrt(  (1/(b*b)) * (da*da) + ((a*a)/(b*b))*db*db);
340	>
341	>	TF1 *gaus_data = new TF1("gaus_data","gaus",0.7,1.3);
342	>	gaus_data->SetParameter(0,hdata->GetMaximum());
343	>	gaus_data->SetParameter(2,hdata->GetMean());
344	>	gaus_data->SetParameter(1,hdata->GetRMS());
345	>	hdata->Fit(gaus_data);
346	>	TF1 *gaus_mc = new TF1("gaus_mc","gaus",0.7,1.3);
347	>	gaus_mc->SetParameter(0,hmc->GetMaximum());
348	>	gaus_mc->SetParameter(1,hmc->GetMean());
349	>	gaus_mc->SetParameter(2,hmc->GetRMS());
350	>	hmc->Fit(gaus_mc);
351	>
352	>
353	>	cn->cd();
354	>	hdata->GetYaxis()->SetTitle("events / 0.1");
355	>	hdata->SetMarkerColor(kRed);
356	>	/*  hdata->Rebin(4);
357	>	hmc->Rebin(4);
358	>
359	>	gaus_mc->SetParameter(0,gaus_mc->GetParameter(0)*4);//rebinning
360	>	gaus_data->SetParameter(0,gaus_data->GetParameter(0)*4);//rebinning
361	>	*/
362	>	gaus_mc->SetLineColor(kBlue);
363	>	gaus_data->SetLineColor(kRed);
364	>
365	>	hdata->GetXaxis()->SetRangeUser(0,2);
366	>	hdata->GetYaxis()->SetTitleOffset(1.3);
367	>	hdata->Draw("e1");
368	>	hmc->Draw("histo,same");
369	>	hdata->Draw("e1,same");
370	>	gaus_data->Draw("same");
371	>	gaus_mc->Draw("same");
372	>	stringstream summary;
373	>	summary << "#splitline{#bf{" << varname << "}}{#splitline{" << pt1 << " < p_{T}^{Z} < " << pt2 << ", " << a1 << " #leq #alpha < " << a2 << "}{#splitline{data: " << std::setprecision(4) << a << " #pm " << std::setprecision(4) << da << "}{";
374	>	summary << "#splitline{mc:   " << std::setprecision(4) << b << " #pm " << std::setprecision(4) << db << "}{#splitline{ratio: " << factor;
375	>	summary << " #pm " << error << "}{#splitline{}{#splitline{data fit: " << std::setprecision(4) << gaus_data->GetParameter(1) << "+/-";
376	>	summary << std::setprecision(4) << gaus_data->GetParError(1) << "}{#splitline{";
377	>	summary << "mc fit:" << std::setprecision(4) << gaus_mc->GetParameter(1) << "+/-" << std::setprecision(4) << gaus_mc->GetParError(1) << "}{ratio: ";
378	>	summary << gaus_data->GetParameter(1)/gaus_mc->GetParameter(1) << "+/-" << gaus_data->GetParameter(1)/gaus_mc->GetParameter(1) * sqrt((gaus_data->GetParError(1) * gaus_data->GetParError(1))/(gaus_data->GetParameter(1)*gaus_data->GetParameter(1)) + (gaus_mc->GetParError(1) * gaus_mc->GetParError(1))/(gaus_mc->GetParameter(1)*gaus_mc->GetParameter(1)) ) << "}}}}}}}}";
379	>
380	>	TText *infobox = write_title(summary.str());
381	>	infobox->SetX(0.75);
382	>	infobox->SetTextSize(0.03);
383	>	infobox->SetY(0.75);
384	>	infobox->Draw();
385	>
386	>
387	>	DrawPrelim();
388	>	CompleteSave(cn,ContainerName+"/ResponsePlots/"+saveas);
389	>
390	>	res.MPF_Rabs_location.push_back(ZbPointResult(variable,pt1,pt2,a1,a2,Value(a,da),Value(b,db)));
391	>	cout << "Data;" << a << ";"<<da<<";MC"<<b<<";"<<db<<endl;
392	>	delete cn;
393	>	delete hdata;
394	>	delete hmc;
395	>
396	>	return Value(factor,error);
397	>	}
398
399	<	vector<int> emptyvector;
400	<	float binX[5]={10,50,100,200,1000};
107	<	float binY[5]={0.0,0.1,0.15,0.2,0.3};
108	<	vector<float>binningx;
109	<	vector<float>binningy;
110	<	binningx.assign(binX,binX+5);
111	<	binningy.assign(binY,binY+5);
399	>	ZbResultContainer new_data_mc_agreement_2d(string AlphaVariable, string ContainerName, TCut BTagCut, TCut BTagWgt) {
400	>	ZbResultContainer results(ContainerName);
401
402	<	TCut cut(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
402	>	cutWeight=TCut("(weight*(weight<1000)*(is_data+(!is_data)))")*BTagWgt;
403	>	LeadingB=BTagCut;
404
405	<	TCanvas *agrcan = new TCanvas("agrcan","agrcan");
405	>	gStyle->SetOptFit(0);
406	>	TCut basecut(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
407	>	//  const int nalphacuts=6;
408	>	//  float alphacuts[nalphacuts] = {0.1,0.15,0.2,0.25,0.3,0.35};
409	>
410	>	//  const int nptcuts=5;
411	>	//  float ptcuts[nptcuts]={30,50,75,125,1000};
412	>
413	>	write_info(__FUNCTION__,"JUST TESTING - PLEASE REMOVE THIS LINE");const int nalphacuts=2;float alphacuts[nalphacuts] = {0.1,0.3};const int nptcuts=3;float ptcuts[nptcuts]={30,75,1000};
414	>
415	>
416	>	float MPF_Results[nalphacuts][nptcuts];
417	>	float MPF_Errors[nalphacuts][nptcuts];
418	>	float RABS_Results[nalphacuts][nptcuts];
419	>	float RABS_Errors[nalphacuts][nptcuts];
420	>
421	>
422	>	for(int ia=0;ia<nalphacuts;ia++) {
423	>	for(int ipt=0;ipt<nptcuts-1;ipt++) {
424	>	stringstream specialcut;
425	>	float alow,ahigh;
426	>	if(ia>0) {
427	>	specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (" << AlphaVariable << "<" << alphacuts[ia] << " && " << AlphaVariable << ">" << alphacuts[ia-1] << "))";
428	>	alow=alphacuts[ia-1];
429	>	ahigh=alphacuts[ia];
430	>	} else {
431	>	specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (" << AlphaVariable << "<" << alphacuts[ia] << "))";
432	>	alow=0;
433	>	ahigh=alphacuts[ia];
434	>	}
435	>	cout << specialcut.str() << endl;
436	>	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);
437	>	Value RABS_data_over_mc = get_Zb_data_over_mc(ContainerName,results,"Zb3010_pfJetGoodPt[0]/pt","",TCut(basecut && specialcut.str().c_str()),"Rabs___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__alpha_"+any2string(alphacuts[ia]), ptcuts[ipt], ptcuts[ipt+1],alow,ahigh);
438	>
439	>	MPF_Results[ia][ipt]=MPF_data_over_mc.getValue();
440	>	MPF_Errors[ia][ipt]=MPF_data_over_mc.getError();
441	>
442	>	RABS_Results[ia][ipt]=RABS_data_over_mc.getValue();
443	>	RABS_Errors[ia][ipt]=RABS_data_over_mc.getError();
444	>
445	>	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;
446	>	}
447	>	}
448
449	<	cout << "Drawing 2d for data " << endl;
450	<	TH2F *alphazptplane_data = allsamples.Draw("alphazptplane_data","alpha:pt",binningx,binningy,"Z p_{T}","#alpha",cut,data,luminosity,emptyvector,false);
119	<	cout << "Drawing 2d for MC " << endl;
120	<	TH2F *alphazptplane_mc   = allsamples.Draw("alphazptplane_mc",  "alpha:pt",binningx,binningy,"Z p_{T}","#alpha",cut,mc,  luminosity,emptyvector,false);
449	>	float MPF_ExtraPolatedResults[nptcuts];
450	>	float RABS_ExtraPolatedResults[nptcuts];
451
452	+	TCanvas *can = new TCanvas("can");
453	+	can->cd(1)->SetLeftMargin(0.15);
454
455	<	TFile *f = new TFile("test.root","RECREATE");
456	<	alphazptplane_data->Write();
457	<	f->Close();
455	>	TGraphErrors *MPF_FinalGraph = new TGraphErrors();
456	>	TGraphErrors *RABS_FinalGraph = new TGraphErrors();
457	>	TGraphErrors *MPF_FaceValueAtPoint3 = new TGraphErrors();
458	>	TGraphErrors *RABS_FaceValueAtPoint3 = new TGraphErrors();
459
460	<	for(int ix=1;ix<binningx.size();ix++) {
461	<	cout << "pt in " << binningx[ix-1] << " , " << binningx[ix] << endl;
462	<	TGraphErrors gr;
463	<	for(int iy=1;iy<binningy.size();iy++) {
464	<	Value mc   = Value((getfrom2Dmap(alphazptplane_mc,ix,iy)).getValue(),0.0);// not taking MC statistics into account
465	<	Value data = getfrom2Dmap(alphazptplane_data,ix,iy);
466	<	cout << " For an alpha cut at " << binningy[iy] << " we get : mc (" << mc << ") , data (" << data << ") --> data/mc = " << data/mc << endl;
467	<	gr.SetPoint(iy-1,binningy[iy],(data/mc).getValue());
468	<	gr.SetPointError(iy-1,0,(data/mc).getError());
460	>	for(int ipt=0;ipt<nptcuts-1;ipt++) {
461	>
462	>	stringstream filename;
463	>	filename << "Extrapolation/Zplusb_data_over_mc___" << ptcuts[ipt] << "_to_" << ptcuts[ipt+1];
464	>	cout << "Going to create histo for " << filename.str() << endl;
465	>	TGraphErrors *mpf_gr =new TGraphErrors();
466	>	TGraphErrors *rabs_gr =new TGraphErrors();
467	>
468	>	int rabs_counter=0;
469	>	int mpf_counter=0;
470	>
471	>	for(int ia=0;ia<nalphacuts;ia++) {
472	>	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;
473	>	if(MPF_Results[ia][ipt]==MPF_Results[ia][ipt] && MPF_Errors[ia][ipt]==MPF_Errors[ia][ipt]) { // remove nan's
474	>	mpf_gr->SetPoint(mpf_counter,alphacuts[ia],MPF_Results[ia][ipt]);
475	>	mpf_gr->SetPointError(mpf_counter,0,MPF_Errors[ia][ipt]);
476	>	mpf_counter++;
477	>	} else {
478	>	dout << "Detected BS for MPF method (bin ignored) : " << endl;
479	>	dout << "       alpha: " << alphacuts[ia] << " , pt range: " << ptcuts[ipt] << " < pt < " << ptcuts[ipt+1] << endl;
480	>	cout << "       value (MPF) : " << MPF_Results[ia][ipt] << " +/- " << MPF_Errors[ia][ipt] << endl;
481	>	}
482	>	if(RABS_Results[ia][ipt]==RABS_Results[ia][ipt] && RABS_Errors[ia][ipt]==RABS_Errors[ia][ipt]) {
483	>	rabs_gr->SetPoint(rabs_counter,alphacuts[ia],RABS_Results[ia][ipt]);
484	>	rabs_gr->SetPointError(rabs_counter,0,RABS_Errors[ia][ipt]);
485	>	rabs_counter++;
486	>	} else {
487	>	dout << "Detected BS for RABS method (bin ignored) : " << endl;
488	>	dout << "       alpha: " << alphacuts[ia] << " , pt range: " << ptcuts[ipt] << " < pt < " << ptcuts[ipt+1] << endl;
489	>	cout << "       value (RABS) : " << RABS_Results[ia][ipt] << " +/- " << RABS_Errors[ia][ipt] << endl;
490	>	}
491		}
492	+
493	+	stringstream specialcut;
494	+	specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (" << AlphaVariable << "<0.3))";
495	+	Value MPF_data_over_mc = get_Zb_data_over_mc(ContainerName,results,"mpf","",TCut(basecut && specialcut.str().c_str()),"MPF___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__INCLUSIVEalpha_0.3",ptcuts[ipt],ptcuts[ipt+1],0,0.3);
496	+	Value RABS_data_over_mc = get_Zb_data_over_mc(ContainerName,results,"Zb3010_pfJetGoodPt[0]/pt","",TCut(basecut && specialcut.str().c_str()),"Rabs___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__INCLUSIVEalpha_0.3",ptcuts[ipt],ptcuts[ipt+1],0,0.3);
497	+	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;
498	+	MPF_FaceValueAtPoint3->SetPoint(ipt,(ptcuts[ipt]+ptcuts[ipt+1])/2,MPF_data_over_mc.getValue());
499	+	MPF_FaceValueAtPoint3->SetPointError(ipt,0,MPF_data_over_mc.getError());
500	+	RABS_FaceValueAtPoint3->SetPoint(ipt,(ptcuts[ipt]+ptcuts[ipt+1])/2,RABS_data_over_mc.getValue());
501	+	RABS_FaceValueAtPoint3->SetPointError(ipt,0,RABS_data_over_mc.getError());
502	+
503	+	can->cd();
504	+	mpf_gr->SetMarkerStyle(21);
505	+	mpf_gr->SetMarkerColor(kBlue);
506	+	mpf_gr->Draw("AP");
507	+	mpf_gr->Fit("pol1","");
508	+	mpf_gr->GetXaxis()->SetTitle("#alpha");
509	+	mpf_gr->GetXaxis()->CenterTitle();
510	+	mpf_gr->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
511	+	mpf_gr->GetYaxis()->CenterTitle();
512	+	TF1 *mpf_pol=(TF1*)mpf_gr->GetFunction("pol1");
513	+	mpf_pol->SetLineWidth(0);
514	+	TF1 *mpf_pol_complete = new TF1("mpf_pol_complete","[0]+[1]*x");
515	+	mpf_pol_complete->SetParameters(mpf_pol->GetParameters());
516	+	mpf_pol_complete->SetLineWidth(1);
517	+	mpf_pol_complete->SetLineColor(kBlue);
518	+	float min,max;
519	+	FindMinMax(mpf_gr,min,max);
520	+	TH1F *histo = new TH1F("histo","histo",1,0,0.4);
521	+	histo->GetXaxis()->SetTitle("#alpha");
522	+	histo->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
523	+	histo->GetXaxis()->CenterTitle();
524	+	histo->GetYaxis()->CenterTitle();
525	+	histo->GetYaxis()->SetTitleOffset(1.3);
526	+	histo->SetStats(0);
527	+	histo->GetXaxis()->SetRangeUser(0,0.4);
528	+	histo->GetYaxis()->SetRangeUser(min,max);
529	+	mpf_gr->GetYaxis()->SetRangeUser(min,max);
530	+
531	+	TPolyLine *mpf_fit_uncert = GetFitUncertaintyShape(mpf_gr, "pol1", min, max,0.0,0.4);
532	+	mpf_pol->SetLineColor(TColor::GetColor("#0101DF"));
533	+	mpf_fit_uncert->SetFillColor(TColor::GetColor("#5353FC"));
534	+	histo->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
535	+	histo->Draw();
536	+	mpf_fit_uncert->Draw("F");
537	+	mpf_fit_uncert->Draw("");
538	+	mpf_gr->Draw("P");
539	+	mpf_pol_complete->Draw("same");
540	+
541	+	float mpf_a=mpf_pol->GetParameter(0);
542	+	float mpf_b=mpf_pol->GetParameter(1);
543	+	MPF_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),mpf_a);
544	+	MPF_FinalGraph->SetPointError(ipt,0,mpf_pol->GetParError(0));
545	+
546	+	MPF_ExtraPolatedResults[ipt]=mpf_a;
547	+
548	+	stringstream mpf_info;
549	+	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() << "}";
550	+
551	+	TText *mpf_mark = write_title(mpf_info.str());
552	+	mpf_mark->SetX(0.75);
553	+	mpf_mark->SetTextSize(0.03);
554	+	mpf_mark->SetY(0.75);
555	+	mpf_mark->Draw();
556	+
557	+	string filenamebkp=filename.str();
558	+	filename << "__MPF";
559	+	DrawPrelim();
560	+	CompleteSave(can,ContainerName+"/"+filename.str());
561	+	cout << "MPF : " << mpf_a << " + " << mpf_b << " * alpha " << endl;
562	+
563	+	filename.str("");
564	+
565	+	rabs_gr->SetMarkerStyle(21);
566	+	rabs_gr->SetMarkerColor(kRed);
567	+	rabs_gr->Draw("AP");
568	+	rabs_gr->Fit("pol1","");
569	+	rabs_gr->GetXaxis()->SetTitle("#alpha");
570	+	rabs_gr->GetXaxis()->CenterTitle();
571	+	rabs_gr->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
572	+	rabs_gr->GetYaxis()->CenterTitle();
573	+	TF1 *rabs_pol=(TF1*)rabs_gr->GetFunction("pol1");
574	+	rabs_pol->SetLineWidth(0);
575	+	TF1 *rabs_pol_complete = new TF1("rabs_pol_complete","[0]+[1]*x");
576	+	rabs_pol_complete->SetParameters(rabs_pol->GetParameters());
577	+	rabs_pol_complete->SetLineColor(kRed);
578	+	rabs_pol_complete->SetLineWidth(1);
579	+	FindMinMax(rabs_gr,min,max);
580	+	rabs_gr->GetYaxis()->SetRangeUser(min,max);
581	+	histo->GetYaxis()->SetRangeUser(min,max);
582	+	TPolyLine *rabs_fit_uncert = GetFitUncertaintyShape(rabs_gr, "pol1", min, max,0.0,0.4);
583	+	rabs_pol->SetLineColor(TColor::GetColor("#FA5858"));
584	+	rabs_pol->SetFillColor(TColor::GetColor("#FA5858"));
585	+	rabs_fit_uncert->SetLineColor(TColor::GetColor("#FA5858"));
586	+	rabs_fit_uncert->SetFillColor(TColor::GetColor("#FA5858"));
587	+	histo->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
588	+	histo->Draw();
589	+	rabs_fit_uncert->Draw("F");
590	+	rabs_fit_uncert->Draw("");
591	+	rabs_gr->Draw("P");
592	+	rabs_pol_complete->Draw("same");
593
594	<	gr.Draw("AP*");
595	<	gr.GetYaxis()->SetTitle("data / mc");
596	<	gr.GetYaxis()->CenterTitle();
597	<	gr.GetXaxis()->SetTitle("#alpha");
598	<	gr.GetXaxis()->CenterTitle();
599	<	gr.GetXaxis()->SetRangeUser(0,0.4);
600	<	gr.GetYaxis()->SetRangeUser(0.5,1.5);
601	<	gr.Draw("AP*");
594	>
595	>	float rabs_a=rabs_pol->GetParameter(0);
596	>	float rabs_b=rabs_pol->GetParameter(1);
597	>	RABS_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),rabs_a);
598	>	RABS_FinalGraph->SetPointError(ipt,0,rabs_pol->GetParError(0));
599	>
600	>	cout << "!+!+!+!+!+!+!+!+!+ Just added a point to the final plot : " << rabs_a << " +/- " << rabs_pol->GetParError(0) << endl;
601	>
602	>	stringstream rabs_info;
603	>	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() << "}";
604	>	TText *rabs_mark = write_title(rabs_info.str());
605	>	rabs_mark->SetX(0.75);
606	>	rabs_mark->SetTextSize(0.03);
607	>	rabs_mark->SetY(0.75);
608	>	rabs_mark->Draw();
609	>
610	>
611	>	RABS_ExtraPolatedResults[ipt]=rabs_a;
612	>	filename << filenamebkp << "__RABS";
613		DrawPrelim();
614	<	stringstream ptbin;
615	<	ptbin << binningx[ix-1] << " < p_{T}^{Z} < " << binningx[ix] << " GeV";
149	<	TText *ptt = write_title(ptbin.str());
150	<	ptt->SetX(0.7);
151	<	ptt->SetY(0.7);
152	<	ptt->SetTextSize(0.03);
153	<	ptt->Draw();
154	<	TLine *l = new TLine(0.0,1.0,0.4,1.0);
155	<	l->SetLineStyle(2);
156	<	l->SetLineColor(kBlue);
157	<	l->Draw();
158	<
159	<	stringstream saveas;
160	<	saveas << "Zplusb_data_over_mc___" << binningx[ix-1] << "to" << binningx[ix];
161	<	CompleteSave(agrcan,saveas.str());
162	<	//    agrcan->SaveAs("tester.png");
614	>	CompleteSave(can,ContainerName+"/"+filename.str());
615	>	cout << "RABS : " << rabs_a << " + " << rabs_b << " * alpha " << endl;
616		}
617
618	<	delete agrcan;
618	>	MPF_FinalGraph->SetMarkerStyle(21);
619	>	MPF_FinalGraph->SetMarkerColor(kBlue);
620	>	MPF_FinalGraph->Fit("pol0",""); // quiet, use minos
621	>
622	>	TF1 *mpf_pol0=(TF1*)MPF_FinalGraph->GetFunction("pol0");
623	>	TF1 *mpf_pol0_complete = new TF1("mpf_pol0_complete","[0]+[1]*x");
624	>	mpf_pol0_complete->SetLineWidth(1);
625	>	mpf_pol0_complete->SetLineColor(kBlue);
626	>	mpf_pol0->SetLineColor(kBlue);
627	>	MPF_FinalGraph->Draw("AP");
628	>	MPF_FinalGraph->GetXaxis()->SetTitle("p_{T}^{Z}");
629	>	MPF_FinalGraph->GetXaxis()->CenterTitle();
630	>	MPF_FinalGraph->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
631	>	MPF_FinalGraph->GetYaxis()->CenterTitle();
632	>	MPF_FinalGraph->Draw("AP");
633	>	mpf_pol0_complete->Draw("same");
634	>
635	>	stringstream mpf_result;
636	>	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)}";
637	>	TText *rmpf_final_mark = write_title(mpf_result.str());
638	>	rmpf_final_mark->SetX(0.75);
639	>	rmpf_final_mark->SetY(0.75);
640	>	rmpf_final_mark->SetTextSize(0.03);
641	>	rmpf_final_mark->Draw();
642	>	DrawPrelim();
643	>	float MPFResult=1/mpf_pol0->GetParameter(0);
644	>	float MPFResultError=mpf_pol0->GetParError(0)/(mpf_pol0->GetParameter(0)*mpf_pol0->GetParameter(0));
645	>
646	>	stringstream filename2;
647	>	filename2 << "Extrapolation/FINAL_RESULT_MPF";
648	>	CompleteSave(can,ContainerName+"/"+filename2.str());
649	>
650	>
651	>	RABS_FinalGraph->SetMarkerStyle(21);
652	>	RABS_FinalGraph->SetMarkerStyle(21);
653	>	RABS_FinalGraph->SetMarkerColor(kRed);
654	>	RABS_FinalGraph->Fit("pol0","QE"); // quiet, use minos
655	>	RABS_FinalGraph->Draw("AP");
656	>	RABS_FinalGraph->GetXaxis()->SetTitle("p_{T}^{Z}");
657	>	RABS_FinalGraph->GetXaxis()->CenterTitle();
658	>	RABS_FinalGraph->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
659	>	RABS_FinalGraph->GetYaxis()->CenterTitle();
660	>	RABS_FinalGraph->Draw("AP");
661	>	TF1 *rabs_pol0=(TF1*)RABS_FinalGraph->GetFunction("pol0");
662	>	stringstream rabs_result;
663	>	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})}";;
664	>	TText *rabs_final_mark = write_title(rabs_result.str());
665	>	rabs_final_mark->SetX(0.75);
666	>	rabs_final_mark->SetY(0.75);
667	>	rabs_final_mark->SetTextSize(0.03);
668	>	rabs_final_mark->Draw();
669	>	DrawPrelim();
670	>
671	>	float RabsResult=1/rabs_pol0->GetParameter(0);
672	>	float RabsResultError=rabs_pol0->GetParError(0)/(rabs_pol0->GetParameter(0)*rabs_pol0->GetParameter(0));
673	>
674	>	filename2.str("");
675	>	filename2 << "Extrapolation/FINAL_RESULT_RABS";
676	>	CompleteSave(can,ContainerName+"/"+filename2.str());
677	>
678	>	can->cd();
679	>	MPF_FaceValueAtPoint3->SetMarkerStyle(21);
680	>	MPF_FaceValueAtPoint3->SetMarkerColor(kBlue);
681	>	MPF_FaceValueAtPoint3->Fit("pol0","QE"); // quiet, use minos
682	>	MPF_FaceValueAtPoint3->Draw("AP");
683	>	MPF_FaceValueAtPoint3->GetXaxis()->SetTitle("p_{T}^{Z}");
684	>	MPF_FaceValueAtPoint3->GetXaxis()->CenterTitle();
685	>	MPF_FaceValueAtPoint3->GetYaxis()->SetTitle("<MPF>__{data}/<MPF>_{mc} for #alpha<0.3");
686	>	MPF_FaceValueAtPoint3->GetYaxis()->CenterTitle();
687	>	MPF_FaceValueAtPoint3->Draw("AP");
688	>	TF1 *faceval_pol0=(TF1*)MPF_FaceValueAtPoint3->GetFunction("pol0");
689	>	faceval_pol0->SetLineColor(kBlue);
690	>	faceval_pol0->SetLineWidth(1);
691	>	faceval_pol0->Draw("same");
692	>	TF1 *faceval_pol0_complete = new TF1("faceval_pol0_complete","[0]+[1]*x");
693	>	faceval_pol0_complete->SetParameters(faceval_pol0->GetParameters());
694	>	stringstream faceval_result;
695	>	faceval_result << "#splitline{C_{abs}= " << std::setprecision(5) << 1/faceval_pol0->GetParameter(0) << " #pm " << std::setprecision(2) << faceval_pol0->GetParError(0)/(faceval_pol0->GetParameter(0)*faceval_pol0->GetParameter(0)) << "}{ (MPF at #alpha<0.3)}";;
696	>	TText *faceval_final_mark = write_title(faceval_result.str());
697	>	faceval_final_mark->SetX(0.75);
698	>	faceval_final_mark->SetY(0.75);
699	>	faceval_final_mark->SetTextSize(0.03);
700	>	faceval_final_mark->Draw();
701	>	DrawPrelim();
702	>
703	>	float FaceValResult=1/faceval_pol0->GetParameter(0);
704	>	float FaceValResultError=faceval_pol0->GetParError(0)/(faceval_pol0->GetParameter(0)*faceval_pol0->GetParameter(0));
705	>
706	>	filename2.str("");
707	>	filename2 << "Extrapolation/FINAL_RESULT_MPF_FaceValp3";
708	>	CompleteSave(can,ContainerName+"/"+filename2.str());
709	>
710	>	can->cd();
711	>	RABS_FaceValueAtPoint3->SetMarkerStyle(21);
712	>	RABS_FaceValueAtPoint3->SetMarkerColor(kRed);
713	>	RABS_FaceValueAtPoint3->Fit("pol0","QE"); // quiet, use minos
714	>	RABS_FaceValueAtPoint3->Draw("AP");
715	>	RABS_FaceValueAtPoint3->GetXaxis()->SetTitle("p_{T}^{Z}");
716	>	RABS_FaceValueAtPoint3->GetXaxis()->CenterTitle();
717	>	RABS_FaceValueAtPoint3->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc} for #alpha<0.3");
718	>	RABS_FaceValueAtPoint3->GetYaxis()->CenterTitle();
719	>	RABS_FaceValueAtPoint3->Draw("AP");
720	>	TF1 *faceval_pol0RABS=(TF1*)RABS_FaceValueAtPoint3->GetFunction("pol0");
721	>	TF1 *faceval_pol0RABS_complete = new TF1("faceval_pol0_complete","[0]+[1]*x");
722	>	faceval_pol0RABS_complete->SetParameters(faceval_pol0RABS->GetParameters());
723	>	faceval_pol0RABS->SetLineWidth(1);
724	>	faceval_pol0RABS->SetLineColor(kRed);
725	>	faceval_pol0RABS_complete->Draw("same");
726	>	stringstream RABSfaceval_result;
727	>	RABSfaceval_result << "#splitline{C_{abs}= " << std::setprecision(5) << 1/faceval_pol0RABS->GetParameter(0) << " #pm " << std::setprecision(2) << faceval_pol0RABS->GetParError(0)/(faceval_pol0RABS->GetParameter(0)*faceval_pol0RABS->GetParameter(0)) << "}{ (R_{abs} at #alpha<0.3)}";;
728	>	TText *RABSfaceval_final_mark = write_title(RABSfaceval_result.str());
729	>	RABSfaceval_final_mark->SetX(0.75);
730	>	RABSfaceval_final_mark->SetY(0.75);
731	>	RABSfaceval_final_mark->SetTextSize(0.03);
732	>	RABSfaceval_final_mark->Draw();
733	>	DrawPrelim();
734
735	+	float RABSFaceValResult=1/faceval_pol0RABS->GetParameter(0);
736	+	float RABSFaceValResultError=faceval_pol0RABS->GetParError(0)/(faceval_pol0RABS->GetParameter(0)*faceval_pol0RABS->GetParameter(0));
737	+
738	+	filename2.str("");
739	+	filename2 << "Extrapolation/FINAL_RESULT_RABS_FaceValp3";
740	+	CompleteSave(can,ContainerName+"/"+filename2.str());
741	+
742	+	cout << "FINAL RESULTS: " << endl;
743	+	cout << "MPF : " << MPFResult  << " +/- " << MPFResultError  << endl;
744	+	cout << "Rabs: " << RabsResult << " +/- " << RabsResultError << endl;
745	+	cout << "Face value at 0.3 :" << FaceValResult << " +/- " << FaceValResultError << "  (MPF) " << endl;
746	+	cout << "Face value at 0.3 :" << RABSFaceValResult << " +/- " << RABSFaceValResultError << "  (RABS) " << endl;
747	+
748	+	delete can;
749	+
750	+	results.MPF_Result_FaceValue=Value(FaceValResult,FaceValResultError);
751	+	results.Rabs_Result_FaceValue=Value(RABSFaceValResult,RABSFaceValResultError);
752	+	results.MPF_Result_Extrapolated=Value(MPFResult,MPFResultError);
753	+	results.Rabs_Result_Extrapolated=Value(RabsResult,RabsResultError);
754	+	cout << results << endl;
755	+	return results;
756		}
757
169	–	void print_all_b_yields() {
170	–	cout << "Basic selection with a b jet" << endl;
171	–	print_yield(ZplusBsel&&"pfJetGoodNumBtag>0");
172	–	cout << "Leading jet is a b " << endl;
173	–	print_yield(ZplusBsel&&LeadingB);
174	–	cout << "|#eta_{b}|<1.3 " << endl;
175	–	print_yield(ZplusBsel&&LeadingB&&EtaB);
176	–	cout << "|#delta#phi(b<Z)|>2.7" << endl;
177	–	print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
178	–	cout << "10<ptZ<100 GeV" << endl;
179	–	print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000");
180	–	cout << "#alpha < 0.3" << endl;
181	–	print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.3");
182	–	cout << "#alpha < 0.2" << endl;
183	–	print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.2");
184	–	cout << "#alpha < 0.15" << endl;
185	–	print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.15");
186	–	cout << "#alpha < 0.1" << endl;
187	–	print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.1");
188	–	cout << "#alpha < 0.05" << endl;
189	–	print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.05");
190	–	}
758
759	<	void draw_Zb_kin_vars() {
760	<	draw_kin_variable("pfBJetDphiZ",ZplusBsel&&LeadingB,100,0,3.2,"#delta#phi (Z,b lead)","DeltaPhiZBlead",1);
761	<	draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",50,0,200,"Z p_{T}","Zpt",1);
762	<	draw_kin_variable("pfJetGoodPt[1]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",50,0,200,"Sub-Leading Jet Pt","Jet2Pt",1);
763	<	draw_kin_variable("pfJetGoodPtlBtag[0]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",50,0,200,"Leading Jet Pt (B)","LeadingJetPt",1);
764	<	draw_kin_variable("alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,2,"#alpha","alpha",1);
759	>	void TEST_DoPUStudy(string identifier) {
760	>
761	>	float numVtxcuts[7]={0,10,15,20};
762	>
763	>	TCanvas *can = new TCanvas("can","can");
764	>	TH1F *malpha[8];
765	>	TH1F *dalpha[8];
766	>	cout << identifier << ";";
767	>	cout << endl;
768	>	for(int i=1;i<5;i++) {
769	>	stringstream sSpecialCut;
770	>	if(i<4) {
771	>	sSpecialCut << "numVtx>" << numVtxcuts[i-1] << " && numVtx<" << numVtxcuts[i];
772	>	cout << numVtxcuts[i-1] << "  < nVtx < " << numVtxcuts[i] << ";";
773	>	} else {
774	>	sSpecialCut << "numVtx>" << numVtxcuts[i-1];
775	>	cout << numVtxcuts[i-1] << ";";
776	>	}
777	>
778	>	TCut SpecialCut(sSpecialCut.str().c_str());
779	>
780	>
781	>	malpha[i] = allsamples.Draw("malpha"+any2string(i),"mpf",1000,0,10,"MPF","events",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&SpecialCut,mc,luminosity);
782	>	dalpha[i] = allsamples.Draw("dalpha"+any2string(i),"mpf",1000,0,10,"MPF","events",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&SpecialCut,data,luminosity);
783	>
784	>
785	>	float a=dalpha[i]->GetMean();
786	>	float b=malpha[i]->GetMean();
787	>	float da=dalpha[i]->GetMeanError();
788	>	float db=malpha[i]->GetMeanError();
789	>	float factor = a / b;
790	>	float error = TMath::Sqrt(  (1/(b*b)) * (da*da) + ((a*a)/(b*b))*db*db);
791	>
792	>	cout << malpha[i]->GetMean() << ";" << malpha[i]->GetMeanError() << ";" << dalpha[i]->GetMean() << ";" << dalpha[i]->GetMeanError() << ";" << malpha[i]->Integral()<<";"<<dalpha[i]->Integral() << ";"<<factor << ";" << error <<  endl;
793	>	malpha[i]->SetLineColor(i+1);
794	>	if(i==0) malpha[i]->Draw("histo");
795	>	else malpha[i]->Draw("histo,same");
796	>	}
797	>
798	>	CompleteSave(can,"AlphaOverview_"+identifier);
799	>	delete can;
800	>
801	>	/*  for(int i=1;i<8;i++) {
802	>	delete malpha[i];
803	>	delete dalpha[i];
804	>	}
805	>	*/
806	>
807		}
808
809	<	void draw_mpf_vars() {
810	<	ZbData::data_over_mc.clear();
811	<	draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.3",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p3",0);
812	<	draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.2",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p2",0);
813	<	draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.15",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p15",0);
814	<	draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.1",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p1",0);
815	<	draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.05",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p05",0);
816	<
817	<	draw_kin_variable("pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.3",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p3",0);
818	<	draw_kin_variable("pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.2",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p2",0);
819	<	draw_kin_variable("pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.15",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p15",0);
820	<	draw_kin_variable("pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.1",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p1",0);
821	<	draw_kin_variable("pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.05",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz__mpf_alpha_smaller_0p05",0);
822	<	}
809	>	void TEST_ScenarioComparison() {
810	>	//very dumb way to do this.
811	>	TGraphAsymmErrors *gr[5];
812	>	TF1 *fit[5];
813	>	bool dofit=false;
814	>
815	>	TCanvas *can = new TCanvas("can","can",500,500);
816	>
817	>	TLegend *leg = new TLegend(0.2,0.2,0.5,0.45);
818	>	leg->SetBorderSize(0);
819	>	leg->SetFillColor(kWhite);
820	>
821	>	for(int i=0;i<5;i++) {
822	>	gr[i] = new TGraphAsymmErrors();
823	>	double x[4]   = {5,12.5,17.5,20};
824	>	double dxl[4] = {5,2.5,2.5,2.5};
825	>	double dxh[4] = {5,2.5,2.5,10};
826	>
827	>	string name="";
828	>	string color="";
829	>
830	>	if(i==0) {
831	>	double y[4]   = {1.008,1.013,1.003,1.007};
832	>	double dyl[4] = {0.011,0.009,0.013,0.020};
833	>	double dyh[4] = {0.011,0.009,0.013,0.020};
834	>	name="30/30";
835	>	gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
836	>	color="#a01d99";
837	>	}
838	>	if(i==1) {
839	>	double y[4]   = {1.008+0.001,1.013+0.001,1.003+0.001,1.007+0.001};
840	>	double dyl[4] = {0.011,0.009,0.013,0.020};
841	>	double dyh[4] = {0.011,0.009,0.013,0.020};
842	>	gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
843	>	name="30/15";
844	>	color="#2464bc";
845	>	}
846	>	if(i==2) {
847	>	double y[4]   = {1.014,1.026,1.014,1.059};
848	>	double dyl[4] = {0.011,0.010,0.014,0.023};
849	>	double dyh[4] = {0.011,0.010,0.014,0.023};
850	>	gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
851	>	name="15/15";
852	>	color="#f49230";
853	>	}
854	>	if(i==3) {
855	>	double y[4]   = {1.015,1.020,1.008,1.050};
856	>	double dyl[4] = {0.011,0.010,0.014,0.022};
857	>	double dyh[4] = {0.011,0.010,0.014,0.022};
858	>	gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
859	>	name="20/15";
860	>	color="#c72f3c";
861	>	}
862	>	if(i==4) {
863	>	double y[4]   = {1.015+0.001,1.020+0.001,1.008+0.001,1.050+0.001};
864	>	double dyl[4] = {0.011,0.010,0.014,0.022};
865	>	double dyh[4] = {0.011,0.010,0.014,0.022};
866	>	gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
867	>	name="20/20";
868	>	color="#7cb433";
869	>	}
870	>
871	>
872	>	gr[i]->SetTitle();
873	>	gr[i]->GetXaxis()->SetTitle("N(Vtx)");
874	>	gr[i]->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
875	>	gr[i]->GetXaxis()->CenterTitle();
876	>	gr[i]->GetYaxis()->CenterTitle();
877	>	gr[i]->SetName(name.c_str());
878	>	gr[i]->SetMarkerSize(0.00001);
879	>	gr[i]->GetYaxis()->SetRangeUser(0.90,1.10);
880	>	gr[i]->SetLineColor(TColor::GetColor(color.c_str()));
881	>	gr[i]->SetMarkerColor(TColor::GetColor(color.c_str()));
882	>	if(dofit) {
883	>	gr[i]->Fit("pol1");
884	>	fit[i] = (TF1*)gr[i]->GetFunction("pol1");
885	>	fit[i]->SetLineColor(TColor::GetColor(color.c_str()));
886	>	}
887	>	leg->AddEntry(gr[i],name.c_str(),"l");
888	>	if(i==0) gr[i]->Draw("AP*");
889	>	else gr[i]->Draw("P");
890	>	}
891	>	leg->Draw();
892	>	DrawPrelim();
893	>
894	>	CompleteSave(can,"ScenarioComparison");
895	>	delete can;
896	>
897	>	}
898	>
899	>	void TEST_ScenarioComparisonInclusive() {
900	>	//dumbest way ever to do this. but ok we only need to do it once.
901	>	TGraphAsymmErrors *gr[5];
902	>	TF1 *fit[5];
903	>	bool dofit=true;
904	>
905	>	TCanvas *can = new TCanvas("can","can",500,500);
906	>
907	>	TLegend *leg = new TLegend(0.2,0.2,0.5,0.45);
908	>	leg->SetBorderSize(0);
909	>	leg->SetFillColor(kWhite);
910	>
911	>	for(int i=0;i<5;i++) {
912	>	gr[i] = new TGraphAsymmErrors();
913	>	double x[4]   = {5,12.5,17.5,20};
914	>	double dxl[4] = {5,2.5,2.5,2.5};
915	>	double dxh[4] = {5,2.5,2.5,10};
916	>
917	>	string name="";
918	>	string color="";
919	>
920	>	if(i==0) {
921	>	double y[4]   = {1.01103,1.01496,1.01486,1.02785};
922	>	double dyl[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
923	>	double dyh[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
924	>	name="30/30";
925	>	gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
926	>	color="#a01d99";
927	>	}
928	>	if(i==1) {
929	>	double y[4]   = {1.01103+0.001,1.01496+0.001,1.01486+0.001,1.02785+0.001};
930	>	double dyl[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
931	>	double dyh[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
932	>	gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
933	>	name="30/15";
934	>	color="#2464bc";
935	>	}
936	>	if(i==2) {
937	>	double y[4]   = {1.01575,1.02879,1.0325,1.05456};
938	>	double dyl[4] = {0.00288281,0.00278536,0.00367451,0.00590421};
939	>	double dyh[4] = {0.00288281,0.00278536,0.00367451,0.00590421};
940	>	gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
941	>	name="15/15";
942	>	color="#f49230";
943	>	}
944	>	if(i==3) {
945	>	double y[4]   = {1.01221,1.01922,1.02396,1.04677};
946	>	double dyl[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
947	>	double dyh[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
948	>	gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
949	>	name="20/15";
950	>	color="#c72f3c";
951	>	}
952	>	if(i==4) {
953	>	double y[4]   = {1.01221+0.001,1.01922+0.001,1.02396+0.001,1.04677+0.001};
954	>	double dyl[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
955	>	double dyh[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
956	>	gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
957	>	name="20/20";
958	>	color="#7cb433";
959	>	}
960	>
961	>
962	>	gr[i]->SetTitle();
963	>	gr[i]->GetXaxis()->SetTitle("N(Vtx)");
964	>	gr[i]->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
965	>	gr[i]->GetXaxis()->CenterTitle();
966	>	gr[i]->GetYaxis()->CenterTitle();
967	>	gr[i]->SetName(name.c_str());
968	>	gr[i]->SetMarkerSize(0.00001);
969	>	gr[i]->GetYaxis()->SetRangeUser(0.90,1.10);
970	>	gr[i]->SetLineColor(TColor::GetColor(color.c_str()));
971	>	gr[i]->SetMarkerColor(TColor::GetColor(color.c_str()));
972	>	if(dofit) {
973	>	gr[i]->Fit("pol1");
974	>	fit[i] = (TF1*)gr[i]->GetFunction("pol1");
975	>	fit[i]->SetLineColor(TColor::GetColor(color.c_str()));
976	>	}
977	>	leg->AddEntry(gr[i],name.c_str(),"l");
978	>	if(i==0) gr[i]->Draw("AP*");
979	>	else gr[i]->Draw("P");
980	>	}
981	>	leg->Draw();
982	>	DrawPrelim();
983	>
984	>	CompleteSave(can,"ScenarioComparisonInclusive");
985	>	delete can;
986	>
987	>	}
988	>
989	>	void TEST_compare_selection(string identifier) {
990	>	bool recognized_scenario=false;
991	>
992	>	cout << "Running with identifier " << identifier << endl;
993	>	if(identifier=="Zb3010") {
994	>	LeadingB=TCut ("Zb3010_bTagProbCSVBP[0]>0.898");
995	>	EtaB=TCut("abs(Zb3010_pfJetGoodEta[0])<1.3");
996	>	PhiZcut=TCut("abs(Zb3010_pfJetDphiZ[0])>2.7");
997	>	recognized_scenario=true;
998	>	}
999	>
1000	>	if(identifier=="Zb3010") {
1001	>	LeadingB=TCut ("Zb3010_bTagProbCSVBP[0]>0.898");
1002	>	EtaB=TCut("abs(Zb3010_pfJetGoodEta[0])<1.3");
1003	>	PhiZcut=TCut("abs(Zb3010_pfJetDphiZ[0])>2.7");
1004	>	recognized_scenario=true;
1005	>	}
1006	>
1007	>	if(identifier=="Zb40") {
1008	>	LeadingB=TCut ("Zb40_bTagProbCSVBP[0]>0.898");
1009	>	EtaB=TCut("abs(Zb40_pfJetGoodEta[0])<1.3");
1010	>	PhiZcut=TCut("abs(Zb40_pfJetDphiZ[0])>2.7");
1011	>	recognized_scenario=true;
1012	>	}
1013	>
1014	>	if(identifier=="Zb301030") {
1015	>	LeadingB=TCut ("Zb301030_bTagProbCSVBP[0]>0.898");
1016	>	EtaB=TCut("abs(Zb301030_pfJetGoodEta[0])<1.3");
1017	>	PhiZcut=TCut("abs(Zb301030_pfJetDphiZ[0])>2.7");
1018	>	recognized_scenario=true;
1019	>	}
1020	>
1021	>	if(identifier=="Zb1530") {
1022	>	LeadingB=TCut ("Zb1530_bTagProbCSVBP[0]>0.898");
1023	>	EtaB=TCut("abs(Zb1530_pfJetGoodEta[0])<1.3");
1024	>	PhiZcut=TCut("abs(Zb1530_pfJetDphiZ[0])>2.7");
1025	>	recognized_scenario=true;
1026	>	}
1027	>
1028	>	if(identifier=="Zb301010") {
1029	>	LeadingB=TCut ("Zb301010_bTagProbCSVBP[0]>0.898");
1030	>	EtaB=TCut("abs(Zb301010_pfJetGoodEta[0])<1.3");
1031	>	PhiZcut=TCut("abs(Zb301010_pfJetDphiZ[0])>2.7");
1032	>	recognized_scenario=true;
1033	>	}
1034	>
1035	>	if(identifier=="Zb1510") {
1036	>	LeadingB=TCut ("Zb1510_bTagProbCSVBP[0]>0.898");
1037	>	EtaB=TCut("abs(Zb1510_pfJetGoodEta[0])<1.3");
1038	>	PhiZcut=TCut("abs(Zb1510_pfJetDphiZ[0])>2.7");
1039	>	recognized_scenario=true;
1040	>	}
1041	>
1042	>	if(identifier=="Zb301010") {
1043	>	LeadingB=TCut ("Zb301010_bTagProbCSVBP[0]>0.898");
1044	>	EtaB=TCut("abs(Zb301010_pfJetGoodEta[0])<1.3");
1045	>	PhiZcut=TCut("abs(Zb301010_pfJetDphiZ[0])>2.7");
1046	>	recognized_scenario=true;
1047	>	}
1048	>
1049	>	if(identifier=="ZbMikko") {
1050	>	LeadingB=TCut ("ZbMikko_bTagProbCSVBP[0]>0.898 && ZbMikko__IsClean");
1051	>	EtaB=TCut("abs(ZbMikko_pfJetGoodEta[0])<1.3 && ZbMikko__IsClean");
1052	>	PhiZcut=TCut("abs(ZbMikko_pfJetDphiZ[0])>2.7 && ZbMikko__IsClean");
1053	>	recognized_scenario=true;
1054	>	}
1055	>
1056	>	if(identifier=="Zb3010_SecEta3") {
1057	>	LeadingB=TCut ("Zb3010_SecEta3_bTagProbCSVBP[0]>0.898");
1058	>	EtaB=TCut("abs(Zb3010_SecEta3_pfJetGoodEta[0])<1.3");
1059	>	PhiZcut=TCut("abs(Zb3010_SecEta3_pfJetDphiZ[0])>2.7");
1060	>	recognized_scenario=true;
1061	>	}
1062	>
1063	>	if(identifier=="Zb3010_SecEta5") {
1064	>	LeadingB=TCut ("Zb3010_SecEta5_bTagProbCSVBP[0]>0.898");
1065	>	EtaB=TCut("abs(Zb3010_SecEta5_pfJetGoodEta[0])<1.3");
1066	>	PhiZcut=TCut("abs(Zb3010_SecEta5_pfJetDphiZ[0])>2.7");
1067	>	recognized_scenario=true;
1068	>	}
1069	>
1070	>	if(identifier=="Zb3010_p5Clean") {
1071	>	LeadingB=TCut ("Zb3010_p5Clean_bTagProbCSVBP[0]>0.898 && Zb3010_p5Clean_IsClean");
1072	>	EtaB=TCut("abs(Zb3010_p5Clean_pfJetGoodEta[0])<1.3 && Zb3010_p5Clean_IsClean");
1073	>	PhiZcut=TCut("abs(Zb3010_p5Clean_pfJetDphiZ[0])>2.7 && Zb3010_p5Clean_IsClean");
1074	>	recognized_scenario=true;
1075	>	}
1076	>
1077	>	if(identifier=="Zb3010CHS") {
1078	>	LeadingB=TCut ("Zb3010CHS_bTagProbCSVBP[0]>0.898");
1079	>	EtaB=TCut("abs(Zb3010CHS_pfJetGoodEta[0])<1.3");
1080	>	PhiZcut=TCut("abs(Zb3010CHS_pfJetDphiZ[0])>2.7");
1081	>	recognized_scenario=true;
1082	>	}
1083	>
1084	>	assert(recognized_scenario);
1085	>
1086	>	cout << "Cuts: " << endl;
1087	>	cout << "   " << (const char*) LeadingB << endl;
1088	>	cout << "   " << (const char*) EtaB << endl;
1089	>	cout << "   " << (const char*) PhiZcut << endl;
1090	>
1091	>	cout << endl << endl;
1092	>
1093	>	//    ScenarioComparison();
1094	>	//    ScenarioComparisonInclusive();
1095	>
1096	>
1097	>	//     DoPUStudy(identifier);
1098	>
1099	>	//   SpecialCutFlow(identifier);
1100	>	//    draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",17000,30,200,"Z p_{T}","Official/"+identifier+"/Zpt__"+identifier,1);
1101	>	//   draw_kin_variable(identifier+"_pfJetGoodPt[0]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",40,0,200,"Leading Jet Pt (B)","Official/"+identifier+"/LeadingJetPt__"+identifier,1);
1102	>	//   draw_kin_variable(identifier+"_pfJetGoodPt[1]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",40,0,200,"Sub-Leading Jet Pt [GeV]","Official/"+identifier+"/Jet2Pt__"+identifier,1);
1103	>	//   draw_kin_variable(identifier+"_pfJetGoodPt[1]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,2,"p_{T}^{J2} / p_{T}^{Z}","Official/"+identifier+"/SubLeadingJetPt_Over_ZPt__"+identifier,1);
1104	>	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);
1105	>	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);
1106	>	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);
1107	>	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);
1108	>	draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut,20,0,2,"#alpha","Official/"+identifier+"/alpha_full__"+identifier,1);
1109	>	//   draw_kin_variable(identifier+"_pfBJetDphiZ[0]",ZplusBsel&&LeadingB&&"pt>10&&pt<1000&&pfBJetDphiZ[0]>0",30,0,3.2,"#delta#phi (Z,b lead)","Official/"+identifier+"/DeltaPhiZBlead__"+identifier,0);
1110	>
1111	>	//   draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&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, 10 GeV < p_{T}^{Z} < 1000 GeV");
1112	>	//   draw_kin_variable(identifier+"_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&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, 10 GeV < p_{T}^{Z} < 1000 GeV");
1113	>
1114	>
1115	>	//   */
1116	>
1117	>	}
1118	>
1119	>	void GetNumberEventsInsideOutsideAlphaWindow(TCut specialcut, string title) {
1120	>
1121	>	//  TCut cut =  ZplusBsel&&LeadingB&&PhiZcut&&specialcut&&TCut("Zb3010_alpha<0.3");
1122	>	TCut cut = specialcut;
1123	>	TCut in  =  EtaB;
1124	>	TCut out =  TCut("abs(Zb3010_pfJetGoodEta[0])>1.3");
1125	>
1126	>	TH1F *data_IN  = allsamples.Draw("data_IN",    "mll",1,0,150, "nothing [GeV]", "events", cut&&in, data,luminosity);
1127	>	TH1F *data_OUT = allsamples.Draw("data_OUT",   "mll",1,0,150, "nothing [GeV]", "events", cut&&out,data,luminosity);
1128	>	TH1F *mc_IN    = allsamples.Draw("mc_IN",      "mll",1,0,150, "nothing [GeV]", "events", cut&&in, mc  ,luminosity);
1129	>	TH1F *mc_OUT   = allsamples.Draw("mc_OUT",     "mll",1,0,150, "nothing [GeV]", "events", cut&&out,mc  ,luminosity);
1130	>
1131	>	dout << title << " : " << endl;
1132	>	dout << "  Data: In " << data_IN->Integral() << " vs out " << data_OUT->Integral() << endl;
1133	>	dout << "  MC  : In " << mc_IN->Integral()   << " vs out " << mc_OUT->Integral()   << endl;
1134	>	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;
1135	>	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;
1136	>
1137	>	delete data_IN;
1138	>	delete data_OUT;
1139	>	delete mc_IN;
1140	>	delete mc_OUT;
1141	>	}
1142	>
1143	>	void TEST_GetNumberEventsInsideOutsideAlphaWindow() {
1144	>	TCanvas *ca = new TCanvas("ca","ca");
1145	>	//   GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==0"), "ee");
1146	>	//   GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==1"), "mm");
1147	>	//   GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2"), "ee/mm");
1148	>
1149	>	cout << "BASE SELECTION" << endl;
1150	>	GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==0"), "ee");
1151	>	GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==1"), "mm");
1152	>	GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2"), "ee/mm");
1153	>
1154	>	cout << "BASE SELECTION: Z+b" << endl;
1155	>	GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
1156	>	GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
1157	>	GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2"), "ee/mm");
1158	>
1159	>	cout << "Leading B" << endl;
1160	>	GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
1161	>	GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
1162	>	GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2"), "ee/mm");
1163	>
1164	>	cout << "PhiZcut" << endl;
1165	>	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
1166	>	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
1167	>	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2"), "ee/mm");
1168	>
1169	>	cout << "alpha cut" << endl;
1170	>	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0")&&TCut("Zb3010_alpha<0.3"), "ee");
1171	>	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1")&&TCut("Zb3010_alpha<0.3"), "mm");
1172	>	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2")&&TCut("Zb3010_alpha<0.3"), "ee/mm");
1173	>
1174	>	delete ca;
1175	>
1176	>
1177	>	}
1178	>
1179	>
1180	>
1181	>
1182	>	void do_basic_ZB_analysis(bool DoCompleteAnalysis) {
1183	>
1184	>	//https://twiki.cern.ch/twiki/bin/viewauth/CMS/BtagPOG#2011_Data_and_MC
1185	>	//https://twiki.cern.ch/twiki/pub/CMS/BtagPOG/SFb-ttbar_payload.txt
1186	>
1187	>	cout << "The lepton requirement is " << (const char*) leptoncut << endl;
1188	>	bool doquick=true;
1189	>
1190	>	TCanvas *zbcanvas = new TCanvas("zbcanvas","zbcanvas");
1191	>
1192	>	//     compare_selection("Zb3010_p5Clean");
1193	>	//     compare_selection("Zb3010CHS");
1194	>	//     compare_selection("Zb3010");
1195	>	//     compare_selection("Zb301030");
1196	>	//     compare_selection("Zb1530");
1197	>	//     compare_selection("Zb301010");
1198	>	//     compare_selection("Zb3010_SecEta3");
1199	>	//     compare_selection("Zb3010_SecEta5");
1200	>	//     compare_selection("Zb1510");
1201	>	//     compare_selection("Zb301010");
1202	>	//     compare_selection("ZbMikko");
1203	>	//     compare_selection("Zb3010");
1204	>	//     compare_selection("Zb40");
1205
1206	<	void do_basic_ZB_analysis() {
1206	>	if(!doquick) {
1207	>	print_all_b_yields();
1208	>	draw_mpf_vars();
1209	>	//  DrawEvilCutFlow();
1210	>
1211	>	draw_Zb_kin_vars();
1212	>
1213	>	}
1214	>
1215	>	//  GetNumberEventsInsideOutsideAlphaWindow();
1216	>
1217	>	dout <<"THIS IS WHERE IT GETS REALLY INTERESTING!!!!" << endl;
1218	>	if(DoCompleteAnalysis) {
1219	>	/* need to run the following scenarios:
1220	>	* - normal
1221	>	* - inclusive
1222	>	* - medium WP
1223	>	* - alpha up
1224	>	* - alpha down
1225	>	*/
1226	>
1227	>	ZbResultContainer inclusive = new_data_mc_agreement_2d("Zb3010_alpha","inclusive",TCut("Zb3010_bTagProbCSVBP[0]>-10000"),TCut("1.0"));
1228	>	ZbResultContainer Reference = new_data_mc_agreement_2d("Zb3010_alpha","reference",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"));
1229	>	ZbResultContainer effmistagUp = new_data_mc_agreement_2d("Zb3010_alpha","effmistagUp",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtTUp"));
1230	>	ZbResultContainer effmistagDown = new_data_mc_agreement_2d("Zb3010_alpha","effmistagDown",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtTDown"));
1231	>	ZbResultContainer medium = new_data_mc_agreement_2d("Zb3010_alpha","medium",TCut("Zb3010_bTagProbCSVBP[0]>0.679"),TCut("ZbCHS3010_BTagWgtM"));
1232	>	ZbResultContainer loose = new_data_mc_agreement_2d("Zb3010_alpha","loose",TCut("Zb3010_bTagProbCSVBP[0]>0.244"),TCut("ZbCHS3010_BTagWgtL"));
1233	>	ZbResultContainer AlphaUp = new_data_mc_agreement_2d("Zb3010_alphaUp","AlphaUp",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"));
1234	>	ZbResultContainer AlphaDown = new_data_mc_agreement_2d("Zb3010_alphaDown","AlphaDown",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"));
1235	>
1236	>	//and now let's compute some systematics!
1237	>	//1 alpha variation
1238	>	float SysMPF_Alpha_down = Reference.MPF_Result_FaceValue.getValue()-AlphaDown.MPF_Result_FaceValue.getValue();
1239	>	float SysMPF_Alpha_up = Reference.MPF_Result_FaceValue.getValue()-AlphaUp.MPF_Result_FaceValue.getValue();
1240	>
1241	>	//2 btagging efficiency/mistag correction
1242	>	float SysMPF_EFF_down = Reference.MPF_Result_FaceValue.getValue()-effmistagDown.MPF_Result_FaceValue.getValue();
1243	>	float SysMPF_EFF_up = Reference.MPF_Result_FaceValue.getValue()-effmistagUp.MPF_Result_FaceValue.getValue();
1244	>
1245	>	//3 purity
1246	>	zbcanvas->cd();
1247	>	write_info(__FUNCTION__,"Don't know the exact purity at the moment, still needs to be determined!");
1248	>	TGraphErrors *gSysMPF_Purity = new TGraphErrors();
1249	>	gSysMPF_Purity->SetPoint(0,0.05,inclusive.MPF_Result_FaceValue.getValue());
1250	>	gSysMPF_Purity->SetPointError(0,0.0,inclusive.MPF_Result_FaceValue.getError());
1251	>	gSysMPF_Purity->SetPointError(1,0.20,loose.MPF_Result_FaceValue.getError());//x value?
1252	>	gSysMPF_Purity->SetPointError(1,0.0,loose.MPF_Result_FaceValue.getError());
1253	>	gSysMPF_Purity->SetPointError(2,0.54,medium.MPF_Result_FaceValue.getError());//x value?
1254	>	gSysMPF_Purity->SetPointError(2,0.0,medium.MPF_Result_FaceValue.getError());
1255	>	gSysMPF_Purity->SetPointError(3,0.82,Reference.MPF_Result_FaceValue.getError());//x value?
1256	>	gSysMPF_Purity->SetPointError(3,0.0,Reference.MPF_Result_FaceValue.getError());
1257	>
1258	>	gSysMPF_Purity->GetXaxis()->SetTitle("Purity");
1259	>	gSysMPF_Purity->GetXaxis()->CenterTitle();
1260	>	gSysMPF_Purity->GetYaxis()->SetTitle("C_{abs}");
1261	>	gSysMPF_Purity->GetYaxis()->CenterTitle();
1262	>
1263	>	gSysMPF_Purity->Fit("pol1");
1264	>	gSysMPF_Purity->Draw();
1265	>	CompleteSave(zbcanvas,"Systematics/Purity");
1266	>
1267	>	TF1 *PurityFit = (TF1*)gSysMPF_Purity->GetFunction("pol1");
1268	>	float ExtrapolatedResult=PurityFit->GetParameter(0)+1.0*PurityFit->GetParameter(1);
1269	>
1270	>	float SysMPF_Purity = abs(ExtrapolatedResult-Reference.MPF_Result_FaceValue.getValue());
1271	>
1272	>	float sys_alpha  = abs(SysMPF_Alpha_down)>abs(SysMPF_Alpha_up)?abs(SysMPF_Alpha_down):abs(SysMPF_Alpha_up);
1273	>	float sys_eff    = abs(SysMPF_EFF_down)>abs(SysMPF_EFF_up)?abs(SysMPF_EFF_down):abs(SysMPF_EFF_up);
1274	>	float sys_purity = abs(SysMPF_Purity);
1275	>	float sys        = sqrt(sys_alpha*sys_alpha+sys_eff*sys_eff+sys_purity*sys_purity);
1276	>
1277	>	dout << "MPF METHOD:" << endl;
1278	>	dout << "   Systematics : " << endl;
1279	>	dout << "      Alpha variation      : " << SysMPF_Alpha_down << "   ,  " << SysMPF_Alpha_up << "   ( " << sys_alpha << ") " << endl;
1280	>	dout << "      eff/mistag variation : " << SysMPF_EFF_down << "   ,  " << SysMPF_EFF_up << "   ( " << sys_eff << ") " << endl;
1281	>	dout << "      purity               : " << SysMPF_Purity << "   ( " << sys_purity << " )" << endl;
1282	>	dout << "      total                : " << sys << endl;
1283	>	dout << endl << endl;
1284	>	dout << "   FINAL RESULT : " << endl;
1285	>	dout << "     " << Reference.MPF_Result_FaceValue/inclusive.MPF_Result_FaceValue << " (stat) +/- " << sys << " (sys) " << endl;
1286	>
1287	>	LeadingB=TCut("Zb3010_bTagProbCSVBP[0]>0.898");
1288	>
1289	>	} else {
1290	>	ZbResultContainer Reference = new_data_mc_agreement_2d("Zb3010_alpha","reference",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgt"));
1291	>	}
1292
1293
1294	+	delete zbcanvas;
1295	+	}
1296	+
1297	+
1298	+
1299	+
1300	+	//*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*  DELETE EVERYTHING BELOW THIS LINE /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*//
1301	+
1302	+
1303	+	const char* concatenate(string bla1, string bla2) {
1304	+	stringstream bla;
1305	+	bla << bla1 << bla2;
1306	+	return bla.str().c_str();
1307	+	}
1308	+
1309	+	void SpecialCutFlow(string identifier) {
1310	+	stringstream MegaCut;
1311	+
1312	+
1313	+
1314	+	MegaCut << "   1*(" << (const char*) (ZplusBsel&&TCut(concatenate(identifier,"_pfJetGoodNumBtag>0"))) << ")";
1315	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB) << ")";
1316	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB) << ")";
1317	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut) << ")";
1318	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000").c_str())) << ")";
1319	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.3").c_str())) << ")";
1320	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.2").c_str())) << ")";
1321	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.15").c_str())) << ")";
1322	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.1").c_str())) << ")";
1323	+	MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.05").c_str())) << ")";
1324	+	MegaCut << " ";
1325	+
1326	+	TCanvas *can = new TCanvas("can","can");
1327	+	can->SetLogy(1);
1328	+	TCut basecut(ZplusBsel);
1329	+
1330	+
1331	+
1332	+	TLegend *leg = allsamples.allbglegend();
1333	+
1334	+
1335	+	TH1F *data  = allsamples.Draw("data",    MegaCut.str(),11,-0.5,10.5, "", "events", basecut,data,luminosity);
1336	+	THStack mcm = allsamples.DrawStack("mc", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,mc,luminosity);
1337	+
1338	+	float runningsum=0;
1339	+	for(int i=data->GetNbinsX();i>0;i--) {
1340	+	runningsum+=data->GetBinContent(i);
1341	+	data->SetBinContent(i,runningsum);
1342	+	}
1343	+
1344	+	float minimum=data->GetMaximum();
1345	+
1346	+	TH1F* h;
1347	+	TIter nextOF(mcm.GetHists());
1348	+
1349	+	float mcsum=0;
1350	+	float zb=0;
1351	+	while ( h = (TH1F*)nextOF() ) {
1352	+	float runningsum=0;
1353	+	minimum=data->GetMaximum();//done deliberately at each step so get the minimum of the last (!) contribution
1354	+	for(int i=h->GetNbinsX();i>0;i--) {
1355	+	runningsum+=h->GetBinContent(i);
1356	+	h->SetBinContent(i,runningsum);
1357	+	if(runningsum<minimum&&runningsum>0) minimum=runningsum;
1358	+	}
1359	+	if(Contains(h->GetName(),"Z_b_")) {
1360	+	zb=h->GetBinContent(h->GetNbinsX());
1361	+	}
1362	+	mcsum+=h->GetBinContent(h->GetNbinsX());
1363	+	}
1364	+
1365	+
1366	+
1367	+	data->SetMinimum(0.05*minimum);
1368	+	can->cd(1)->SetBottomMargin(0.25);
1369	+	data->GetXaxis()->SetBinLabel(1,"Pre-selection");
1370	+	data->GetXaxis()->SetBinLabel(2,"Contains b jet");
1371	+	data->GetXaxis()->SetBinLabel(3,"Leading jet is b-jet");
1372	+	data->GetXaxis()->SetBinLabel(4,"|#eta_{b}|<1.3 ");
1373	+	data->GetXaxis()->SetBinLabel(5,"|#delta#phi(b,Z)|>2.7");
1374	+	data->GetXaxis()->SetBinLabel(6,"10<p_{T}^{Z}<1000 GeV");
1375	+	data->GetXaxis()->SetBinLabel(7,"#alpha < 0.3");
1376	+	data->GetXaxis()->SetBinLabel(8,"#alpha < 0.2");
1377	+	data->GetXaxis()->SetBinLabel(9,"#alpha < 0.15");
1378	+	data->GetXaxis()->SetBinLabel(10,"#alpha < 0.1");
1379	+	data->GetXaxis()->SetBinLabel(11,"#alpha < 0.05");
1380	+	data->GetXaxis()->LabelsOption("v");
1381	+
1382	+	stringstream purityinfo;
1383	+	purityinfo << "Purity: " << std::setprecision(3) << 100*zb/mcsum << " %";
1384	+	stringstream neventsinfo;
1385	+	neventsinfo << "Nevents: " << data->GetBinContent(data->GetNbinsX());
1386	+	TH1F *crap = new TH1F("crap","",1,0,1);
1387	+	crap->SetLineColor(kWhite);
1388	+	leg->AddEntry(crap,purityinfo.str().c_str(),"l");
1389	+	leg->AddEntry(crap,neventsinfo.str().c_str(),"l");
1390	+
1391	+	data->Draw("e1");
1392	+	mcm.Draw("same");
1393	+	data->Draw("same,e1,axis");
1394	+	data->Draw("same,e1");
1395	+	//  data->Draw("same,TEXT");
1396	+
1397	+	leg->Draw();
1398	+
1399	+	DrawPrelim();
1400	+
1401	+	CompleteSave(can,"CutFlow___"+identifier);
1402	+
1403	+	delete crap;
1404	+	delete data;
1405
219	–	data_mc_agreement_2d();
1406
221	–	write_warning(__FUNCTION__,"Replaced B jet dphi cut with jet dphi cut - need to put this back!");
1407		}

Diff Legend

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