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.22 by buchmann, Mon Jan 21 14:43:36 2013 UTC vs.
Revision 1.29 by buchmann, Fri Apr 26 07:18:47 2013 UTC

#	Line 17 | Line 17
17		#include <TProfile.h>
18		#include <TLegendEntry.h>
19
20	+	#include <RooRealVar.h>
21	+	#include <RooArgSet.h>
22	+	#include <RooDataSet.h>
23	+	#include <RooDataHist.h>
24	+	#include <RooMCStudy.h>
25	+	#include <RooCategory.h>
26	+
27	+	#include <RooPlot.h>
28	+	#include <RooSimultaneous.h>
29	+	#include <RooAddPdf.h>
30	+	#include <RooFitResult.h>
31	+	#include <RooVoigtian.h>
32	+	#include <RooChebychev.h>
33	+	#include <RooMsgService.h>
34	+
35	+	#include <RooStats/ModelConfig.h>
36	+	#include "RooStats/ProfileLikelihoodCalculator.h"
37	+	#include "RooStats/LikelihoodInterval.h"
38	+	#include "RooStats/HypoTestResult.h"
39	+	#include "RooStats/SimpleLikelihoodRatioTestStat.h"
40	+	#include "RooStats/ProfileLikelihoodTestStat.h"
41	+	#include "RooStats/HybridCalculatorOriginal.h"
42	+	#include "RooStats/HypoTestInverterOriginal.h"
43	+
44	+
45	+
46		using namespace std;
47
48		using namespace PlottingSetup;
#	Line 26 | Line 52 | namespace ZbData {
52		TCut LeadingB("ZbCHS3010_bTagProbCSVBP[0]>0.898");
53		TCut EtaB("abs(ZbCHS3010_pfJetGoodEta[0])<1.3");
54		TCut PhiZcut("abs(ZbCHS3010_pfJetDphiZ[0])>2.7");
55	<
55	>	float LowerHistoBoundary=0.5; // where to cut on MPF (or Rabs)
56	>	float UpperHistoBoundary=1.5; // where to cut on MPF (or Rabs)
57	>
58		}
59
60		TCut STDWEIGHT;
#	Line 60 | Line 88 | class ZbResultContainer {
88		public:
89		string name;
90		vector<ZbPointResult> MPF_Rabs_location;
91	+
92		Value MPF_Result_FaceValue;
93		Value Rabs_Result_FaceValue;
94		Value MPF_Result_Extrapolated;
95		Value Rabs_Result_Extrapolated;
96
97	+	vector<Value> MPF_Result_FaceValue_Bins;
98	+	vector<Value> Rabs_Result_FaceValue_Bins;
99	+	vector<Value> MPF_Result_Extrapolated_Bins;
100	+	vector<Value> Rabs_Result_Extrapolated_Bins;
101	+
102	+	void Print();
103	+
104		ZbResultContainer(string name);
69	–	ZbResultContainer(const ZbResultContainer &old);
105		};
106
107
108		ZbResultContainer::ZbResultContainer(string _name) {
109		name=_name;
110	<	cout<< "Result Container for \"" << name << "\" has been initialized" << endl;
76	<
110	>	dout<< "Result Container for \"" << name << "\" has been initialized" << endl;
111		}
112
113	<	std::ostream &operator<<(std::ostream &ostr, ZbResultContainer &b)
113	>	void ZbResultContainer::Print()
114		{
115	<	ostr << "Results for " << b.name << endl;
116	<	ostr << " MPF and Rabs locations: " << endl;
117	<	for(int i=0;i<b.MPF_Rabs_location.size();i++) ostr << "   " << b.MPF_Rabs_location[i] << endl;
118	<	ostr << " Results: " << endl;
119	<	ostr << "   Face Value: " << endl;
120	<	ostr << "     MPF " << b.MPF_Result_FaceValue << endl;
121	<	ostr << "     Rabs " << b.Rabs_Result_FaceValue << endl;
122	<	ostr << "   Extrapolated: " << endl;
123	<	ostr << "     MPF " << b.MPF_Result_Extrapolated << endl;
124	<	ostr << "             Rabs " << b.Rabs_Result_Extrapolated << endl;
125	<
126	<	return ostr;
115	>	dout << "Results for " << this->name << endl;
116	>	dout << " MPF and Rabs locations: " << endl;
117	>	for(int i=0;i<this->MPF_Rabs_location.size();i++) dout << "   " << this->MPF_Rabs_location[i] << endl;
118	>	dout << " Results: " << endl;
119	>	dout << "   Face Value: " << endl;
120	>	dout << "     MPF " << this->MPF_Result_FaceValue << endl;
121	>	dout << "        Bins: [ ";
122	>	if(MPF_Result_FaceValue_Bins.size()>0) {
123	>	for(int i=0;i<MPF_Result_FaceValue_Bins.size()-1;i++) dout << MPF_Result_FaceValue_Bins[i] << " , ";
124	>	if(MPF_Result_FaceValue_Bins.size()>0) dout << MPF_Result_FaceValue_Bins[MPF_Result_FaceValue_Bins.size()-1] << " ]" << endl;
125	>	}
126	>
127	>	dout << "     Rabs " << this->Rabs_Result_FaceValue << endl;
128	>	if(Rabs_Result_FaceValue_Bins.size()>0) {
129	>	dout << "        Bins: [ ";
130	>	for(int i=0;i<Rabs_Result_FaceValue_Bins.size()-1;i++) dout << Rabs_Result_FaceValue_Bins[i] << " , ";
131	>	if(Rabs_Result_FaceValue_Bins.size()>0) dout << Rabs_Result_FaceValue_Bins[Rabs_Result_FaceValue_Bins.size()-1] << " ]" << endl;
132	>	}
133	>
134	>	dout << "   Extrapolated: " << endl;
135	>	dout << "     MPF " << this->MPF_Result_Extrapolated << endl;
136	>	if(MPF_Result_Extrapolated_Bins.size()>0) {
137	>	dout << "        Bins: [ ";
138	>	for(int i=0;i<MPF_Result_Extrapolated_Bins.size()-1;i++) dout << MPF_Result_Extrapolated_Bins[i] << " , ";
139	>	if(MPF_Result_Extrapolated_Bins.size()>0) dout << MPF_Result_Extrapolated_Bins[MPF_Result_Extrapolated_Bins.size()-1] << " ]" << endl;
140	>	}
141	>
142	>	dout << "             Rabs " << this->Rabs_Result_Extrapolated << endl;
143	>	if(Rabs_Result_Extrapolated_Bins.size()>0) {
144	>	dout << "        Bins: [ ";
145	>	for(int i=0;i<Rabs_Result_Extrapolated_Bins.size()-1;i++) dout << Rabs_Result_Extrapolated_Bins[i] << " , ";
146	>	if(Rabs_Result_Extrapolated_Bins.size()>0) dout << Rabs_Result_Extrapolated_Bins[Rabs_Result_Extrapolated_Bins.size()-1] << " ]" << endl;
147	>	}
148		}
149
150
96	–
151		using namespace ZbData;
152
153
#	Line 161 | Line 215 | CutYields print_yield(TCut cut) {
215		if(Contains(name,"Z_b_")) Yield.Zb+=hist->Integral();
216		}
217
218	<	cout << Yield << endl;
219	<	cout << endl << endl;
218	>	dout << Yield << endl;
219	>	dout << endl << endl;
220
221		delete data;
222		CleanLegends();
#	Line 225 | Line 279 | void draw_kin_variable(string variable,
279		speciall->Draw();
280		}
281
282	<	save_with_ratio(datahisto,mcstack,kinpad->cd(),saveas);
282	>	Save_With_Ratio(datahisto,mcstack,kinpad->cd(),saveas);
283
284		datahisto->Delete();
285		delete kinleg;
#	Line 236 | Line 290 | void draw_kin_variable(string variable,
290
291		void print_all_b_yields() {
292		vector<CutYields> yields;
293	<	cout << "Basic selection with a b jet" << endl;
293	>	dout << "Basic selection with a b jet" << endl;
294		yields.push_back(print_yield(ZplusBsel&&"ZbCHS3010_pfJetGoodNumBtag>0"));
295	<	cout << "Leading jet is a b-jet " << endl;
295	>	dout << "Leading jet is a b-jet " << endl;
296		yields.push_back(print_yield(ZplusBsel&&LeadingB));
297	<	cout << "|#eta_{b}|<1.3 " << endl;
297	>	dout << "|#eta_{b}|<1.3 " << endl;
298		yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB));
299	<	cout << "|#delta#phi(b<Z)|>2.7" << endl;
299	>	dout << "|#delta#phi(b<Z)|>2.7" << endl;
300		yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut));
301	<	cout << "30<ptZ<1000 GeV" << endl;
301	>	dout << "30<ptZ<1000 GeV" << endl;
302		yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"));
303	<	cout << "#alpha < 0.35" << endl;
303	>	dout << "#alpha < 0.35" << endl;
304		yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.35"));
305	<	cout << "#alpha < 0.3" << endl;
305	>	dout << "#alpha < 0.3" << endl;
306		yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.3"));
307	<	/*  cout << "#alpha < 0.2" << endl;
307	>	/*  dout << "#alpha < 0.2" << endl;
308		yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.2"));
309	<	cout << "#alpha < 0.15" << endl;
309	>	dout << "#alpha < 0.15" << endl;
310		yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.15"));
311	<	cout << "#alpha < 0.1" << endl;
311	>	dout << "#alpha < 0.1" << endl;
312		yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.1"));*/
313
314		for(int iy=0;iy<yields.size();iy++) {
315	<	cout << yields[iy] << endl;
315	>	dout << yields[iy] << endl;
316		}
317		}
318
#	Line 345 | Line 399 | void DrawEvilCutFlow() {
399		void draw_Zb_kin_vars() {
400
401		draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000"),25,0,200,"Z p_{T}","Official/Zpt",1);
402	+	draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==0"),25,0,200,"Z p_{T}","Official/Zpt_ee",1);
403	+	draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==1"),25,0,200,"Z p_{T}","Official/Zpt_mm",1);
404	+	draw_kin_variable("pt",ZplusBsel&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000"),25,0,200,"Z p_{T}","Official/Zpt__nonBtagged",1);
405	+	draw_kin_variable("pt",ZplusBsel&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==0"),25,0,200,"Z p_{T}","Official/Zpt__nonBtagged_ee",1);
406	+	draw_kin_variable("pt",ZplusBsel&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==1"),25,0,200,"Z p_{T}","Official/Zpt__nonBtagged_mm",1);
407	+
408		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);
409	+	draw_kin_variable("ZbCHS3010_pfJetGoodPt[1]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==0"),25,0,200,"Sub-Leading Jet Pt","Official/Jet2Pt_ee",1);
410	+	draw_kin_variable("ZbCHS3010_pfJetGoodPt[1]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==1"),25,0,200,"Sub-Leading Jet Pt","Official/Jet2Pt_mm",1);
411	+	draw_kin_variable("ZbCHS3010_pfJetGoodPt[1]",ZplusBsel&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000"),25,0,200,"Sub-Leading Jet Pt","Official/Jet2Pt__nonBtagged",1);
412	+	draw_kin_variable("ZbCHS3010_pfJetGoodPt[1]",ZplusBsel&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==0"),25,0,200,"Sub-Leading Jet Pt","Official/Jet2Pt__nonBtagged_ee",1);
413	+	draw_kin_variable("ZbCHS3010_pfJetGoodPt[1]",ZplusBsel&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==1"),25,0,200,"Sub-Leading Jet Pt","Official/Jet2Pt__nonBtagged_mm",1);
414	+
415		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);
416	+	draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==0"),25,0,200,"Leading Jet Pt (B)","Official/LeadingJetPt_ee",1);
417	+	draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==1"),25,0,200,"Leading Jet Pt (B)","Official/LeadingJetPt_mm",1);
418	+	draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]",ZplusBsel&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000"),25,0,200,"Leading Jet Pt (B)","Official/LeadingJetPt__nonBtagged",1);
419	+	draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]",ZplusBsel&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==0"),25,0,200,"Leading Jet Pt (B)","Official/LeadingJetPt__nonBtagged_ee",1);
420	+	draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]",ZplusBsel&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==1"),25,0,200,"Leading Jet Pt (B)","Official/LeadingJetPt__nonBtagged_mm",1);
421	+
422	+
423		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);
424		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);
425		draw_kin_variable("ZbCHS3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>50&&pt<75"),  40,0,2,"#alpha","Official/alpha_50_to_75",1);
#	Line 371 | Line 444 | void draw_mpf_vars() {
444		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");
445		}
446
447	+	void make_RooFit(string ContainerName, TH1F *h, float &peak, float &error_peak, string saveas) {
448	+	TCanvas *fcan = new TCanvas("fcan","fcan");
449	+	RooRealVar mpf("mpf","mpf", -5, 5, "");
450	+	RooDataHist AllData("AllData", "CMS Real Dataset", mpf , h);
451	+
452	+	RooRealVar fraction("fraction", "fraction", 0.8, 0, 1);
453	+	RooRealVar mean("mean", "mean", 1, 0, 3);
454	+	RooRealVar sigma("sigma", "sigma", 1, 0, 5);
455	+	RooRealVar width("width", "width", 2.4, 0, 5);
456	+	RooRealVar chevvar1("chevvar1", "chevvar1", 0.2,-1.0, 1.0);
457	+	RooRealVar chevvar2("chevvar2", "chevvar2", 0.2,-1.0, 1.0);
458	+
459	+	RooVoigtian signal("signal", "signal", mpf, mean, width, sigma);
460	+	RooChebychev background("background_", "background", mpf, RooArgList(chevvar1, chevvar2));
461	+
462	+	RooAddPdf pdf("pdf",  "pdf", RooArgList(signal, background), fraction);
463	+
464	+	RooFitResult *result = pdf.fitTo(AllData, RooFit::Save());
465	+	RooPlot* frame = mpf.frame(RooFit::Bins(30), RooFit::Title("MPF")) ;
466	+	AllData.plotOn(frame);
467	+	pdf.plotOn(frame, RooFit::LineColor(kBlue)) ;
468	+	pdf.plotOn(frame, RooFit::Components(RooArgSet(background)), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed)) ;
469	+	pdf.plotOn(frame, RooFit::Components(RooArgSet(signal)), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen)) ;
470	+
471	+	RooRealVar* thePeak = (RooRealVar*) result->floatParsFinal().find("mean");
472	+	peak = thePeak->getVal();
473	+	error_peak = thePeak->getError();
474	+
475	+	///******************* just making it beautiful ....
476	+	RooPlot* frame2 = mpf.frame(RooFit::Bins(30), RooFit::Title("MPF")) ;
477	+	TH1F *hcopy = (TH1F*)h->Clone("hcopy");
478	+	hcopy->Rebin(40);
479	+	RooDataHist AllRData("AllRData", "CMS Dataset", mpf , hcopy);
480	+
481	+	AllRData.plotOn(frame2);
482	+	pdf.plotOn(frame2, RooFit::LineColor(kBlue)) ;
483	+	pdf.plotOn(frame2, RooFit::Components(RooArgSet(background)), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed)) ;
484	+	pdf.plotOn(frame2, RooFit::Components(RooArgSet(signal)), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen)) ;
485	+
486	+	string saveasFinal = saveas + "_Fit";
487	+	fcan->cd() ; gPad->SetLeftMargin(0.15) ; frame->GetYaxis()->SetTitleOffset(1.4) ; frame2->Draw();
488	+
489	+	if(Contains(saveas,"_data")) DrawPrelim();
490	+	else DrawMCPrelim();
491	+
492	+	CompleteSave(fcan, ContainerName+"/Fit/"+saveasFinal);
493	+
494	+	delete hcopy;
495	+	delete fcan;
496	+	delete frame;
497	+	delete thePeak;
498	+	}
499	+
500	+	float bother=0.1;
501	+	float bother2=1;
502		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) {
375	–	//  write_warning(__FUNCTION__,"Debugging this function, therefore always returning 3.1415+/-1.010101"); return Value(3.1415,1.010101);
503		TCanvas *cn = new TCanvas("cn","cn");
504		string varname="MPF";
505
506	+	if(a2<15.0/pt2 && a1>0) dout << "this bin ( " << a1 << " < alpha < " << a2 <<" and " << pt1 << " < p_{T}^{Z} < " << pt2 <<" ) is not expected to do a lot." << endl;
507	+	//note that "15.0" here is the subleading jet cut.
508	+
509		bool DoFit=false;
510		bool UseFit=false;
511	+	bool DoRooFit=false;
512
513
514		if(Contains(variable,"pfJetGood")) varname="R_{abs}";
515		//  TH1F *hdata  = allsamples.Draw("hdata",variable,1200,0,30, varname, "events",  cut,data,luminosity);
516		//  TH1F *hmc    = allsamples.Draw("hmc" , variable,1200,0,30, varname, "events", cut, mc,  luminosity);
517	<	TH1F *hdata  = allsamples.Draw("hdata",variable,1200,0.5,1.5, varname, "events",  cut,data,luminosity); // making it more precise
518	<	TH1F *hmc    = allsamples.Draw("hmc" , variable,1200,0.5,1.5, varname, "events", cut, mc,  luminosity); // making it more precise
517	>	TH1F *hdata, *hmc;
518	>	if(!DoRooFit) {
519	>	hdata = allsamples.Draw("hdata",variable,1200,ZbData::LowerHistoBoundary,ZbData::UpperHistoBoundary, varname, "events",  cut,data,luminosity); // making it more precise
520	>	hmc    = allsamples.Draw("hmc" , variable,1200,ZbData::LowerHistoBoundary,ZbData::UpperHistoBoundary, varname, "events", cut, mc,  luminosity); // making it more precise
521	>	} else {
522	>	hdata  = allsamples.Draw("hdata",variable,7200,0.0,6.0, varname, "events",  cut,data,luminosity); // making it more precise
523	>	hmc    = allsamples.Draw("hmc" , variable,7200,0.0,6.0, varname, "events", cut, mc,  luminosity); // making it more precise
524	>	}
525	>
526		hmc->SetLineColor(kBlue);
527
528		float a=hdata->GetMean();
#	Line 407 | Line 545 | Value get_Zb_data_over_mc(string Contain
545		gaus_mc->SetParameter(2,hmc->GetRMS());
546		hmc->Fit(gaus_mc);
547		}
548	<
548	>	if(DoRooFit) {
549	>
550	>	if(hdata->Integral()>10) {
551	>	make_RooFit(ContainerName, hdata, a, da, saveas+"_data");
552	>	make_RooFit(ContainerName, hmc, b, db, saveas+"_mc");
553	>	factor = a / b;
554	>	error = TMath::Sqrt(  (1/(b*b)) * (da*da) + ((a*a)/(b*b))*db*db);
555	>	} else {
556	>	factor=NAN;//set to nan
557	>	error=NAN;//set to nan
558	>	}
559	>	}
560	>
561		cn->cd();
562	<	hdata->Rebin(0.05*1200/(1.5-0.5));
563	<	hmc->Rebin(0.05*1200/(1.5-0.5));
562	>	hdata->Rebin(int(0.05*1200/(ZbData::UpperHistoBoundary-ZbData::LowerHistoBoundary)));
563	>	hmc->Rebin(int(0.05*1200/(ZbData::UpperHistoBoundary-ZbData::LowerHistoBoundary)));
564		hdata->GetYaxis()->SetTitle("events / 0.05");
565		hdata->SetMarkerColor(kRed);
566		if(DoFit) {
#	Line 463 | Line 613 | Value get_Zb_data_over_mc(string Contain
613		CompleteSave(cn,ContainerName+"/ResponsePlots/"+saveas);
614
615		res.MPF_Rabs_location.push_back(ZbPointResult(variable,pt1,pt2,a1,a2,Value(a,da),Value(b,db)));
616	<	cout << "Data;" << a << ";"<<da<<";MC;"<<b<<";"<<db<<endl;
616	>	dout << "Data;" << a << ";"<<da<<";MC;"<<b<<";"<<db<<endl;
617		delete cn;
618		delete hdata;
619		delete hmc;
#	Line 479 | Line 629 | float PurityMedium=0.54;
629		float PurityTight=0.82;
630
631
632	+	void GetPtPointPosition(float AvPos[], const int nptcuts,string AlphaVariable, TCut basecut, float FaceValueToBeTakenAt, float ptcuts[]) {
633	+	for(int i=0;i<nptcuts-1;i++) {
634	+	stringstream specialcut;
635	+	specialcut << "((pt>30&& pt< 1000) && (" << AlphaVariable << "<" << FaceValueToBeTakenAt << "))";
636	+	TH1F *hdata  = allsamples.Draw("hada", "pt",100,ptcuts[i],ptcuts[i+1], "MPF", "events", TCut(basecut && specialcut.str().c_str()), data,  luminosity);
637	+	AvPos[i]=hdata->GetMean();
638	+	dout << AvPos[i] << " : ";
639	+	delete hdata;
640	+	}
641	+	dout << endl;
642	+	}
643	+
644	+
645		void DeterminePurity(TCut basecut, int nptcuts, float ptcuts[], string AlphaVariable, string ContainerName, float FaceValueToBeTakenAt) {
646
647		bool ComputePurity=false;
#	Line 515 | Line 678 | void DeterminePurity(TCut basecut, int n
678		}
679
680		float purity=Zb/total;
681	<	dout << "Determined a purity of " << 100*purity << " % for the container " << ContainerName << endl;
681	>	dout << "Determined a purity of " << 100*purity << " % for the container \"" << ContainerName << "\"" << endl;
682
683		if(ContainerName=="inclusive") PurityInclusive=purity;
684		if(ContainerName=="loose") PurityLoose=purity;
#	Line 534 | Line 697 | string NiceAlphaPrint(float alpha) {
697
698		ZbResultContainer new_data_mc_agreement_2d(bool gofast, string AlphaVariable, string ContainerName, TCut BTagCut, TCut BTagWgt, float FaceValueToBeTakenAt) {
699
537	–
700		if(gofast) write_info(__FUNCTION__,"Fast mode activated for "+ContainerName);
701		ZbResultContainer results(ContainerName);
702
#	Line 550 | Line 712 | ZbResultContainer new_data_mc_agreement_
712
713		const int nptcuts=5;
714		float ptcuts[nptcuts]={30,50,75,125,1000};
715	+	float ptPointPosition[nptcuts];
716	+	GetPtPointPosition(ptPointPosition,nptcuts,AlphaVariable,basecut,FaceValueToBeTakenAt,ptcuts);
717	+
718
719		DeterminePurity(basecut, nptcuts,ptcuts, AlphaVariable, ContainerName, FaceValueToBeTakenAt);
720
#	Line 572 | Line 737 | ZbResultContainer new_data_mc_agreement_
737		alow=0;
738		ahigh=alphacuts[ia];
739		}
740	<	cout << specialcut.str() << endl;
740	>	dout << specialcut.str() << endl;
741		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(NiceAlphaPrint(alphacuts[ia])), ptcuts[ipt], ptcuts[ipt+1],alow,ahigh);
742		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(NiceAlphaPrint(alphacuts[ia])), ptcuts[ipt], ptcuts[ipt+1],alow,ahigh);
743
#	Line 582 | Line 747 | ZbResultContainer new_data_mc_agreement_
747		RABS_Results[ia][ipt]=RABS_data_over_mc.getValue();
748		RABS_Errors[ia][ipt]=RABS_data_over_mc.getError();
749
750	<	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;
750	>	dout << "alpha in [" << alow << " , " << ahigh << " 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;
751		}
752		}
753
#	Line 601 | Line 766 | ZbResultContainer new_data_mc_agreement_
766
767		stringstream filename;
768		filename << "Extrapolation/Zplusb_data_over_mc___" << ptcuts[ipt] << "_to_" << ptcuts[ipt+1];
769	<	cout << "Going to create histo for " << filename.str() << endl;
769	>	dout << "Going to create histo for " << filename.str() << endl;
770		TGraphErrors *mpf_gr =new TGraphErrors();
771		TGraphErrors *rabs_gr =new TGraphErrors();
772	+	mpf_gr->SetName("MPF_Graph");
773	+	rabs_gr->SetName("Rabs_Graph");
774
775		int rabs_counter=0;
776		int mpf_counter=0;
777
778		for(int ia=0;ia<nalphacuts && !gofast;ia++) {
779	<	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;
779	>	dout << "Setting a point for an alpha cut at " << alphacuts[ia] << " with value " << MPF_Results[ia][ipt] << " (MPF) and " << RABS_Results[ia][ipt] << " (RABS)" << endl;
780		if(MPF_Results[ia][ipt]==MPF_Results[ia][ipt] && MPF_Errors[ia][ipt]==MPF_Errors[ia][ipt]) { // remove nan's
781		mpf_gr->SetPoint(mpf_counter,alphacuts[ia],MPF_Results[ia][ipt]);
782		mpf_gr->SetPointError(mpf_counter,0,MPF_Errors[ia][ipt]);
#	Line 617 | Line 784 | ZbResultContainer new_data_mc_agreement_
784		} else {
785		dout << "Detected BS for MPF method (bin ignored) : " << endl;
786		dout << "       alpha: " << alphacuts[ia] << " , pt range: " << ptcuts[ipt] << " < pt < " << ptcuts[ipt+1] << endl;
787	<	cout << "       value (MPF) : " << MPF_Results[ia][ipt] << " +/- " << MPF_Errors[ia][ipt] << endl;
787	>	dout << "       value (MPF) : " << MPF_Results[ia][ipt] << " +/- " << MPF_Errors[ia][ipt] << endl;
788		}
789		if(RABS_Results[ia][ipt]==RABS_Results[ia][ipt] && RABS_Errors[ia][ipt]==RABS_Errors[ia][ipt]) {
790		rabs_gr->SetPoint(rabs_counter,alphacuts[ia],RABS_Results[ia][ipt]);
#	Line 626 | Line 793 | ZbResultContainer new_data_mc_agreement_
793		} else {
794		dout << "Detected BS for RABS method (bin ignored) : " << endl;
795		dout << "       alpha: " << alphacuts[ia] << " , pt range: " << ptcuts[ipt] << " < pt < " << ptcuts[ipt+1] << endl;
796	<	cout << "       value (RABS) : " << RABS_Results[ia][ipt] << " +/- " << RABS_Errors[ia][ipt] << endl;
796	>	dout << "       value (RABS) : " << RABS_Results[ia][ipt] << " +/- " << RABS_Errors[ia][ipt] << endl;
797		}
798		}
799
#	Line 634 | Line 801 | ZbResultContainer new_data_mc_agreement_
801		specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (" << AlphaVariable << "<" << FaceValueToBeTakenAt << "))";
802		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);
803		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);
804	<	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;
804	>	dout << "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;
805		MPF_FaceValueAtPoint3->SetPoint(ipt,(ptcuts[ipt]+ptcuts[ipt+1])/2,MPF_data_over_mc.getValue());
806		MPF_FaceValueAtPoint3->SetPointError(ipt,0,MPF_data_over_mc.getError());
807		RABS_FaceValueAtPoint3->SetPoint(ipt,(ptcuts[ipt]+ptcuts[ipt+1])/2,RABS_data_over_mc.getValue());
808		RABS_FaceValueAtPoint3->SetPointError(ipt,0,RABS_data_over_mc.getError());
809
810
811	+	TH1D *mpf_band, *rabs_band;
812
813		if(!gofast) {
814		can->cd();
815		mpf_gr->SetMarkerStyle(21);
816		mpf_gr->SetMarkerColor(kBlue);
817		mpf_gr->Draw("AP");
818	<	mpf_gr->Fit("pol1","");
818	>
819	>	TF1 *mpf_pol = new TF1("mpf_pol","[0]+[1]*x",0.0,0.4);
820	>
821	>	mpf_gr->Fit(mpf_pol); // better error estimation using MINOS (E), quiet (Q) and using loglikelihood method
822	>	mpf_band = getBand(mpf_pol,"MPF_band");
823		mpf_gr->GetXaxis()->SetTitle("#alpha");
824		mpf_gr->GetXaxis()->CenterTitle();
825		mpf_gr->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
826		mpf_gr->GetYaxis()->CenterTitle();
827	<	TF1 *mpf_pol=(TF1*)mpf_gr->GetFunction("pol1");
827	>	mpf_band->SetMarkerSize(0);
828		mpf_pol->SetLineWidth(0);
829	+	mpf_pol->SetLineColor(mpf_band->GetFillColor());
830		TF1 *mpf_pol_complete = new TF1("mpf_pol_complete","[0]+[1]*x");
831		mpf_pol_complete->SetParameters(mpf_pol->GetParameters());
832	+
833		mpf_pol_complete->SetLineWidth(1);
834		mpf_pol_complete->SetLineColor(kBlue);
835
#	Line 672 | Line 846 | ZbResultContainer new_data_mc_agreement_
846		histo->GetYaxis()->SetRangeUser(min,max);
847		mpf_gr->GetYaxis()->SetRangeUser(min,max);
848
849	<	TPolyLine *mpf_fit_uncert = GetFitUncertaintyShape(mpf_gr, "pol1", min, max,0.0,0.4);
849	>	//TPolyLine *mpf_fit_uncert = GetFitUncertaintyShape(mpf_gr, "pol1", min, max,0.0,0.4);
850	>
851		mpf_pol->SetLineColor(TColor::GetColor("#0101DF"));
852	<	mpf_fit_uncert->SetFillColor(TColor::GetColor("#5353FC"));
852	>	//mpf_fit_uncert->SetFillColor(TColor::GetColor("#A2A2FA"));
853	>	mpf_band->SetFillColor(TColor::GetColor("#A2A2FA"));
854		histo->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
855		histo->Draw();
856	<	mpf_fit_uncert->Draw("F");
857	<	mpf_fit_uncert->Draw("");
856	>	//mpf_fit_uncert->Draw("F");
857	>	//mpf_fit_uncert->Draw("");
858	>	mpf_gr->Draw("P");
859	>	mpf_band->Draw("C E3 SAME");
860		mpf_gr->Draw("P");
861		mpf_pol_complete->Draw("same");
862
863		float mpf_a=mpf_pol->GetParameter(0);
864		float mpf_b=mpf_pol->GetParameter(1);
865	<	MPF_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),mpf_a);
866	<	MPF_FinalGraph->SetPointError(ipt,0,mpf_pol->GetParError(0));
865	>	MPF_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),mpf_pol->GetParameter(0));
866	>	MPF_FinalGraph->SetPointError(ipt,0,mpf_band->GetBinError(mpf_band->FindBin(0)));
867
868		MPF_ExtraPolatedResults[ipt]=mpf_a;
869
#	Line 702 | Line 880 | ZbResultContainer new_data_mc_agreement_
880		filename << "__MPF";
881		DrawPrelim();
882		CompleteSave(can,ContainerName+"/"+filename.str());
883	<	cout << "MPF : " << mpf_a << " + " << mpf_b << " * alpha " << endl;
883	>	dout << "MPF : " << mpf_a << " + " << mpf_b << " * alpha " << endl;
884
885		filename.str("");
886
887	+	TF1 *rabs_pol = new TF1("rabs_pol","[0]+[1]*x",0.0,0.4);
888	+	rabs_gr->Fit(rabs_pol,"E");
889	+	rabs_band = getBand(rabs_pol,"Rabs_band");
890	+	rabs_pol->SetLineColor(rabs_band->GetFillColor());
891		rabs_gr->SetMarkerStyle(21);
892		rabs_gr->SetMarkerColor(kRed);
893		rabs_gr->Draw("AP");
894	<	rabs_gr->Fit("pol1","");
894	>	//rabs_gr->Fit("pol1","E");
895	>	rabs_band->SetMarkerSize(0);
896	>
897	>
898		rabs_gr->GetXaxis()->SetTitle("#alpha");
899		rabs_gr->GetXaxis()->CenterTitle();
900		rabs_gr->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
901		rabs_gr->GetYaxis()->CenterTitle();
902	<	TF1 *rabs_pol=(TF1*)rabs_gr->GetFunction("pol1");
902	>	//TF1 *rabs_pol=(TF1*)rabs_gr->GetFunction("pol1");
903		rabs_pol->SetLineWidth(0);
904		TF1 *rabs_pol_complete = new TF1("rabs_pol_complete","[0]+[1]*x");
905		rabs_pol_complete->SetParameters(rabs_pol->GetParameters());
#	Line 723 | Line 908 | ZbResultContainer new_data_mc_agreement_
908		FindMinMax(rabs_gr,min,max);
909		rabs_gr->GetYaxis()->SetRangeUser(min,max);
910		histo->GetYaxis()->SetRangeUser(min,max);
911	<	TPolyLine *rabs_fit_uncert = GetFitUncertaintyShape(rabs_gr, "pol1", min, max,0.0,0.4);
912	<	rabs_pol->SetLineColor(TColor::GetColor("#FA5858"));
913	<	rabs_pol->SetFillColor(TColor::GetColor("#FA5858"));
914	<	rabs_fit_uncert->SetLineColor(TColor::GetColor("#FA5858"));
915	<	rabs_fit_uncert->SetFillColor(TColor::GetColor("#FA5858"));
911	>	//TPolyLine *rabs_fit_uncert = GetFitUncertaintyShape(rabs_gr, "pol1", min, max,0.0,0.4);
912	>	rabs_pol->SetLineColor(TColor::GetColor("#F99C9C"));
913	>	rabs_pol->SetFillColor(TColor::GetColor("#F99C9C"));
914	>	//      rabs_band->SetLineColor(TColor::GetColor("#F99C9C"));
915	>	rabs_band->SetFillColor(TColor::GetColor("#F99C9C"));
916	>	//rabs_fit_uncert->SetLineColor(TColor::GetColor("#F99C9C"));
917	>	//rabs_fit_uncert->SetFillColor(TColor::GetColor("#F99C9C"));
918		histo->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
919		histo->Draw();
920	<	rabs_fit_uncert->Draw("F");
921	<	rabs_fit_uncert->Draw("");
920	>	rabs_band->Draw("C E3 SAME");
921	>	//rabs_fit_uncert->Draw("F");
922	>	//rabs_fit_uncert->Draw("");
923		rabs_gr->Draw("P");
924		rabs_pol_complete->Draw("same");
925
926		float rabs_a=rabs_pol->GetParameter(0);
927		float rabs_b=rabs_pol->GetParameter(1);
928	+
929		RABS_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),rabs_a);
930	<	RABS_FinalGraph->SetPointError(ipt,0,rabs_pol->GetParError(0));
930	>	RABS_FinalGraph->SetPointError(ipt,0,rabs_band->GetBinError(rabs_band->FindBin(0)));
931
932	<	cout << "!+!+!+!+!+!+!+!+!+ Just added a point to the final plot : " << rabs_a << " +/- " << rabs_pol->GetParError(0) << endl;
932	>	dout << "!+!+!+!+!+!+!+!+!+ Just added a point to the final plot : " << rabs_a << " +/- " << rabs_band->GetBinError(rabs_band->FindBin(0)) << endl;
933
934		stringstream rabs_info;
935		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() << "}";
#	Line 754 | Line 943 | ZbResultContainer new_data_mc_agreement_
943		filename << filenamebkp << "__RABS";
944		DrawPrelim();
945		CompleteSave(can,ContainerName+"/"+filename.str());
946	<	cout << "RABS : " << rabs_a << " + " << rabs_b << " * alpha " << endl;
946	>	dout << "RABS : " << rabs_a << " + " << rabs_b << " * alpha " << endl;
947	>
948	>	delete mpf_band;
949	>	delete rabs_band;
950		}
951		}
952
#	Line 890 | Line 1082 | ZbResultContainer new_data_mc_agreement_
1082		filename2 << "Extrapolation/FINAL_RESULT_RABS_FaceValp3";
1083		CompleteSave(can,ContainerName+"/"+filename2.str());
1084
1085	<	cout << "FINAL RESULTS: " << endl;
1086	<	if(!gofast) cout << "MPF : " << MPFResult  << " +/- " << MPFResultError  << endl;
1087	<	if(!gofast) cout << "Rabs: " << RabsResult << " +/- " << RabsResultError << endl;
1088	<	cout << "Face value at " << FaceValueToBeTakenAt << " :" << FaceValResult << " +/- " << FaceValResultError << "  (MPF) " << endl;
1089	<	cout << "Face value at " << FaceValueToBeTakenAt << " :" << RABSFaceValResult << " +/- " << RABSFaceValResultError << "  (RABS) " << endl;
1085	>	dout << "FINAL RESULTS: " << endl;
1086	>	if(!gofast) dout << "MPF : " << MPFResult  << " +/- " << MPFResultError  << endl;
1087	>	if(!gofast) dout << "Rabs: " << RabsResult << " +/- " << RabsResultError << endl;
1088	>	dout << "Face value at " << FaceValueToBeTakenAt << " :" << FaceValResult << " +/- " << FaceValResultError << "  (MPF) " << endl;
1089	>	dout << "Face value at " << FaceValueToBeTakenAt << " :" << RABSFaceValResult << " +/- " << RABSFaceValResultError << "  (RABS) " << endl;
1090
1091		delete can;
1092
1093	+
1094	+
1095	+	for(unsigned int i=0;i<MPF_FinalGraph->GetN();i++) {
1096	+	double x,y,dy;
1097	+	MPF_FinalGraph->GetPoint(i,x,y);
1098	+	dy = MPF_FinalGraph->GetErrorY(i);
1099	+	results.MPF_Result_Extrapolated_Bins.push_back(Value(y,dy));
1100	+	}
1101	+
1102	+	for(unsigned int i=0;i<MPF_FaceValueAtPoint3->GetN();i++) {
1103	+	double x,y,dy;
1104	+	MPF_FaceValueAtPoint3->GetPoint(i,x,y);
1105	+	dy = MPF_FaceValueAtPoint3->GetErrorY(i);
1106	+	results.MPF_Result_FaceValue_Bins.push_back(Value(y,dy));
1107	+	}
1108	+
1109	+	for(unsigned int i=0;i<RABS_FinalGraph->GetN();i++) {
1110	+	double x,y,dy;
1111	+	RABS_FinalGraph->GetPoint(i,x,y);
1112	+	dy = RABS_FinalGraph->GetErrorY(i);
1113	+	results.Rabs_Result_Extrapolated_Bins.push_back(Value(y,dy));
1114	+	}
1115	+
1116	+	for(unsigned int i=0;i<RABS_FaceValueAtPoint3->GetN();i++) {
1117	+	double x,y,dy;
1118	+	RABS_FaceValueAtPoint3->GetPoint(i,x,y);
1119	+	dy = RABS_FaceValueAtPoint3->GetErrorY(i);
1120	+	results.Rabs_Result_FaceValue_Bins.push_back(Value(y,dy));
1121	+	}
1122	+
1123	+
1124		results.MPF_Result_FaceValue=Value(FaceValResult,FaceValResultError);
1125		results.Rabs_Result_FaceValue=Value(RABSFaceValResult,RABSFaceValResultError);
1126		if(!gofast) results.MPF_Result_Extrapolated=Value(MPFResult,MPFResultError);
1127		if(!gofast) results.Rabs_Result_Extrapolated=Value(RabsResult,RabsResultError);
1128	<	cout << results << endl;
1128	>	results.Print();
1129		return results;
1130		}
1131
#	Line 914 | Line 1137 | void DoPUStudy(string identifier) {
1137		TCanvas *can = new TCanvas("can","can");
1138		TH1F *malpha[8];
1139		TH1F *dalpha[8];
1140	<	cout << identifier << ";";
1141	<	cout << endl;
1140	>	dout << identifier << ";";
1141	>	dout << endl;
1142		for(int i=1;i<5;i++) {
1143		stringstream sSpecialCut;
1144		if(i<4) {
1145		sSpecialCut << "numVtx>" << numVtxcuts[i-1] << " && numVtx<" << numVtxcuts[i];
1146	<	cout << numVtxcuts[i-1] << "  < nVtx < " << numVtxcuts[i] << ";";
1146	>	dout << numVtxcuts[i-1] << "  < nVtx < " << numVtxcuts[i] << ";";
1147		} else {
1148		sSpecialCut << "numVtx>" << numVtxcuts[i-1];
1149	<	cout << numVtxcuts[i-1] << ";";
1149	>	dout << numVtxcuts[i-1] << ";";
1150		}
1151
1152		sSpecialCut << " && " << identifier << "_alpha<0.3";
#	Line 941 | Line 1164 | void DoPUStudy(string identifier) {
1164		float factor = a / b;
1165		float error = TMath::Sqrt(  (1/(b*b)) * (da*da) + ((a*a)/(b*b))*db*db);
1166
1167	<	cout << malpha[i]->GetMean() << ";" << malpha[i]->GetMeanError() << ";" << dalpha[i]->GetMean() << ";" << dalpha[i]->GetMeanError() << ";" << malpha[i]->Integral()<<";"<<dalpha[i]->Integral() << ";"<<factor << ";" << error <<  endl;
1167	>	dout << malpha[i]->GetMean() << ";" << malpha[i]->GetMeanError() << ";" << dalpha[i]->GetMean() << ";" << dalpha[i]->GetMeanError() << ";" << malpha[i]->Integral()<<";"<<dalpha[i]->Integral() << ";"<<factor << ";" << error <<  endl;
1168		delete malpha[i];
1169		delete dalpha[i];
1170		}
#	Line 1022 | Line 1245 | void ScenarioComparison() {
1245		}
1246
1247
1248	+	gStyle->SetOptFit(0);
1249		gr[i]->SetTitle();
1250		gr[i]->GetXaxis()->SetTitle("N(Vtx)");
1251		gr[i]->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
#	Line 1154 | Line 1378 | void ScenarioComparisonInclusive() {
1378		void compare_selection(string identifier) {
1379		bool recognized_scenario=false;
1380
1381	<	cout << "Running with identifier " << identifier << endl;
1381	>	dout << "Running with identifier " << identifier << endl;
1382		if(identifier=="Zb30") {
1383		LeadingB=TCut ("Zb30_bTagProbCSVBP[0]>0.898");
1384		EtaB=TCut("abs(Zb30_pfJetGoodEta[0])<1.3");
#	Line 1266 | Line 1490 | void compare_selection(string identifier
1490
1491		assert(recognized_scenario);
1492
1493	<	cout << "Cuts: " << endl;
1494	<	cout << "   " << (const char*) LeadingB << endl;
1495	<	cout << "   " << (const char*) EtaB << endl;
1496	<	cout << "   " << (const char*) PhiZcut << endl;
1493	>	dout << "Cuts: " << endl;
1494	>	dout << "   " << (const char*) LeadingB << endl;
1495	>	dout << "   " << (const char*) EtaB << endl;
1496	>	dout << "   " << (const char*) PhiZcut << endl;
1497
1498	<	cout << endl << endl;
1498	>	dout << endl << endl;
1499
1500		//  ScenarioComparison();
1501		//  ScenarioComparisonInclusive();
#	Line 1326 | Line 1550 | void TEST_GetNumberEventsInsideOutsideAl
1550		//   GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==1"), "mm");
1551		//   GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2"), "ee/mm");
1552
1553	<	cout << "BASE SELECTION" << endl;
1553	>	dout << "BASE SELECTION" << endl;
1554		GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==0"), "ee");
1555		GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==1"), "mm");
1556		GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2"), "ee/mm");
1557
1558	<	cout << "BASE SELECTION: Z+b" << endl;
1558	>	dout << "BASE SELECTION: Z+b" << endl;
1559		GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
1560		GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
1561		GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2"), "ee/mm");
1562
1563	<	cout << "Leading B" << endl;
1563	>	dout << "Leading B" << endl;
1564		GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
1565		GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
1566		GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2"), "ee/mm");
1567
1568	<	cout << "PhiZcut" << endl;
1568	>	dout << "PhiZcut" << endl;
1569		GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
1570		GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
1571		GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2"), "ee/mm");
1572
1573	<	cout << "alpha cut" << endl;
1573	>	dout << "alpha cut" << endl;
1574		GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0")&&TCut("ZbCHS3010_alpha<0.3"), "ee");
1575		GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1")&&TCut("ZbCHS3010_alpha<0.3"), "mm");
1576		GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2")&&TCut("ZbCHS3010_alpha<0.3"), "ee/mm");
#	Line 1357 | Line 1581 | void TEST_GetNumberEventsInsideOutsideAl
1581		}
1582
1583
1584	<
1584	>	void MetSpectrumDepletionIllustration(int isdata) {
1585	>
1586	>	TCanvas *can = new TCanvas("can","can");
1587	>	TCut basecut(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
1588	>
1589	>	cutWeight=STDWEIGHT*TCut("ZbCHS3010_BTagWgtTUp*(1.0)*(ZbCHS3010_bTagProbCSVBP[0]>0.898)");
1590	>	TH1F *standard  = allsamples.Draw("standard","met[1]",80,0,400, "MET [GeV]", "events", basecut,isdata,luminosity);
1591	>
1592	>	cutWeight=STDWEIGHT*TCut("ZbCHS3010_BTagWgtTUp*(1.0*softMuon)*(ZbCHS3010_bTagProbCSVBP[0]>0.898)");
1593	>	TH1F *Neutrinos  = allsamples.Draw("Neutrinos","met[1]",80,0,400, "MET [GeV]", "events", basecut,isdata,luminosity);
1594	>
1595	>	cutWeight=STDWEIGHT*TCut("ZbCHS3010_BTagWgtTUp*(1.0*(1.0-softMuon))*(ZbCHS3010_bTagProbCSVBP[0]>0.898)");
1596	>	TH1F *NoNeutrinos  = allsamples.Draw("NoNeutrinos","met[1]",80,0,400, "MET [GeV]", "events", basecut,isdata,luminosity);
1597	>
1598	>	standard->Scale(1.0/standard->Integral());
1599	>	Neutrinos->Scale(1.0/Neutrinos->Integral());
1600	>	NoNeutrinos->Scale(1.0/NoNeutrinos->Integral());
1601	>
1602	>	Neutrinos->SetLineColor(kGreen);
1603	>	NoNeutrinos->SetLineColor(kRed);
1604	>
1605	>	can->SetLogy(1);
1606	>	TLegend *leg = make_legend();
1607	>	leg->AddEntry(standard,"Standard sample","ep");
1608	>	leg->AddEntry(Neutrinos,"#nu enriched sample","l");
1609	>	leg->AddEntry(NoNeutrinos,"#nu depleted sample","l");
1610	>	leg->SetX1(0.9*leg->GetX1());
1611	>	standard->Draw("e1");
1612	>	Neutrinos->Draw("histo,same");
1613	>	NoNeutrinos->Draw("histo,same");
1614	>	if(isdata==data) leg->SetHeader("Data:");
1615	>	else leg->SetHeader("MC:");
1616	>	leg->Draw();
1617	>	DrawPrelim();
1618	>
1619	>	if(isdata==data) CompleteSave(can,"NeutrinoMetSpectrum__Data");
1620	>	else CompleteSave(can,"NeutrinoMetSpectrum__MC");
1621	>	delete standard;
1622	>	delete Neutrinos;
1623	>	delete NoNeutrinos;
1624	>	delete can;
1625	>	CleanLegends();
1626	>	}
1627	>
1628	>	string PrintVector(vector<float> Systematics) {
1629	>	stringstream fullanswer;
1630	>	fullanswer << " [ ";
1631	>	if(Systematics.size()>0) fullanswer << Systematics[0];
1632	>	for(unsigned int is=1;is<Systematics.size();is++) fullanswer << " , " << Systematics[is];
1633	>	fullanswer << " ] ";
1634	>	return fullanswer.str();
1635	>	}
1636	>
1637
1638		void do_basic_ZB_analysis(bool DoCompleteAnalysis) {
1639		STDWEIGHT=TCut(cutWeight);
1640		cutWeight=TCut(STDWEIGHT*TCut("ZbCHS3010_BTagWgtT"));
1641	<	cout << "The lepton requirement is " << (const char*) leptoncut << endl;
1641	>	dout << "The lepton requirement is " << (const char*) leptoncut << endl;
1642		bool doquick=false;
1643
1644		TCanvas *zbcanvas = new TCanvas("zbcanvas","zbcanvas");
1645
1646	<	//   ScenarioComparison();
1647	<	//   ScenarioComparisonInclusive();
1648	<
1646	>
1647	>	//   MetSpectrumDepletionIllustration(data);
1648	>	//   MetSpectrumDepletionIllustration(mc);
1649	>
1650		//   compare_selection("ZbCHS3010_p5Clean");
1651		//   compare_selection("ZbCHS3010");
1652		//   compare_selection("Zb30");
#	Line 1388 | Line 1665 | void do_basic_ZB_analysis(bool DoComplet
1665		//   compare_selection("Zb40");
1666
1667		if(!doquick) {
1668	<	//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.)
1669	<	//print_all_b_yields();
1668	>	//    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.)
1669	>	//    print_all_b_yields();
1670		//     draw_mpf_vars();
1671	<	//DrawEvilCutFlow();
1671	>	//    DrawEvilCutFlow();
1672
1673		//     draw_Zb_kin_vars();
1674
#	Line 1410 | Line 1687 | void do_basic_ZB_analysis(bool DoComplet
1687
1688		bool fast=true;bool slow=false;//don't change these
1689
1690	<	ZbResultContainer inclusive = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","inclusive",TCut("mll>0"),TCut("1.0"),0.3);
1690	>	ZbResultContainer inclusive = new_data_mc_agreement_2d(slow,"ZbCHS3010_alpha","inclusive",TCut("mll>0"),TCut("1.0"),0.3);
1691		ZbResultContainer Reference = new_data_mc_agreement_2d(slow,"ZbCHS3010_alpha","reference",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"),0.3);
1692
1693	+
1694		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);
1695		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);
1696
#	Line 1434 | Line 1712 | void do_basic_ZB_analysis(bool DoComplet
1712		ZbResultContainer AlphaThresholdVariationUp10inc =   new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationUp10inc",  TCut("mll>0"),                              TCut("1.0"),0.3+0.1*0.3);
1713		ZbResultContainer AlphaThresholdVariationDown30inc = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationDown30inc",TCut("mll>0"),                              TCut("1.0"),0.3-0.3*0.3);
1714		ZbResultContainer AlphaThresholdVariationUp30inc =   new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationUp30inc",  TCut("mll>0"),                              TCut("1.0"),0.3+0.3*0.3);
1715	+
1716	+
1717	+	ZbResultContainer OnlyOneJet   = new_data_mc_agreement_2d(slow,"ZbCHS3010_alpha","OnlyOneJet",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898&&ZbCHS3010_pfJetGoodNum==1"),TCut("ZbCHS3010_BTagWgtT"),0.3);
1718	+	ZbResultContainer MultipleJets = new_data_mc_agreement_2d(slow,"ZbCHS3010_alpha","MultipleJets",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898&&ZbCHS3010_pfJetGoodNum>1"),TCut("ZbCHS3010_BTagWgtT"),0.3);
1719	+
1720	+	ZbResultContainer OnlyOneJetinc   = new_data_mc_agreement_2d(slow,"ZbCHS3010_alpha","OnlyOneJetinc",TCut("ZbCHS3010_pfJetGoodNum==1"),TCut("1.0"),0.3);
1721	+	ZbResultContainer MultipleJetsinc = new_data_mc_agreement_2d(slow,"ZbCHS3010_alpha","MultipleJetsinc",TCut("ZbCHS3010_pfJetGoodNum>1"),TCut("1.0"),0.3);
1722
1723		//and now let's compute some systematics!
1724		//1 alpha variation
1725		float SysMPF_Alpha_down = Reference.MPF_Result_FaceValue.getValue()-AlphaDown.MPF_Result_FaceValue.getValue();
1726		float SysMPF_Alpha_up = Reference.MPF_Result_FaceValue.getValue()-AlphaUp.MPF_Result_FaceValue.getValue();
1727
1728	+	vector<float> SysMPF_Alpha_down_Bins;
1729	+	vector<float> SysMPF_Alpha_up_Bins;
1730	+	for(unsigned int iFV=0;iFV<Reference.MPF_Result_FaceValue_Bins.size();iFV++) SysMPF_Alpha_down_Bins.push_back(Reference.MPF_Result_FaceValue_Bins[iFV].getValue()-AlphaDown.MPF_Result_FaceValue_Bins[iFV].getValue());
1731	+	for(unsigned int iFV=0;iFV<Reference.MPF_Result_FaceValue_Bins.size();iFV++) SysMPF_Alpha_up_Bins.push_back(Reference.MPF_Result_FaceValue_Bins[iFV].getValue()-AlphaUp.MPF_Result_FaceValue_Bins[iFV].getValue());
1732	+
1733	+
1734		//2 btagging efficiency/mistag correction
1735		float SysMPF_EFF_down = Reference.MPF_Result_FaceValue.getValue()-effmistagDown.MPF_Result_FaceValue.getValue();
1736		float SysMPF_EFF_up = Reference.MPF_Result_FaceValue.getValue()-effmistagUp.MPF_Result_FaceValue.getValue();
1737
1738	+	vector<float> SysMPF_EFF_down_Bins;
1739	+	vector<float> SysMPF_EFF_up_Bins;
1740	+	for(unsigned int iFV=0;iFV<Reference.MPF_Result_FaceValue_Bins.size();iFV++) SysMPF_EFF_down_Bins.push_back(Reference.MPF_Result_FaceValue_Bins[iFV].getValue()-effmistagDown.MPF_Result_FaceValue_Bins[iFV].getValue());
1741	+	for(unsigned int iFV=0;iFV<Reference.MPF_Result_FaceValue_Bins.size();iFV++) SysMPF_EFF_up_Bins.push_back(Reference.MPF_Result_FaceValue_Bins[iFV].getValue()-effmistagUp.MPF_Result_FaceValue_Bins[iFV].getValue());
1742	+
1743		//3 Presence of soft muons (->neutrino question)
1744		float SysMPF_Neutrino_up = Reference.MPF_Result_FaceValue.getValue()-NeutrinoQ.MPF_Result_FaceValue.getValue();
1745		float SysMPF_Neutrino_down = Reference.MPF_Result_FaceValue.getValue()-ANeutrinoQ.MPF_Result_FaceValue.getValue();
1746
1747	+	vector<float> SysMPF_Neutrino_up_Bins;
1748	+	vector<float> SysMPF_Neutrino_down_Bins;
1749	+	for(unsigned int iFV=0;iFV<Reference.MPF_Result_FaceValue_Bins.size();iFV++) SysMPF_Neutrino_up_Bins.push_back(Reference.MPF_Result_FaceValue_Bins[iFV].getValue()-NeutrinoQ.MPF_Result_FaceValue_Bins[iFV].getValue());
1750	+	for(unsigned int iFV=0;iFV<Reference.MPF_Result_FaceValue_Bins.size();iFV++) SysMPF_Neutrino_down_Bins.push_back(Reference.MPF_Result_FaceValue_Bins[iFV].getValue()-ANeutrinoQ.MPF_Result_FaceValue_Bins[iFV].getValue());
1751	+
1752	+
1753		//4 purity
1754		zbcanvas->cd();
1755		gStyle->SetOptFit(0);
#	Line 1510 | Line 1812 | void do_basic_ZB_analysis(bool DoComplet
1812
1813		CompleteSave(zbcanvas,"Systematics/Purity");
1814
1815	+	dout << "Bin-by-bin systematics : " << endl;
1816	+	dout << "    Oversmearing : " << endl;
1817	+	dout << "      down: " << PrintVector(SysMPF_Alpha_down_Bins) << endl;
1818	+	dout << "      up  : " << PrintVector(SysMPF_Alpha_up_Bins) << endl;
1819	+	dout << "    eff/mistag variation: " << endl;
1820	+	dout << "      down: " << PrintVector(SysMPF_EFF_down_Bins) << endl;
1821	+	dout << "      up  : " << PrintVector(SysMPF_EFF_up_Bins) << endl;
1822	+	dout << "    neutrinos: " << endl;
1823	+	dout << "      down: " << PrintVector(SysMPF_Neutrino_down_Bins) << endl;
1824	+	dout << "      up  : " << PrintVector(SysMPF_Neutrino_up_Bins) << endl;
1825
1826
1827		float SysMPF_Purity = abs(ExtrapolatedResult-Reference.MPF_Result_FaceValue.getValue());
#	Line 1554 | Line 1866 | void do_basic_ZB_analysis(bool DoComplet
1866		dout << "             +30%: " << AlphaThresholdVariationUp30.Rabs_Result_FaceValue/AlphaThresholdVariationUp30inc.Rabs_Result_FaceValue << " (stat)" << endl;
1867
1868
1869	+	dout << "  Also checked what happens when you only use events with one jet and only with more than one jets: " << endl;
1870	+	dout << "      Only 1      : " << OnlyOneJet.MPF_Result_FaceValue/OnlyOneJetinc.MPF_Result_FaceValue << " (stat) " << endl;
1871	+	dout << "      More than 1 : " << MultipleJets.MPF_Result_FaceValue/MultipleJetsinc.MPF_Result_FaceValue << " (stat) " << endl;
1872	+
1873		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());
1874		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());
1875		float sys_alphavar = sys_alphavar_up>sys_alphavar_down?sys_alphavar_up:sys_alphavar_down;

Diff Legend

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