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.13 by buchmann, Fri Nov 23 12:10:58 2012 UTC vs.
Revision 1.14 by buchmann, Mon Nov 26 12:06:26 2012 UTC

#	Line 274 | Line 274 | Value get_Zb_data_over_mc(string variabl
274		if(Contains(variable,"pfJetGood")) varname="R_{abs}";
275		TH1F *hdata  = allsamples.Draw("hdata",variable,1200,0,30, varname, "events",  cut,data,luminosity);
276		TH1F *hmc    = allsamples.Draw("hmc" , variable,1200,0,30, varname, "events", cut, mc,  luminosity);
277	+	hmc->SetLineColor(kBlue);
278	+
279	+
280	+
281
282		float a=hdata->GetMean();
283		float b=hmc->GetMean();
#	Line 282 | Line 286 | Value get_Zb_data_over_mc(string variabl
286		float factor = a / b;
287		float error = TMath::Sqrt(  (1/(b*b)) * (da*da) + ((a*a)/(b*b))*db*db);
288
289	+	TF1 *gaus_data = new TF1("gaus_data","gaus",0.7,1.3);
290	+	gaus_data->SetParameter(0,hdata->GetMaximum());
291	+	gaus_data->SetParameter(2,hdata->GetMean());
292	+	gaus_data->SetParameter(1,hdata->GetRMS());
293	+	hdata->Fit(gaus_data);
294	+	TF1 *gaus_mc = new TF1("gaus_mc","gaus",0.7,1.3);
295	+	gaus_mc->SetParameter(0,hmc->GetMaximum());
296	+	gaus_mc->SetParameter(1,hmc->GetMean());
297	+	gaus_mc->SetParameter(2,hmc->GetRMS());
298	+	hmc->Fit(gaus_mc);
299	+
300	+
301		cn->cd();
302		hdata->GetYaxis()->SetTitle("events / 0.1");
303	<	hdata->Rebin(4);
303	>	hdata->SetMarkerColor(kRed);
304	>	/*  hdata->Rebin(4);
305		hmc->Rebin(4);
306	+
307	+	gaus_mc->SetParameter(0,gaus_mc->GetParameter(0)*4);//rebinning
308	+	gaus_data->SetParameter(0,gaus_data->GetParameter(0)*4);//rebinning
309	+	*/
310	+	gaus_mc->SetLineColor(kBlue);
311	+	gaus_data->SetLineColor(kRed);
312	+
313		hdata->GetXaxis()->SetRangeUser(0,2);
314		hdata->Draw("e1");
315		hmc->Draw("histo,same");
316		hdata->Draw("e1,same");
317	+	gaus_data->Draw("same");
318	+	gaus_mc->Draw("same");
319		stringstream summary;
320	<	summary << "#splitline{" << varname << "}{#splitline{data: " << std::setprecision(4) << a << " #pm " << std::setprecision(4) << da << "}{#splitline{mc:   " << std::setprecision(4) << b << " #pm " << std::setprecision(4) << db << "}{ratio: " << factor;
321	<	summary << " #pm " << error << "}}}";
320	>	summary << "#splitline{" << varname << "}{#splitline{data: " << std::setprecision(4) << a << " #pm " << std::setprecision(4) << da << "}{";
321	>	summary << "#splitline{mc:   " << std::setprecision(4) << b << " #pm " << std::setprecision(4) << db << "}{#splitline{ratio: " << factor;
322	>	summary << " #pm " << error << "}{#splitline{}{#splitline{data fit: " << std::setprecision(4) << gaus_data->GetParameter(1) << "+/-";
323	>	summary << std::setprecision(4) << gaus_data->GetParError(1) << "}{#splitline{";
324	>	summary << "mc fit:" << std::setprecision(4) << gaus_mc->GetParameter(1) << "+/-" << std::setprecision(4) << gaus_mc->GetParError(1) << "}{ratio: ";
325	>	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)) ) << "}}}}}}}";
326		TText *infobox = write_title(summary.str());
327		infobox->SetX(0.75);
328		infobox->SetTextSize(0.03);
#	Line 317 | Line 347 | Value get_Zb_data_over_mc(string variabl
347
348		void new_data_mc_agreement_2d() {
349
320	–	write_error(__FUNCTION__,"need to make alpha<0.3 inclusive again!!!! IMPORTANT SINCE THAT'S THE FINAL RESULT!!!");
321	–
350		gStyle->SetOptFit(0);
351		TCut basecut(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
352		const int nalphacuts=6;
353		float alphacuts[nalphacuts] = {0.1,0.15,0.2,0.25,0.3,0.35};
326	–	int PointThree=-1;
327	–	for(int i=0;i<nalphacuts;i++) {
328	–	if((alphacuts[i]-0.3)<0.01) PointThree=i-1;
329	–	}
330	–	if(PointThree<0) {
331	–	write_warning(__FUNCTION__,"alpha=0.3 MUST BE ONE OF THE ALPHA VALUES!");
332	–	assert(PointThree<0);
333	–	}
354
355		const int nptcuts=5;
356		float ptcuts[nptcuts]={30,50,75,125,1000};
#	Line 364 | Line 384 | void new_data_mc_agreement_2d() {
384		float RABS_ExtraPolatedResults[nptcuts];
385
386		TCanvas *can = new TCanvas("can");
387	+	can->cd(1)->SetLeftMargin(0.15);
388
389		TGraphErrors *MPF_FinalGraph = new TGraphErrors();
390		TGraphErrors *RABS_FinalGraph = new TGraphErrors();
#	Line 402 | Line 423 | void new_data_mc_agreement_2d() {
423		cout << "       value (RABS) : " << RABS_Results[ia][ipt] << " +/- " << RABS_Errors[ia][ipt] << endl;
424		}
425		}
405	–	cout << "About to add a point in pt (" << (ptcuts[ipt]+ptcuts[ipt+1])/2 << " ) to face value : " << MPF_Results[PointThree][ipt] << endl;
426
427		stringstream specialcut;
428		specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (Zb3010_alpha<0.3))";
429		Value MPF_data_over_mc = get_Zb_data_over_mc("mpf","",TCut(basecut && specialcut.str().c_str()),"MPF___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__INCLUSIVEalpha_0.3");
430		Value RABS_data_over_mc = get_Zb_data_over_mc("Zb3010_pfJetGoodPt[0]/pt","",TCut(basecut && specialcut.str().c_str()),"Rabs___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__INCLUSIVEalpha_0.3");
431	+	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;
432		MPF_FaceValueAtPoint3->SetPoint(ipt,(ptcuts[ipt]+ptcuts[ipt+1])/2,MPF_data_over_mc.getValue());
433		MPF_FaceValueAtPoint3->SetPointError(ipt,0,MPF_data_over_mc.getError());
434		RABS_FaceValueAtPoint3->SetPoint(ipt,(ptcuts[ipt]+ptcuts[ipt+1])/2,RABS_data_over_mc.getValue());
#	Line 423 | Line 444 | void new_data_mc_agreement_2d() {
444		mpf_gr->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
445		mpf_gr->GetYaxis()->CenterTitle();
446		TF1 *mpf_pol=(TF1*)mpf_gr->GetFunction("pol1");
447	+	mpf_pol->SetLineWidth(0);
448		TF1 *mpf_pol_complete = new TF1("mpf_pol_complete","[0]+[1]*x");
449		mpf_pol_complete->SetParameters(mpf_pol->GetParameters());
450		mpf_pol_complete->SetLineWidth(1);
#	Line 430 | Line 452 | void new_data_mc_agreement_2d() {
452		float min,max;
453		FindMinMax(mpf_gr,min,max);
454		TH1F *histo = new TH1F("histo","histo",1,0,0.4);
455	+	histo->GetXaxis()->SetTitle("#alpha");
456	+	histo->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
457	+	histo->GetXaxis()->CenterTitle();
458	+	histo->GetYaxis()->CenterTitle();
459		histo->SetStats(0);
460		histo->GetXaxis()->SetRangeUser(0,0.4);
461		histo->GetYaxis()->SetRangeUser(min,max);
#	Line 437 | Line 463 | void new_data_mc_agreement_2d() {
463
464		TPolyLine *mpf_fit_uncert = GetFitUncertaintyShape(mpf_gr, "pol1", min, max,0.0,0.4);
465		mpf_pol->SetLineColor(TColor::GetColor("#0101DF"));
466	<	mpf_fit_uncert->SetFillColor(TColor::GetColor("#0101DF"));
466	>	mpf_fit_uncert->SetFillColor(TColor::GetColor("#5353FC"));
467	>	histo->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
468		histo->Draw();
469	<	mpf_fit_uncert->Draw("same,f");
469	>	mpf_fit_uncert->Draw("F");
470	>	mpf_fit_uncert->Draw("");
471		mpf_gr->Draw("P");
472		mpf_pol_complete->Draw("same");
473
#	Line 451 | Line 479 | void new_data_mc_agreement_2d() {
479		MPF_ExtraPolatedResults[ipt]=mpf_a;
480
481		stringstream mpf_info;
482	<	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() << "}";
482	>	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() << "}";
483
484		TText *mpf_mark = write_title(mpf_info.str());
485		mpf_mark->SetX(0.75);
#	Line 473 | Line 501 | void new_data_mc_agreement_2d() {
501		rabs_gr->Fit("pol1","");
502		rabs_gr->GetXaxis()->SetTitle("#alpha");
503		rabs_gr->GetXaxis()->CenterTitle();
504	<	rabs_gr->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
504	>	rabs_gr->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
505		rabs_gr->GetYaxis()->CenterTitle();
506		TF1 *rabs_pol=(TF1*)rabs_gr->GetFunction("pol1");
507	+	rabs_pol->SetLineWidth(0);
508		TF1 *rabs_pol_complete = new TF1("rabs_pol_complete","[0]+[1]*x");
509		rabs_pol_complete->SetParameters(rabs_pol->GetParameters());
510		rabs_pol_complete->SetLineColor(kRed);
#	Line 484 | Line 513 | void new_data_mc_agreement_2d() {
513		rabs_gr->GetYaxis()->SetRangeUser(min,max);
514		histo->GetYaxis()->SetRangeUser(min,max);
515		TPolyLine *rabs_fit_uncert = GetFitUncertaintyShape(rabs_gr, "pol1", min, max,0.0,0.4);
516	<	rabs_pol->SetLineColor(TColor::GetColor("#0101DF"));
517	<	rabs_pol->SetFillColor(TColor::GetColor("#0101DF"));
516	>	rabs_pol->SetLineColor(TColor::GetColor("#FA5858"));
517	>	rabs_pol->SetFillColor(TColor::GetColor("#FA5858"));
518	>	rabs_fit_uncert->SetLineColor(TColor::GetColor("#FA5858"));
519	>	rabs_fit_uncert->SetFillColor(TColor::GetColor("#FA5858"));
520	>	histo->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
521		histo->Draw();
522	<	rabs_fit_uncert->Draw("same,f");
522	>	rabs_fit_uncert->Draw("F");
523	>	rabs_fit_uncert->Draw("");
524		rabs_gr->Draw("P");
525		rabs_pol_complete->Draw("same");
526
#	Line 500 | Line 533 | void new_data_mc_agreement_2d() {
533		cout << "!+!+!+!+!+!+!+!+!+ Just added a point to the final plot : " << rabs_a << " +/- " << rabs_pol->GetParError(0) << endl;
534
535		stringstream rabs_info;
536	<	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() << "}";
536	>	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() << "}";
537		TText *rabs_mark = write_title(rabs_info.str());
538		rabs_mark->SetX(0.75);
539		rabs_mark->SetTextSize(0.03);
#	Line 1011 | Line 1044 | void compare_selection(string identifier
1044
1045		}
1046
1047	<	void GetNumberEventsInsideOutsideAlphaWindow() {
1047	>	void GetNumberEventsInsideOutsideAlphaWindow(TCut specialcut, string title) {
1048
1049	<	TCut cut =  ZplusBsel&&LeadingB&&EtaB&&PhiZcut;
1050	<	TCut in  =  TCut("Zb3010_alpha<=0.3");
1051	<	TCut out =  TCut("Zb3010_alpha>0.3");
1049	>	//  TCut cut =  ZplusBsel&&LeadingB&&PhiZcut&&specialcut&&TCut("Zb3010_alpha<0.3");
1050	>	TCut cut = specialcut;
1051	>	TCut in  =  EtaB;
1052	>	TCut out =  TCut("abs(Zb3010_pfJetGoodEta[0])>1.3");
1053
1054		TH1F *data_IN  = allsamples.Draw("data_IN",    "mll",1,0,150, "nothing [GeV]", "events", cut&&in, data,luminosity);
1055		TH1F *data_OUT = allsamples.Draw("data_OUT",   "mll",1,0,150, "nothing [GeV]", "events", cut&&out,data,luminosity);
1056		TH1F *mc_IN    = allsamples.Draw("mc_IN",      "mll",1,0,150, "nothing [GeV]", "events", cut&&in, mc  ,luminosity);
1057		TH1F *mc_OUT   = allsamples.Draw("mc_OUT",     "mll",1,0,150, "nothing [GeV]", "events", cut&&out,mc  ,luminosity);
1058
1059	+	dout << title << " : " << endl;
1060		dout << "  Data: In " << data_IN->Integral() << " vs out " << data_OUT->Integral() << endl;
1061		dout << "  MC  : In " << mc_IN->Integral()   << " vs out " << mc_OUT->Integral()   << endl;
1062	+	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;
1063	+	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;
1064	+
1065	+	delete data_IN;
1066	+	delete data_OUT;
1067	+	delete mc_IN;
1068	+	delete mc_OUT;
1069		}
1070
1071	+	void GetNumberEventsInsideOutsideAlphaWindow() {
1072	+	TCanvas *ca = new TCanvas("ca","ca");
1073	+	//   GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==0"), "ee");
1074	+	//   GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==1"), "mm");
1075	+	//   GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2"), "ee/mm");
1076	+
1077	+	cout << "BASE SELECTION" << endl;
1078	+	GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==0"), "ee");
1079	+	GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==1"), "mm");
1080	+	GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2"), "ee/mm");
1081	+
1082	+	cout << "BASE SELECTION: Z+b" << endl;
1083	+	GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
1084	+	GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
1085	+	GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2"), "ee/mm");
1086	+
1087	+	cout << "Leading B" << endl;
1088	+	GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
1089	+	GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
1090	+	GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2"), "ee/mm");
1091	+
1092	+	cout << "PhiZcut" << endl;
1093	+	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
1094	+	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
1095	+	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2"), "ee/mm");
1096	+
1097	+	cout << "alpha cut" << endl;
1098	+	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0")&&TCut("Zb3010_alpha<0.3"), "ee");
1099	+	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1")&&TCut("Zb3010_alpha<0.3"), "mm");
1100	+	GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2")&&TCut("Zb3010_alpha<0.3"), "ee/mm");
1101	+
1102	+	delete ca;
1103	+
1104	+
1105	+	}
1106	+
1107	+
1108	+
1109
1110		void do_basic_ZB_analysis(int do_inclusive) {
1111
1112		//https://twiki.cern.ch/twiki/bin/viewauth/CMS/BtagPOG#2011_Data_and_MC
1113		//https://twiki.cern.ch/twiki/pub/CMS/BtagPOG/SFb-ttbar_payload.txt
1114
1115	+	cout << "The lepton requirement is " << (const char*) leptoncut << endl;
1116		if(do_inclusive) {
1117		LeadingB=TCut("Zb3010_bTagProbCSVBP[0]>=-10000");
1118	<	cutWeight=TCut("(weight*(weight<1000)*(is_data+(!is_data)*((id1==id2&&id1==0)*0.95+(id1==id2&&id1==1)*0.88+(id1!=id2)*0.92)))");
1118	>	cutWeight=TCut("(weight*(weight<1000)*(is_data+(!is_data)))");
1119		write_info(__FUNCTION__,"Inclusive selection - no btag requirement, no btag efficiency/mistag correction (weight)");
1120		} else {
1121		LeadingB=TCut("Zb3010_bTagProbCSVBP[0]>0.898");
#	Line 1069 | Line 1150 | void do_basic_ZB_analysis(int do_inclusi
1150
1151		}
1152
1153	<	GetNumberEventsInsideOutsideAlphaWindow();
1154	<	//  new_data_mc_agreement_2d();
1153	>	//  GetNumberEventsInsideOutsideAlphaWindow();
1154	>	new_data_mc_agreement_2d();
1155
1075	–	write_error(__FUNCTION__,"Need to implement switch for inclusive result; inclusive result needs to have btag efficiency reweighting deactivated!");
1156		delete zbcanvas;
1157		}
1158

Diff Legend

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